Genetic mapping and identification of a QTL determining tolerance to freezing stress in Fragaria vesca L.

Extreme cold and frost cause significant stress to plants which can potentially be lethal. Low temperature freezing stress can cause significant and irreversible damage to plant cells and can induce physiological and metabolic changes that impact on growth and development. Low temperatures cause physiological responses including winter dormancy and autumn cold hardening in strawberry (Fragaria) species, and some diploid F. vesca accessions have been shown to have adapted to low-temperature stresses. To study the genetics of freezing tolerance, a F. vesca mapping population of 143 seedlings segregating for differential responses to freezing stress was raised. The progeny was mapped using ‘Genotyping-by-Sequencing’ and a linkage map of 2,918 markers at 851 loci was resolved. The mapping population was phenotyped for freezing tolerance response under controlled and replicated laboratory conditions and subsequent quantitative trait loci analysis using interval mapping revealed a single significant quantitative trait locus on Fvb2 in the physical interval 10.6 Mb and 15.73 Mb on the F. vesca v4.0 genome sequence. This physical interval contained 896 predicted genes, several of which had putative roles associated with tolerance to abiotic stresses including freezing. Differential expression analysis of the 896 QTL-associated gene predictions in the leaves and crowns from ‘Alta’ and ‘NCGR1363’ parental genotypes revealed genotype-specific changes in transcript accumulation in response to low temperature treatment as well as expression differences between genotypes prior to treatment for many of the genes. The putative roles, and significant interparental differential expression levels of several of the genes reported here identified them as good candidates for the control of the effects of freezing tolerance at the QTL identified in this investigation and the possible role of these candidate genes in response to freezing stress is discussed.

Introduction a peat-based compost (90% peat, 10% clay), with the addition of 1:5 (v/v) of granulated perlite, and were raised in a glasshouse for five weeks at 20 ± 2˚C and an 18-hour photoperiod. The plants were watered twice a week with a balanced nutrient solution containing 7.8 mmol N, 1 mmol P, and 4.6 mmol K per litre (used in 1:100 dilution).

Freezing tolerance phenotyping
Prior to low temperature testing, the plants were acclimated at 2˚C and a 12-hour photoperiod for six weeks. Plants were watered with cold water as needed. After six weeks, plants were transferred to three programmable freezers where they were first kept at -2˚C for 12 hours. Subsequently the temperature was lowered by 2˚C/h until the target temperature was reached. The target temperature was held for 4 hours before raising it by 2˚C/h to 2˚C where they were maintained for a further 10 hours. Subsequently, the plants were kept overnight at room temperature, following which they were transferred to a greenhouse and maintained at 18 ± 2˚C with an 18-hour photoperiod for five weeks before survival (plants were observed to be dead or alive) was scored.

Determination of suitable freezing temperatures for mapping population screen
The parental lines 'Alta' and 'NCGR1363,' selected for their differential response to freezing stress in a previous study [7] and an F 1 hybrid line from the resultant cross 'NCGR1363' × 'Alta' was phenotyped to determine suitable temperatures at which to screen an F 2 test population. The lines were subjected to freezing stress temperatures of -18˚C, -15˚C, -12˚C, -9˚C, -6˚C, and 0˚C following the procedure described above. For each temperature, 13 test plants of each genotype were screened.

Freezing tolerance phenotyping of a segregating mapping population
The grandparental lines, F 1 parent and an F 2 mapping population consisting of 143 plants from a selfing of the cross 'NCGR1363' × 'Alta' were propagated and subjected to cold tolerance screening. The optimal stress temperatures calculated from the progenitor germplasm screen (described above), to which the genotypes were subjected, were -5˚C, -8.5˚C, and -12˚C. Each F 2 genotype was replicated nine times at each temperature, whilst the grandparental and F 1 hybrid genotypes were replicated 18 times at each temperature. Five replicates of the entire experiment were performed at each temperature.

Statistical data analysis
The analysis of the survival data (alive/dead) from both the experiment to determine optimal stress temperatures and the subsequent screening of the F 2 mapping population, employed the following logistic model: where π ijkt is the observation [alive(1)/dead(0)] made on plant k from genotype i, in replicate j, exposed to temperature t, β 0 is an unknown constant, α i is the main effect of the genotype (i = 1, . . ..,n), βt is the coefficient that estimates the effect temperature (t) has on plant survival, E j is the effect of replicate j (j = 1,. . ., 5), (Eα) ij is the interaction between the genotype i and replicate j.
The LT 50 for genotype i was estimated as : ÊðLT 50 Þ i ¼ À ðb 0^þ a î Þ=b 1 ð2Þ The Glimmix procedure in SAS1 was used to implement this model for the F 2 -population while the similar calculations and the corresponding survival plot of the initial experiment (the parents and the hybrids) was drawn using R [13] following the code of Luke Miller (https:// lukemiller.org/index.php/archive/).

Mapping population development
The cross 'NCGR1363' (susceptible to low temperature stress) × 'Alta' (tolerant to low temperature stress) was performed at the NIBIO outstation in Kvithamar in Stjørdal, Norway. The hybrid nature of the plants were confirmed using a set of microsatellites [14]. A confirmed hybrid F 1 plant from the cross was self-pollinated and a population of 143 segregating F 2 seedlings ('NCGR1363×Alta') was raised for low-temperature stress tolerance phenotyping as detailed above and subsequent genetic map construction.

Genotyping by sequencing (GBS) and SNP calling
DNA was extracted from emerging fresh leaf tissue of 12-week old plants of the selfed 'NCGR1363' × 'Alta' mapping population, the F 1 parental plant and the grandparental plants 'NCGR1363' and 'Alta' using the DNeasy Plant Minikit (Qiagen) and sample quality was determined using a QiAgility spectrophotometer (Qiagen). Samples were considered suitable for genotyping if they had a 260/280 ratio in the 1.8 to 2.0 range. DNA quantification was performed with a Qubit fluorometer against manufacturer-supplied standards (Thermo Scientific) and was normalised to 10 ng/ul. Genotyping data were generated from the grandparents, the F 1 parent and the 143 progeny of the mapping population following the 'Genotyping-by-Sequencing' (GBS) protocol described by Elshire et al. [12]. Briefly, DNA was digested with the enzyme ApeKI and multiplexed fragment libraries were sequenced on an Illumina NextSeq 500 v2 sequencing machine, generating, on average, 1.5 million 75 bp single reads from each sample.
Demultiplexed raw reads from each sample were quality trimmed and aligned to the F. vesca v4.0 genome sequence using BWA-MEM version 0.7.12 [15] to create BAM files from which SNP variants were called using FreeBayes v1.0.2-16 [16] using the following specific parameters (-min-base-quality 10-min-supporting-allele-qsum 10-read-mismatch-limit 3 -min-coverage 5-no-indels-min-alternate-count 4-exclude-unobserved-genotypesgenotype-qualities-ploidy 2 or 3-no-mnps-no-complex-mismatch-base-quality-threshold 10). Filtering of variants was performed with a GBS-specific rule set where read counts for a locus must exceed 8, minimum allele frequency across all samples must exceed 5% and genotypes must have been observed in at least 66% of samples.

Linkage map construction and quantitative trait loci analysis (QTL) analysis
The resultant SNP data were used for mapping linkage map construction using JOINMAP 4.1 (Kyasma, NL). Following grouping, initial marker placement was determined using Maximum Likelihood with a minimum logarithm of odds (LOD) score threshold of 3.0, a recombination fraction threshold of 0.35, ripple value of 1.0, jump threshold of 3.0 and a triplet threshold of 5.0, and mapping distances were calculated using the Kosambi mapping function to produce individual linkage groups. Imputation was then performed following the protocol described by [17] and a second round of mapping using the parameters described above was performed to produce the final linkage map of the selfed 'NCGR1363' × 'Alta' mapping progeny. The linkage map presented was plotted using MapChart 2.3 [18], and the physical positions of the SNP sequence tags on the F. vesca genome sequence were used to plot MareyMaps of genetic vs. physical position of all mapped genetic markers.
QTL analysis was performed for the LT 50 phenotypic data for the mapping progeny using MAPQTL 6.0 (Kyazma, NL). The non-parametric Kruskal-Wallis test was used to identify significant associations between markers and traits individually, subsequently followed by interval mapping employing a step size of 1.0 cM, and the percentage phenotypic variance explained and associated LOD values were calculated. A LOD significance threshold of 3.2 was calculated following a permutation test with 10,000 reps, and was used to determine significance. The calculated LOD values were plotted with MapChart 2.3 [18] using the chart function.

Functional variant identification and candidate gene analysis
Gene sequences for the predicted genes within the interval spanning the QTL identified on Fvb2 of the F. vesca v4.0 genome sequence were extracted from the sequence data repository on the Genome Database for Rosaceae [19] and functionally annotated using OmicsBox (https:// www.biobam.com/omicsbox/) running default parameters. Candidate genes within the QTL interval were identified based on the relevance to cold tolerance of their functional annotation.

RNASeq analysis of QTL interval genes
Crown and leaf tissue from parental genotypes 'Alta' and 'NCGR1363' exposed to cold temperature treatment (2˚C for 42 d) and untreated plants (0 h controls) were ground to a powder under liquid nitrogen and RNA was extracted using the Spectrum Plant Total RNA Kit (Sigma) according to manufacturer's instructions. This cold treatment corresponds precisely to the acclimation conditions employed for all plants prior to low temperature stress testing as described above. Twenty-four libraries (two tissues types × two parental genotypes × two timepoints [0 hours control and 42 days] × three biological replicates) were prepared from RNA samples with a RIN (RNA integrity number) above 8 using the strand-specific TruSeq™ RNAseq library (Illumina). Paired-end 150 bp read sequencing from the libraries was performed over three lanes of the Illumina HiSeq4000 sequencing platform at the Norwegian Sequencing Centre, at the University of Oslo, Norway.
The fastq files generated were analyzed using FastQC (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/) and TrimGalore (https://www.bioinformatics.babraham.ac. uk/projects/trim_galore/) was used for trimming of adapter sequences and low quality bases. Using the OmicsBox platform, trimmed reads were mapped to the predicted mRNA sequences from the F. vesca v.4.a2 reference genome annotation on the Genome Database for Rosaceae [19] using RSEM [20] and Bowtie2 [21]. All RNASeq data have been deposited in the NCBI database under accession number GSE166374.
Differentially expressed genes were identified using the edgeR (version 3.28.0;Robinson, McCarthy [22]) Bioconductor package that is integrated in the OmicsBox platform utilizing the following parameters: A false discovery rate (FDR) cut off of 0.05; generalized linear model (GLM) quasi-likelihood F-test; counts per million reads (CPM) cutoff of 0.5 in a minimum of 2 of 3 biological replicates; and sample normalization based on TMM (weighted trimmed mean of M-values) as recommended by the package user guide. Differentially expressed genes were then evaluated according to their functional annotation and those with a potential role in freezing tolerance as well as others observed to be highly differentially expressed were considered as potential candidates contributing to the identified QTL.

Freezing tolerance phenotyping
The estimates of the temperatures at which 50% of the cohort of clones of a given genotype survived (LT 50 ) along with their standard errors were -12.3 ± 0.25˚C for 'Alta' and -7.9 ± 0.28˚C for 'NCGR1363' which agreed with the results previously reported by Davik et al. [7]. The F 1 hybrid, 'NCGR1363×Alta' had an LT 50 estimate of -10.7 ± 0.38˚C (Fig 1) and estimated LT 50 values were calculated for 141 of the 143 F 2 progeny of the selfed 'NCGR1363' × 'Alta' mapping population (S1 Table). Data were plotted as frequency histograms (Fig 2) and the LT 50 estimates were used as phenotypes for quantitative trait locus detection.

Genotyping and linkage map construction
Four of the genotyped seedlings revealed genotyping data that was either of poor quality or suggested contamination in the sequencing due to an abundance of heterozygous genotypes. As such, theses genotypes were removed from further analysis and subsequent genetic analysis was performed with data from the remaining 139 seedlings (S2 Table). A total of 16,551 putatively polymorphic sequence variants were identified between the grandparental genotypes that were heterozygous in the F 1 parent of the selfed 'NCGR1363' × 'Alta' mapping progeny after data were analysed with the criteria described in the materials and methods. Of these, 3,294 clustered into one of seven discrete linkage groups corresponding to the seven

PLOS ONE
pseudomolecules of the F. vesca v4.0 genome sequence following an initial round of linkage mapping. Following imputation and map construction, the seven resolved linkage groups contained a total of 2,918 markers at 851 loci and spanned a total genetic distance of 593.7 cM (Fig 3; S2 Table), equating to a total physical distance on the F. vesca v4.0 genome sequence of 217.1 Mb. Linkage group 6 (Fvb3) was the longest, spanning 108.5 cM (36.7 Mb), whilst LG1 (Fvb1) was the shortest, spanning 57.

Quantitative trait loci analysis, functional variant identification, and candidate gene analysis
Kruskal-Wallis analysis for the seedlings (n = 134) for which a phenotype and a genotype were available suggested a QTL for freezing tolerance, and revealed significant marker trait associations on linkage group Fvb2 of the selfed 'NCGR1363' × 'Alta' mapping population, with both marker Fvb2_15730261 and Fvb2_10601614 significant to P�0.001 (test statistic 13.279 and 13.846 respectively). Following interval mapping implemented in MAPQTL 6.0 (Kyazma, NL), a single significant QTL was identified on Fvb2 with a peak LOD score of 3.73 explaining 11.5% of the observed trait variance (Fig 5). The most significant associations were with three

Differential expression of candidate genes
Differential expression analysis of the 896 QTL-associated candidate genes in the leaves and crowns from 'Alta' and 'NCGR1363' parental genotypes treated for 0 H or 42 D at 2˚C revealed genotype-specific changes (fold change, FC) in transcript accumulation in response to low temperature treatment as well as expression differences (fold difference, FD) between genotypes prior to-(i.e. basal transcript levels), and after treatment ( Table 1). Fourteen of the candidate genes identified in the QTL interval that were putatively associated with plant response to temperature and/or osmotic stresses, exhibited parental genotype differences in transcript accumulation either prior to treatment or in response to cold treatment, sometimes both (S4 Table). Candidate genes showing genotype-specific cold-responsive expression differences primarily in crown tissue were DRP3A (FvH4_2g11510), SWEET1-like (FvH4_2g14860.1), hAT (FvH4_2g12511.1), BYPASS1-like (FvH4_2g13680.1), ASMT (FvH4_2g15840.1) and NAC029 (FvH4_2g16180.1). The hAT, BYPASS1-like, ASMT, and NAC029 homologues exhibited a greater cold-responsive down regulation in 'Alta' crowns when compared to 'NCGR1363,' while SWEET1-like showed a greater cold-induced increase in 'Alta' when compared to 'NCGR1363', and the t3 isoform of FvH4_2g11510.1 (DRP3A) displayed a FD of greater than -1600 in cold-treated crowns of 'NCGR1363' when compared to 'Alta.' Additional differentially-expressed genes were identified with a role in abiotic plant stress, including those encoding a vacuolar iron transporter (FvH4_2g13960.1), an acyl-[acyl-carrier-protein] (FvH4_2g14690.1), two MAD3/HGMR homologs (FvH4_2g17051.1 and FvH4_2g17060.1), a  Table, section B). With the exception of DRP3A, the directions of LTS-induced FC effects appear to nearly offset disparate basal expression levels between genotypes.
Candidate genes displaying genotype-specific cold-responsive expression differences primarily in leaf tissue, and basal pre-LT-treatment levels in crowns were a PIP2;7 homolog

Discussion
One of the aims of a previous investigation was to identify parental genotypes that would be suitable for studying the genetics of freezing tolerance in F. vesca [7]. Here, the genotypes displaying the greatest differences in freezing tolerance were crossed and from the selfing of the resultant F 1 progeny, a mapping population segregating for freezing tolerance was raised and phenotyped. A freezing tolerance QTL identified on chromosome Fvb2 of the F. vesca genome spanned a physical interval of 5.1 Mb and contained 896 gene predictions encoding proteins including ADH1, ERD10, PIP2 aquaporins, ascorbate oxidase, a hAT dimerization domaincontaining protein, CTR1 kinase, B1L, ASMT, EXPLA2, transcription factor RDUF2, a ninjafamily AFP2-like transcription factor and two NAC transcription factors, all of which have been reported to have putative roles in plant temperature or osmotic stress response and, for some, freezing tolerance. Here we discuss the candidate genes identified in the context of freezing tolerance in the selfed 'NCGR1363' × 'Alta' F. vesca mapping population.

Alcohol dehydrogenase
An increase in the production of C1 to C9 alcohols in plants enhances membrane fluidity which prevents phase transition occurring in plant cell membranes, thus promoting greater tolerance to freezing stress in plants [23]. Alcohol dehydrogenase (ADH) genes play a role in the production of C1 to C9 alcohols and have been shown to be cold-induced in Arabidopsis and cereal crops [24]. More recently, Song [25] reported that in Arabidopsis, ADH1 was significantly upregulated in response to cold treatment and the ADH1 knockout mutants they screened showed lower basal freezing tolerance than wild-type plants, and a higher percentage of ion leakage after freezing treatment, suggesting a pivotal role for ADH1 in the protection of plasma membranes and thus in freezing stress tolerance. An ADH1 homologue was the first protein-encoding gene to be completely sequenced in cultivated strawberry [26], was the first gene sequence to be genetically mapped in Fragaria [27], and has been used to infer phylogenetic relationships in the genus [28]. Koehler et al. [29] reported a strong correlation between ADH levels and cold tolerance in the cultivated strawberry F. ×ananassa and similarly, Davik et al. [7] demonstrated that LT 50 was strongly correlated (r = -0.86) with ADH protein levels. Davik et al. [7] also reported that ADH levels were very low in F. vesca control crowns, but strongly induced in cold-treated crowns, with up to a 200-fold increase in ADH protein levels observed after 42 days of cold treatment in accessions that were shown to be highly tolerant to freezing stress. The authors concluded that ADH likely contributes to cold hardiness in F. vesca.
The ADH1 homologue first mapped by Davis and Yu [27] is located at 12,948,939 bp on chromosome Fvb2 of the F. vesca genome, placing it within the mapping interval of the QTL identified in this investigation for freezing tolerance in F. vesca 'Alta.' Both ADH homologs were expressed at higher basal levels in 'NCGR1363' crowns but appeared to show a lower cold-induced increase, especially for FvH4_2g14760.1, in leaves in comparison with 'Alta'. Neither of these homologs displayed the cold-induced increases at the transcript level in crowns that was previously observed for immunoreactive ADH proteins [7]. The higher expression of ADH transcripts in leaves of 'Alta' however correlates well with its lower LT 50 . FvH4_2g16030.1, FH4_2g09610.1) The expression of dehydrins has previously been reported to be highly correlated to cold-stress tolerance in cultivated strawberry [30]. More recently, Koehler et al [29] performed gene expression and proteomic profiling of the commercial cultivars 'Jonsok' and 'Frida' following cold exposure and demonstrated that the transcript levels of two dehydrin-like genes, a COR47-like and a XERO2-like gene were strongly correlated with cold stress. The authors speculated that the strong increase in observed levels of a dehydrin protein identified through a one-dimensional electrophoresis western-blot that used an anti-K peptide diagnostic for dehydrin was the XERO2-like dehydrin. In a proteomics study of F. vesca, dehydrin accumulation was observed following 14 days of cold treatment, with higher levels of seven distinct dehydrins accumulating after 42 days cold treatment [7]. Further, an examination of the natural variation in cold/freezing tolerance in diploid Fragaria genotypes showed a strong correlation of plant survival with the expression of total dehydrins [7]. However, due to nonspecificity of Arabidopsis dehydrin antibodies, the authors were unable to determine which specific dehydrins accumulated in the study.

Dehydrins (
The physical region spanning the QTL identified in this investigation contains a predicted dehydrin gene (FvH4_2g16030.1) with strong homology to the acidic class of dehydrins (exemplified by the Arabidopsis ERD10, ERD14, and COR47). The ERD (early responsive to dehydration stress) genes, while first identified as rapidly upregulated in response to dehydration in Arabidopsis, were subsequently observed to be upregulated by cold and other stresses [31]. ERD10, an ABA-dependent dehydrin, was characterised in Brassica napus [32] where it was shown to be induced in leaf tissue by cold stress, and subsequently it was also reported to be induced in response to cold stress in Arabidopsis [33,34]. Strong evidence in strawberry [30] and other plants demonstrate that over-or trans-expression of dehydrins increase cold and other stress tolerances [35][36][37][38]. However, the ERD-encoding homolog FvH4_2g16030.1 did not show expression in the material utilized for RNASeq analysis of QTL candidate genes, and further characterisation of gene expression for these genes was not performed in this investigation. FvH4_2g15440.1, FvH4_2g15450.

1)
Aquaporins are a highly conserved group of membrane proteins which help transport water across biological membranes and are known as major intrinsic proteins. The plasma-membrane intrinsic proteins (PIPs) are a class of aquaporins that are highly responsive to environmental stimuli and have roles in various physiological functions including response to drought stress [39]. The PIP gene family, comprising 13 genes in Arabidopsis thaliana have been shown to be expressed under various abiotic conditions including drought, cold, and high salinity stress, as well as abscisic acid (ABA) treatment [40]. In the study of [40], PIP2;5 was shown to be up-regulated by cold stress, while most of the other members of the family were down-regulated. Similarly, in a proteomics study of cold stress in banana species (Musa spp. 'Dajiao' and 'Cavendish'), the abundance of aquaporins significantly increased after 3 hours of cold stress in 'Dajiao' seedlings [41] and the authors concluded that the aquaporins MaPIP1;1, MaPIP1;2, MaPIP2;4, MaPIP2;6, and MaTIP1;3 were all involved in decreasing lipid peroxidation and maintaining leaf cell water potential in cold stressed seedlings, which were likely the cellular adaptations responsible for increased cold tolerance of 'Dajiao' over 'Cavendish' seedlings. A total of ten PIP aquaporins have previously been reported in the genome of F. vesca [39], where diurnal expression was observed in the transcript levels of three of the characterised genes. More recently, substrate-specific expression profiles were shown for aquaporins in the cultivated strawberry F. ×ananassa [42], suggesting functional specialisation of aquaporins within the same class. Two predicted genes identified within the QTL interval displayed high homology to PIP2 aquaporins and could thus play a role in freezing-stress tolerance in F. vesca. The RNASeq analyses performed here revealed the t3 mRNA isoform of FvH4_2g15440 was expressed at higher levels in pre-LT-treated crowns and displayed a greater LT stress induction in leaves of 'Alta' when compared with 'NCGR1363'. As with the ADH homologs in the QTL interval, involvement of this single PIP2;7 isoform, one of five from this gene, in contributing to freezing tolerance would be predicated with such a role being exercised primarily in leaves. FvH4_2g16180.1, FvH4_2g12690.1,  FvH4_2g13330.1, FvH4_2g13320.1) NAC transcription factors are one of the largest families of transcription factors in plants and have been implicated in enhancing tolerance to various abiotic stresses including drought, high salinity and cold, in a number of plants [43,44]. In apple (Malus pumila) a close relative of Fragaria in the Rosaceae family, a NAC transcription factor MdNAC029 was shown to be a negative modulator of the cold stress response, directly repressing the expression of two Crepeat binding factors, MdCBF1 and MdCBF4, which are regarded key regulators of the plant response to cold stress [45]. Similarly, the role of NAC transcription factors in cold-stress response was studied in Prunus mume another member of the Rosaceae family, and 113 PmNAC genes were identified and characterised [46]. Seventeen of the genes identified were highly up-regulated in stem tissue during cold temperature stress during winter. Further analysis of a subset of 15 NAC genes showed that they were up and down-regulated in response to low-temperature treatment and were suggested to be putative candidates for regulating freezing resistance in the species. Within the freezing tolerance QTL identified in this investigation, candidate genes were identified with homology to three NAC transcription factors, NAC017, shown to negatively regulate drought-stress responses in Arabidopsis [47], NAC082, reported to be a ribosomal stress response mediator [48] and a homologue of NAC029, involved in cold-stress in apple [45] and upregulated in response to cold stress in Gossypium barbadense [49]. In this investigation, NAC029 showed a greater down-regulation in 'Alta' crowns when compared to 'NCGR1363,' however the latter showed a much lower pre-LT basal crown expression of this gene, and no significant expression differences between parents were detectable in cold-treated crowns. Despite this, the downregulation of NAC029 in response to cold is consistent with a role of this protein as a negative regulator of cold stress response in Fragaria as was observed in apple [45].

Ascorbate oxidase (FvH4_2g16000.1)
Abiotic stress induces excess reactive oxygen species (ROS) which cause oxidative stress in plants resulting in damage to lipids, DNA, RNA and proteins. ROS detoxification systems are needed to protect plant cells against the toxic effect of these species [50,51]. The ascorbate/ glutathione pathway can ameliorate the oxidative stress. Ascorbate redox status in the cell wall is regulated by the apoplastic ascorbate oxidase (AO) where it catalyses oxidation of ascorbate to monodehydroascorbate (MDHA). The short-lived MDHA may then be reduced by a membrane-associated cytochrome B or disproportionate to ascorbate and dehydroascorbate (DHA). The increased transport of DHA into the cell would be expected to lead to an alteration of the overall redox status of ascorbate, decreasing its ability to provide antioxidative support. This possibility is consistent with the observation that AO-deficient RNAi antisense mutants are more tolerant to salt and oxidative stresses than WT while overexpressing plants are susceptible to these treatments [52][53][54]. 'NCGR1363', with an LT 50 of -7.9˚C, exhibited much lower non-LT (-11.1 FD) and LT levels (-6.5 FD) of AO crown expression compared to 'Alta' (LT 50 of -12.3˚C), contrasting with these previously observed effects relating AO expression to tolerance to salt and osmotic stresses.

Other QTL-related candidate genes showing interparental differential expression
In addition to candidate genes with an identifiable role in freezing tolerance from previous studies, several genes within the QTL interval, whilst lacking clear association with low temperature stress tolerance in the Rosaceae, have been previously connected with temperature and/or osmotic stress response in other plants species.

RDUF2 homolog (FvH4_2g16170.1)
RDUF2 is an E3 ubiquitin ligase whose expression in Arabidopsis was shown to be enhanced by salt, drought and ABA-treatment, and a knock-out mutant of this AtRDUF2 exhibited markedly reduced tolerance to drought stress [55]. RDUF2 is likely part of ABA-mediated positive regulation of drought responses in plants. In both pre-and post-LT treated crowns and leaves, 'NCGR1363' expressed the RDUF2 homolog (FvH4_2g16170.1) at levels 200-to nearly 900-fold lower than in 'Alta', making it a strong candidate gene in the identified QTL interval.

CONSTITUTIVE TRIPLE RESPONSE1 (CTR1; FvH4_2g15800.1)
The CONSTITUTIVE TRIPLE RESPONSE1 (CTR1), a Raf-like Ser/Thr protein kinase, is a negative regulator that inhibits ethylene signal transduction [56,57] which functions as an essential upstream positive regulator of EIN3 in ethylene signalling [58]. Shi et al. [59] demonstrated that both a ctr1 mutant and an EIN3-over-expressing line displayed enhanced freezing tolerance in Arabidopsis. While the CTR1 homologue FvH4_2g15800.1 showed no coldinduced expression changes in parental crown tissue, the basal level in 'NCGR1363' was over 5-fold lower than 'Alta.' The data presented here showed a cold-induced fold change increase in CTR1 transcript accumulation in 'Alta' leaves while changes in 'NCGR1363' leaves were not significant.

EXPANSIN-like A2 (EXLA2; FvH4_2g16110.1)
The Arabidopsis EXPANSIN-like A2 (EXLA2) gene was first characterised by its regulation and role in responses to biotic stress, namely infections with the necrotrophic pathogen Botrytis cinerea, Pseudomonas syringae pv. tomato, and the necrotrophic fungus Alternaria brassicicola [60]. Expansins cause loosening and extension of the cell wall, possibly by disruption of noncovalent bonding between cellulose microfibrils and matrix glucans [61]. The exla2 mutant described by Abuqamar et al. [60] exhibited hypersensitivity to salt and cold stress.
The exla2 homologue in this study (FvH4_2g16110.1) displayed lower cold-responsive changes in transcript accumulation in both crowns and leaves of 'NCGR1363' compared with 'Alta.'

BYPASS1-like (FvH4_2g13680.1)
BYPASS1-like is a DUF793 family protein rapidly induced under cold treatment in Arabidopsis and is thought to enhance freezing tolerance in plants through stabilizing CBF3 and ensuring normal CBF and CBF target gene expression [62]. While 'NCGR1363' exhibited a 6.4-fold lower basal expression of a B1L homolog (FvH4_2g13680.1) in crowns compared to 'Alta,' transcript levels did not change in response to 42 d cold treatment, whereas 'Alta' reduced B1L transcript levels 7.4-fold to levels closely matching 'NCGR1363' expression following cold treatment.

SWEET1-like homologue (FvH4_2g14860.1)
Expression of the SWEET1-like homologue FvH4_2g14860.1 was 9-fold lower in crown tissue of cold-treated 'NCGR1363' than in 'Alta.' SWEET proteins are a family of oligomeric sugar transporters in plants, and in Arabidopsis, disruption of AtSWEET11 and AtSWEET12, normally down-regulated in response to cold stress, display increased freezing tolerance in an AtS-WEET11 AtSWEET12 double mutant [63]. It is conceivable that the higher cold-responsive expression of SWEET1-like in 'Alta' crowns therefore could sequester the Fragaria homolog of SWEET11 in a complex to further lower the functional levels of this protein during cold stress.

hAT dimerization domain-containing protein (FvH4_2g12511.1)
The hAT transposon superfamily encodes transposase proteins harbouring dimerisation domains [64]. The gene FvH4_2g12511.1 in F. vesca encodes a homolog hAT dimerization domain-containing protein. Transcript accumulation from this gene was over 300-fold lower in cold-treated crown tissue in 'Alta' compared with the less cold-tolerant 'NCGR1363'. Although it is not known whether this gene harbours functions unrelated to transposon activity, it has been shown that conserved genes derived from transposable elements are associated with abiotic stress phenotypes in Arabidopsis [65]. If it functions to increase transposition, elevated expression in 'NCGR1363' in response to cold stress would likely be detrimental.

AFP2-like homolog (FvH4_2g18440.1)
AFP2 is one of four members of a family of ABI FIVE binding proteins in Arabidopsis. Knockdown afp2-1 mutant plants were shown to be hypersensitive to salt, glucose and osmotic stress, but only mildly hypersensitive to ABA [66]. In addition to induction of stomatal closure and tolerance of drought, and salt stress, vegetative responses to ABA include cold stress tolerance (reviewed in Leung and Giraudat [67]). Significantly lower expression of the ninja-family AFP2-like homolog (FvH4_2g18440.1) in 'NCGR1363' leaves compared with 'Alta' suggests a possible role in the freezing tolerance observed here.

ASMT homolog (FvH4_2g15840.1)
Phyto-melatonin, synthesized by ASMT, is postulated to mediate plant stress responses by counteracting stress-induced ROS [68]. Direct evidence of the cold-tolerance promoting properties of melatonin, also in the Rosaceae, stem from effects of exogenous application (e.g., Gao, Lu [69]). The ASMT homolog FvH4_2g15840.1 contained in the QTL interval displayed a 2.5 lower basal expression in 'NCGR1363' crowns compared to 'Alta', and both exhibited a downregulation of ASMT transcripts in cold-treated crowns, the extent of which was over 100 times greater in 'Alta'. LT-induced ASMT downregulation appears inconsistent for melatonin playing a role in enhanced cold tolerance in 'Alta', but may indicate that higher basal, pre-LT expression is the mechanism through which melatonin prepares plants for improved tolerance to low temperatures; certainly, the effects of melatonin on cold tolerance in the Rosaceae [69] and Arabidopsis [70] are based solely on pre-chilling treatments.

Conclusions
Freezing tolerance is a quantitative complex trait with numerous genetic factors and a strong environmental component contributing to its expression. In this investigation, we identified a significant QTL that explained 10.4% of the phenotypic variance observed in a F. vesca mapping population which was located in a wide physical interval on Fvb2 of the F. vesca v4.0 genome sequence. The physical interval was relatively large, spanning 5.1 Mb, and gene expression studies of the crowns and leaves of parental cultivars during cold-stress highlighted several potential candidate genes within the interval that could be responsible for the variation observed in freezing tolerance of the selfed 'NCGR1363' × 'Alta' progeny.
Significant interparental differential expression levels of several of the genes reported here, along with previous evidence for roles for many of them in cold-and freezing-temperature stress responses, identified them as good candidates for the control of the effects of freezing tolerance at the QTL identified in this investigation. In order to determine the causal genetic factor for the freezing tolerance observed, further functional annotation and characterisation of the candidate genes identified will need to be performed, including the identification of causal genetic variants in the grand-parental, parental and progeny lines of the selfed 'NCGR1363' × 'Alta' mapping population and additional studies of candidate genes expression at further time-points during challenge with low-temperature stress which were beyond the scope of this current investigation. A greater knowledge of the genetic elements influencing tolerance to low-temperature stress and freezing could help develop new strawberry varieties adapted to growing environments at higher latitudes and capable of surviving in extreme winter conditions in years with no snow cover.
Supporting information S1 Table. The