Genetic loci associated with freezing tolerance in a European rapeseed (Brassica napus L.) diversity panel identified by genome-wide association mapping.

Winter biotypes of rapeseed (Brassica napus L.) require a vernalization treatment to enter the reproductive phase and generally produce greater yields than spring rapeseed. To find genetic loci associated with freezing tolerance in rapeseed, we first performed genotyping-by-sequencing (GBS) on a diversity panel consisting of 222 rapeseed accessions originating primarily from Europe, which identified 69,554 high-quality single-nucleotide polymorphisms (SNPs). Model-based cluster analysis suggested that there were eight subgroups. The diversity panel was then phenotyped for freezing survival (visual damage and Fv/Fo and Fv/Fm) after 2 months of cold acclimation (5°C) and a freezing treatment (-15°C for 4 h). The genotypic and phenotypic data for each accession in the rapeseed diversity panel was then used to conduct a genome-wide association study (GWAS). GWAS results showed that 14 significant markers were mapped to seven chromosomes for the phenotypes scored. Twenty-four candidate genes located within the mapped loci were identified as previously associated with lipid, photosynthesis, flowering, ubiquitination, and cytochrome P450 in rapeseed or other plant species. Published 2022. This article is a U.S. Government work and is in the public domain in the USA. Plant Direct published by American Society of Plant Biologists and the Society for Experimental Biology and John Wiley & Sons Ltd.

INTRODUCTION

Canola (Brassica napus L.) is the second major oilseed crop cultivated worldwide after soybeans with a grain production of 70.4 million tons in the marketing year of 2020/2021 (UFOP, 2021). Canola is an allotetraploid (also called amphidiploid, AACC; 2n = 4x = 38) and has a genome size of 1129–1235 Mb (Chalhoub et al., 2014) and evolved through multiple hybridization events from two diploid progenitor species, Brassica rapa L. (AA, 2n = 2x = 20) and Brassica oleracea L. (CC, 2n = 2x = 18) (Allender & King, 2010; Song et al., 1988; Song & Osborn, 1992).

Two biotypes (spring and winter) of canola are planted in the United States, depending upon their geographic adaptation and winter survival. Canola biotypes that require an extended period of cold to flower (vernalization) are referred as winter canola and are generally planted in the fall for grain harvest the following year. Canola biotypes that are planted in early spring and can flower without a vernalization period are referred as spring canola. Winter canola commonly gets 20%–30% better yield than spring canola (Bushong et al., 2018). Grain yields of winter canola normally surpass 2000 kg ha−1 and can achieve a yield of up to 5000 kg ha−1 (Stamm et al., 2019). However, the winter survival of canola is generally poor and unpredictable in the northern USA (Gusta & Wisniewski, 2013), limiting growth of yearly winter canola acreage.

Winter hardiness of canola is a complex trait requiring both cold acclimation and freezing tolerance (Rife & Zeinali, 2003). The process of cold acclimation comprises a series of coordinated physiological and biochemical alterations (Miura & Furumoto, 2013) and the reprogramming of gene expression networks (Fowler & Thomashow, 2002; Kreps et al., 2002; Seki et al., 2001, 2002). Cold acclimation allows plants to survive periods of freeze/thaw cycles and to resist osmotic dehydration generated by extracellular ice nucleation under freezing conditions (Hincha & Zuther, 2014). Cold acclimation is in part activated through C‐repeat binding factors (CBFs) that bind to the C‐repeat cis‐element present in some COLD RESPONSIVE (COR) genes to induce their expression (Lee et al., 2005). Many of the COR proteins play critical roles in triggering cold acclimation and subsequent freezing tolerance (Liu et al., 2019; Zhao, Zhang, et al., 2016). Besides the CBF‐mediated cold‐response pathway, some studies have suggested that CBF‐independent pathways also play various roles in controlling cold acclimation and freezing tolerance in Arabidopsis (Chinnusamy et al., 2006; Liu et al., 2019; Park et al., 2015). Furthermore, some responses to cold acclimation are gated by circadian rhythm (Fowler et al., 2005) as expression of some diurnal‐ and circadian‐regulated genes were significantly altered during cold acclimation in Arabidopsis (Arabidopsis thaliana Heyn) (Bieniawska et al., 2008; Espinoza et al., 2010). Photosynthetic cold acclimation is also regarded as an important factor of cold hardiness and freezing tolerance; the process involves the interactions of light, low temperature, and chloroplast redox responses (Gray et al., 1997; Wanner & Junttila, 1999). Hence, chlorophyll fluorescence is commonly applied as an assessment of the light‐processing efficiency of plant photosystems and the level of plant stress (Cessna et al., 2010).

In this study, we used GWAS to locate markers or genes associated with target traits for freezing tolerance and photosynthetic responses. GWAS is normally applied to evaluate and locate genomic regions that affect polygenic traits (Visscher et al., 2017). Several canola populations were used previously for GWAS studies of cold responses in rapeseed (Bus et al., 2011; Hatzig et al., 2015; Rahman et al., 2016; Schiessl et al., 2017; Wrucke et al., 2019; Xu et al., 2016; Zhou et al., 2017). Most of these diversity panels were a combination of winter, semi‐winter and spring canola plants. Recently, we used a diversity panel of 399 accessions to conduct GWAS (Chao et al., 2021); this diversity panel was assembled from a rapeseed breeding population at Kansas State University that represents primarily winter biotypes from mostly North America and Europe. The current study evaluated freezing tolerance from another diversity panel of 222 accessions consisting of winter rapeseed obtained from the Genebank of the Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) in Gatersleben, Germany. The accessions making up this new diversity panel were assembled from wild winter rapeseed collected primarily from Europe including Belgium, Canada, Denmark, Germany, Ireland, Lithuania, Netherlands, Poland, Sweden, Russia, and United Kingdom (Table S1, spreadsheet “IPK Accessions”).

MATERIALS AND METHODS

Plant material

Seeds of 222 winter rapeseed accessions were obtained from IPK in Gatersleben, Germany. Common cultivar names are listed in Table S1 (Spreadsheet “IPK Accessions”). The growth conditions were the same as cold‐acclimation treatment of winter canola described previously (Chao et al., 2021). Briefly, seeds were planted into potting soil in Deepot Cells and put in support trays. Plants were grown for 4 weeks under greenhouse conditions with daily watering and weekly fertilization.

Phenotypic data collection

Freezing experiments were performed using custom‐built Percival walk‐in chambers (Chao et al., 2021). Each chamber can hold up to 30 flats, and each flat has 15 cells. For each experiment, plants were germinated in the greenhouse in a complete randomized block design with three technical replications per accession; randomization was done using an Excel sub‐routine. The flats containing randomized potted plants were also randomly placed inside the chambers using the Excel sub‐routine. Four‐week‐old plants were acclimated for 2 months at 5°C under 12‐h photoperiods, and subsequently frozen in walk‐in freezing chambers. The growth conditions in walk‐in chambers, freezing treatments (−15°C for 4 h), and phenotypic data collection (chlorophyll fluorescence in Fv/Fm and Fv/Fo values at 3 and 7 days and damage scores at 2 and 4 weeks) were the same as cold‐acclimation treatment of winter canola described previously (Chao et al., 2021). The phenotypic data were collected from six independent experiments (18 plants total per accession) for GWAS analysis (see below).

Phenotypic data analysis

Phenotypic data of visual damage scale (2 and 4 weeks) and chlorophyll fluorescence (Fv/Fo and Fv/Fm at 3 and 7 days) were analyzed with a linear mixed model using SAS (SAS Institute, 2015). The statistical model was ; where y is the observed phenotypic value, μ is overall mean, E is the effect of the i experiment, r(E) is the effect of the j replication within the i experiment, C is the effect of the k growth chamber, and G is the effect of the l genotype. Experiment, chamber, and replication were taken into account as random factors and genotype was considered as fixed factor. Best linear unbiased estimators (BLUEs) were assessed for each trait and used for further genome wide association analysis. Variance components for genotype (V ), genotype by experiment interaction (V ), genotype by chamber interaction (V ), and random error (V ) were assessed by considering all factors as random. Broad sense heritability (H ) was assessed as , where i is the number of experiment, j is number of replication, and k is number of chamber.

SNP marker discovery and genotyping

The 222 accessions were genotyped using genotyping‐by‐sequencing (GBS). DNA was isolated with the Wizard Genomic DNA Purification Kit (A1125; Promega) per the manufacturer's instructions and quantified with a Quant‐iT PicoGreen dsDNA assay kit (P7589; Thermo Fisher Scientific). The GBS library was constructed using the same protocol described in Horvath et al. (2020). The GBS library was sequenced on an Illumina (San Diego, CA) HiSeq 4000 to generate single‐end, 100‐bp reads at the Genomic Sequencing and Analysis Facility at the University of Texas Southwestern Medical Center at Dallas, TX. The SNP discovery and genotype calling was performed using the TASSEL‐GBS pipeline (Glaubitz et al., 2014) with the B. napus v4.1 (Chalhoub et al., 2014) as reference genome.

Population structure and linkage disequilibrium

A principle component (PC) analysis was performed using TASSEL v5.2 (Bradbury et al., 2007). Based on a screen plot (Figure S1), the first four PCs were selected for a model‐based cluster analysis using R function mclust (Scrucca et al., 2016). Linkage disequilibrium (LD) was estimated as r 2 between pairs of SNP markers using TASSEL v5.2 (Bradbury et al., 2007) with a sliding window of 50. Figure S2 shows the distribution of r 2 measure of pairwise linkage disequilibrium depending on the physical distance between SNPs.

Genome wide association analysis

Genome wide association analysis was conducted for each trait using TASSEL v5.2 (Bradbury et al., 2007). A linear mixed model including population structure and kinship was used. The first four PCs were used as covariates and included in the mixed model as fixed factors. Kinship matrix was calculated with TASSEL and included in the mixed model as a random factor. Significance of marker‐trait association was defined by false discovery rate (FDR) as a q value < .1 (Table S2, spreadsheet “List of all markers”).

Selection of potential candidate genes

Local LD decay rates for 500 Kb with a 10 Kb sliding window was determined based on the method done previously (Horvath et al., 2020) to produce an LD decay table established for this population (under column “Local LD decay interval” in Table S2, spreadsheet “List of all markers”). The local LD decay rate for regions centered on each significant marker was determined, and genes linked to those loci were located using the Canola Genome Browser (http://www.genoscope.cns.fr/brassicanapus/cgi-bin/gbrowse/colza/). Candidate genes inside the associated loci were determined based on gene expression data from an earlier RNAseq study (Accession # PRJNA560411, Horvath et al., 2020) and/or known or expected gene function. Cold‐regulated AGI orthologues are identified based on “Electronic Fluorescent Pictograph” (Winter et al., 2007; https://www.arabidopsis.org/).

RESULTS

Phenotyping rapeseed accessions for freezing‐induced stress

Six phenotypic data sets generated from a winter rapeseed diversity panel (222 accessions) were used to conduct GWAS (Figure 1 and Table S3). At 3 and 7 days following the freezing treatment, chlorophyll fluorescence values (Fv/Fm and Fv/Fo) were measured, and at 2 and 4 weeks, visual damage scores were determined. Correlation coefficients (r) among the six phenotypic data sets are provided in Table 1 (also Table S3). The range of r values for damage score versus chlorophyll fluorescence is .52–.67, and the range of r values for chlorophyll fluorescence (3 days) versus chlorophyll fluorescence (7 days) is .71–.74. Three pairs of phenotypic data, damage score (2 weeks) versus damage score (4 weeks), Fv/Fo (3 days) versus Fv/Fm (3 days), and Fv/Fo (7 days) versus Fv/Fm (7 days), were closely associated, and their r values are .96, .96, and .97, respectively (Table 1). Thus, damage scores are only moderately associated with chlorophyll fluorescence scores.

FIGURE 1

Alteration in the winter rapeseed diversity panel (222 accessions) with respect to freezing survival after acclimation for 2 months at 5°C. (a and b) The value of each accession is the average damage score of six independent experiments 2 and 4 weeks following freezing treatment. (c and d) The value of each accession is the average chlorophyll fluorescence (Fv/Fo and Fv/Fm, respectively) of six independent experiments 3 day after freezing treatment. (e and f) The value of each accession is the average chlorophyll fluorescence (Fv/Fo and Fv/Fm) of six independent experiments 7 day after freezing treatment. Vertical bars refer to standard error

TABLE 1

Correlation coefficients (r) among the six phenotypic data sets

	Damage score, 2 weeks	Damage score, 4 weeks	Fv/Fo, 3 days	Fv/Fm, 3 days	Fv/Fo, 7 days	Fv/Fm, 7 days
Damage score, 2 weeks	1.00
Damage score, 4 weeks	.96	1.00
Fv/Fo, 3 days	.60	.58	1.00
Fv/Fm, 3 days	.52	.49	.96	1.00
Fv/Fo, 7 days	.67	.66	.73	.71	1.00
Fv/Fm, 7 days	.62	.61	.74	.74	.97	1.00

Heritability, population structure, and linkage disequilibrium analysis

The ranges of average values, means and medians of the estimated BLUEs for 2‐ and 4‐week freezing damage score and 3‐ and 7‐day fluorometer for 222 accessions are listed in Table 2 and Table S1 (spreadsheet “Estimated BLUEs”). The heritability for 2‐week damage score, 4‐week damage score, 3‐day Fv/Fm, 3‐day Fv/Fo, 7‐day Fv/Fm, and 7‐day Fv/Fo were .39, .42, .71, .67, .74, and .74, respectively.

TABLE 2

Ranges of average damage scores and fluorescence readings, as well as means and medians for 222 accessions

	Range of avg value	Mean	Median
Damage score, 2 weeks	.23–2.35	1.41	1.40
Damage score, 4 weeks	.17–2.76	1.46	1.51
Fv/Fm, 3 days	.23–.79	.63	.64
Fv/Fo, 3 days	.72–3.87	2.81	2.85
Fv/Fm, 7 days	.12–.79	.57	.58
Fv/Fo, 7 days	.59–4.18	2.94	2.96

In total, 69,554 SNP markers were identified for the 222 accessions with missing values less than 50%. Model‐based cluster analysis suggested that there were eight subgroups (Figure 2 and Table S1, spreadsheet “Class and Country of origin”); there was no clear relationship between geographic distribution and the subgroups. The reason for this result may be due to that the number of accessions from each region is swayed heavily (50%) toward one location (Sweden—113 accessions), with Ireland only having one accession (.45%) and Russia with 15 (6.8%) in the IPK collection. However, the large collection of accessions from Sweden do not all cluster together, indicating that the Swedish collection may have originated from different regions of Sweden, which can be very diverse in a south to north gradient.

FIGURE 2

Population structure of rapeseed accessions. Scatter plots of PC1 and PC2 derived from a principal component analysis

Pairwise LD was calculated using all 69,554 SNP markers across all 222 accessions. The average r value of all the SNP pairs was .17. The mean r value was decayed to .26 between markers with distance of <.1 Mb and .09 with distance <1.0 Mb (Figure S2).

Markers associated with freezing tolerance based on GWAS

GWAS analyses were performed using the estimated BLUEs from six phenotypic data sets for 3‐day Fv/Fm, 3‐day Fv/Fo, 7‐day Fv/Fm, 7‐day Fv/Fo, 2‐week damage score, and 4‐week damage score. A total of 11 loci on six chromosomes were significantly associated with 7‐day Fv/Fm (Table S2, spreadsheet “List of all markers”). The region on chromosome A01 contained four markers and the most significant marker SA01_22176724 explained 12.4% of total variation. It should be noted that these four markers are tightly linked; three markers (SA01_22143868, SA01_22143872, and SA01_22143879) shared the same gene models inside 50 Kb to both sides of the marker, whereas marker SA01_22176724 partially shared gene models with the above three markers inside 50 Kb to both sides of the marker. The genomic region on chromosome A02 contained only one marker and explained 9% of total variation. The region on chromosome A03 contained four markers and the most significant marker SA03_15761546 explained 10% of total variation. Here, marker SA03_15761546 and marker SA03_15761556 were tightly linked, and they share the same gene models inside 50 Kb to both sides of the marker. The region on chromosome A07 contained four markers and the most significant marker SA07_23526647 explained 11.4% of total variation. These four markers were also tightly linked; marker SA07_23415428 partially shared gene models with marker SA07_23433209, whereas marker SA07_23433209 shared the same gene models with marker SA07_23433215 but partially shared gene models with marker SA07_23526647 inside 50 Kb to both sides of the marker. The genomic region on chromosome A08 contained three markers and the most significant marker SA08_14970416 explained 10.2% of total variation. These three markers (SA08_14970416, SA08_14970445, and SA08_14970446) were tightly linked and shared the same gene models inside 50 Kb to both sides of the marker. The genomic region on chromosome C01 contained only one marker and explained 9.8% of total variation. The same six genomic regions were identified for 7‐day Fv/Fo with 7‐day Fv/Fm as noted above.

A total of three loci on two chromosomes were significantly associated with 3‐day Fv/Fm at FDR of q‐value smaller than .05 (Table S2, spreadsheet “List of all markers”). The first region was located on chromosome A06 containing two markers and the second region on chromosome C04 contained two markers. The most significant markers at each locus individually explained about 11% of total variation.

No significant markers were associated with visual damage scores at either time point nor for 3‐day Fv/Fo.

Underlying candidate genes associated with freezing tolerance

A total of 32 significant markers representing 14 significant loci (Table 3) were associated with 3‐day Fv/Fm, 7‐day Fv/Fm, and 7‐day Fv/Fo (Table S2, spreadsheet “List of all markers” and “List of markers‐no duplicates”). Representative Manhattan plots are shown in Figure 3. The accurate borders of the associated loci according to the local LD decay rate in the area encompassing each marker was also established (Table S4, spreadsheet “List of surrounding gene models” and under the column “ROI based on local LD‐decay rate”). The LD decay interval ranged from 5 to 129 Kb to either side of the markers, and 13 loci had recognized gene models (1 to 24 gene models) inside the LD decay intervals (highlighted in pink under column “Locus”). In 11 of the 13 loci, the marker fell inside the transcribed segment of a gene (highlighted in red under column “Locus”). The local LD decay interval for marker SC04_16855437 cannot be determined using the method mentioned above; however, the marker fell close to a non‐annotated expressed sequence tag (CDMBras.unigeneT00161455001). Genes inside the LD decay interval might represent candidates associated with freezing tolerance especially when there are only a small number of genes. Genes identified to be differentially expressed in winter rapeseed after a 1‐month cold acclimation (Horvath et al., 2020) are also specified (up‐regulated genes are high‐lighted in green and down‐regulated genes are highlighted in orange under column “Canola Gene Model”).

TABLE 3

Fourteen significant loci were identified based on GWAS

Trait	Marker	Chromosome	Position	p value	q value	Marker R 2
7‐day Fv/Fm	SA01_22143868	SA01	22143868	.000022359	.089036723	.08602
7‐day Fv/Fm	SA01_22176724	SA01	22176724	4.3784E−07	.02950954	.12362
7‐day Fv/Fm	SA02_8739425	SA02	8739425	.000016457	.079226349	.09014
7‐day Fv/Fm	SA03_15761546	SA03	15761546	5.0346E−06	.037702441	.1
7‐day Fv/Fm	SA03_16201884	SA03	16201884	.00001038	.05829927	.09556
7‐day Fv/Fm	SA03_16566033	SA03	16566033	.000013196	.068414154	.09052
7‐day Fv/Fm	SA07_23415428	SA07	23415428	2.1637E−06	.036457263	.11075
7‐day Fv/Fm	SA07_23433209	SA07	23433209	8.6714E−06	.053130456	.10328
7‐day Fv/Fm	SA07_23526647	SA07	23526647	1.5619E−06	.036457263	.11357
7‐day Fv/Fm	SA08_14970416	SA08	14970416	3.9429E−06	.037702441	.10243
7‐day Fv/Fm	SC01_9232921	SC01	9232921	5.8459E−06	.036457263	.09804
3‐day Fv/Fm	SA06_4056139	SA06	4056139	1.4549E−06	.024514338	.11016
3‐day Fv/Fm	SC04_13335664	SC04	13335664	1.3523E−06	.024514338	.11377
3‐day Fv/Fm	SC04_16855437	SC04	16855437	1.3673E−06	.024514338	.1128

FIGURE 3

Manhattan plots. Manhattan plot of marker–trait association study for the six phenotypic data sets of 222 rapeseed accessions. The X axis is the chromosomes, and the Y axis is the −log10(FDR)

There are 19 gene models inside the 50 Kb around both sides of the significant marker SA01_22143868 on chromosome A01. This marker had significant association with Fv/Fm readings at 7 days after freezing treatment. The LD decay interval is short (<5 Kb to either side of the marker) and only two of these gene models fall inside the LD decay interval. However, the marker SA01_22143868 falls into the exon of a gene encoding SENSITIVE TO FREEZING 2 (SFR2). This protein plays a role in membrane lipid remodeling and freezing tolerance (Barnes et al., 2016; Thorlby et al., 2004). Another gene known to encode RAB ESCORT PROTEIN (REP) in Arabidopsis is within the LD decay region. Other genes of interest including RSA1 INTERACTING TRANSCRIPTION FACTOR 1 (RITF1), ELO HOMOLOG 1&2, and ACYL‐ACYL CARRIER PROTEIN DESATURASE 1‐6 (AAD1‐6), fall just outside the LD decay interval. These genes are involved in detoxification of reactive oxygen species, fatty acid elongation, and fatty acid metabolic process, respectively.

A nearby significant marker SA01_22176724 on chromosome A01 also has a small LD decay region (<5 Kb), and only two gene models fall within the LD decay interval. This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. There are 21 gene models inside the 50 Kb around both sides of this marker; however, 14 gene models overlap between this marker and the marker SA01_22143868 mentioned above (highlighted in blue and light green under column “Protein match” in Table S4, spreadsheet “List of surrounding gene models”). The marker SA01_22176724 falls into the exon of the gene encoding HYDROXYPROLINE O‐GALACTOSYLTRANSFERASE (GALT3). Another gene encoding EMBRYO DEFECTIVE (EMB) in Arabidopsis is within the LD decay region. Several genes of interest (mentioned also for marker SA01_22143868) fall just outside the LD decay interval. There is also an ABA‐responsive (ABR) gene located around 46 Kb to the marker.

There are 16 gene models inside the 50 Kb around both sides of the significant marker SA02_8739425 on chromosome A02. This marker had significant association with Fv/Fm readings at 7 days after freezing treatment. Four of these gene models fall within the LD decay interval. However, only one (BnaA02g15190D) of them encodes an orthologue of Arabidopsis MLP‐LIKE PROTEIN 31, and this gene was down regulated following a 1‐month cold acclimation treatment; the other three have no homology with any Arabidopsis genes. Although outside the area of likely LD decay, a photosynthetic gene CHLORORESPIRATORY REDUCTION 23 (CRR23) is physically linked to this marker.

There are 24 gene models inside the 50 Kb around both sides of the significant marker SA03_15761546 on chromosome A03, and all of these gene models fall within the LD decay interval. This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. The marker SA03_15761546 falls into the intron of a gene encoding NICOTINATE TRANSPORTER (NIAP). Other genes of interest include DNA J PROTEIN C66 (DJC66) and polyketide cyclase/dehydrase (BnaA03g32620D). DJC66 encodes DNA J PROTEIN C66 (BnaA03g32750D). This protein is postulated to play a role in cell‐cycle regulation and plant growth (Eloy et al., 2015; Schulz et al., 2014). Polyketide cyclase/dehydrase appears to play a role in seed dormancy and plant defense (Chibani et al., 2006; Jones et al., 2006), and the transcript of this gene is cold‐regulated in Arabidopsis (Winter et al., 2007).

There are 19 gene models inside the 50 Kb around both sides of the significant marker SA03_16201884 on chromosome A03, and 17 of these gene models fall within the LD decay interval. This marker had significant association with Fv/Fm readings at 7 days after freezing treatment. The marker SA03_16201884 falls into the exon of a gene encoding TERPENE SYNTHASE 20 (TPS20); TPS enzymes are the gatekeepers in generating diversity of bioactive terpenoids across the plant kingdom (Karunanithi & Zerbe, 2019). Other genes of interest within the LD decay interval include NINE‐CIS‐EPOXYCAROTENOID DIOXYGENASE 3 (NCED3), CYTOCHROME P450 (CYP), and LYSM RLK1‐INTERACTING KINASE 1 (LIK1). NCED3 encodes a vital enzyme in the biosynthesis of abscisic acid. CYPs play a crucial role in plant development and defense, and one (BnaA03g33530D) of the CYP gene models is up regulated following a 1‐month cold acclimation treatment in canola. LIK1 encodes LRR‐RLK protein and appears to be involved in jasmonic acid and ethylene signaling pathways (Le et al., 2014) and is down‐regulated after a one‐month cold acclimation in canola (Table S4, spreadsheet “List of surrounding gene models” and under the column “Canola Gene Model”).

There are 19 gene models inside the 50 Kb around both sides of the significant marker SA03_16566033 on chromosome A03, and eight of these gene models fall within the LD decay interval. This marker had significant association with Fv/Fm readings at 7 days after freezing treatment. The marker SA03_16566033 falls into the intron of a gene encoding 60S ribosomal protein L37‐2. There is a gene encoding cytochrome b6‐f complex subunit 5 (PETG) that is a little down‐regulated by cold in Arabidopsis (Winter et al., 2007) and is nearby to the gene model in which the marker falls. PETG is required for the stability of the cytochrome b6‐f complex (Cyt b6/f), which facilitates electron transfer between photosystem II and photosystem I (Schwenkert et al., 2007).

There are 30 gene models inside the 50 Kb around both sides of the significant marker SA07_23415428 on chromosome A07, and only six of these gene models fall inside the 13‐Kb LD decay interval. This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. The marker SA07_23415428 falls into the intron of a gene encoding UBIQUITIN CONJUGATING ENZYME 13A/35, which may not affect the function of this gene. Other genes inside the LD decay interval encode three hypothetical proteins, a curculin‐like (mannose‐binding) lectin family protein, and a ubiquitin‐specific protease family C19‐related protein (SHOU4L). SHOU4L is cold induced (12–24 h) in Arabidopsis (Winter et al., 2007). Although outside the area of LD decay, there are some genes that are differentially expressed in canola after a 1‐month cold acclimation. These include VITAMIN E DEFICIENT6 (BnaA07g34270), EMBRYO DEFECTIVE 1473 (BnaA07g34280D), GAMMA‐GLUTAMYL HYDROLASE 3 (BnaA07g34300D), and a gene encoding D‐mannose binding lectin protein (BnaA07g34420D).

A nearby locus surrounds the significant marker SA07_23433209 on chromosome A07. There are 27 gene models inside the 50 Kb around both sides of this marker; however, 24 gene models overlap between this marker and the marker SA07_23415428 mentioned above (highlighted in blue and light green under column “Protein match” in Table S4, spreadsheet “List of surrounding gene models”). This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. There are six gene models within the 26‐Kb LD decay interval, and the marker SA07_23433209 falls into the exon of a gene encoding PSEUDOURIDINE SYNTHASE 3 (PUS3), which is involved in transforming uridine into pseudouridine (Rintala‐Dempsey & Kothe, 2017). Other genes within the LD decay interval encode a hypothetical protein, a VACUOLAR ATP SYNTHASE SUBUNIT A (VHA‐A), a VACUOLAR‐PYROPHOSPHATASE LIKE PROTEIN 1 (VPL1), and a mitochondrial transcription termination factor family protein (MTERF16). VPL1 plays a role in proton transport in the plant vacuolar membrane; overexpression OVP1 enhances cold tolerance in rice (Zhang et al., 2011) and salt and drought tolerance in Arabidopsis (Gaxiola et al., 2001).

Another nearby locus surrounds the significant marker SA07_23526647 on chromosome A07. There are 23 gene models inside the 50 Kb around both sides of this marker, and 3 gene models overlap with genes of neighboring marker SA07_23433209 (highlighted in blue and light green under column “Protein match” in Table S4, spreadsheet “List of surrounding gene models”). This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. Five gene models are localized within the 26.5‐Kb LD decay interval, and this marker falls into the exon of a gene encoding NUCLEOLAR COMPLEX ASSOCIATED 3 (NOC3), which interacts with a floral regulator REBELOTE to control specific developmental processes, for example, floral meristem termination (de Bossoreille et al., 2018). NOC3 is also up regulated following a 1‐month cold acclimation treatment (Table S4, spreadsheet “List of surrounding gene models” and under the column “Canola Gene Model”). Other genes within the LD decay interval includes BOI‐RELATED GENE 2 (BRG2), UBIQUITIN CARBOXYL‐TERMINAL HYDROLASE (UCH), SMALL AUXIN UPREGULATED 40 (SAUR40), and MYB DOMAIN PROTEIN 63 (MYB63). MYB63 is a transcription factor that is involved in the lignin biosynthetic pathway during secondary wall formation in Arabidopsis (Zhou et al., 2009).

There are 19 gene models inside the 50 Kb around both sides of the significant marker SA08_14970416 on chromosome A08, and 9 of these gene models fall within the 43.5‐Kb LD decay interval. This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. The marker SA08_14970416 falls into the intron of a gene encoding LYSINE‐HISTIDINE‐LIKE TRANSPORTER 6 (LHT6). Other genes within the LD decay interval includes CESTA (CES), CONSTANS‐LIKE 16 (COL16), ALUMINIUM ACTIVATED MALATE TRANSPORTER 4 (ALMT4), ENHANCED ETHYLENE RESPONSE 1 (EER1), CTL‐LIKE PROTEIN DDB, PHOTOSYNTHESIS‐AFFECTED MUTANT 71 LIKE 5 (PML5), PHYTOCHROME AND FLOWERING TIME 1 (PFT1), and ETHYLENE RESPONSE DNA BINDING FACTOR 1 (EDF1). CES is a positive regulator of brassinosteroid biosynthesis (Poppenberger et al., 2011). COL16 plays a role in the regulation of chlorophyll biosynthesis and accumulation; overexpression of COL16 enhances chlorophyll accumulation in petunia corollas (Ohmiya et al., 2019). PHOTOSYNTHESIS AFFECTED MUTANT71 is an integral thylakoid membrane protein that functions in manganese uptake across the chloroplast envelope and the thylakoid membrane (Schneider et al., 2016). PFT1 acts downstream of phyB to control the expression of FLOWERING LOCUS T, and it stimulates flowering in response to suboptimal light conditions (Cerdán & Chory, 2003). EDFs serve as downstream genes in the ethylene response to regulate flower senescence and abscission (Chen et al., 2015).

There are 9 gene models inside the 50 Kb around both sides of the significant marker SC01_9232921 on chromosome C01. This marker had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. Only one gene model (BnaC01g14040D) falls within the 24‐Kb LD decay interval, and the marker SC01_9232921 falls into the exon of this gene, which encodes an INAPERTURATE POLLEN1 (INP1) with a role in the formation of pollen surface apertures (Dobritsa & Coerper, 2012). Although outside the area of LD decay, a HSP70‐INTERACTING PROTEIN 1 (HIP1) is associated with this marker; HIP1 protein binds to Hsp70, stabilizing it and promoting binding of target polypeptides (Webb et al., 2001).

Two markers were identified to have significant association with Fv/Fm readings at 3 days after freezing treatment; however, no markers were found to have significant association with Fv/Fo readings at 3 days after freezing treatment. There are 23 gene models inside the 50 Kb around both sides of the significant marker SA06_4056139 on chromosome A06, and five of these gene models fall within the 10.2‐Kb LD decay interval. The marker SA06_4056139 falls into the exon of a gene, which encodes a RNA‐binding (RRM/RBD/RNP motifs) family protein. Other genes within the LD decay interval encode a RING/U‐box superfamily protein, a RNA‐binding family protein, an ACTIN‐RELATED PROTEIN 3 (ARP3), and a succinate‐‐CoA ligase subunit alpha. Although outside the area of LD decay, three genes of interest (BnaA06g07590D, BnaA06g07620D, and BnaA06g07630D) encode different class of CYTOCHROME P450 (i.e., CYP86C3, CYP71B7, and CYP71B23, respectively). There are seven gene models inside the 50 Kb around both sides of the significant marker SC04_13335664 on chromosome C04, and all fall within the LD decay interval. These genes encode three unnamed proteins, a hypothetical protein, a C2H2‐type zinc finger family protein, a NAC DOMAIN CONTAINING PROTEIN 95 (ANAC095), and a TOX high mobility group box protein.

DISCUSSION

Relationship among six phenotypic data sets

We used a rapeseed diversity panel assembled from European winter varieties to conduct GWAS. Differing from previous reports (Chao et al., 2021), we applied six phenotypic data to conduct GWAS, including an additional 4‐week damage score and fluorescence readings (Fv/Fm at 3 and 7 days after a freezing treatment). Though these rapeseed accessions were collected from Europe (Table S1, spreadsheet “IPK Accessions”), the phenotypic data were similar to the rapeseed diversity panel used previously (Chao et al., 2021); for example, the range of the average varietal damage scores (2 weeks) were .23–2.35 for this diversity panel as opposed to .11–2.44 for the previous one, and the range of chlorophyll fluorescence readings were .72–3.87 (Fv/Fo, 3 days) and .59–4.18 (Fv/Fo, 7 days) for this diversity panel as opposed to .87–4.29 (Fv/Fo, 3 days) and .87–4.39 (Fv/Fo, 7 days) for the previous one. These results indicated that the scores for freezing damage were similar between these two rapeseed populations.

A total of 14 significant loci were identified in association with 3‐day Fv/Fm, 7‐day Fv/Fm, and 7‐day Fv/Fo based on GWAS after excluding duplicates; however, no significant markers were associated with 3‐day Fv/Fo and 2‐ and 4‐week freezing damage score (Table 3 and Table S2, spreadsheet “List of all markers”). The identification of significant markers by GWAS appears to be positively correlated with heritability as the heritability values were .71, .74, and .74 for 3‐day Fv/Fm, 7‐day Fv/Fm, and 7‐day Fv/Fo and were .39, .42, .67 for 2‐week damage score, 4‐week damage score, and 3‐day Fv/Fo, respectively. Although there were no significant markers associated with freezing damage score, it should be noted that there were reasonable correlations (r range = .49–.67, Table 1) observed between visual damage and chlorophyll fluorescence ratings, indicating that some loci may contribute to both phenotypes (i.e., freezing damage/tolerance and photosynthetic efficiency).

Candidate genes associated with cold acclimation and freezing tolerance

When comparing the data of this study with our previous GWAS work (Chao et al., 2021), no overlapping significant loci were identified between these two studies. On chromosome A07, significant markers were identified in both studies; however, the distance between a cluster of three markers (23415428, 23433209, 23526647) on A07 from this study and the one (10641317) from our previous study is too large (over 12.7 Mb) to be relevant. The reasons for these disparate results may be (1) different diversity panels were used to conduct GWAS, (2) different approaches were used to estimate marker significance, and (3) different methods were used to determine statistical significance (FDR as a q value < .1 in this study versus p values < .005 in previous study).

All significant markers in this study were associated with chlorophyll fluorescence (3‐day Fv/Fm, 7‐day Fv/Fm, and 7‐day Fv/Fo). Since photosynthetic efficiency correlated with freezing damage scores, this phenotype acts as an acceptable indicator of freezing damage. Chlorophyll fluorescence has been used as an indicator of cellular damage, predominantly damage to photosystems, following abiotic stress in other plant systems (Cessna et al., 2010), and indeed several photosynthetic and light signaling genes were identified within the associated loci. Other genes associated with lipid and fatty acids, photosynthesis, flowering, ubiquitination, and cytochrome were also identified. Our confidence in the associations of these markers with phenotypes is supported by many candidate genes with known or expected roles in cold acclimation and freezing tolerance. The prospective roles of these candidate genes are discussed below (Table 4 and Table S4, spreadsheet “Candidate genes”).

TABLE 4

Candidate genes involved in cold acclimation and freezing tolerance

Marker	Chromosome	Position	ROI based on local LD‐decay rate	Canola gene model	AGI orthologue	Other names
SA01_22143868	SA01	22143868	22139102..22148634	BnaA01g32790D	AT3G06510.2	SENSITIVE TO FREEZING 2 (SFR2)
SA01_22143868	SA01	22143868	22139102..22148634	BnaA01g32810D	AT3G06470.1	ELONGATION OF FATTY ACIDS PROTEIN 2 (ELO2)
SA01_22143868	SA01	22143868	22139102..22148634	BnaA01g32820D	AT3G06460.1	ELO HOMOLOG 1 (ELO1)
SA02_8739425	SA02	8739425	8725309..8753541	BnaA02g15140D	AT1G70760.1	CHLORORESPIRATORY REDUCTION 23 (CRR23)
SA02_8739425	SA02	8739425	8725309..8753541	BnaA02g15190D	AT1G70840.1	MLP‐LIKE PROTEIN 31 (MLP31)
SA02_8739425	SA02	8739425	8725309..8753541	BnaA02g15220D	AT1G70840.1	MLP‐LIKE PROTEIN 31 (MLP31)
SA02_8739425	SA02	8739425	8725309..8753541	BnaA02g15230D	AT1G70830.1	MLP‐LIKE PROTEIN 28 (MLP28)
SA02_8739425	SA02	8739425	8725309..8753541	BnaA02g15240D	AT1G70830.1	MLP‐LIKE PROTEIN 28 (MLP28)
SA03_15761546	SA03	15761546	15657085..15866007	BnaA03g32620D	AT3G13062.2	Polyketide cyclase/dehydrase
SA03_15761546	SA03	15761546	15657085..15866007	BnaA03g32770D	AT3G13330.1	PROTEASOME ACTIVATING PROTEIN 200 (PA200)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33380D	AT3G14415.2	GLYCOLATE OXIDASE 2 (GOX2)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33480D	AT3G14630.1	CYTOCHROME P450, 72A9 (CYP72A9)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33490D		CYTOCHROME P450, 72A15 (CYP72A15)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33500D	AT3G14640.1	CYTOCHROME P450, 72A10 (CYP72A10)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33510D		CYTOCHROME P450, 72A15 (CYP72A15)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33520D		CYTOCHROME P450, 72A15 (CYP72A15)
SA03_16201884	SA03	16201884	16161988..16241780	BnaA03g33530D	AT3G14640.1	CYTOCHROME P450, 72A10 (CYP72A10)
SA03_16566033	SA03	16566033	16551949..16580117	BnaA03g34060D	ATCG00600.1	Cytochrome b6‐f complex subunit 5 (PETG)
SA03_16566033	SA03	16566033	16551949..16580117	BnaA03g34070D	AT2G16005.1	INTERACTOR OF SYNAPTOTAGMIN1 (ROSY1)
SA07_23415428	SA07	23415428	23406287..23424569	BnaA07g34450D	AT1G78870.1	UBIQUITIN CONJUGATING ENZYME 13A (UBC13A)
SA07_23526647	SA07	23526647	23515656..23537638	BnaA07g34660D	AT1G79110.1	BOI‐RELATED GENE 2 (BRG2)
SA07_23526647	SA07	23526647	23515656..23537638	BnaA07g34670D	AT1G79120.1	Ubiquitin carboxyl‐terminal hydrolase (UCH)
SA08_14970416	SA08	14970416	14945662..14995170	BnaA08g19480D	AT1G25540.1	PHYTOCHROME AND FLOWERING TIME 1 (PFT1)
SA08_14970416	SA08	14970416	14945662..14995170	BnaA08g19490D	AT1G25560.1	ETHYLENE RESPONSE DNA BINDING FACTOR 1 (EDF1)

Note: Different markers are in gray.

Lipid and fatty acid associated genes

SENSITIVE TO FREEZING2 (SFR2): SFR2 (BnaA01g32790D) is associated with significant marker SA01_22143868 and SA01_22176724 on chromosome A01. SFR2 encodes a protein homologous to glycosyl hydrolases (β‐glycosidases), and Arabidopsis SFR2 exhibited a hydrolytic activity against β‐d‐glucosides. SFR2 is essential for freezing tolerance because mutations to this gene causes freezing sensitivity in Arabidopsis (Thorlby et al., 2004). The sfr2 mutation exhibited severe leaf damage in the whole‐plant freezing assays from the Columbia ecotype of Arabidopsis (Warren et al., 1996). Some β‐glucosidases alter the properties of the cell wall (Gerardi et al., 2001; Li et al., 2001), a structural layer that is often associated with freezing tolerance. Barnes et al. (2016) showed that SFR2 was activated posttranslationally by pH (acidification) and Mg2+and proposed that SFR2 activation is a molecular freezing sensing mechanism causing chloroplast membrane lipid remodeling and subsequent freezing tolerance. Although SFR2 is correlated with posttranslational activation, this gene is not cold regulated in Arabidopsis (Thorlby et al., 2004) and is not differentially expressed in canola following a one‐month cold acclimation treatment. The marker SA01_22143868 falls into the exon of SFR2, and SFR2 is one of the two gene models falling within the LD decay interval; thus, assuming the association was not spurious, this SNP could be the sequence variations responsible for the phenotypic effect.

MAJOR LATEX PROTEIN (MLP): MLP‐LIKE PROTEIN 31 (BnaA02g15190D) is one of the four gene models within the LD decay interval and is about 6.7 Kb away from the marker SA02_8739425; in addition, this gene was down regulated following a one‐month cold acclimation treatment. Interestingly, three other gene models (BnaA02g15220D, BnaA02g15230D, and BnaA02g15240D) also encode MLP‐LIKE PROTEIN in the same locus but are outside the linkage decay range of the significant marker, SA02_8739425, and one (BnaA02g15230D) of the three was down regulated following a 1‐month cold acclimation treatment. MLP proteins are part of the Bet v 1 family, also known as the pathogenesis related 10 (PR10)‐like protein family (Radauer et al., 2008). Two AtMLPs (At2g01520 and At2g01530) could be induced by cis‐cinnamic acid and stimulate vegetative growth and prolong flowering (Guo et al., 2011). MLPs are also known to have binding affinity toward various hydrophobic compounds, such as long‐chain fatty acids and steroids, as well as involving in drought and salt tolerance and resistance against pathogens (Fujita & Inui, 2021). Considering various roles MLPs play in response to biotic and abiotic stresses and our cold acclimation treatments, the identification of MLPs in GWAS analysis upholds the assumption that they may be involved in weakening cold stress responses; however, the decisive role of MLPs in rapeseed remains to be clarified.

INTERACTOR OF SYNAPTOTAGMIN1 (ROSY1): ROSY1 is a lipid‐binding protein which plays a role in root gravitropism. Knockout mutants reveal reduced basipetal auxin transport, a faster root gravitropic response, and an upsurge in salt stress tolerance (Dalal et al., 2016). In this study, ROSY1 (BnaA03g34070D) is mapped within the LD decay interval and is about 3 Kb away from the marker SA03_16566033 on chromosome A03. Giving ROSY1's function in response to salt stress, the regulation of ROSY1 may directly or indirectly contribute to other stress signals such as cold acclimation enhanced freezing tolerance.

Several other lipid‐associated genes are linked to different significant markers; for example, a gene model (BnaA03g32620D) encoding polyketide cyclase/dehydrase is within the LD decay interval and is about 1 Kb away from the marker SA03_15761546. Many PYRABACTIN RESISTANCE 1‐LIKE (PYL) proteins contain the polyketide cyclase domain. They are part of the steroidogenic acute regulatory related lipid transfer (START) protein superfamily (Park et al., 2009) and play a role in ABA binding and drought tolerance (Dittrich et al., 2019; Santiago et al., 2009; Zhao, Chan, et al., 2016). Also, two ELONGATION OF FATTY ACIDS PROTEIN (ELO) genes (BnaA01g32810D and BnaA01g32820D) are associated with marker SA01_22143868 and SA01_22176724 on chromosome A01 but outside the LD decay interval. ELO2 is a component of the membrane‐bound fatty acid elongation system. It is involved in the production of 26‐carbon very long‐chain fatty acids and play a role in ceramide and inositol sphingolipid biosynthesis (Oh et al., 1997); ceramides are sphingolipids that exhibit unique physical properties in the lipid bilayer (Alonso & Goñi, 2018) and play critical roles in maintaining membrane function and integrity (Hussain et al., 2012). The observation of many lipid and fatty acid associated genes linked to the significant markers strengthens the hypothesis that lipid and fatty acid metabolism play significant roles in freezing tolerance.

Photosynthetic and light signaling genes

Cytochrome b6‐f complex subunit 5 (PetG): Photosynthesis is sensitive to changes in environmental stimuli, and several lines of evidence suggested that it is vital for plants to guard the photosynthetic apparatus from damage at low temperature (Sung et al., 2003). Cold acclimation modifies the photosynthetic machinery and allows plants to survive at freezing temperatures (Crosatti et al., 2013). One of the mechanisms may be that cold acclimation improves resistance to low‐temperature‐induced photoinhibition and promotes recovery after returning to above freezing temperatures. The recovery of PSII at low temperature in winter‐hardy ecotype of Colobanthus quitensis is associated with its capability to sustain PsaA, Cyt b6/f and D1 protein after photoinhibitory conditions (Bascuñán‐Godoy et al., 2012). Our study identified a PetG gene (BnaA03g34060D) within the linkage decay interval to marker SA03_16566033 (around 10 Kb) on chromosome A03. PetG encodes a bitopic plastid protein (i.e., subunit V) of the Cyt b6/f complex, which resides in the thylakoid membranes and facilitates electron transfer between photosystem II (PSII) and photosystem I (PSI) (Schwenkert et al., 2007). PetG protein plays a role in the assembly and stability of Cyt b6/f complex (Schwenkert et al., 2007). The linkage of PetG to SA03_16566033 suggests that this gene may be involved in cold acclimation process, and it supports a link between cold and photosynthetic processes to alleviate damage of photosynthetic apparatus during cold or freezing conditions.

CHLORORESPIRATORY REDUCTION 23 (CRR23) and GLYCOLATE OXIDASE 2 (GOX2): Two photosynthetic‐ and light associated genes (CRR23 and GOX2) falling outside the LD decay range were cold‐regulated in our earlier RNAseq study (Table S4, spreadsheet “List of surrounding gene models” and under the column “Canola Gene Model”). CRR23 (BnaA02g15140D) is less than 34 Kb from significant marker, SA02_8739425 on chromosome A02. CRR23 is essential for the stabilization and accumulation of the chloroplast NAD(P)H dehydrogenase complex in Arabidopsis and is involved in PSI cyclic electron transport (Shimizu et al., 2008). This gene was down regulated following a one‐month cold acclimation treatment. Since CRR23 is important for the integrity of chloroplast, the linkage of CRR23 to SA02_8739425 suggests that it may play a role in freezing tolerance. GOX2 (BnaA03g33380D) is about 55 Kb from significant marker, SA03_16201884, and it was down regulated following a 1‐month cold acclimation treatment. GOX2 protein catalyzes the conversion of glycolate into glyoxylate in peroxisomes during photorespiration with the generation of H2O2, which are known as reactive oxygen species (ROS) in plant defense responses against pathogens (Rojas et al., 2012). Since there is some evidence that photorespiration can prevent light‐induced damage to the molecules involved in photosynthesis (Eisenhut et al., 2017), the linkage of GOX2 to SA03_16201884 suggests that this gene may be involved in cold acclimation and freezing tolerance.

Flowering associated genes

PHYTOCHROME AND FLOWERING TIME 1 (PFT1) and ETHYLENE RESPONSE DNA BINDING FACTOR 1 (EDF1): PFT1 (BnaA08g19480D) is within the LD decay interval and is about 5 Kb away from the marker SA08_14970416 on chromosome A08. PFT1 is a nuclear protein, and it induces flowering during shade avoidance (far‐red enrichment/suboptimal light conditions) (Cerdán & Chory, 2003; Wollenberg et al., 2008). Far‐red enrichment is analogous to vernalization because both processes can stimulate rapid flowering (Bagnall, 1992; Martinez‐Zapater & Somerville, 1990). PFT1 also functions as a positive regulator of JA signaling that controls plant defense responses to leaf‐infecting fungal pathogens (Kidd et al., 2009). Moreover, it regulates plant response to abiotic stress as mutant Arabidopsis exhibits a heightened sensitivity to salt stress but a decreased sensitivity to drought stress (Elfving et al., 2011). Giving the significant link between the response to vernalization and far‐red enrichment and other stress responses mentioned above, PFT1 function could be considered as one of the many cellular mechanisms occurred during cold acclimation.

Another gene EDF1 (BnaA08g19490D) is adjacent to PFT1. EDF1 protein belongs to the EDF4/RELATED TO ABI3/VP1 (RAV) subfamily of ERF/AP2‐related proteins, and it functions as a repressor to activate the senescence and abscission in flower organs (Chen et al., 2015). If there are sequence variations in PFT1 and EDF1 genes, they could have a joint effect at this locus on cold acclimation and/or vernalization.

Ubiquitination and freezing tolerance

UBIQUITIN CONJUGATING ENZYME 13A (UBC13A), PROTEASOME ACTIVATING PROTEIN 200 (PA200), BOI‐RELATED GENE 2 (BRG2), and UBIQUITIN CARBOXYL‐TERMINAL HYDROLASE (UCH): UBC13A (BnaA07g34450D) are associated with significant marker SA07_23415428 and SA07_23433209 on chromosome A07, and the marker SA07_23415428 falls into the exon of UBC13A. Protein ubiquitination is responsible for protein turn over via the ubiquitin proteasome system (UPS) and is a fundamental mechanism of post‐translational modifications (Varshavsky, 2012). UBC13A belongs to the subfamily V of Arabidopsis E2s for the conjugation of ubiquitin to a substrate and is involved in the response to pathogen infection and low temperature stress (Wang et al., 2018). Loss of UBC13 blocks the expression of many cold‐responsive genes and causes hypersensitivity to low‐temperature stress in Arabidopsis (Wang et al., 2018). Also, a PA200 gene (BnaA03g32770D) is associated with marker SA03_15761546 within the LD decay interval on chromosome A03. PA200 is involved in DNA repair through proteasomal hydrolysis of peptides (Ustrell et al., 2002) and partakes in genomic stability that may be mediated through postglutamyl cleavage by proteasomes (Blickwedehl et al., 2008).

BRG2 (BnaA07g34660D) is one of the five gene models within the LD decay interval and is less than 13 Kb from the marker SA07_23526647; in addition, this gene was down regulated following a one‐month cold acclimation treatment. BOTRYTIS SUSCEPTIBLE1 INTERACTOR (BOI) and BRGs represent a subclass of RING E3 ligases that are critical in recognizing and binding poly‐ubiquitins on target proteins for degradation, and they are known to impact plant disease resistance and abiotic stress tolerance (Cho et al., 2017; Dong et al., 2006; Luo et al., 2010). Some Arabidopsis RING E3 ligases such as ARABIDOPSIS TÓXICOS EN LEVADURA78 and HIGH EXPRESSION OF OSMOTICALLY RESPONSIVE GENES1 are known to be negative regulators of cold stress in Arabidopsis (Ishitani et al., 1998; Suh & Kim, 2015). We also identified a RING E3 Ligase (BnaA04g27970D) that linked to marker S1_656249548 on chromosome A04 from our previous GWAS analysis (Chao et al., 2021). Thus, the BRG2 gene identified in this study could be involved in rapeseed stress tolerance. Interestingly, another ubiquitin‐related gene UCH (BnaA07g34670D) is adjacent to BRG2. UCHs belong to a subclass of the de‐ubiquitylating enzyme (DUB) family, which plays a role to regenerate the protein in a non‐modified form by cleaving ubiquitin from ubiquitin‐conjugated proteins (D'Andrea & Pellman, 1998; Wilkinson, 1997). UCHs are involved in the control of shoot architecture (Yang et al., 2007) and period maintenance of the circadian clock (Hayama et al., 2019) in Arabidopsis. Thus, ubiquitination (through UBC13A and BRG2) and de‐ubiquitination (through UCH) pathways may act together to sustain rapeseed freezing tolerance.

CYTOCHROME P450 in stress response

Six CYTOCHROME P450 genes (CYPs) from the subfamily CYP72A are mapped within the LD decay interval, and they are located 8 to 21 Kb away from the significant marker SA03_16201884 on chromosome A03. These gene models BnaA03g33480D, BnaA03g33490D, BnaA03g33500D, BnaA03g33510D, BnaA03g33520D, and BnaA03g33530D encode CYP72A9, CYP72A15, CYP72A10, CYP72A15, CYP72A15, CYP72A10, respectively. In addition, one (BnaA03g33530D) of them was up regulated following a one‐month cold acclimation treatment. CYPs belong to the oxidoreductases class of enzymes, and they catalyze NADPH‐ and/or O2‐dependent hydroxylation reactions in many organisms (Pandian et al., 2020). CYPs play important roles in heat and/or cold acclimation and are involved in the biosynthesis of secondary metabolites, antioxidants, and phytohormones as well as in the detoxification of xenobiotics (Pandian et al., 2020; Tao et al., 2017). A wheat CYP gene from the subfamily CYP72A, TaCYP72A, contributes to host resistance to mycotoxic Fusarium virulence factor deoxynivalenol as well as a positive impact on grain development (Gunupuru et al., 2018). It is therefore that some of these CYPs may play a role in freezing‐injury repair processes.

CONCLUSION

We have identified 24 candidate genes in the 14 loci associated with the chlorophyll fluorescence ratings used to measure freezing tolerance. Most markers had significant association with Fv/Fm and Fv/Fo readings at 7 days after freezing treatment. Thus, the candidate genes linked to those markers may play a general role in protecting photosynthetic machinery. Interestingly, a strong candidate gene encoding SFR2 is known to be involved in freezing sensing mechanism. In this study, we did not find CBFs, genes playing a major role in plant cold acclimation, in our GWAS results. Thus, the CBF regulon appears not to be a common target for evolution of variations in freezing tolerance of rapeseed, regardless of its well documented involvement in the process. The fact that many loci had genes associated with lipid, photosynthesis, flowering, ubiquitination, and cytochrome P450 strengthen the confidence in these observed associations. Although GWAS can offer invaluable first insight into genetic architecture, validation bias is a known limitation in GWAS. Therefore, it is critical to examine the functions of these candidate genes in rapeseed cold acclimation and freezing tolerance.

CONFLICT OF INTEREST

The authors declare that there is no conflict of interest.

AUTHOR CONTRIBUTIONS

Conceptualization, W.S.C., D.P.H., and J.V.A.; data analysis, W.S.C., X.L. and D.P.H.; investigation, W.S.C., D.P.H., J.V.A., and X.L.; methodology, W.S.C., X.L., D.P.H., and J.V.A.; writing—original draft preparation, W.S.C. and X.L.; writing—review and editing, W.S.C., D.P.H., J.V.A., and X.L. All authors have read and given approval to the published version of the manuscript.

Supplemental Figure S1. Scree plot of PCs.

Click here for additional data file.

Supplemental Figure S2. Distribution of the pairwise linkage disequilibrium measure r2 depending on the physical distance between SNPs.

Click here for additional data file.

Supplemental Table S1. Rapeseed accessions, estimated BLUEs, and class of 222 rapeseed accessions.

Click here for additional data file.

Supplemental Table S2. List of significant markers.

Click here for additional data file.

Supplemental Table S3. Six phenotypic data.

Click here for additional data file.

Supplemental Table S4. Significant markers and surrounding genes.

Click here for additional data file.

Supplemental Table S5. FASTQ key file.

Click here for additional data file.