Genome-Wide Association Mapping of Freezing Tolerance Loci in Canola (Brassica napus L.)

Abstract

Winter canola generally produces greater yields than spring canola. However, its range is limited due to its inability to withstand the harsh winter conditions that occur in many northern regions of the U.S.A. To identify loci associated with freezing tolerance in canola, we conducted a genome-wide association study (GWAS) using a genotyped diversity panel containing 399 accessions consisting primarily of winter canola. One-month-old greenhouse grown plants were subsequently cold-acclimated for two months in an environmental growth chamber prior to phenotyping for freezing survival using a visual damage scale and chlorophyll fluorescence (Fv/Fo). There was reasonable correlation observed between visual damage and chlorophyll fluorescence ratings among the top associated loci; the results indicated that some loci contributed to both freezing damage/tolerance and photosynthetic efficiency. The resulting numerical values for phenotypes were used for association analyses with the identified SNPs. Thirteen significant markers were identified on nine chromosomes for the phenotypes scored, with several showing significance for multiple phenotypes. Twenty-five candidate genes were identified as previously associated with freezing tolerance, photosynthesis, or cold-responsive in canola or Arabidopsis.

1. Introduction

Canola is an important oilseed crop contributing to the global demand for oil production. Winter canola generally has 20–30% greater yield potential than spring canola because its earlier spring flowering trait avoids later season heat that can impact flowering and yield [1]. Yields of winter canola typically exceed 2000 kg ha⁻¹ and can reach 5000 kg ha⁻¹ [2]. However, the winter survival of canola is in general low and inconstant in many regions of the northern U.S.A. [3]; thus, limiting yearly winter canola acreage.

Winter hardiness in canola requires both cold acclimation and freezing tolerance (intensity and duration) [4]. The process of cold acclimation involves environmentally-induced alterations in cellular constituents [5]. This includes the ability of plants to withstand osmotic dehydration triggered by extracellular ice nucleation that causes water to diffuse from cells [6] and reprogramming of gene expression networks [7,8,9,10], allowing plant tissues to survive periods of freeze/thaw cycles.

Cold acclimation is partially mediated by C-repeat (CRT) binding factors (CBFs), which function as transcription factors. CBFs bind to the CRT cis-element in the promoters of many cold responsive (COR) genes to regulate their expression [11]. Many products of the COR genes play vital roles in regulating cold acclimation and subsequent freezing tolerance [12,13]. Besides the CBF-mediated pathway, CBF-independent pathways have also been implicated in regulating cold acclimation and freezing tolerance in Arabidopsis [13,14,15]. Moreover, cold acclimated gene expression is gated by circadian responses [16], and expression of diurnally regulated genes was greatly altered during cold acclimation in Arabidopsis [17,18]. Photosynthetic cold acclimation is also considered an essential component of cold hardiness and freezing tolerance; it requires the complex interactions of low temperature, light and chloroplast redox responses [19,20]. Consequently, chlorophyll fluorescence is often used as a measure of the light-processing efficiency of plant photosystems [21]. The ratio of variable to maximum fluorescence (Fv/Fm) has been used to evaluate low temperature stress in plants [19,22]; however, the ratio of variable to minimal fluorescence (Fv/Fo) is considered a more sensitive stress detector that allows comparison of samples in the same dark-adapted state [23].

Although some genetic components of cold acclimation and freezing tolerance have been identified in model systems and confirmed in canola [24,25], it is clear that freezing tolerance is a complex trait under the control of many genes [3,12,13]. The advantages of next generation sequencing and high-throughput phenotyping offer an opportunity to identify novel genetic markers associated with freezing tolerance. Canola is an allo-tetraploid (AACC; 4x = 38) with a genome size of 1129–1235 MBp, and next generation sequencing technology has been used to produce a reference genome for canola [26]. Reference genomes allow researchers to re-sequence the genome to identify genetic differences between any individuals and the reference sequence [27], or to use techniques such as genotype-by-sequencing (GBS) to inexpensively obtain and map hundreds of thousands of single nucleotide polymorphisms (SNPs) from populations [28], facilitating the implementation of genome wide association studies (GWAS). GWAS is commonly used to assess and mark genomic regions that influence polygenic traits [29]. It is a method complementary to traditional QTL mapping in biparental crosses and has the capability to map multiple traits simultaneously using publicly available computer programs such as TASSEL [30]. However, GWAS requires a diverse population of a given species to utilize the genetic diversity that has built up through evolutionary time [29]. This genetic diversity provides the large number of polymorphic genetic markers needed to associate specific loci to the phenotypic diversity for any given trait in that population. Because the scale of diversity within the populations often results in high crossover rates between markers and the various mutations that might impact a trait of interest, GWAS requires exceptionally high marker density when the linkage disequilibrium decay rate is high [31].

Previously genotyped populations of canola have been used for GWAS studies [32,33,34,35,36,37,38]. However, the majority of these diversity panels were a mix of winter, semi-winter and spring lines, and none had more than 200 winter varieties included in the population. Here, we use a characterized GBS mapped population of primarily winter canola biotypes containing a greater number of high-quality SNPs and larger population size [39] than previously associated with survival after freezing stress in cold acclimated canola. A collection of 399 accessions consisting primarily of winter canola were used to conduct GWAS in this study, and freezing tests were performed under controlled environment conditions to limit unpredictable stochastic events common to field conditions.

2. Materials and Methods

2.1. Plant Material

Seeds from a well characterized winter canola diversity panel (399 accessions including 15 likely spring accessions that were previously incorrectly designated as winter biotypes) were obtained from Kansas State University. A list of 429 common cultivar names associated with the accession numbers are available [39]; however, 20 spring canola accessions were not included in freezing experiments because these plants would not survive the freezing treatments of −15 °C. An additional 10 winter accessions were also excluded from the freezing experiments due to lack of seeds. Plants were grown for one month under greenhouse conditions (16 h photoperiod with temperatures of 22–28 °C). The growth conditions were similar to cold-acclimation treatment of winter and summer-annual biotypes of Camelina described previously [40,41]. Briefly, seeds were sown into potting soil (Pro-Mix BX; Premier Horticulture Inc., Quakertown, PA, USA) in Deepot Cells (D60L: 6.4 cm × 35.6 cm, 983 mL volume; Stuewe & Sons, Inc., Tangent, OR, USA), placed in support trays, and grown in a greenhouse for 4 weeks with daily watering and weekly fertilization (20-20-20 General Purpose, Jack’s Professional Water-Soluble Fertilizer, J.R. Peters, Inc., Allentown, PA).

2.2. Phenotype Canola for Freezing Tolerance

Two freezing experiments were performed. Each experiment included three replications, and each replication contained 399 accessions with three individual plants per accession. Three Percival walk-in chambers (two freezing chambers and one vernalization chamber, see Figure 1) were used for these experiments. Each chamber can accommodate 30 flats, and each flat contains 15 cells. All flats and cells were numbered. For each replication, plants were arranged at random in the greenhouse using Excel. The flats, of randomized potted plants, were also randomized using Excel for their placement in the chambers. For vernalization, two replications were vernalized in each of the three chambers, and for freezing treatments, three replications were frozen in each of the two freezing chambers.

Optimal conditions for conducting freezing tests were determined by subjecting 2-month cold-acclimated plants to various temperatures (−10, −15, −17, and −20 °C) for 4 h and phenotyping for damage ratings 3 times/week over a 4-week period (data not shown). All plants survived at −10 °C; in contrast, all plants died at −20 °C and, occasionally, all died at −17 °C. Survival rates were around 50% after plants were subjected to a −15 °C freezing treatment followed by four weeks of recovery (Figure 2); thus, all freezing treatments used for this study were conducted at −15 °C.

Four weeks old plants were acclimated for two months at 5 °C under 12 h photoperiods in above mentioned walk-in chambers, and subsequently frozen in walk-in freezing chambers. All walk-in chambers were equipped with LED grow lights with an average intensity of 4500 lux at plant level (Lumibar pro, Lumigrow Inc. 1480 64th Street, Suite 150, Emeryville, CA 94608, USA) (Figure 1). Watering was applied as needed (normally once every two-three weeks) during cold acclimation. Freezing treatments consisted of a linear ramp-down in temperature from 4 °C to −15 °C over 8 h, held at −15 °C for 4 h, followed by a linear ramp-up in temperature to 4 °C over 8 h under a 12:12 h dark:light regime. Subsequently, the plants were transferred to a greenhouse under 16 h photoperiod with supplemental halogen lighting and photosynthetic efficiency was evaluated by chlorophyll fluorescence (Fv/Fo) at 3 days and 7 days after freezing as described by Horvath et al. [41]. Freeze damage was also evaluated by scoring the number of youngest leaves with or without visible injury (0 = no viable leaves; 1 = less than 2 unaffected leaves; 2 = 2 unaffected leaves; 3 = more than 2 unaffected leaves) 2 weeks following freezing treatment (see Supplemental Figure S1 for damage scores). The phenotypic data were collected from two separate experiments. Each experiment included three runs (with three technical reps per run) for a total of 9 plants assessed per accession. The data from each of two experimental runs were maintained and averaged separately for GWAS analysis (see below).

2.3. Genome-Wide Association Mapping of Winter Hardiness

These data averages were fed into a model to estimate the means taking into account the variation in the experimental design using the LS means method. Briefly, LSmean (least squares mean) estimates or fixed effects parameter estimates were calculated for the three phenotype response variables (Fv/Fo at 3 days, Fv/Vo at 7 days, and visual damage score) via the LSMEANS statement in the GLIMMIX Procedure in SAS/STAT 14.3/SAS ^® 9.4 (© 2016 SAS Institute Inc., Cary, NC, USA). For the normally distributed Fv/Fo values, the replicated accessions were fixed effects and hierarchical components of the experimental design (potted plant position in flat, flat position in chamber, and chamber within replicate) were modeled as random effects. The visual damage score was multinomial, and the cumulative logit link was applied in the model, and the fixed effects parameter estimates were generated with the SOLUTION option in the MODEL statement. These means were then used to identify SNPs associated with freezing damage using the program TASSEL5 using the 251,576 SNPs previously characterized for this population [39]. Twelve different models were assessed for all sets of phenotypic data. Six of the models were: a naïve General Linear Model (GLM) that provided permutated p-values as false discovery controls with no kinship or principle components (PC), a GLM with a single PC, a GLM with a PCA of 17 (controlling for about ½ of the available variation), and these same conditions with the addition of population structure added as defined by the program Structure indicating support for three primary populations (K = 3). Additionally, a Mixed Linear Model (MLM) with the same conditions were used for the GLM analysis plus a kinship matrix derived using a TASSEL5 sub-routine. The best models were identified using the least mean square difference (MSD) method [42], and they were a MLM with kinship and a PCA of 0, 1 and 17 (Supplemental Tables S1–S6 are examples of a GLM without kinship but with false discovery statistics and MLM with kinship and a PCA of 0). Q-Q plots were also generated to provide a graphic representation of the model fit (12 Q-Q plots are shown in Supplemental Figure S2 to reflect 12 different accessed models using “Fv/Fo 3d run2” as an example). In all but one case, the MLM models that did or did not include PCs plus kinship were the top three best models. The p-values less than 0.004 (representing the top 0.4% or roughly the top 1000 markers) of the SNPs were identified for each experimental run for all phenotypic measures (3 d Fv/Fo, 7 d Fv/Fo, and 2-week damage). Markers that were significant (p-values < 0.005) in both experimental runs for all three top models were considered significant, and 13 significant loci were identified (Supplemental Table S7).

2.4. Selection of Candidate Genes

Local LD decay rates for 500 kb region centered on each significant marker was determined from the LD decay table established for this population [39]. Genes within these associated loci were identified using the Canola Genome Browser (http://www.genoscope.cns.fr/brassicanapus/cgi-bin/gbrowse/colza/). Candidate genes within the associated loci were identified based on known or suspected gene function (Supplemental Table S8) and/or gene expression data from a related RNAseq study—Bioproject accession # PRJNA560411 in the NCBI SRA database [43].

3. Results

3.1. Phenotyping Canola Accessions Freezing Treatments

One-month-old plants from a winter canola diversity panel (399 accessions) were cold-acclimated for two months and phenotyped for freezing survival using a visual damage scale and chlorophyll fluorescence (Fv/Fo). Three phenotypic data sets were collected for conducting GWAS (Figure 3 and Supplemental Table S9). At 3 and 7 days after the freezing treatment, plants were phenotyped for chlorophyll fluorescence (Fv/Fo), and at 2 weeks, plants were phenotyped for visual damage. At 3 days, the average fluorescence readings of six separate experiments (18 plants total) were between 0.9 and 4.3, and the mean, medium, and mode values were 3.2, 3.3, and 3.2, respectively. At 7 days, the average fluorescence readings were between 0.9 and 4.4, and the mean, medium, and mode values were 3.0, 3.0, and 3.3, respectively. At 2 weeks, the average damage scores were between 0.1 and 2.4, and the mean, medium, and mode values were 1.2, 1.2, and 1.3, respectively. Among the three phenotypic data sets evaluated in this study, correlation coefficients (R) between fluorometer readings taken at three- and seven-days was greater (0.84) than coefficients observed between the two-week visual damage scores and the three- and seven-day fluorometer readings (0.6 and 0.72, respectively;

Table 1. Correlation coefficient between fluorometer readings (Fv/Fo) and visual damage scores (Frz).

Correlation of Visual Damage Scores at 2 Weeks Post Freezing with 3- and 7-Day Fv/Fo Data
	Frz Avg	3 Day Fv/Fo Avg	7 Day Fv/Fo Avg
Frz Avg	1
3 day Fv/Fo Avg	0.5975	1
7 day Fv/Fo Avg	0.7208	0.8409	1

).

3.2. Assessing Model Fit

We ran 12 different models to determine the impact of population structure, kinship, and principle component of 0, 1, and 17 (the later explain ~25% of the variation) on model fit. In almost every case, for individual experimental runs with LSMeans-estimated averages, the best three models were a MLM with kinship and a PCA of 0, 1 and 17 based on both the MSD analysis (

Table 2. Mean Squared Distance (MSD) analysis to assess the fit of the various models for each of the experimental runs for the three phenotypic analyses. The table is ordered based on the average MSD for all six analyses. Highlighted values indicate the best model for each analysis.

Model	MSD 2week_run1_Lsmeans	MSD 2week_run2_Lsmeans	MSD FvFo3d_run1_Lsmeans	MSD FvFo3d_run2_Lsmeans	MSD FvFo7d_run1_LSmeans	MSD FvFo7d_run2_LSmeans	Average
MLM Naïve kinship	2.36 × 10−5	3.10 × 10−4	7.98 × 10−5	1.02 × 10−5	1.18 × 10−5	2.71 × 10−6	7.30 × 10−5
MLM PC17 kinship	4.40 × 10−5	3.72 × 10−4	1.24 × 10−5	1.54 × 10−5	1.56 × 10−5	2.92 × 10−6	7.71 × 10−5
MLM PC1 kinship	4.48 × 10−4	4.33 × 10−4	8.47 × 10−5	1.23 × 10−5	1.03 × 10−5	2.46 × 10−6	1.65 × 10−4
GLM PC17 population	1.38 × 10−4	1.51 × 10−4	1.32 × 10−4	8.13 × 10−4	2.57 × 10−5	4.80 × 10−4	2.90 × 10−4
MLM PC17 kinship population	3.03 × 10−5	2.26 × 10−4	7.19 × 10−4	1.70 × 10−5	4.90 × 10−4	6.02 × 10−4	3.47 × 10−4
GLM PC17	3.25 × 10−4	3.10 × 10−4	1.76 × 10−4	8.75 × 10−4	3.01 × 10−4	5.09 × 10−4	4.16 × 10−4
MLM Naïve kinship population	1.08 × 10−3	1.19 × 10−4	5.90 × 10−4	2.49 × 10−3	1.78 × 10−3	2.09 × 10−3	1.36 × 10−3
GLM Naïve population	1.41 × 10−3	1.11 × 10−4	1.24 × 10−3	3.19 × 10−3	1.22 × 10−3	2.02 × 10−3	1.53 × 10−3
GLM PC1 population	1.48 × 10−3	1.09 × 10−4	1.21 × 10−3	3.20 × 10−3	1.24 × 10−3	2.3 × 10−3	1.54 × 10−3
GLM PC1	2.54 × 10−3	1.47 × 10−4	1.75 × 10−3	3.26 × 10−3	2.96 × 10−3	2.23 × 10−3	2.15 × 10−3
MLM PC1 kinship population	6.24 × 10−4	2.83 × 10−4	1.23 × 10−3	5.16 × 10−3	1.13 × 10−3	4.77 × 10−3	2.20 × 10−3
GLM Naïve	1.59 × 10−2	1.10 × 10−3	3.78 × 10−3	3.45 × 10−3	2.47 × 10−3	2.18 × 10−3	4.82 × 10−3

) and the representative Q-Q plots (Supplemental Figure S2: used Fv/Fo 3d run2 as an example). Although the addition of population structure and principle components improved the model fit for the GLM, population structure degraded model fit for most of the MLMs. MLM with kinship, with or without PCA generally provided the best fit in all but one case where a GLM with a PCA of 17 produced the best fit. This suggests that family structure rather than population structure has a greater impact on the fit of the models for this population, and kinship should be incorporated in any association models using this population. Representative Manhattan and interactive Q-Q plots are shown in Figure 4 and Figure 5.

3.3. Loci Associated with Freezing Tolerance Based on GWAS

The estimated means from the phenotypic data collected from two freezing experiments were entered into Tassel5 for association with SNP markers. For each phenotype evaluated, the top 1000 markers were identified (p-values were all less than 0.005 for any given analysis). Only thirteen markers were identified as being significant in both experiments in all three of the top models using the estimated phenotypic averages (

Table 3. p-Value of SNP–phenotype associations based on mixed linear model. Thirteen significant loci were identified among both experimental runs for any one of the phenotypes scored.

		MLM PC1						MLM PC17						MLM Naïv
		2 Week Damage		3 Day Fv/Fo		7 Day Fv/Fo		2 Week Damage		3 Day Fv/Fo		7 Day Fv/Fo		2 Week Damage		3 Day Fv/Fo		7 Day Fv/Fo
Chromosom	Mark	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2	Exp 1	Exp 2
A04_RANDOM	S1_656249548	1.53 × 10−4	1.20 × 10−4	1.01 × 10−4	2.97 × 10−3	9.39 × 10−6	1.05 × 10−3	2.82 × 10−3	1.14 × 10−3	9.42 × 10−4	6.70 × 10−3	5.77 × 10−5	1.05 × 10−3	8.04 × 10−4	3.24 × 10−4	1.13 × 10−4	3.04 × 10−3	8.67 × 10−6	1.07 × 10−3
A07	S1_155086644	8.06 × 10−3	6.93 × 10−5	9.51 × 10−3	1.00 × 10−3	1.15 × 10−4	1.19 × 10−3	8.94 × 10−4	1.42 × 10−4	8.86 × 10−3	2.78 × 10−3	6.15 × 10−5	1.19 × 10−3	1.20 × 10−3	7.78 × 10−5	9.88 × 10−3	1.04 × 10−3	1.15 × 10−4	1.23 × 10−3
A09	S1_209198432	5.74 × 10−2	2.50 × 10−3	2.09 × 10−3	3.16 × 10−4	2.51 × 10−3	3.94 × 10−4	9.03 × 10−2	6.76 × 10−3	2.39 × 10−3	8.64 × 10−4	1.42 × 10−3	3.94 × 10−4	8.12 × 10−2	4.21 × 10−3	2.17 × 10−3	3.23 × 10−4	2.53 × 10−3	3.98 × 10−4
ANN_RANDOM	S1_739997351	1.45 × 10−2	2.86 × 10−1	1.24 × 10−1	1.94 × 10−3	1.44 × 10−3	1.56 × 10−3	1.55 × 10−2	3.45 × 10−1	1.56 × 10−1	1.62 × 10−3	7.17 × 10−4	1.56 × 10−3	1.01 × 10−2	2.86 × 10−1	1.19 × 10−1	1.83 × 10−3	1.42 × 10−3	1.47 × 10−3
ANN_RANDOM	S1_738423768	4.77 × 10−2	7.73 × 10−3	9.94 × 10−3	9.17 × 10−3	2.19 × 10−3	2.45 × 10−3	1.08 × 10−1	1.37 × 10−2	1.11 × 10−2	2.38 × 10−3	1.08 × 10−3	2.45 × 10−3	4.43 × 10−2	1.32 × 10−2	1.09 × 10−2	3.75 × 10−3	1.82 × 10−3	4.86 × 10−4
C02	S1_316742524	2.05 × 10−3	1.60 × 10−3	2.49 × 10−1	1.45 × 10−1	1.64 × 10−1	9.41 × 10−2	6.35 × 10−4	2.87 × 10−3	2.50 × 10−1	1.74 × 10−1	1.50 × 10−1	9.41 × 10−2	6.03 × 10−4	6.85 × 10−4	2.59 × 10−1	1.35 × 10−1	1.64 × 10−1	8.81 × 10−2
C03	S1_359625016	1.97 × 10−1	1.55 × 10−2	2.10 × 10−2	3.33 × 10−3	1.44 × 10−3	2.60 × 10−3	2.10 × 10−1	2.96 × 10−2	2.45 × 10−2	1.09 × 10−2	2.07 × 10−3	2.60 × 10−3	2.40 × 10−1	1.81 × 10−2	2.20 × 10−2	3.45 × 10−3	1.46 × 10−3	2.69 × 10−3
C03	S1_376295769	2.37 × 10−3	1.69 × 10−4	3.09 × 10−1	3.23 × 10−3	4.36 × 10−1	1.36 × 10−3	7.00 × 10−4	7.86 × 10−4	3.80 × 10−1	7.00 × 10−3	4.13 × 10−1	1.36 × 10−3	6.28 × 10−4	2.05 × 10−4	3.14 × 10−1	3.30 × 10−3	4.31 × 10−1	1.38 × 10−3
C07	S1_528342773	5.40 × 10−2	1.11 × 10−2	1.04 × 10−2	3.63 × 10−3	1.81 × 10−3	4.74 × 10−4	4.10 × 10−2	2.55 × 10−2	9.65 × 10−3	5.54 × 10−4	1.01 × 10−3	4.74 × 10−4	1.06 × 10−1	1.54 × 10−2	9.41 × 10−3	9.52 × 10−3	2.22 × 10−3	2.52 × 10−3
C07	S1_515244307	1.59 × 10−3	4.79 × 10−4	6.92 × 10−1	2.58 × 10−2	3.90 × 10−1	4.00 × 10−3	4.89 × 10−4	1.66 × 10−4	6.57 × 10−1	4.67 × 10−3	3.91 × 10−1	4.00 × 10−3	4.62 × 10−4	6.20 × 10−4	6.82 × 10−1	2.61 × 10−2	3.85 × 10−1	4.06 × 10−3
C08	S1_567184601	7.64 × 10−3	4.51 × 10−2	4.29 × 10−1	2.65 × 10−3	6.63 × 10−4	2.37 × 10−3	8.74 × 10−4	7.76 × 10−2	3.86 × 10−1	7.53 × 10−3	5.02 × 10−4	2.37 × 10−3	1.17 × 10−3	4.19 × 10−2	4.32 × 10−1	2.59 × 10−3	6.56 × 10−4	2.31 × 10−3
C09	S1_641455172	1.05 × 10−1	8.96 × 10−1	1.22 × 10−4	1.49 × 10−3	6.94 × 10−3	1.34 × 10−2	3.97 × 10−1	8.76 × 10−1	2.35 × 10−3	4.82 × 10−3	1.99 × 10−2	1.34 × 10−2	4.23 × 10−1	9.65 × 10−1	1.29 × 10−4	1.86 × 1003	6.53 × 10−3	1.55 × 10−2
C09	S1_641455176	1.05 × 10−1	8.96 × 10−1	1.22 × 10−4	1.49 × 10−3	6.94 × 10−3	1.34 × 10−2	3.97 × 10−1	8.76 × 10−1	2.35 × 10−3	4.82 × 10−3	1.99 × 10−2	1.34 × 10−2	4.23 × 10−1	9.65 × 10−1	1.29 × 10−4	1.86 × 10−3	6.53 × 10−3	1.55 × 10−2

and Supplemental Table S7). Although there was poor correlation between the two experiments, the contaminating spring lines were consistently over-represented in both experiments among the lines with low freezing tolerance, suggesting that the phenotypic values were not random and thus subsequent associations are likely valid (Supplemental Table S10). Several markers were indicated in multiple phenotypes. There were reasonable correlations observed between chlorophyll fluorescence ratings and visual damage among the top associated loci, indicating that some loci contributed to both phenotypes (i.e., photosynthetic efficiency and freezing damage/tolerance).

Thirteen markers mapping to twelve loci were identified as significant (Table 3 and Supplemental Table S7); two markers on C09 (S1_641455172 and S1_641455176) are 3 base pair apart and therefore are represented by S1_641455172 from hereon. Associated loci are located on Chromosomes A04_RANDOM, A07, A09, ANN_RANDOM, C02, two on C03, two on C07, C08, and 2 on C09. A total of 6 Manhattan plots and 6 Q-Q plots were presented, representing 3 phenotypic data (2-week damage, 3 days Fv/Fo, and 7 days Fv/Fo) from two separate experiments ( Figure 4 and Figure 5).

3.4. Underlying Candidate Genes Associated with Freezing Tolerance

Gene models within 50 kb to either side of the twelve markers were identified using the Canola Genome Browser (Supplemental Table S8). Additionally, the precise boundaries of the associated loci based on the local LD decay rate in the region surrounding each marker was determined. Eight of the twelve loci had recognized gene models mapping within the LD decay intervals. In the remaining four cases, the markers (S1_656249548, S1_739997351, S1_376295769, and S1_515244307) fell within gene poor regions, and/or the LD decay rate was high and thus no strong candidate genes were identified for these loci. The LD decay interval was very small for seven of the twelve loci (highlighted in blue under column “Region of Interest” in Supplemental Table S8), and thus, assuming the association was not spurious, the markers from these loci were expected to map within a specific gene. In four of these cases, the marker itself fell within the transcribed portion of a gene (S1_738423768, S1_316742524, S1_567184601, and S1_641455172; highlighted in red and orange under column “Locus” in Supplemental Table S8), and two of these mapped to the coding region of the candidate gene (highlighted in red). In five remaining cases, the LD decay rate was large enough to cover several genes; however, two markers (S1_209198432 and S1_528342773) still fell into the transcribed region of a gene (highlighted in orange and red under column “Locus” in Supplemental Table S8). In these cases, any of the genes within the LD decay interval could be good candidates (highlighted in yellow, orange and red under column “Locus” in Supplemental Table S8).

Cold-regulated AGI orthologues are identified based on “Electronic Fluorescent Pictograph” [44]. Candidate genes that were cold-regulated in Arabidopsis (highlighted in blue under column “AGI orthologue” in Supplemental Table S8), or that played an obvious role in photosynthesis/light signaling (highlighted in green) are indicated. Genes observed to be differentially expressed in canola following a one-month cold acclimation treatment are also indicated (up-regulated genes are highlighted in green and down-regulated genes are highlighted in orange under column “Canola Gene Model” in Supplemental Table S8).

Six genes with various functions were found within 50 kb to either side of the marker S1_656249548 on chromosome A04_RANDOM. However, none were within the very narrow LD decay interval. Yet, a RING E3 LIGASE gene associated with E3 ligases play pivotal roles in protein degradation [45], and this particular gene was shown to increase biomass under drought stress when over expressed in Arabidopsis [46].

There are 29 gene models within the 50 kb surrounding either side of marker S1_155086644 on chromosome A07. LOW TEMPERATURE INDUCED 29 (LTI29) is very close to the marker as was an orthologue of SMALL AUXIN UPREGULATED RNA 60. Although outside the LD decay range, a CATALASE 1 (CAT1) gene that is cold regulated in Arabidopsis is also nearby; LTI29 is a dehydrin protein and CAT1 is an antioxidant enzyme. The auxin-regulated gene was previously associated with adaptation to high elevation [47].

There are 24 genes within the 50 kb surrounding either side of the significant marker S1_209198432 on chromosome A09. The marker S1_209198432 falls into an intron of the gene encoding tetratricopeptide repeat (TPR)-like superfamily protein, and another cold-regulated gene known to encode the 14-3-3 protein GENERAL REGULATORY FACTOR 12 (GRF12) protein in Arabidopsis is within the LD decay region. Another gene encoding a UREIDE PERMEASE 5 is also present within the LD decay interval that was shown to be induced by cold in canola (NCBI-SRA database accession # PRJNA560411). Two other genes falling just outside the LD decay interval were also cold-regulated in our parallel RNAseq study; they include an up-regulated orthologue of the Arabidopsis gene ROOT INITIATION DEFECTIVE 1 during cold acclimation and a down-regulated orthologue of UDP-GLUCOSE DEHYDROGENASE 1. Another gene of interest is a CYTOCHROME B5-A (CB5-A); CB5s are involved in electron-transfer in the desaturation of fatty acid [48].

There are eight genes within the 50 kb surrounding either side of the marker S1_739997351 on chromosome ANN_RANDOM. However, the LD decay interval is very short, and this marker does not fall within any of these genes. Thus, although there are several interesting gene models including an orthologue of Arabidopsis CYCLIN B2;3, CHLOROPLAST J-LIKE DOMAIN 1, JUMONJI DOMAIN-CONTAINING PROTEIN 16, and BASIC LEUCINE-ZIPPER 70, there is no strong evidence that any of these genes are involved in cold acclimation or freezing tolerance.

Another nearby locus surrounds the significant marker S1_738423768 on ANN_RANDOM. Again, the LD decay region is very small; however, there are 15 gene models within the 50 kb on either side of the marker. The marker itself, falls into the noncoding region of a gene model (BnaAnng22960D) that has no homology with any Arabidopsis genes. There is a carboxyl-terminal peptidase encoding gene that is slightly down-regulated by cold in Arabidopsis and is adjacent to the gene model in which the marker falls. Moreover, there is a gene encoding an orthologue of Arabidopsis PHOTOSYSTEM II LIGHT HARVESTING COMPLEX GENE B1B2 within 50 kb of the marker, as well as a RIBOSOMAL PROTEIN L18 and several F-box protein encoding genes.

There are only six gene models within the 50 kb of marker S1_316742524 on chromosome C02. Only three of these have homology to known Arabidopsis gene models. The LD decay intervals are exceptionally small in this area of the genome, and thus the only solid candidate gene is BnaC02g36140D which encodes an orthologue of a cold-regulated hypothetical protein in Arabidopsis.

There are eight genes within 50 kb of the significant marker S1_359625016 on Chromosome C03. Six of these fall within the LD decay interval, and five of these six are cold-regulated in Arabidopsis. These include orthologues of Arabidopsis genes encoding a SGNH hydrolase-type esterase superfamily protein, ribonuclease H-like superfamily protein, RESPONSE REGULATOR 6, RING/U-box superfamily protein, and basic-leucine zipper transcription factor K. The xanthine/uracil permease family protein falls within the LD decay interval but is not cold-regulated in Arabidopsis. Likewise, there are seven additional genes on chromosome C03 near the significant marker S1_376295769 located at position 52566692. Although none fall within the 23 kb LD decay interval on either side of this marker, there are a set of three pectinacetylesterase encoding genes and a disease resistance protein (TIR-NBS-LRR class) whose Arabidopsis orthologues are all cold-regulated.

Only three gene models fall into the LD decay interval near S1_528342773 on chromosome C07. However, the marker itself falls into the coding region of an orthologue of Arabidopsis CHLOROPHYLL A/B PROTEIN 140 (CAB140). This gene plays a role in photosynthesis and is up-regulated during cold acclimation in canola. The other two genes encode a NUCLEAR RNA POLYMERASE A2 and SU(VAR)3-9-RELATED PROTEIN 5 with a likely role in chromatin modifications. Either of these genes could conceivably be important in cold acclimation process. There is also a basic helix-loop-helix (bHLH) DNA-binding superfamily protein located within 50 kb of the marker that is cold-regulated in Arabidopsis. However, the strongest candidate is the CAB140 gene. Looking further out from this marker, there is an AMINOPHOSPHOLIPID ATPASE 12 (ALA12). Arabidopsis has 12 putative aminophospholipid translocases, ALA1 to ALA12, and ALA1 is involved in phospholipid transportation across the cell membrane and low-temperature tolerance [49,50]. There were no genes within 50 kb of our significant marker S1_515244307 on chromosome C07, and the narrow LD decay window suggest that this marker may be a false positive or is associated with a non-transcribed regulator of cold acclimation processes.

The LD decay rate around the significant marker S1_567184601 on chromosome C08 is quite small; however, this marker does fall into the noncoding region of an orthologue of Arabidopsis KNOTTED1-LIKE HOMEOBOX GENE 6 (AT1G23380). Although there is no indication that this gene plays a role in cold acclimation or freezing tolerance, it likely regulates meristem growth and development [51]. As such, it could be involved in the growth suppression and/or vernalization processes that generally accompanies long term cold treatments. Although outside the region of likely LD decay, there is a cluster of genes that are cold regulated in Arabidopsis and canola (Supplemental Table S8). These include orthologues of AT1G23310, AT1G23340 AT1G23350, AT1G23360 and AT1G23490 encoding ALANINE-2-OXOGLUTARATE AMINOTRANSFERASE 1, a carboxyl-terminal proteinase, a plant invertase/pectin methylesterase inhibitor superfamily protein, a 2-phytyl-1,4-naphthoquinone methyltransferase that catalyzes the final step in phylloquinone (vitamin K1) biosynthesis, and ADP-RIBOSYLATION FACTOR 1, respectively.

There are 19 gene models within the 50 kb surrounding either side of marker S1_641455172 (and S1_641455176) on C09. The LD decay rate around the S1_641455172 marker is also quite small; however, this marker falls into the coding region of an orthologue of Arabidopsis gene (AT5G13950) encoding a nuclear factor related to kappa-B-binding protein. Although outside the region of likely LD decay, two orthologues of AT5G13970 and AT5G13890 are cold regulated in Arabidopsis. In addition, a CHALCONE SYNTHASE gene is linked to this marker; CHS is a key enzyme in the flavonoid biosynthesis pathway [52]. A MITOCHONDRIAL PHOSPHATE TRANSPORTER (MPT) is also associated with this marker; MPT protein regulates phosphate transport and ATP synthesis to maintain cellular energy homeostasis in plant cells [53].

4. Discussion

In this study, significant variation in freezing tolerance and photosynthetic responses were observed within our population of canola (Figure 3), which made it possible to identify alleles associated with these traits by GWAS. The GWAS analyses identified thirteen markers with p-values < 0.005 (The p value equivalent of the top 1000 or 0.4% of the markers) in both experimental runs in all three of the top association models. Although the p value cutoff for any given run was fairly relaxed compared to other such studies, this method was stringent in the requirement for significance across multiple independent experiments and multiple association models. Thus, our confidence that few false positives were identified is bolstered by the observation that areas with rapid linkage decay, in most cases, had identified markers that fell within a recognized gene. We do have to note that although we phenotyped our population for freezing tolerance, the fact that the contaminating spring lines were over-represented among the varieties with consistently low freezing tolerance opens the possibility that some of the markers we identified may be associated with loci that control both freezing tolerance and spring biotype growth habits.

4.1. Comparison to Other Association Mapping Studies Identify Additional Regions of Interest

There have been many other genetic analyses of freezing tolerance in canola, but most have not associated the markers with the known physical map of the canola genome. However, one such study [54] did identify a single quantitative trait locus on chromosome A08 between 1,069,000 to 15,652,000 bp. Although no markers in this region met our criteria for significance, a region on A08 between 12,330,000 and 12,370,000 bp has multiple markers with low p-values across all three phenotyping studies in the second run of our experiments.

Several GWAS studies on freezing tolerance in canola, under various conditions and phenotyping strategies, have been published [32,33,34,35,36,37,38]. However, none of the loci associated with these studies were identified as significant in both experimental runs in our study. There are, however, multiple makers with low p-values near the marker on chromosome A06 near 21,440,000 bp that were significantly associated with differences in freezing tolerance as defined by electrolyte leakage analysis [55], but these markers are only in association with Fv/Fo at day 3 and are not consistent in both experimental runs. There are two markers that have fairly low p-values in experimental run 2 for visual damage scores and for experimental run 1 for Fv/Fo at 7 days that are within 50 kb of a marker on A03 at position 20,587,761 bp associated with freezing tolerance under simulated conditions [38]. There is also a cluster of markers within 150 kb of the significant marker on chromosome A07 at 13,505,270 bp, another cluster near the marker on chromosome A09 at 8,152,488 bp, and one more near the marker on A09 at position 8,867,241 bp that were significant in a GWAS study of freezing tolerance under growth chamber freezing conditions [56]. There is one marker (S1_306342432) very near the significant locus on C02 for freezing tolerance under field conditions that has low p-values in experiment one for both 2-week damage scores and Fv/Fo at 7 days. Likewise, there is a cluster of markers mapping to chromosome C04 between 23,260,477 and 23,260,521 bp that are closely linked to a marker associated with freezing tolerance under growth chamber conditions that mapped to 23,256,454 bp in a previous study [56]. Thus, although we did not find any highly significant loci previously identified by earlier studies, our study provides some additional support for freezing tolerance loci identified in those studies.

4.2. Phenotypic Measures of Visual Damage and Photosynthetic Damage are Correlated

Different phenotypic data evaluated in this study complemented one another, as indicated by their correlation coefficients. However, the two fluorometer readings had a stronger positive correlation with each other compared to their correlation with the visual damage scores. Because these correlations have a similar trend, photosynthetic efficiency appears to be a reasonable early indicator of freezing damage. Chlorophyll fluorescence can be a useful indicator of cellular damage, particularly damage to photosystems, following plant stress in other systems [21]. Indeed, several markers (loci) had strong associations across all three phenotypes, indicating common mechanisms are likely involved. However, these related phenotypic measures were not completely correlated, suggesting that other factors may be involved as freezing-induced cellular degradation progresses over time. For example, freezing-induced damage to photosystems impacts the dissipation of light energy (fluorescence or heat), often resulting in the formation of free radical oxygen or reactive oxygen species (ROS), leading to downstream damage to cellular homeostasis. The long-term accumulation and consequences of exposure to ROS can lead to progressive damage over time, unless mitigated by cellular repair mechanisms [21]. Our results suggest that different physiological and molecular mechanisms are involved in cold acclimation and freezing tolerance; the potential roles of these candidate genes are discussed below. Our confidence in the associations of these markers with phenotypes is strengthened by the multiple candidate genes with known or suspected roles in freezing tolerance identified within the associated loci.

4.3. Candidate Genes Associated with Freezing Tolerance

RING E3 Ligases: Although the gene sequence for RING E3 Ligase (BnaA04g27970D) was not associated with marker S1_656249548, based on the local linkage decay rate (

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	Names
S1_656249548	A04_random	522883	521883-523883	BnaA04g27970D	AT5G37890	RING E3 LIGASE (SINAL7)
S1_155086644	A07	10641317	10625317-10657317	BnaA07g11360D	AT1G20630	CATALASE 1 (CAT1)
S1_155086644	A07	10641317	10625317-10657317	BnaA07g11370D	AT1G20620	CATALASE 3 (CAT3)
S1_155086644	A07	10641317	10625317-10657317	BnaA07g11440D	AT1G20470	SMALL AUXIN UPREGULATED RNA 60 (SAUR60)
S1_155086644	A07	10641317	10625317-10657317	BnaA07g11450D	AT1G20450	LOW TEMPERATURE INDUCED 29 (LTI29)
S1_209198432	A09	21784443	21771443-21797443	BnaA09g28870D	AT1G26340	CYTOCHROME B5-A (CB5-A)
S1_209198432	A09	21784443	21771443-21797443	BnaA09g28880D	AT1G26370	ROOT INITIATION DEFECTIVE 1 (RID1)
S1_209198432	A09	21784443	21771443-21797443	BnaA09g28910D	AT1G26440	UREIDE PERMEASE 5 (UPS5)
S1_209198432	A09	21784443	21771443-21797443	BnaA09g28950D	AT1G26480	GENERAL REGULATORY FACTOR 12 (GRF12)
S1_209198432	A09	21784443	21771443-21797443	BnaA09g28990D	AT1G26560	BETA GLUCOSIDASE 40 (BGLU40)
S1_209198432	A09	21784443	21771443-21797443	BnaA09g29000D	AT1G26570	UDP-GLUCOSE DEHYDROGENASE 1 (UGD1)
S1_739997351	ANN-RANDOM	27349546	27348546-27350546	BnaAnng23820D	AT1G08640	CHLOROPLAST J-LIKE DOMAIN 1, CJD1
S1_738423768	ANN_RANDOM	25775963	25774963-25776963	BnaAnng22920D	AT2G34420	PHOTOSYSTEM II LIGHT HARVESTING COMPLEX B1B2 (LHB1B2)
S1_738423768	ANN_RANDOM	25775963	25774963-25776963	BnaAnng22950D	AT5G18460	Carboxyl-terminal peptidase
S1_316742524	C02	39235351	39234351-39236351	BnaC02g36140D	AT2G19890	Hypothetical protein
S1_359625016	C03	35895939	35877939-35913939	BnaC03g51340D	AT5G62920	RESPONSE REGULATOR 6 (ARR6)
S1_376295769	C03	52566692	52543692-52589692	BnaC03g63190D	AT4G19420	PECTIN ACETYLESTERASE 8 (PAE8)
S1_528342773	C07	14698486	14397486-14999486	BnaC07g09400D	AT1G29930	CHLOROPHYLL A/B PROTEIN 140 (CAB140)
S1_567184601	C08	8769737	8770737-8768737	BnaC08g06270D	AT1G23310	ALANINE-2-OXOGLUTARATE AMINOTRANSFERASE 1 (AOAT1)
S1_567184601	C08	8769737	8770737-8768737	BnaC08g06310D	AT1G23360	MENG
S1_567184601	C08	8769737	8770737-8768737	BnaC08g06320D	AT1G23380	KNOTTED1-LIKE HOMEOBOX GENE 6 (KNAT6)
S1_641455172	C09	44563121	44562121-44564121	BnaC09g43140D	AT5G14040	MITOCHONDRIAL PHOSPHATE TRANSPORTER (MPT)
S1_641455172	C09	44563121	44562121-44564121	BnaC09g43220D	AT5G13970	Midasin-like protein
S1_641455172	C09	44563121	44562121-44564121	BnaC09g43230D	AT5G13950	Nuclear factor kappa-B-binding protein
S1_641455172	C09	44563121	44562121-44564121	BnaC09g43250D	AT5G13930	CHALCONE SYNTHASE (CHS)

and Supplemental Table S8), it was less than 30 kb from the most highly significant marker on chromosome A04. This marker had significant association with freezing damage scores at 2-weeks and Fv/Fo readings at 3 and 7 days after freezing treatment. Thus, the candidate gene associated with this marker may have a common mechanism in protecting photosynthetic machinery and mitigating tissue damage. Ubiquitin proteasome system, a mechanism for protein turn over, is involved in the stress responsive pathways, and E3 ligases play pivotal roles in recognizing and tethering poly-ubiquitins on target proteins for degradation [45]. Some Arabidopsis RING E3 ligases such as HIGH EXPRESSION OF OSMOTICALLY RESPONSIVE GENES1 (HOS1) and ARABIDOPSIS TÓXICOS EN LEVADURA78 (AtATL78) were found to be negative regulators of cold stress in Arabidopsis [11,57,58]. For example, the interaction of HOS1 and INDUCER OF CBF EXPRESSION 1 (ICE1) mediates the ubiquitination of ICE1, resulting in cold induced degradation of ICE1, and the RNAi mediated suppression of AtATL78 increases cold tolerance. In contrast, the Arabidopsis SEVEN IN ABSENTIA LIKE 7 (SINAL7) protein also possesses a RING motif and E3 ubiquitin-ligase activity, and over-expression of SINAL7 protein resulted in an increase in plant biomass and drought tolerance [46]. Thus, the RING E3 Ligases in canola could be involved in drought tolerance or play similar roles as that of HOS1 and AtATL78 in weakening cold stress responses; however, the definitive role of RING E3 Ligases in canola remains to be elucidated.

Dehydrin Protein LTI29: BnaA07g11450D encodes a LOW TEMPERATURE INDUCED 29 (LTI29) and is within the linkage decay interval and is the closest gene (less than 500 bases) to marker S1_155086644 (Table 4 and Supplemental Table S8) on chromosome A07. This marker was associated with all three phenotypes evaluated in this study. LTI29 can be induced by low temperature, ABA, and dehydration in Arabidopsis [59], and was strongly down regulated in a parallel transcriptomic analysis of acclimation responses of canola. LTI29 belongs to a family of acidic SK-type dehydrins [60], and mass spectrometry studies indicate LTI29 (AT1G20450) was rapidly synthesized during cold acclimation [61]. Immunolocalization studies showed that LTI29 is located in the cytoplasm and intracellular membranes (envelopes of plastids, plasma membrane, tonoplast) in Arabidopsis, and LTI29 labeling increased dramatically in the membrane after cold acclimation [62]. The overexpression of LTI29 together with LTI30 in Arabidopsis confers freezing tolerance both in nonacclimated and cold-acclimated plants [62]. How dehydrins contribute to plant survival during dehydration and freezing stress is not well understood, although considerable evidence points to a role in stabilizing membrane structure and sustaining the activity of essential enzymes [63]. Thus, LTI29 may be involved in early protection and mitigation of cell membranes from freezing damage. It seems highly likely that sequence variation in this gene may be responsible for some of the phenotypic difference we observed in our study, and identifying the causative mutations associated with enhanced freezing tolerance could allow for targeted modification to increase freezing tolerance of less tolerant but otherwise better adapted lines.

Antioxidant Enzyme Catalases CAT1 and CAT3: BnaA07g11360D and BnaA07g11370D are near the same locus as LTI29 but are outside the linkage decay range of the significant marker, S1_155086644. These two genes are orthologues of Arabidopsis CATALASE 1 and 3 (CAT1 & 3), which encode antioxidant enzymes (Table 4 and Supplemental Table S8). Freezing stress increases the production of various reactive oxygen molecules via free radicals, which lead to plant oxidative stress through the oxidation of biotic molecules [64]. Hydrogen peroxide (H₂O₂) damages membranes through the peroxidation of unsaturated fatty acids in membrane lipids [65]. Catalase is an important antioxidant enzyme that plays a role in mitigating effects of H₂O₂. Several studies have shown that cold tolerant genotypes had more antioxidant activity under cold stress [66,67,68]. Likewise, Valizadeh-Kamran et al. [68] found that H₂O₂ increase was more evident in the sensitive genotypes of barley (Hordeum vulgare) than in tolerant ones. The identification of CAT1 and CAT3 in our GWAS analysis supports the hypothesis that they may play a role in canola freezing tolerance.

CYTOCHROME B5-A (CB5-A): CB5-A is a membrane protein and is involved in the desaturation of fatty acid. The Arabidopsis homolog, CB5-A, is similar to BnaA09g28870D that, although not in the linkage decay intervals, is within 50 kb of marker S1_209198432 on chromosome A09 (Table 4 and Supplemental Table S8). CB5-A is a heme binding membrane protein that localizes to the endoplasmic reticulum in plants [69]. Five cytochrome b5 proteins (CB5-A to CB5-E) are identified in Arabidopsis [70,71]. However, the specific functions of individual CB5 proteins are unknown; there is likely to be some degree of functional redundancy. CB5 plays a role in electron-transfer in the desaturation of fatty acid [48], and plants can exhibit enhanced freezing tolerance by increasing the levels of unsaturated fatty acids [72]. CB5-A is linked to a marker (S1_209198432) associated with chlorophyll fluorescence (Fv/Fo) readings at 3 and 7 days after freezing treatment. Thus, CB5-A may be involved in protecting cells against cold stress by increasing membrane fluidity.

CHALCONE SYNTHASE (CHS): CHS is a key enzyme that catalyzes the first committed step in the flavonoid biosynthesis pathway [52]. Flavonoids play many vital roles in plant development and defense [73]. These metabolites also accumulate after cold acclimation and are thought to be important for freezing tolerance. Several Arabidopsis mutants with reduced flavonoid content exhibited impaired freezing tolerance though the involvement of flavonoids to freezing tolerance appeared to be genotype-dependent [74]. A chalcone synthase gene CHS (BnaC09g43250D) is linked to markers S1_641455172 and S1_641455176 on C09 (Table 4 and Supplemental Table S8). These two markers showed significant association primarily with chlorophyll fluorescence (Fv/Fo) readings at 3 day after freezing treatment. Thus, this CHS gene could play a role in flavonoid biosynthesis and subsequent involvement in freezing tolerance through protection of early phase of cellular damage.

MITOCHONDRIAL PHOSPHATE TRANSPORTER (MPT): The mitochondrial membrane localized phosphate transporter (MPT) regulates phosphate transport and ATP synthesis to maintain cellular energy homeostasis in plant cells [53]. A canola MPT gene (BnaC09g43140D) was found to be associated with marker S1_641455172 and S1_641455176 on C09 (Table 4 and Supplemental Table S8). These markers showed significant association primarily with chlorophyll fluorescence (Fv/Fo) readings at 3 day after freezing treatment, and thus the associated candidate MPT may play a role in protecting the early phase of cellular damage, as aforementioned CHS. Transcript abundance of MPT genes (AtMPTs) can be induced by high salinity stress in Arabidopsis seedlings [75] and by chilling treatment in tree peony (Paeonia suVruticosa) [76]. In addition, both transcript and protein levels increased dramatically during freezing in the Wood Frog Rana sylvatica, a vertebrate animal [77]. This canola MPT gene may play a role in modulation of cellular energy homeostasis during early phase of damage recovery.

4.4. Other Genes of Interest

Any of the genes near our significant markers could be responsible for variation in freezing tolerance observed in our population of canola. In addition to the genes discussed above, other linked genes of interest include those also cold-regulated in canola and/or other systems (Table 4 and Supplemental Table S8). For example, orthologues of Arabidopsis genes SMALL AUXIN UPREGULATED RNA 60, on chromosome A07 within the linkage decay region has been associated with adaptation to high altitudes and low temperatures in Arabidopsis [47]. This gene is adjacent to the LTI29 gene, and if there is sequence variation in these two genes, they could have a combined influence at this locus on freezing tolerance.

On chromosome A09, there are also two additional genes orthologous to Arabidopsis ROOT INITIATION DEFECTIVE 1—just outside the LD decay interval, and UREIDE PERMEASE 5—inside the LD decay interval, that were up-regulated in the parallel RNAseq study. A 14-3-3 protein encoding gene GENERAL REGULATORY FACTOR 12, located within the LD decay interval, was also up-regulated by cold in Arabidopsis, as is a cold responsive orthologue of BETA GLUCOSIDASE 40—mapping just outside the LD decay interval. Finally, a UDP-GLUCOSE DEHYDROGENASE 1—encoding gene located just outside the LD decay interval was strongly down-regulated in canola after four weeks of acclimation.

On Chromosome ANN_RANDOM, a gene encoding a carboxyl-terminal peptidase—an orthologue of AT5G18460, which is slightly down-regulated by cold in Arabidopsis could have potential role in freezing survival. However, the gene encoded by BnaAnng22960D (Supplemental Table S8) that has no homologue in Arabidopsis is the most likely candidate gene for this locus since the LD decay interval is very small and the marker fall into the transcribed region of this gene.

Only one likely candidate gene is present in the gene-poor region surrounding the significant marker on chromosome C02. The LD decay interval in this area is very small, but this marker does fall into the coding region of an orthologue of AT2G19890. This gene encodes a putative protein involved in gene regulation that is down-regulated by cold in Arabidopsis.

There are five out of six orthologous genes mapping to chromosome C03 within the LD decay interval that are all cold regulated in Arabidopsis. This suggests that perhaps this whole region may constitute a cold-responsive chromosomal segment. Although most of these genes are non-regulatory in nature, the one closest to the significant marker S1_359625016 encodes an orthologue of RESPONSE REGULATOR 6. This gene is known to play a role in cytokinin signaling in Arabidopsis [78]. Likewise, there are a series of pectinacetylesterase-encoding genes near the significant marker S1_376295769 at position 52,566,692 on chromosome C03 that are all cold regulated in Arabidopsis (Table 4 and Supplemental Table S8). These three genes all are 3′-untranslated region of the gene and adjacent to this locus, but outside the LD decay interval around the significant marker. The only other potential candidate gene in this region (outside the narrow LD decay interval) is BnaC03g63230D that encodes a disease resistance protein (TIR-NBS-LRR class). This gene is repressed by cold in Arabidopsis (Supplemental Table S8).

Two genes are located close to the significant marker S1_567184601 on chromosome C08—encoding orthologues of ALANINE-2-OXOGLUTARATE AMINOTRANSFERASE 1 and MENG (a 2-phytyl-1,4-naphthoquinone methyltransferase that catalyzes the final step in vitamin K1 biosynthesis); both are downregulated in canola after four weeks of acclimation. Due to their responsiveness to prolonged cold treatment, they may play a functional role in cold acclimation and freezing tolerance. Likewise, two genes encoding a carboxyl-terminal proteinase and a putative plant invertase/pectin methylesterase inhibitor are located close to the significant marker and are both differentially expressed in response to cold in Arabidopsis (Supplemental Table S8). However, as was noted in other cases, the significant marker falls into a LD decay interval that is quite small but is mapped to the transcribed region of an orthologue of KNOTTED1-LIKE HOMEOBOX GENE 6 from Arabidopsis, and thus is best candidate gene for this region.

Similarly, BnaC09g43230D encodes an orthologue of a nuclear factor kappa-B-binding protein gene, and the significant marker S1_641455172 maps into the coding region of this gene. No other genes fall within the very narrow LD decay interval at this locus. However, there is an adjacent gene encoding a midasin-like protein that is slightly cold-regulated in Arabidopsis. Likewise, the previously discussed CHALCONE SYNTHASE gene also maps near the significant marker on chromosome C09. Further study will be needed to determine if sequence variation in any of these other genes play a significant role in the diversity of the freezing tolerance of canola.

4.5. Photosynthetic and Light Signaling Processes Are Correlated with Cold Tolerance

Three photosynthetic and light signaling genes, CHLOROPLAST J-LIKE DOMAIN 1 (CJD1), PHOTOSYSTEM II LIGHT HARVESTING COMPLEX GENE B1B2 (LHB1B2), and CHLOROPHYLL A/B PROTEIN 140 (CAB140) were located outside the region of the LD decay intervals. However, they did fall within 50 kb of two markers on chromosomes ANN_RANDOM, and one which fell into the coding region of a photosynthetic gene on chromosome C07 that was also significantly differentially expressed in our parallel RNAseq analysis (Table 4 and Supplemental Table S8). Down-regulation of photosynthetic processes has been associated with survival after freezing in multiple species, and several lines of evidence suggested that it is important for plants to protect the photosynthetic apparatus from damage during low-temperature light exposure [79]. Indeed, differences between the growth and freezing tolerance of cold acclimated winter and spring cereals have been reported to involve different levels of photoinhibition [80,81,82]. The linkage of photosynthetic/light signaling genes to these markers strengthens an association between cold and photosynthetic processes to mitigate damage of photosynthetic apparatus during cold or freezing conditions.

Transcriptomics analysis comparing naturally occurring variants with differences in freezing tolerance have provided additional evidence for photosynthetic processes playing a significant role in freezing tolerance differences within species. For example, microarray studies on two near isogenic lines of wheat selected for variation in freezing tolerance identified genes involved in photosynthetic processes as over-represented in a gene set enrichment analysis (GSEA) [83]. Likewise, Liu et al. [84] observed over-representation of photosynthetic genes associated with genes down-regulated in the northern biotype after differential expression analysis in alligator weed (Alternanthera philoxeroides) that compared populations from southern China to a population that rapidly increased its range northward into China. Finally, a transcriptome analysis of two camelina (C. sativa) varieties with strong differences in freezing tolerance also implicated photosynthetic processes as a major difference between the two lines. Specifically, genes associated with chloroplasts were correlated with differential survival after low temperature treatments between these two camelina biotypes. Additionally, a sequence ‘TGGCCCTCGCTCAC’ was over-represented among genes that were down-regulated in the freezing tolerant camelina biotype Joelle [41]. This sequence is over-represented in genes with the ontology of “chloroplast” in Arabidopsis and may be the binding site for light-induced promoters [85]. This result suggested that some sequence-specific DNA-binding factors were less functional in freezing tolerant camelina biotype Joelle relative to the sensitive biotype CO46. In fact, these observations may imply that freezing-tolerant biotypes of camelina, such as Joelle, have undergone greater modifications to reduce oxidative stress during cold-acclimating conditions. Thus, the observation of many photosynthetic and light signaling genes linked to markers associated with freezing tolerance in canola strengthens the hypothesis that plants have evolved effective mechanisms to protect the photosynthetic apparatus from damage and/or induction of ROS during low-temperature light exposure.

4.6. CBF and Freezing Tolerance

Unexpectedly, we did not see the conventional freezing tolerance gene CBF in our GWAS results, although it is known to enhance freezing survival in several Brassica species [86,87]. However, in all of the transcriptomics studies noted above, expression of CBF was only moderately indicated as a major factor differing between the tolerant and susceptible biotypes. Similar phenomena were also observed in soybean (Glycine max) [88], potato (Solanum ssp) [89], and tomato (Solanum lycopersicum) [79], which all appeared to have cold-induced CBF regulons, but these regulons were inadequate for significantly improving freezing survival. Therefore, the CBF regulon may not be a strong target for evolution of differences in freezing tolerance, despite its well documented importance in the process. However, a QTL identified in a cross between Swedish and Italian Arabidopsis varieties with a known difference in freezing tolerance did highlight a locus that included CBF and in which variation in the CBF gene itself was indicated [90]. Given the importance of CBF in multiple stress responses, it is possible that changes directed at genes encoding this small gene family are highly deleterious on an evolutionary scale, and that allelic variation are either rapidly fixed in the population or lost. However, there are some cases where CBF variants were associated with freezing tolerance when variants of CBF were specifically collected and examined for their impact on freezing tolerance (see for example [91]). Thus, it is possible that CBF variants in our diversity panel may have been missed.

5. Conclusions

Twenty-five candidate genes were identified in the 13 loci associated with the three phenotypes used to quantify freezing tolerance in this study. Many of these genes fall within physiological and molecular processes previously shown to be involved in protecting plant photosystems, cellular membranes, and RNA integrity, among others, from the consequences of oxidative stress associated with freezing damage. Interestingly, a candidate gene encoding RING E3 Ligase is known to weaken cold tolerance, and three significant markers (S1_316742524, S1_528342773, and S1_641455172) fell into the coding region of the candidate genes. The fact that many loci contained cold tolerance associated genes strengthen the confidence in these observed associations. However, confirmation bias is a recognized hazard in GWAS. Thus, it would be important to functionally test these candidate genes for their role in canola freezing tolerance. However, these data do provide a starting point for breeders to apply genomic selection for improving freezing tolerance in canola.