Genome wide association study of frost tolerance in wheat

Winter wheat growing areas in the Northern hemisphere are regularly exposed to heavy frost. Due to the negative impact on yield, the identification of genetic factors controlling frost tolerance (FroT) and development of tools for breeding is of prime importance. Here, we detected QTL associated with FroT by genome wide association studies (GWAS) using a diverse panel of 276 winter wheat genotypes that was phenotyped at five locations in Germany and Russia in three years. The panel was genotyped using the 90 K iSelect array and SNPs in FroT candidate genes. In total, 17,566 SNPs were used for GWAS resulting in the identification of 53 markers significantly associated (LOD ≥ 4) to FroT, corresponding to 23 QTL regions located on 11 chromosomes (1A, 1B, 2A, 2B, 2D, 3A, 3D, 4A, 5A, 5B and 7D). The strongest QTL effect confirmed the importance of chromosome 5A for FroT. In addition, to our best knowledge, eight FroT QTLs were discovered for the first time in this study comprising one QTL on chromosomes 3A, 3D, 4A, 7D and two on chromosomes 1B and 2D. Identification of novel FroT candidate genes will help to better understand the FroT mechanism in wheat and to develop more effective combating strategies.

Bread wheat (Triticum aestivum L.) is an allohexaploid (2n = 6x = 42, AABBDD) species derived from two hybridization events in the region of the Near Eastern Fertile Crescent 1 . The spread of domesticated hexaploid wheat from the Near Eastern Fertile Crescent to today's growing regions required phenological adaption to different environments, e.g., selection of spring/winter types or phenotypes with reduced photoperiod sensitivity 2, 3 . In general, three different types of bread wheat adapted to specific environments according to the vernalization requirement, which is necessary for the transition from the vegetative to the generative stage, are known (spring, facultative and winter type ) 4 . The exposure to cold temperature needed for the vernalization response varies between 2-4 weeks and 4-6 weeks in semi and strong winter types, respectively 5 . Winter wheat types are higher yielding compared to spring wheat types, because they can take advantages of autumn rainfall 6 and have a longer vegetation period 7 . Today, wheat is cultivated on around 220 M ha worldwide, which resulted in an annual production of 757 M tones in 2018 8 . It is estimated that approximately one-third of the wheat growing area is cultivated with winter or facultative wheat types. Winter hardy wheat varieties are mostly needed in the Great Plains of North America, the Russian Federation, as well as Eastern Turkey, Northwestern Iran and China 9 .
Low temperatures significantly reduce yield performance in Australia 10 , Europe 11 and United States 12 . Frost damage is observed when sensitive tissue of plants is faced with low temperature in different growth stages 13 . Economic damage of frost events on crop performance depends on the time point of occurrence 14 . Due to the cold acclimation phenomenon, winter cereals survive frost by regulating their metabolism at low temperatures and protecting critical structures of cells against freezing temperatures 15,16 .
In order to reduce the negative effects of frost on crop production, it is necessary to identify genes or genomic regions involved in FroT 17 . The mechanism of plants that describes the response to low temperature by increasing the freezing tolerance is called cold acclimation. Some physiological and biochemical changes occur during cold Scientific Reports | (2022) 12:5275 | https://doi.org/10.1038/s41598-022-08706-y www.nature.com/scientificreports/ acclimation, e.g., soluble sugars, proline and cold-resistance proteins are synthesized to protect proteins at the physiological level 18 . These substances play a role in cold stress response of plants through the regulation of the osmotic potential, ice crystal formation, stability of cell membranes and reactive oxygen species. Several components encompassing messenger molecules, protein kinases, phosphatases and transcription factors assumed to be involved in cold-stress signaling pathways have been reported during the last decades 19 . Two important FroT loci namely FROST RESISTANCE 1 (FR-A1) and FROST RESISTANCE 2 (FR-A2) are located on the long arm of wheat chromosome 5A. These loci influence freezing tolerance and winter hardiness. The first locus is closely located to the VRN-A1 gene, but no information on the effect of this gene on the FR1 locus is known. Vrn-B1 and Vrn-D1 genes are mapped on the long arm of chromosome 5B and 5D, respectively. VRN-A1 has a major effect in the determination of spring/winter habit 20 and plays a major role in FroT 15 . Regarding other genes involved in FroT, the ICE (inducer of CBF expression)-CBF (C-repeat binding factor)-COR (cold-responsive or cold-regulated) pathway has been known as the main cold signaling pathway in many plant species [21][22][23] . Under low temperatures, DELLA releases ICE1 from its JASMONATE ZIM-DOMAIN (JAZs) enabling the induction of CBF genes, which are members of the AP2/ERF (APETALA2/ETHLENE RESPONSIVE FACTOR) family 24 . CBF genes bind to C-repeat/ dehydration-responsive elements (CRT/DRE) and regulate the expression of cold-responsive/late embryogenesis-abundant (COR/LEA) genes 25 .
In addition to the ICE-CBF-COR pathway, vernalization genes e.g., VRN1, VRN2 and VRN3 respond to low temperature through the flowering pathway. Therefore, changes in the regulatory regions of vernalization genes cause a delay in flowering time in plants 26 . For instance, VRN1 reduces FroT by decreasing the transcript level of CBF and COR genes 27 . Therefore, a delay in flowering time increases FroT, which indicates a connection between the flowering and cold response pathway by the interaction of VRN1 with CBF and COR genes.
Genome wide association studies (GWAS) are widely applied in many crop plants to identify quantitative trait loci (QTLs) associated with traits of interest [28][29][30][31] . Development of high-throughput single nucleotide polymorphism (SNP) genotyping platforms, e.g., Illumina 32 and Affymetrix 33 enabled the conduction of GWAS in plants and became a useful approach to detect QTL and allelic variation for complex traits 34,35 . GWAS was successfully applied in wheat to identify QTL regions associated with abiotic stress tolerance (e.g. 36 ), yield components (e.g. [37,387]), grain quality (e.g. 39 ) or diseases resistance (e.g. 40 ). However, up to now, only a few studies were published dealing with FroT candidate identification by GWAS 41 .Nevertheless, during the last decade, several FroT QTL mapping analyses using bi-parental populations were conducted and QTL on different wheat chromosomes (with the exception of chromosome 4D) were identified [42][43][44][45][46] . However, the majority of these QTL regions was identified by bi-parental QTL studies. Vagujfalvi et al. 15 pointed out that ten wheat chromosomes were assumed to be involved in the regulatory gene networks associated with FroT. However, until now, the majority of genes which are assumed to be involved in FroT have been identified on chromosomes 5A, 5B and 5D 20,25,47 .
Therefore, the aims of the present study were (i) to conduct GWAS to identify genome regions associated with FroT, (ii) to investigate potential candidate genes from QTL regions using the wheat reference genome and (iii) to compare our results with previously published FroT regions and genes in winter wheat.

Material and methods
Phenotypic data. A panel of 276 bread wheat genotypes from 31 countries was evaluated for FroT. This panel comprised 83, 4, 143 and 46 genotypes from Asia, Australia, Europe and USA, respectively (Supplementary Table S1). Out of these 216, 48, and 12 were cultivars, breeding lines and doubled haploid lines, respectively. In a previous study, Babben et al. 25 used 235 out of the 276 genotypes under investigation to identify polymorphism in known FroT genes and to conduct a candidate gene based association study (CGAS).
All genotypes were tested in four environments during 2012 and 2013 (Gatersleben, Germany; Ranzin, Germany; Puskin, Russia; Roshchinskiy, Russia) and one environment in 2012 and 2014 (Novosibirsk, Russia), according to Babben et al. 25 . Genotypes were tested in a random design in double rows and two replications per genotype. However, in Roshchinskiy, genotypes were tested as a miniplot (2.5m 2 ) trial with only one replication. FroT was assessed as percentage winter survival of plants per plot for each genotype after winter (0% = all plants died, 100% = no plant died; for further information see Babben et al. 25 ). The quality check of phenotypic data was done as described by Babben et al. 25 . Least Square means (LSmeans) per genotype were estimated by fitting a mixed linear model in SAS 9.4 48 (for further information please see Babben et al. 25 ).

Genome wide association studies (GWAS). Genotyping of the 276 wheat genotypes was conducted at
Trait Genetics, Gatersleben (Germany), by using the 90 K iSelect array (Illumina Inc., San Diego, USA). Flanking sequences of the 90 K array were mapped against the reference genome of Chinese Spring RefSeqv1.0 49 . All mapped markers with more than 30% missing values were excluded. The resulted marker set was imputed by using the Beagle 4.1 software 50 . Next, the imputed marker data set was filtered for minor allele frequency (MAF) ≥ 3% and heterozygosity ≤ 12.5%. The filtered marker data set was combined with 182 SNP markers indicating polymorphisms in 15 candidate genes for FroT known from a previously published study 25 . This final marker data set consisting of 17,566 SNP markers was used for LD decay and GWAS.
To estimate linkage disequilibrium (LD) and to determine LD decay the software package R was used (R Core Team. 2014, packages "genetics" and "LDheatmap" 51,52 ). The LD was estimated as squared allelic correlation (r 2 ) between all pairs of markers within a chromosome. For graphical display, the genetic distances between markers in base pairs were plotted against the estimated r 2 . The critical value of r 2 was set to r 2 = 0.2 as described by Voss-Fels et al. 53 . Furthermore, a smooth locally weighted polynomial regression (LOESS) curve was fitted to calculate the LD decay . Finally, the LD decay was determined as intersection point of the LOESS curve and the critical r 2 value 54 . LD decay was estimated for each chromosome and across all 21 wheat chromosomes. www.nature.com/scientificreports/ In order to get comparable results with previously published analysis on the same material, a reduced marker set (249 markers, for further information see 25 ) was used to determine population structure and to calculate kinship matrix.
Kinship matrix was calculated based on Roger's distance for each pairwise genotype -genotype combination 55 . Population structure was investigated by using Bayesian cluster analysis implemented in Structure 56 and principal coordinates analysis (PCoA) implemented in the DARwin 6 software 57 . To determine the population structure by using the software package STRU CTU RE (source), ten independent runs were performed setting the number of populations (k) from 1 to 10. Furthermore, the number of burn-in and Markov Chain Monte Carlo (MCMC) iterations was set to 100,000. To determine the optimal number of subpopulations, the Evanno method (ΔK method) implemented in the software package STRU CTU TRE HARVESTER version 2.3.4 58 was used.
GWAS was conducted by using the software package TASSEL 5 59 . A compressed Mixed Linear Model (CMLM) was used to examine associations between SNP markers and FroT data. Two association models were tested: 1) Q + K model (CMLM with Q-matrix and K-matrix as correction for population structure and kinship relationship); 2) K model (CMLM with K-matrix as correction for kinship relationship). All marker trait associations with LOD ≥ 4 (-log 10 of P value) were assumed to be significantly associated with FroT according to Babben et al. 25 and Zhao et al. 41 . All markers, which were significantly associated with FroT, were assigned to QTL regions according to their chromosomal position and the estimated LD decay (3.5 million base pairs). The peak marker of each QTL region is defined by the highest LOD value. Markers within a distance of ± 3.5 million base pairs to the QTL peak marker were assigned to one QTL region. Genes located within the QTL regions were identified based on their position on the reference genome of Chinese Spring 49 . All high and low confidential (HC and LC) genes located within a QTL region were identified. Additionally, published functional gene annotations 49 were used to identify gene onthology (GO) terms associated with FroT. All GO terms associated with frost or cold tolerance were downloaded from the QuickGO website (https:// www. ebi. ac. uk/ Quick GO/). Genes within QTL regions were filtered based on GO terms associated with frost or cold tolerance (GO, Supplementary Table S2).
Identification of candidate genes via BLASTn. The sequences of associated candidate genes to FroT were used to identify gene IDs. For this purpose, a BLASTn analysis (nucleotide Basic Local Alignment Search Tool) from National Center of Biotechnology Information (NCBI; https:// www. ncbi. nlm. nih. gov/) based on the whole coding sequence (CDS) of candidate genes was used. These sequences were aligned to Triticum species with default settings. Perfect match or high similarity was identified based on 100% query coverage, an Expect (E) value of 0 and an identity higher than 99%.

Results
Genotyping of the 276 wheat accessions resulted in a raw data set of 81,587 SNP markers. In total, 54,340 markers were excluded from further analyses, due to the absence of a hit or no unique map position according to the reference genome sequence of Chinese Spring 49 . The 27,247 uniquely mapped markers were filtered for minor allele frequency (≥ 3%) and heterozygosity (≤ 12.5%), resulting in a set of 17,566 (17,384 90 K SNPs and 182 CG polymorphic sites) informative markers.
The LD decay ranged between 609,686 bp (chromosome 6B) and 8,434,298 bp (chromosome 1D). However, it was not possible to estimate LD decay for chromosome 4D. The LD decay across all chromosomes was 3,377,883 bp (Supplementary Table S3). The calculated LD decay was used to define QTL regions.
Population structure was determined by a Bayesian cluster analysis implemented in in the software package Structure and a PCoA implemented in DARwin 6. The Bayesian cluster analysis revealed an optimal number of K = 2 or K = 3 subpopulations. Due to the origin of the genotypes (North America, Asia and Europe), K = 3 was used as optimal number of subpopulations ( Figure S1). Genotypes were assigned to one of the three subpopulations based on their membership coefficients. However, genotypes with membership coefficients < 0.7 to a subpopulation were considered as admixture 60 . In total, 23, 78, 88 and 87 genotypes were assigned to subpopulation 1, 2, 3 or the admixed group. Additionally, a PCoA plot was used to visualize the results of the Bayesian cluster analysis. The first and the second PCoA explained 8.4% and 4.7% of the whole variance. Results of both analyses pointed to a weak to moderate population structure ( Figure S2).
The compressed mixed linear model with correction for kinship relatedness (CMLM with K) turned out to be the most appropriate GWAS model ( Figure S3). Therefore, GWAS was conducted based on the CMLM K model.
In total, 53 markers were found to be significantly associated with FroT at a significance threshold of LOD ≥ 4 (Supplementary Table S4). Out of these, 16 SNPs were associated with polymorphisms in a candidate gene for FroT (CBF-14 on chromosome 5A) previously published by Babben et al. 25 (Supplementary Table S4). The significantly associated markers were assigned to 23 QTL regions, which explained with a range of 6.1% to 16.2% from the total phenotypic variance. These QTLs were located on 11 chromosomes (1A, 1B, 2A, 2B, 2D, 3A, 3D, 4A, 5A, 5B and 7D) ( Table 1, Fig. 1 and Supplementary Table S4), whereby the majority of significantly associated markers was detected on chromosome 2B (13 markers, Supplementary Table S4 and Fig. 1).
All identified QTL regions were screened for known functional genes located within the QTL regions. In total, 3,112 HC and 2,004 LC genes were identified within the QTL regions (Supplementary Table S2). The number of low confidence (LC) and high confidence (HC) genes per interval ranged between 111 and 589. As mentioned above, 16 of the markers significantly associated with FroT indicate polymorphisms in the CBF-A14 gene located within the QTL region QTL_5A_2 on chromosome 5A. The CBF 14 gene is known as a putative candidate gene for FroT in wheat 25 and in barley 63 . Additionally, all LC and HC genes located within the QTL regions were screened for GO terms associated with cold stress or FroT. In total, in 23 QTLs from the present study a set of 5116 genes was annotated according to the Chinese Spring reference genome (Supplementary  Table S5). At chromosome 2A four QTLs containing 1317 Chinese Spring annotated genes, of which two HC and three LC are related to FroT, were identified. Interestingly, two of these genes code for a cold-responsive protein (WCOR15, QTL_2A_3) and the other one was a cold shock domain protein 1 (QTL_2A_2). In addition, a gene coding for a low temperature and salt responsive protein was identified within QTL_2A_1 on chromosome 2A. Furthermore, genes coding for 16 basic helix-loop-helix (bHLH) transcription factors were detected on chromosome 1B, 2A, 2B, 3A, 5A and 5B. Three genes coding for flowering locus T-like proteins (FT-like) and four genes for flowering time control proteins (FPA) were identified on chromosome 2A and 5A, respectively. An important gene involved in the FroT pathway, ICE1, was observed on chromosome 3A. In addition, 25 CBF genes have been identified. Fifteen of them were located on chromosome 5A (specifically on the locus of frost resistance A2 (FR-A2)) and 10 on chromosome 5B. Seventeen specific gene names out of these 25 CBF candidate sequences were identified via BLASTn (Supplementary Table S6). Twelve and five out of these 17 identified genes were located at the FR-A2 locus on chromosome 5A and 5B, respectively.

Discussion
Bread wheat is grown worldwide in temperate latitudes and subtropical regions 64 and constitutes the main source of proteins and calories in human diets 65 . However, many wheat-growing areas are regularly exposed to heavy low temperature events during the early stage of wheat development 9 causing severe yield losses 66 . Therefore, FroT is an important trait in breeding in order to improve winter hardiness of wheat 16 .
In the recent years, high-throughput sequencing technologies fostered the availability of large SNP data sets and therefore the conduction of population genetic and GWAS studies 32 . Furthermore, these technologies made the first fully annotated reference genome sequence of wheat available 49 . Altogether, this progress in plant genetics and genomics helps to increase the understanding of wheat biology and the molecular basis of important agronomic traits 67,68 . Based on this progress, this study aimed to identify QTL regions and candidate genes associated with FroT in wheat.
Recently, QTL mapping studies or GWAS identified several QTL regions associated with FroT in wheat on all wheat chromosomes except chromosome 4D 25, 41, 43-46, 69, 70 . In general, it is difficult to compare QTL regions identified by different studies using different marker systems and different genetic maps. Therefore, to compare the records from literature with findings of this study, known flanking sequences of markers associated with QTL for FroT in literature were mapped to the reference genome sequence of Chinese Spring 49 . However, flanking markers were not available for all published QTL regions. Hence, these QTLs could not be anchored on the reference genome sequence of Chinese Spring. Furthermore, for some flanking marker sequences no unique position on the reference genome sequence of Chinese Spring 49 could be identified. Therefore, for markers and QTL regions that could be not uniquely mapped to the reference genome, comparison with the results of this study was conducted based on the chromosome (Supplementary Table S7).
To get comparable results, identified candidate genes in the current study with previous study 25 , LOD ≥ 4 was subjected as threshold for significantly associated markers with FroT. In addition, Zhao et al. 41 reported associated significant markers with FroT with LOD ≥ 4. The P value was also adjusted by Bonferroni-Holm (LOD ≥ 5.5) correction in the present study. Finally, since only three markers were identified at Bonferroni-Holm threshold at LOD ≥ 5.5, we considered markers with LOD ≥ 4 for further analysis. In total, 53 markers were found to be significantly associated (LOD ≥ 4) with FroT in this study. These markers were assigned to 23 QTL regions on 11 chromosomes (1A, 1B, 2A, 2B, 2D, 3A, 3D, 4A, 5A, 5B and 7D).
In the present study, three out of 23 identified QTL regions namely QTL_5A_1, QTL_5A_2 and QTL_5B are overlapping with previously reported QTL regions 41,[44][45][46] . It is known that important genes associated with FroT are located on chromosome 5A, e.g., CBF genes and VRN genes. In this study, 17 CBF genes, three bHLH family transcription factors and one FPA gene were identified within QTL regions associated with FroT on chromosome 5A and 5B. Polymorphisms within two out of 17 identified CBF genes (CBF-A14 and CBF-A15) were previously detected by CGAS 25 . In total, 16 markers associated with polymorphic sites for CBF-A14 were significantly associated with FroT in the present study. www.nature.com/scientificreports/ In plants, Hormones and photoreceptors, such as phytochromes, regulate the expression of the CBF regulon as a response to the modulated spectrum of the incident light, to enhance cold acclimation and freezing tolerance in plants 71 . The positive role of PHYA on the transcription levels of CBF pathway genes is reported in tomato 72 , wheat and barley 73,74 . In the present study, PHYA was associated with a QTL for FroT on chromosome 5A. Furthermore, twelve QTL regions, i.e., QTL_1A, QTL_2A-1, QTL_2A-2, QTL_2A-3, QTL_2A-4, QTL_2B_1, QTL_2B_2, QTL_2B_3, QTL_2B_4, QTL_2B_5, QTL_2B_6 and QTL_2B_7 are potentially overlapping with previously reported QTLs on the same chromosomes 42,61,62 . However, due to the unavailability of flanking markers or due to the fact that the flanking sequences of available markers could not be anchored on the reference genome 49 ; it was not possible to confirm this assumption. Three candidate genes associated with cold or low temperature tolerance were identified on chromosome 2A. The QTL_2A_1, QTL_2A_2 and QTL_2A_3 were co-localized with genes coding for Low temperature and salt responsive protein, cold shock domain protein 1, WCOR15 and LEA 19-like protein, respectively.
CORs are referred to as proteins encoded by cold-responsive or cold-regulated genes, which are involved in the cold tolerance acquisition and subsequent freezing tolerance. These genes, i.e., LEA, stress responsive protein (SRP), cold induced (KIN) and low temperature induced (LTI) 21 , increase cold tolerance in plants. For instance, accumulation of COR/LEA proteins during cold acclimation protects cell structures and functions from freezing damage 62 . The Wcor15 is expressed under low temperature 75 and encodes a chloroplast-targeted protein in wheat and barley.
QTL_2A_4 includes a gene coding for Flowering Locus T-like protein. It has been shown that flowering time genes are not only responsible for the transition from the vegetative to the reproductive phase, but also involved in various environmental stress responses. The relation between flowering and cold response is well known 76 .
For eight of the QTLs associated with FroT in this study, we did not find any evidence that these QTL were previously reported in literature. These QTLs are located on chromosome 1B, 2D, 3A, 3D, 4D and 7D and will be discussed in the following.
One out of two identified QTL regions on chromosome 1B contains a gene encoding a bHLH transcription factor. bHLH transcription factors play diverse roles in different physiological processes 77 . Several studies have been shown that bHLH is involved in different responses, which are provoked by cold and other abiotic stresses in Arabidopsis and rice [78][79][80][81][82] . Wang et al. 83 identified 159 bHLH -encoding genes in wheat, which are involved in abiotic and biotic stress response. Furthermore, they pointed out that 98.7% of these genes are associated with more than one stress. In addition, the expression of these genes under different stresses was evaluated. In total, 38.44% of these genes were upregulated under cold stress in wheat 83 .
The identified MYC2-like transcription factor is also called JAM (JASMONATE ASSOCIATED MYC2-LIKE) and bHLH14 is called JAM1. Both transcription factors are members of the IIId bHLH subfamily, which is phylogenetically closely related to MYC proteins interacting with JAZ proteins. These bHLH subfamily acts as transcription repressor of MYC2 and so as a negative regulator of jasmonate mediated response 84,85 . Furthermore, Xiang et al. 86 described the role of cold induced transcription factor bHLH112, which promotes a positive regulation of AP2/ERF transcription factor in Artemisia annua and Jiang et al. 87 identified that bHLH35 is involved in cold tolerance in Anthurium andraeanum.
In addition, we identified one QTL region on chromosome 3A. This region comprises transcription factor ICE1 and bHLH transcription factor genes. The transcription factor ICE1 is known as an important gene involved in freezing tolerance (ICE -CBF -COR) pathway. ICE1 genes are known in wheat, but until now, no ICE1 gene was found to be located within a QTL region associated with FroT. Two ICE homologs, i.e., TaICE41 (accession no. EU562183) and TaICE87 (accession no. EU562184) have been identified in wheat 21 . The identified AP2/B3-like transcriptional factor protein (REM18) on chromosome 5A is also a member of the DREB/ERF subfamily and it is accordingly maybe involved in FroT 88,89 .
The bHLH35-like gene was identified on chromosome 3A. Less knowledge is available for this gene. Jiang et al. 87 have reported the positive role of bHLH35 in response to abiotic stresses in Arabidopsis. They reported that bHLH35 from Anthurium andraeanum (AabHLH35) increases stress tolerance to cold and drought in Arabidopsis. The expression of CBF1 and COR15A in wild type (WT) and AabHLH35 transgenic lines of Arabidopsis was significantly increased under cold stress compared to control plants. Expression of COR15A was threefold higher in AabHLH35 transgenic lines relative to WT lines. Therefore, they assumed that AabHLH35 might promote COR15A expression in response to cold stress. Furthermore, OsbHLH35 increased salinity tolerance in rice 90 and PebHLH35 from Populus euphratica increases drought tolerance in Arabidopsis 91 .
As mentioned above, several genes are involved in the enhancement of FroT in plants. In addition to the ICE-CBF-COR pathway and flowering time genes, we identified nine genes of Jasmonates (JA) in the present study, which play a major role in the ICE-CBF-COR pathway (Fig. 2) by activating transcription factors. Activated transcription factors bind to the cis-acting element in the promoter of target genes to increase FroT in plants 24 .

Conclusion
This study dealt with the identification of QTL regions and putative candidate genes associated with FroT in wheat. GWAS resulted in the identification of 23 QTL regions associated with FroT. The identified QTL regions on chromosome 5A and 5B are in accordance with known genomic regions and candidate genes previously described for FroT in wheat. Moreover, the findings reported here, confirm the results of the previous study of Babben et al. 25 in regard to polymorphisms in candidate genes for FroT (CBF-14) on chromosome 5A. To the best of our knowledge, eight of the detected QTL regions can be assumed to be novel, as these regions were not described in literature before. Furthermore, within the QTL regions on chromosome 1B, 2A, 2B, 3A, 5A and 5B genes with GO terms associated with cold stress response or FroT were identified. The findings reported here give hints to known and previously unknown genome regions and candidate genes, which are putatively Scientific Reports | (2022) 12:5275 | https://doi.org/10.1038/s41598-022-08706-y www.nature.com/scientificreports/ associated with FroT in wheat and therefore mark the starting point for further research. Prospectively, these findings will help to develop diagnostic markers for FroT in wheat and to fine map putative candidate genes. Both will foster the better understanding of FroT in wheat and the improvement of winter hardiness and FroT in bread wheat elite breeding pools.