Non-escaping frost tolerant QTL linked genetic loci at reproductive stage in six wheat DH populations

https://doi.org/10.1016/j.cj.2021.02.015 2214-5141/ 2021 Crop Science Society of China and Institute of Crop Science, CAAS. Publishing services by Elsevier B.V. on behalf of KeAi Communications Co. L This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author. E-mail address: w.ma@murdoch.edu.au (W. Ma). 1 These authors contributed equally to this work.


Introduction
Frost can cause significant grain yield and quality losses in wheat crops [1]. In spring, when plants sense the gradual increase of temperature and their development proceeds beyond the jointing stage, both winter and spring types show considerable sensitivity to low temperatures (0-12°C) and frost (<0°C), particularly to short chilling and frost events at night [2]. The early-flowering phenotype in modern wheat and barley cultivars has resulted in significant grain yield and quality losses from frost damage [3]. Frost damage during reproductive stage can lead to a multiplicity of symptoms, including dead stems, floret and spikelet abortion, and empty shells along the spikes, thus significantly reducing the seed number per spike. In spring wheat cultivation in Australia, late sowing plants and long season varieties can escape the low temperatures and get less impact, which is called as a frost escaping mechanism. However, those plants most likely will face drought and heat stresses during the grain filling period, resulting in significant yield losses. Farmers and breeders are looking forward to frost tolerant varieties with early maturity. Therefore, short season varieties with non-escaping frost mechanism are desirable in Australia.
In a previous abiotic research, during the reproductive stage, male organs increase the sensitivity dramatically from the start of meiosis to the break-up of the tetrad. A single anther needs approximately 24 h in that event [4]. The significant grain number loss was happened during the period of 8-17 days before anthesis when exposing to drought stress [5]. Meiosis (10 days before anthesis) is the most sensible stage to abiotic stress [6]. At meiosis, male sterility occurs under non-freezing temperatures below 10°C in cereals [7,8].
Many previous wheat frost tolerance studies focused on the vegetative development stages [9][10][11][12][13][14]. Genetic segregation for vegetative frost tolerance (or susceptibility) have been reported on chromosomes 5A, 5B, 5D, and 7B in wheat, and 5H in barley. The frost tolerance QTL Fr-A1 and Fr-B1 on chromosomes 5A and 5B were closely linked with the vernalisation genes Vrn-A1a and Vrn-B1a [11,15]. Frost resistance 2 (FR2) genes are in control of delayed heading. In combination with VRN1 showed reproductive frost tolerance [16], it is not clear yet whether these frost tolerance loci are identical with Vrn-A1a and Vrn-B1a [17]. Similarly, it remains unclear whether the frost locus on chromosome 7B is influenced by VRN3 (VrnB4) [9]. In barley, three doubled haploid (DH) populations were used to identify reproductive stage frost tolerance QTL [3]. Two major QTL were identified on chromosomes 2H and 5H. The QTL for frost-induced floret sterility and grain damage overlapped with the anthesis QTL on the Vrn-H1 locus in all three DH populations, whereas in two of the populations the floret sterility QTL on 2H was not close to the growth QTL or Ppd-H1, but close to the earliness per se gene (Eps 2), the cold-regulated gene (Cor14b) and the barley low-temperature gene (Blt14) loci. It seems that the frost tolerance QTL are closely associated with the vernalisation genes in both vegetative and reproductive stages, although the effects kick in at different growth stages.
Regarding the anthesis related genes, three vernalization gene groups including VRN1, VRN2, and VRN3 have been well studied [18]. Three homologous copies of the VRN1 gene, known as Vrn-A1a, Vrn-B1a and Vrn-D1a were mapped on chromosome 5A, 5B and 5D, respectively [19][20][21][22][23][24]. The VRN2 gene was mapped on the distal region of 5AL which was a repressor of flowering. When VRN2 is down-regulated by vernalization, the gene expression of VRN1 was promoted [19,25]. The mutated and dysfunctional VRN2 resulted in spring lines [25]. VRN3, similar to the Flowering Locus T (FT) gene, encodes a Rapidly Accelerated Fibrosarcoma (RAF) kinase inhibitor-like protein [18]. The mutated lines with Vrn-B3 flowered considerably earlier than the recessive vrn-B3 allele [26]. One of the Vrn-B3 genes was mapped on 7BS (VRN-B3) [26]. The results of Yan et al. [26] implied that the VRN2 modulated the quantitative levels of FT (directly or indirectly) and the absence of VRN2 function showed no effect to VRN1 and VRN3 mutations. Lately, another gene for developing spring growth habit VRN-D4 was identified in the short arm and close to the centromeric region of chromosome 5D [27,28]. VRN-D4 is a homologous gene of Vrn-A1.
Photoperiod insensitivity provides wheat with the ability to flower in short day as well as in long day conditions. The genes (PPD1) involved in this process are Ppd-A1, Ppd-B1, and Ppd-D1 (formerly Ppd3, Ppd2, and Ppd1), located on 2A, 2B, and 2D, respectively [29][30][31]. The homologous genes Ppd-A1 and Ppd-B1 showed less effect on flowering in short days than the Ppd-D1 [29]. The insensitivity of Ppd-A1 is greater than Ppd-B1 [32], and Ppd-A1 tends to increase thousand grain weight and yield while Ppd-B1 seems associated with high kernel number through increasing spikelet number [33]. In recent study, a heading time QTL detected on 2B in durum wheat and explain 26.2% of the phenotypic variations. The early heading QTL correspond to higher copy number of Ppd-B1 [34].
The influence of earliness per se (eps) genes tend to influence developmental rate at a much lower level as compared to vernalization and photoperiod. Numerous eps and the related flowering-time QTL in wheat have been mapped to chromosomes 1DL, 2B, 3A, 4A, 4B, and 6B [35][36][37]. Eps loci are associated with spikelet number and size, thereby affecting wheat yield [38]. An ortholog to the Arabidopsis thaliana LUX ARRHYTHMO/PHYTO-CLOCK1 (LUX/PCL1) gene was identified as Eps-3A m in einkorn wheat (Triticum monococcum L.). Lines containing Eps-3A m showed a distorted circadian clock, spikelet number variation and temperature sensitivity [39]. In wheat, homologs to Arabidopsis Early flowering 3 (ELF3) gene have been identified on chromosome group 1 [40], namely TaELF3-1AL, TaELF3-1BL, and TaELF3-1DL. The gene ELF3 was identified as a candidate gene of Eps-A m 1. Wheat lines harbouring ELF3 showed flowered earlier and less spikelets per spike, and stronger photoperiod sensitivity, which indicate the significant epistatic interaction with Ppd1 [41,42]. Eps-D1 deletion reduced the total expression of TaELF3 indicating TaELF3-1DL is the major isoform of gene TaELF3 [43]. On average, lines harbouring allele TaELF3-1DLb headed two days earlier compared with those holding TaELF3-1DLa [40].
Apart from the well-known anthesis-related genes mentioned above, on 3A, 3B and 3D, there is a set of T. aestivum GIGANTEA (TaGI) encoding genes, whose products interact with FLAVIN-BINDING, KELCH REPEAT, and F-BOX 1 (FKF1) domains to form a complex regulating photoperiod-dependent flowering by regulating CONSTANS (CO) expression [44].
A SOC1 (Suppressor of Overexpression of CO 1)-like gene on chromosome 4DL, WSOC1, was reported to influence flowering time in wheat [30]. A gene for wheat vegetative to reproductive transition on the chromosome 7 group, TaVRT-2, interacts with VRN1 and VRN2 and regulates the floral transition [45]. Three short-day flowering-time genes on 1B, including flowering locus T3 (TaFT3-B1), WUSCHEL-like (TaWUSCHELL-B1) and TARGET OF EAT1 (TaTOE1-B1) have been cloned [46], with the early-flowering function for TaFT3-B1 having been validated. A set of heading-date genes (TaHD1) identified on 6A, 6B and 6D, are regulated by long-day conditions and the circadian clock, directly affecting vernalisation genes under long-day conditions. Its mutants showed a delayed flowering response in a long-day environment [47,48].
It is difficult to screen frost tolerance in the field at the reproductive stage, since trials need to be hit by a natural frost event at the right developmental stage, which is purely a matter of chance. To increase the chance, trials consisting of a wide range of genotypes usually need to be planted on multiple sowing dates [3,49]. The type of damage to the affected plants will also be influenced by the weather on days leading to the frost event and after. This goes a long way in explaining why untangling the complex genetic basis of frost tolerance under controlled conditions in glasshouses and cold chambers has been of limited utility to breeders.
Our 2018 large-scale field trials encountered such a chance event exactly at the right stage, which provides a valuable resource for QTL detection for frost damage. Six DH populations, namely, Bethlehem/Westonia (BW); Gregory/Bethlehem substitution line 7AS (G7A); Spitfire/Bethlehem substitution line 7AS (Sp7A). Spitfire/Bethlehem (SpB); Spitfire/Mace (SpM); and Suntop/Bethlehem substitution line 3BL (St3B), representing genetically divergent origins were impacted by frost in two distinct environments. All populations suffered considerable frost damage. Genes and markers contributing to the frost tolerance or susceptibility phenotype were identified. Frost tolerance segregation patterns in the current study provide valuable genetic information that can be used in wheat frost tolerance breeding.

Field and glasshouse experiments
In 2017, individual lines from the six DH populations were planted as 1-m 2 plots in Katanning (Kat), Wongan Hills (WH) and South Perth, while four DH populations, namely BW, SpM, SpB, and Sp7A were also planted in 4-litre pots in a glasshouse at Murdoch University.
In 2018, the field experiments with the six DH populations were conducted at three locations across Australia representing distinct environments including 1716 plots at Narrabri in New South Wales, 1884 plots at Muresk and 2316 plots at Williams in Western Australia. The majority of the DH lines (greater than 96%) were replicated two times at one or more than one of the three locations. Eight parental lines and three control varieties were also utilised in each of the field experiments. The partially replicated experiments were designed using DiGGer in R [50]. Since no frost were recognized in Narrabri, the Zadok data recorded in Narrabri were used to validate the flowering QTL in this study.
To identify the young microspore stage, the parental lines, together with other varieties, were grown in Muresk in 2019.
Field-grown plants were rainfed under standard agronomic management practices, whereas adequate water and fertilizer were provided to plants grown in pots in the glasshouse [51].

Growth stage measurement
Given the importance of establishing the plant growth stage during the frost events, anthesis date and plant growth stages during ear emergence were recorded for all populations in the glasshouse, in Wongan Hills and Katanning in 2017, and in Williams and Muresk during the 2018 trials (Fig. S1). The days to anthesis were calculated by the anthesis date minuses the sowing date. In the Narrabri trials in 2018, the Zadok stages around 50-60 were recorded. In Williams trials 2018, days to maturity were estimated as follows for each plot except the St3B population: on October 31, based on the plant maturity performance, the plant maturity scoring was entered as 0 (mature) to 15 (based on an estimate of 15 days to maturity). The days to maturity were the estimated maturity date minuses the sowing date (June 4, 2018).

Temperature recording
In 2018, probably because of high rainfall and low soil temperatures in August, the average anthesis time was delayed by 15 to 20 days in the WA field trials. Flowering time ranged from September 25 to October 12 in Williams, and from September 13 to October 6 in Muresk. The lowest temperatures (À1.1 and 0.3°C, respectively) were recorded from whether station beside the field trials on September 14, 15 and 16 in Williams, and September 15 and 16 in Muresk ( Fig. 1a; Table S5). The low temperatures ( 2°C) lasted 975 and 300 min in Williams and Muresk, respectively. The low temperatures during heading and anthesis caused frost damage in all six DH populations. The days between frost and anthesis were calculated by the average days to anthesis of each line minus the days to frost (September 15, 2018).
To evaluate the reliability of field scoring of frost damage, a set of randomly selected lines were assessed for the floret sterility and the data was compared with the field visual scoring results. For validating the visual scoring, ten heads were picked from each plot around 20 days after anthesis for frost assessment. Three parental lines (Bethlehem, Mace, and Westonia) and two control lines (Wyalkatchem, Yitpi) from three different anthesis time windows were used, and 7 to 14 plots for each line were randomly selected for floret sterility assessment (Fig. S3a). The two outside florets of each spikelet in a spike were used for sterility calculation. The sterility percentage equals to the sterile florets divided by the total florets in each spike, then times 100. The average sterility of the ten heads of each plot formed the whole plot frost damage level. The correlation coefficient between visual frost damage and floret sterility was 0.79 (P < 0.01) in randomly selected parental and control lines (Fig. S3b).

Plant height measurement
Plant height was measured using three main tillers of three representable plants in each plot. The plant height of G7A and St3B populations was used for QTL analysis in this study as the semidwarf genes were segregating in these two populations (Fig. S4).

Meiosis stage measurement
In 2019 field trials in Muresk, plants were sampled in 2-3 days interval since the auricle distance was 1 cm and the main tillers with the same auricle distance as sampled plants were tagged. The sampled plants were stored above ice and the spikes were dissected and photographed with a ruler. The anthers stained by acetocarmine solution (45%) were placed on a light microscopy (400 Â magnification) and the pollen developmental stages were captured. The days between meiosis and anthesis were calculated by the days to anthesis minus the meiosis date.

Linkage map construction
Genomic DNA was extracted from a single plant for each DH line and their parents [52]. SNP linkage maps were constructed for six DH populations, including two linkage maps of 90K SNP and four linkage maps of 12K Targeted Genotyping-By-Sequencing (tGBS). SNP genotyping was performed using an Infinium iSelect assay on an Illumina iScan instrument according to the manufacturer's protocols (Illumina, San Diego, CA). SNP clustering and genotype calling was performed using GenomeStudio v2011.1 software (Illumina, San Diego, CA) with the custom genotype-calling algorithm described by Cavanagh group [53]. Identical lines were detected and removed using non-metric multidimensional scaling (MDS) of genetic dissimilarity using software from Numerical Taxonomy System (NTsys) v2.2 and Plymouth Routines in Multivariate Ecological Research (PRIMER v6) [54,55]. Lines with large proportion of missing values on SNP genotyping were also removed, together with distorted markers and doublecross markers. Most co-segregating markers were made redundant and removed from the genetic map. As results, a fine-map of BW population was constructed using 77 lines and 2387 SNP markers; and so as the maps of G7A by 218 DH lines and 3592 SNP markers, Sp7A by 191 DH lines and 2367 SNP markers, SpB by 94 DH lines and 2570 SNP markers, SpM by 188 DH lines and 2235 SNP markers, and St3B by 185 DH lines and 1924 SNP markers. (Table S3).

QTL mapping
As proposed for additive, dominant and epistatic QTL mapping in biparental populations [59], we employed inclusive composite interval mapping (ICIM) for QTL detection, using IciMapping V4.1 (http://www.isbreeding.net). A LOD score of 2.5 was used as significance of QTL detection. Permutations were set to 1000 at a significant level of 0.05. The inclusive composite interval mapping addition (ICIM-ADD) method was selected for QTL mapping [60].

TaqMan assays for Ppd1-2B copy number determination
Significant anthesis and frost QTL were detected in the Ppd1-2B region. As variable copy number of the Ppd1-2B gene among parental lines was suspected, TaqMan assays were conducted to determine copy number in those lines. Based on a published protocol [61], 20 lL PCR reactions were set up including 10 lL ddPCR Supermix for Probes (Bio-Rad), 0.4 lL probe plus forward and reverse primers (10 lmol L À1 ), 5 lL DNA (10 ng lL À1 ), and 4.6 lL RNase/ DNase-free water. The primer and probe sequences for Ppd1-2B gene and TaCO2 internal control were as those published by Díaz et al. [61]. PCR cycling parameters were 95°C for 15 min; 40 cycles of 95°C for 15 s; and 60°C for 60 s. Ppd1-2B copy number was analysed based on the ratios of absolute copy numbers against the TaCO2 control.

Statistical analysis
Linear mixed models (LMM) were fitted with ASReml-R [62] in the analysis of the frost and growth stage phenotypic traits, where the variance parameters in the mixed model are estimated using the residual maximum likelihood (REML) procedure [63]. Residual diagnostics were performed to examine the validity of the model assumption (normality and homogeneity of variance). The best linear unbiased predictions (BLUPs) were used for the phenotypic traits.
Phenotypic data were analyzed by multivariate analysis of variance (MANOVA) using the general linear model implemented in IBM SPSS statistics 24 (https://www.ibm.com/au-en/products/ spss-statistics). Wilks Lambda was used as the multivariate test statistic. Post-hoc Tukey's Multiple Range tests were used to identify significant groupings. Pearson correlations of the parameters investigated were analysed by SPSS software using the BLUP values across environments. The broad-sense heritability was calculated through R-studio.

Frost impact levels associated with growth stages in DH parental and control lines
In the 2018 field trials sown at Williams and Muresk (WA), low temperatures during heading and anthesis caused frost damage to all six DH populations (Fig. 1a). The most visible frost damage was on spikes (Fig. 1b). Probably because of the different developmental stages between main spikes and tiller spikes, the levels of frost damage to spikes of the same wheat variety in the whole plot show large variation, such as ''white spikelets"-spikelets appear white and dead; ''half-cut" and ''bald-pointed" spikes-only the lower part of the ear has grains; ''toothed ears"-grains set in several spikelets. Since flowering times were different not only between field trial locations but also between sowing areas within the same location, the representative anthesis data for the parental and control lines in the areas sown to Spitfire/Bethlehem (SpB) and Bethlehem/Westonia (BW) are presented in Table S2 and Fig. 1c, d. On average, anthesis times for cultivars Gregory, Tungsten and Yitpi were significantly later than those for the other parental cultivars at both locations. The cultivar Suntop flowered two and five days later than the rest of cultivars in Williams and Muresk, respectively. The cultivar Bethlehem and its substitution lines B_3B and B_7A were the earliest flowering lines, while Westonia, Spitfire and Mace flowered two to three days later. The cultivar Wyalkatchem was one and two days later than Mace and Westonia. In Muresk, the anthesis time of Suntop showed large differences (three days) between the SpB and BW sown areas (119 and 122), which may lead to significantly different levels of frost damage.
In Williams, the earlier-flowering lines of Bethlehem, B_3B and B_7A, were impacted the most by frost whereas the late-flowering lines Gregory, Tungsten and Yitpi were not affected (Fig. 2a, b; Table S7). Among the cultivars with a similar flowering window, such as Spitfire, Mace and Wyalkatchem, Westonia was the most susceptible to frost, whereas Mace showed the highest level of tolerance. The same trend was observed in Muresk. Spitfire displayed frost tolerance in Williams but was susceptible in Muresk. Its anthesis window was very close to that of Mace in the SpB sown area and to Wyalkatchem in the BW area. Suntop showed high frost tolerance in Williams but was susceptible in Muresk, where its flowering window was shorter. Floret sterility also showed the same levels of frost impact in selected parental and control lines belonging to the three anthesis windows (Fig. S3a, b).
According to the date of the frost event (September 14, 15 and 16, 2018), the time of impact was between 9 and 18 days before the average anthesis time. The 2019 Muresk field trial data for Suntop at 16 days before anthesis (DBA) are presented in Fig. S3c-g. The data show that wheat plants were at the flower development and anther differentiation stages at 16 DBA, according to Bonnett's description [64], and some flowers were at meiosis stage. Based on Koonjul et al. [4], lines are most vulnerable to abiotic stress during spikelet development and meiosis stages. The 2018 frost events incidentally hit the most susceptible meiosis stages of our trial.

Phenotype differentiations across DH populations and environments
Different frost impact levels were observed in six DH populations. On average, the frost impact levels in BW (8.0), Sp7A (8.0), SpB (8.1) and St3B (6.2) were higher than that in G7A (5.0) and SpM (3.7) populations. The impact levels in Muresk (8.6) were almost doubled in Williams (4.4) ( Table S8). The coefficients of variation were significantly higher in Williams compared with Muresk. The high heritability in frost impact was observed (0.71-0.84) in all six populations across two environments. Apart from the frost impact, plant height showed the highest heritability (0.96), followed by Zadok stages (0.79-0.95), anthesis (0.52-0.82) and maturity (0.4-0.87). The significantly positive correlation levels of frost impact (P < 0.01) between replicates and environments further demonstrated that the traits were under genetic control (Table S9).  only Suntop and Gregory carry both the Rht-B1b and the wild-type Rht-D1a alleles. The other six parental lines have the opposite genotype for these alleles. Therefore, Rht-B1b and Rht-D1b segregated in two populations of Suntop Â B_3B (St3B) and Gregory Â B_7A (G7A). No amplification was detected by primer markers of gene TaELF3-1A and TaELF3-1B. Fragments of 709 bp amplified by primers of TaELF-1DL-F5/R4 presented in Gregory, Mace, Suntop, Westonia and Chinese Spring which hold allele TaELF3-1DLa (delaying heading two days) whereas no amplification was in Bethlehem, B_3B, B_7A and Spitfire (Fig. S5d). The results indicate the TaELF3-1D segregations in BW, SpM, G7A and St3B population (Fig. 2c). Fragments of 1.4 kb and 1.28 kb were amplified by primer pairs of VRND4-ins.F4/R3 and VRND4-ins2. F1/R1 in both upstream and downstream of VRND4 insertion, respectively, and showed that Bethlehem, B_3B, B_7A, Gregory, Mace and Spitfire harbour VRN-D4 while the Suntop and Westonia were negative (Fig. S5d), indicating the VRN-D4 segregations in the populations of St3B and BW (Fig. 2c). Copy number differences for the Ppd-B1 allele were detected in parental lines and the segregations in SpB, Sp7A, G7A, and St3B populations were expected (Fig. S5e).

Cross population verified frost QTL and related potential candidate genes
Thirty individual frost QTL were detected in the six DH populations. Out of these, there were 18 major QTL that were each responsible for a phenotypic variation greater than 9.5% and were distributed across 13 chromosomes (2A, 2B, 2D, 3A, 4A, 4B, 4D, 5A, 5D, 6D, 7A, 7B, and 7D). Most frost QTL were closely linked to the QTL for anthesis and physiological maturity Zadok stages and to anthesis-related genes. The QTL results for individual DH populations are included in Table S10.
A large proportion (83%) of the detected frost QTL were showed up consistently in two or more DH populations, except for the QTL on SpM_1A, SpM_1D, SpM_3D, Sp7A_6D, and G7A_3B (Fig. S6). A frost QTL was detected on the homologous region on the short arm of chromosome 2A in the SpB and SpM populations, while the other two frost QTL distant from each other in SpM were closely linked to the anthesis and physiological maturity loci on the long arm of chromosome 2A (Table S10; Fig. S7). The photoperiod (Ppd) allele Ppd-A1 is located on the short arm of 2A (Table S11) A highly significant frost QTL (the highest LOD score 9.2) was detected on homologous regions of the 2B chromosome in four populations of Sp7A, SpB, G7A and St3B, and the additive effects were all contributed by Bethlehem and its substitution lines, while another QTL (LOD score 3.9) in SpM was also closely associated with the homologous region, with the phenotype contributed by Mace (Table S10; Fig. 3). A Ppd1-2B gene (gene bank number DQ885765) was used to identify the physical map location. The potential location of Ppd1-2B was on the short arm of 2B at 56.2 Mb, which was close the marker 2BM46469 (2B: 58.8 Mb) in Sp7A, SpB and G7A, and between markers of 2BM29812 (2B: 53.4 Mb) and 2BM73250 (2B: 77.2 Mb) in St3B (Table S11). Strikingly, the location of Ppd1-2B was in the frost and anthesis QTL region in those four populations, and the QTL results matched the potential segregation of Ppd-B1 gene copy number (Fig. 2c;  Fig. S5e), which points toward a contribution of Ppd-B1 copy number to the phenotype. This QTL corresponds to the newly identified Ppd1-2B location in this study.
On linkage groups 2D3 and 2D2, frost QTL were detected in both SpM and Sp7A, contributed by Mace and Spitfire, respectively. These QTL were associated with anthesis time, as anthesis QTL were closely linked to the homologous region of 2D (Fig. S8). Gene Ppd1-2D was located at 2D: 33.9 Mb from the upstream on the short arm of 2D. SNP M22365 (2D: 18.2 Mb) on Sp7A_2D1 and SNP M15795 (2D: 17.3 Mb) on SpM_2D1 were close to Ppd1-2D on the physical map. The frost and anthesis QTL on SpM_2D3 and Sp7A_2D2 tend not to be associated with the Ppd1-2D locus, since the Ppd-D1a locus does not expect to segregate in those populations (Fig. S5b). This frost QTL may imply a new phenology gene on 2D, associated with anthesis time.
On chromosome 3A, the frost QTL in SpB population was located on the short arm, whereas in St3B it was located on the long arm. The QTL on the short arm overlapped with anthesis QTL (Fig. S9). TaGI gene sequences (AF543844) were isolated [65] and were located on 3A (3A: 84.1 Mb), 3B (3B: 117.9 Mb), and 3D (3D: 71.9 Mb). One frost QTL in SpB_3A was located between SNP markers 3AM63098 (3A: 56.4 Mb) and 3AM34922 (3A: 107.7 Mb) (Table S11). TaGI-3A is most likely the candidate gene contributing to this frost QTL and the associated anthesis QTL. The 3A long arm QTL in the St3B population was located between 3AM9626 (3A: 711.0 Mb) and 3AM39004 (3A: 729.7 Mb), which are not associated with the TaGI gene. The Genbank number of KF769443 was used to identify the location of Eps-3A m on the physical map. Eps-3A m (3A: 740.1 Mb) was located between the markers 3AM10770 (3A: 739.3 Mb) and 3A73079 (3A: 741.2 Mb), about 10 cM away from the frost QTL (Table S11).
Significant frost QTL were consistently detected on the long arm of chromosome 4A in the SpM, G7A, and St3B populations (Fig. S10). Interestingly, the phenotypes were all attributable to the male parents of Mace, B_7A and B_3B. The frost QTL in SpM and G7A populations were located on the anthesis QTL region, and the late anthesis QTL were contributed by the female parents, Spitfire and Gregory. In other words, the early-flowering phenotype contributed by male parents led to frost susceptibility. The homologous gene of the DELLA protein (rht1-D1a; AJ242531) on 4A is TraesCS4A02G271000 (4A: 582.4 Mb), which is located between markers 4AM76744 (4A: 575.0 Mb) and 4AM43375 (4A: 597.6 Mb), and above the common marker 4AM77169 (Table S11). Interestingly, the newly identified rht1-4A locus in the current study was located 20 cM away from the minor frost QTL regions on the short arm of 4A in SpM and G7A populations. The homologous gene sequence of WSOC1-4D is TraesC-S4A02G320300, on chromosome 4A: 608.8 Mb, near SNP marker 4AM77381(4A: 612.1 Mb), which is 30-40 cM away from the frost QTL on the long arm regions of 4A in the G7A and St3B populations (Table S11). The potential WSOC1-4A locus is further away from the frost QTL in the SpM population.
Significant frost QTL were detected on the Rht-B1b region on chromosome 4B in G7A and St3B populations segregating for Rht-B1b ( Fig. 4a; Table S11), which were contributed by the female parents, Gregory and Suntop, and both harbour the Rht-B1b allele. No flowering QTL were detected in those regions. The plant height QTL with the LOD values over 14 were overlapped with the frost QTL in the semi-dwarf genes Rht-B1b region and were contributed by male parents B_7A and B_3B containing the wild-type Rht-B1a on 4B. Wide-type of Rht-B1a attributed the average height phenotype of 20%-30% and the semi-dwarf gene Rht-B1b contributed to frost damage QTL on 4B in G7A and St3B populations (Table S10; Fig. 4a).
A significant frost QTL with a LOD score of 16.0, attributable to Mace, was repeatedly detected on the distal region of chromosome 4B short arm in the SpM population. The QTL region was about 80 cM away from the recessive Rht-B1a gene loci. The contributing gene for this QTL is unknown. A minor frost QTL was  detected in the BW and G7A populations on the homologous region of the long arm of chromosome 4B. In the BW population, the frost QTL were closely linked to the anthesis QTL. Interestingly, frost QTL in the BW population was in the 4BM4887 (4B: 646.6 Mb) region, which is only 6.4 Mb away from the WSOC1-4B and the anthesis QTL, close to 4BM8859 (4B: 660.7 Mb) (Table S11). These results suggest that WSOC1-4B might be a new gene influencing flowering time and contributing to segregation for frost.
Likewise, significant frost QTL were detected on the Rht-D1b gene region on 4D in the G7A and St3B populations segregating for the Rht-D1b gene, contributed by B_7A and B_3B, which harbour the Rht-D1b allele (Fig. 4b; Table S11). The plant height QTL with high LOD value over 14 were located on Rht-D1b region and contributed 20% À 44% of the total phenotype by Gregory and Suntop containing Rht-D1a on 4D. The wide-type of Rht-D1a contributed to plant height whereas the semi-dwarf Rht-D1b attributed to the frost damage QTL on 4D in G7A and St3B populations (Table S10; Fig. 4b).
Frost QTL were consistently detected on the homologous region of chromosome 5A in SpM, Sp7A and SpB populations, with the highest LOD score of 7.7 found in Sp7A (Fig. 5). The frost QTL were in the Vrn-A1a region, attributed to B_7A and Bethlehem in the Sp7A and SpB populations (Table S11), respectively, whereas it was contributed by Spitfire in the SpM population. The frost QTL in the Sp7A population overlapped with the anthesis QTL contributed by Spitfire. Both B_7A and Bethlehem host a Vrn-A1a mutant and flower early. The results further indicate that frost damage was caused by the early flowering genotypes of B_7A and Bethlehem. A minor frost QTL in the St3B population was on the distal region of the short arm 5A, away from the anthesis QTL, and not associated with the flowering genes or their related phenotypes.
In the case of chromosome 5D, frost QTL were detected in SpM, BW and St3B populations (Fig. 6). The largest frost QTL (the highest LOD score 14.2) on 5D was in the Vrn-D1a region in the SpM population while in the BW population it was close to Vrn-D1a (Table S11). The frost phenotype in the SpM and BW populations was contributed by Spitfire and Bethlehem, which harbour Vrn-D1a, whereas for the St3B population in the Williams trial, the contributed phenotype was attributed to Suntop. Several anthesis and maturity QTL were detected in the Vrn-D1a region, indicating that frost damage is closely associated with plant growth stages. Another vernalisation related gene, VRN-D4, which is a homologue of Vrn-A1a, originated from a large segment of chromosome 5A inserted into the short arm of 5D [27,28]. According to the physical map of VRN-D4 (5D: 193 Mb), the gene location should be close to marker 5DM45442 (5D: 278 Mb) on the St3B _5D map, while it is next to 5DM62708 (5D: 154 Mb) on the BW map (Table S11). The frost QTL may stem from the VRN-D4 locus as the VRN-D4 was segregating between the parental lines of Bethlehem and Westonia, Suntop and B_3B, respectively.
Significant frost QTL were detected on chromosome 7A in SpM, SpB and BW populations (Fig. 7), with the highest LOD score of 10. The phenotype in the SpB and BW populations was contributed by Bethlehem. Frost QTL in the SpM population were located in three positions, mainly attributing to Mace. According to the physical map of TaVRT-2-7A (7A: 128.8 Mb), the closest SNP marker is 7AM1849 (7A:128.5 Mb), close to 7AM75587, which was tightly linked to the frost QTL in the SpB and BW populations, while in the SpM population the location was closely linked to the anthesis QTL on 7A (Table S11). One frost QTL was detected in the SpM and SpB populations on the homologous distal region on chromosome 7A. This corresponds to a new segregating locus for frost impact on 7A, which was not associated with the flowering genes or their related phenotypes.
For chromosome 7B, frost QTL were detected in the SpM and BW populations (Fig. S11a). However, the frost QTL regions were distinct from each other. The frost QTL in the SpM population, contributed by Mace, was on the terminal region of the short arm of chromosome 7B, away from the anthesis QTL, whereas in the BW population, the frost QTL were closely linked to the QTL for anthesis and maturity, and were contributed by Westonia. The anthesis QTL in the SpM and BW populations were located on the homologous region on 7B. Another anthesis QTL detected exclusively in the SpM population was located on the distal region of the long arm of 7B. In the SpM population, the physical map of VRN-B3 (7B: 9.7 Mb) was close to marker 7BM76084 (7B: 6.7 Mb), which may contribute the frost phenotype in SpM (Table S11) (Table S11), closely linked to the QTL for anthesis and frost, suggesting the involvement of this locus in frost damage in the BW population. In the SpM population, TaVRT-2-7B is not close to the frost QTL.
A QTL with a LOD score of 4.6 also appeared on the homologous regions on chromosome 7D in the Sp7A and SpB populations (Fig. S11b). The frost phenotype was contributed by B_7A and Bethlehem and was highly associated with anthesis QTL. TaVRT  (Table S11). This gene tends to contribute the significant frost QTL on 7D.

Anthesis QTL not segregating for frost tolerance
No frost damage QTL were detected on chromosome 1B. However, anthesis and maturity QTL were present on the central chromosome regions in the G7A and SpM populations and on the distal region of 1B short arm in the SpB population (Fig. 8). Gene TaFT3-1B (1B: 581.4 Mb) was tightly linked to SNP marker 1BM77588 (1B: 581.6 Mb), which is next to the anthesis QTL region (1BM54518) in the G7A population ( Fig. 8; Table S11). In the SpM population, the TaFT3-1B locus was possibly located above 1BM73284 (1B: 662.9 Mb), on the anthesis QTL region (Table S11). On the physical map, TaELF3-1B was located on the distal region of the long arm in these three populations. Two other anthesis-related genes, TaWUSCHELL-1B (1B: 53.3 Mb) and TaTOE1-1B (1B: 59.1 Mb), were close to markers 1BM42781 (1B: 56.8 Mb) and 1BM47932 (1B: 58.7 Mb) on the Chinese Spring physical map (Table S11). These two genes were located 25 cM away from the maturity QTL in the SpB population. It is not clear whether the QTL was influenced by them.
Another chromosome without frost damage QTL was 5B, even though anthesis QTL are present in the St3B, SpB, Sp7A and SpM populations, and Vrn-B1a was segregating in the St3B, SpB and Sp7A populations (Fig. 9). The anthesis QTL was on or closely linked to the Vrn-B1a allele in the St3B and SpB populations (Table S11), and to the distal regions of either the long arm in the Sp7A population or the short arm in the SpM population.

Allele frequency of VRN1 and Rht genes in historical lines
According to the frost trial results, except for Vrn-B1a, the mutated Vrn-A1a, Vrn-D1a and Rht alleles were highly associated with the frost QTL and contributed to the frost damage phenotype. Since global warming accelerates terminal drought severity, early flowering genes with early maturity phenotype have been gradually introduced into modern varieties through wheat breeding for avoiding late dry season. In the meantime, the proportions of Rht genes were also increased. The utilization of these gene may increase the frost sensitivity in modern varieties. Therefore, it is interesting to see the allele frequency of these mutated genes in historical varieties. In the current study, a set of 171 commercial wheat cultivars, released between 1890 and 2015, were surveyed for the VRN1 and Rht allele frequencies. For VRN1 and Rht, recessive alleles made up 95%-100% of varieties released before 1967, except for vrn-B1 (15%) (Fig. S12). After 1968, Vrn-A1a was maintained at around 61% whereas Vrn-B1a dropped from 85% to 50% and remained at that level thereafter. The Vrn-D1a allele appeared mainly in the varieties released after 2000. The Rht-B1b allele was presented in a large proportion of varieties released after 1968 (56%-69% of lines) while the Rht-D1b allele gradually increased to near 41% in recent years.

The most sensitive growth stage to frost
In previous drought studies in wheat, the young microspore stage of pollen development, before anthesis, appears to be the most sensitive to mild water stress [4,66]. It has been reported that in rice cold-induced pollen sterility at the young microspore stage had effects comparable to those of drought stress [2]. Cold-tolerant lines at the microspore stage are also tolerant to drought stress. The same mechanism was observed in sorghum [67,68].
The young microspore stage in wheat is the time when the auricle distance (AD) between flag leaf and penultimate leaf is 5-8 cm, around 10 days before anthesis [5,66]. In our 2019 wheat trials at the Muresk site, the young microspore stage of Suntop was 16 days before anthesis and the AD was 6 cm. In 2018, the frost events at the Williams and Muresk locations occurred 9-18 days before the average anthesis time, which was the most vulnerable growth stage for wheat to endure low temperatures. Although the lowest temperature in Muresk remained above 0°C, plants still suffered from the sudden temperature drop. According to previous studies [69,70], cold stress induces ABA accumulation in rice anthers, which represses anther cell wall invertase activity. This in turn hampers sugar transport from the tapetum to the pollen and pollen sterility occurs. It is interesting that cold and drought stresses share the same pathways in inducing pollen sterility.

Non-escaping mechanism related frost tolerance
Since the reproductive stage is the period most sensitive to frost, flowering time becomes very critical. Ideally, late flowering is helpful in avoiding frost events, which happens mostly during early spring. This is, however, contradictory to the needs of avoiding terminal drought, which requires early maturity [71]. Being able to combine early flowering with reproductive frost tolerance is thus highly desirable for wheat breeding. In the previous study, frost QTL were detected in the same region as the Vrn-A1a allele on chromosome 5A [10,72]. It has been reported that plants with VRN1 copy showed normal flowering but reducing frost tolerance [73]. In our study, significant frost QTL were detected close to the Vrn-A1a allele in both the Sp7A and SpB populations at the Williams location. In the SpM population grown in Muresk, significant frost QTL were detected close to the Vrn-D1a allele, where the QTL for anthesis and maturity were located for all three environments,   which further validated the tight connections between frost damage and the early anthesis alleles Vrn-A1a and Vrn-D1a.
On the other hand, the Vrn-B1a allele segregated in the BW, SpB, Sp7A, and St3B populations, with anthesis QTL closely linked to the Vrn-B1a gene in the SpB and St3B populations. Strikingly, no frost damage QTL was detected on chromosome 5B in any of the six pop-ulations studied at the two locations. The link between frost and anthesis was not observed on 5B. In our previous study on days to anthesis, the contribution of winter allele vrn-B1 was less than that of vrn-D1 and vrn-A1 [56]. In the historical lines studied, the frequency of spring alleles Vrn-A1a and Vrn-D1a in recent varieties amounted to 0.62 and 0.40, respectively, while the frequency of Vrn-B1a diminished from 0.85 to 0.52 over time. An increased proportion of Vrn-A1a and Vrn-D1a may increase the risk of frost damage. The most popular variety, Mace, only possesses the Vrn-B1a allele and shows good frost tolerance. This clearly demonstrates that the early flowering phenotype induced by Vrn-B1a is associated with higher frost tolerance. In wheat, VRN1 determines the most natural variation in flowering. The functional VRN1 proteins are responsible for the apical meristem transition from the vegetative to reproductive phase. However, VRN1 protein are not essential for wheat flowering [74]. Although the underlying mechanism is unknown, it can be speculated that the metabolites associated with the Vrn-B1a functional network may contribute to the increased tolerance [75]. Because of the global warming, the earlier flower lines without frost damage are highly demanded [76]. Our results revealed that the Vrn-B1a allele can be utilised in breeding for frost tolerance. Another potential anthesis gene associated with the nonescaping frost mechanism is TaFT3-B1 on 1B. Gene TaFT3-1B was very tightly linked to the anthesis QTL on 1B in G7A and SpM populations. The study of Zikhali et al. [46] showed that TaFT3-B1 deletion lines was associated with late flowering while increased the gene copy number was related to early flowering. TaFT3-B1 gene was suggested to promote flowering. Two maturity QTL were detected on 1B in SpB populations, where two potential anthesis genes TaWUSCHELL-B1 and TaTOE1-B1 were closely linked to. In Maize, TARGET OF EAT1 was also named as TaSRR1-B1, TaWUSCHELL-B1 and TaTOE1-B1, respectively. Similar to TaFT3-B1, TaTOE1-B1 had higher gene expression during short days and showed early flowering function [46]. Nevertheless, no frost QTL were detected on 1B in those three populations in two locations. It implies that the markers or the copy numbers for the TaFT3-1B are able to be used for the purposes of shortening the days to anthesis and reducing frost damage risk.

Recessive Rht-B1a and Rht-D1a genes associated with frost tolerance
Significant frost QTL were detected in the proximity of the Rht-B1b and Rht-D1b gene regions in the G7A and St3B populations, in which the two genes were segregating. In both populations, the frost QTL were contributed by the mutated type of Rht-B1b and Rht-D1b, the semi-dwarf alleles that are associated with short plant height as the plant height QTL were contributed by the wide-type of Rht-B1a and Rht-D1a (recessive rht genes). This indicates that the mutated or short plant genotypes are more prone to frost. On the other hand, the wild types Rht-B1a and Rht-D1a are relatively frost tolerant.
Wheat life cycle can be divided into two phases: one is the stem elongation from jointing to anthesis and another is the grain filling stage from anthesis to maturity. The former stage is more important for yield as the number of fertile florets at anthesis will be determined during this stage [77]. The pre-anthesis phase may be more sensitive to photoperiod or temperature, and different levels of hormones (gibberellin, auxin or cytokinin) [78][79][80][81]. Since several dwarfing genes are associated with GA biosynthesis or signalling, the possible effect on spike development may exist [82,83]. It is well known that lines with Rht-B1b and Rht-D1b alleles have reduced sensitivity to gibberellic acid (GA), and GA-responsive growth is therefore repressed [84]. This may affect GA levels among the different combinations of dwarfing genes, and thus lead to varying flowering times. It has been reported that GA can increase the transcription level of SOC1 including a MADS-box gene, promoting flowering [47]. In a study of ethylene effects in the GA-GID1-DELLA signalling pathway, ethylene leads to a reduction of GA, which delays floral induction [84]. Nevertheless, previous research has shown that GA-insensitive dwarfing genes have no effect on spikelet primordia at the shoot apex nor on the number of leaves and internodes [85]. In the current study, no flowering QTL were detected in the Rht-B1b and Rht-D1b gene regions. This suggests that the regulatory function of the Rht-B1b and Rht-D1b genes on anthesis is minor. On the other hand, the variations in plant canopy architecture linked to dwarfing genes, for example, height, tillering, ear structure etc. may affect the frost tolerant levels. Selecting taller or wild-type Rht-B1a and Rht-D1a alleles for frost tolerance in breeding programs will not cause delayed flowering time or maturity, thus these Rht-B1a and Rht-D1a alleles can be recommended for frost tolerance breeding.
In previous studies, the semi-dwarfing alleles Rht-B1b and Rht-D1b did not lead to higher grain yields in drought environments [51,86]. The QTL analysis showed that the recessive Rht-B1a and Rht-D1a alleles are associated with higher grain yields while the semi-dwarfing alleles are associated with lower grain yields [51]. Our results further demonstrated the disadvantages of the Rht-B1b and Rht-D1b alleles in frost-prone environments. From 1967 to 2015, the frequencies of the Rht-B1b and Rht-D1b alleles have increased from 0 to 0.6 and 0.4 in Australian wheat varieties, respectively. This is in line with the first green revolution in that the dwarfing genes are utilised widely. However, a recent study has revealed that the short plant phenotype has become a bottleneck for wheat grain yield improvement [51,86]. Achieving high grain yield by increasing plant height has become a newly established breeding approach. Our study indicates that eliminating the dwarfing alleles Rht-B1b and Rht-D1b in breeding can increase frost tolerance in wheat.

High copy number of Ppd1-2B increase the sensitivity to frost
Day length (photoperiod) affects plant flowering time. QTL for frost and the time of anthesis and maturity were repeatedly detected on chromosomes 2A, 2B and 2D. The Ppd1-2A allele was tightly linked to the frost QTL on 2A. Remarkably, four out of five frost QTL (next to anthesis QTL) were detected in the proximity of Ppd1-2B in five populations. Due to the Ppd1-2B copy number differences between parental lines of these four DH populations, it is likely that the significant frost QTL were contributed by the segregation of Ppd1-2B gene copy numbers, with a high copy number of Ppd1-2B shortening the days to anthesis and inducing early flowering [34], thus contributing to the frost QTL. Ppd1-2D did not segregate in our six populations and no strong linkage to the frost QTL on 2D was detected.

Escaping mechanism related frost tolerance
Besides the Vrn-A1a and Vrn-D1a, the following early flowering related genes are associated with the frost damage QTL and the recessive flowering genes (flower late genes) correspond to frost tolerance, which led a frost escaping mechanism. Gene VRN-B3 was mapped to the short arm of chromosome 7B (7B: 9.7 Mb) [26]. In the SpM population, VRN-B3 is tightly linked to a 7B short arm frost QTL that mapped on the top region (7BM55522) of the chromosome. This implies that VRN-B3 contributes to the frost phenotype in the SpM population. In addition, the TaVRT genes on group 7 appeared to be tightly linked to the QTL of anthesis and frost. The TaVRT gene is independently modulated by photoperiod and vernalization. It can inhibit VRN1 activity through binding to the CArG box of the VRN1 promoter in vivo. After vernalization, both TaVRT2 and VRN2 functions are repressed, which causes VRN1 accumulation for the transition from vegetative to reproductive phase [87]. The functions of TaVRT genes seem to be critical in the SpM, Sp7A, SpB, and BW populations.
Numbers of genes in photoperiod pathway also regulate anthesis. TaGI located on 3A, 3B and 3D, regulate CO which then mediates photoperiodic flowering [44,65]. In previous studies, TaGI was found to be initiated by photoperiod and then expressed in both vegetative and reproductive tissues, suggesting a function in anthesis regulation [44,65]. Coincidentally, in our study, the frost and anthesis QTL, and maturity QTL were localized on the TaGI regions on 3A and 3B in SpB and G7A populations, respectively. This indicates that TaGI functions are both in anthesis and frost.
Apart from the vernalisation and photoperiod genes, eps genes are the third factor regulating heading and anthesis. One eps like QTL identified on the long arm of 1DL, namely, TaELF3, held the function to the T. monococcum Eps-Am1 locus [37,43]. The homologs of TaELF3 were also identified in 1AL and 1DL [40]. One significant frost damage QTL contributed by Spitfire holding TaELF3-1DLb (heading two days earlier) on the long arm of 1D in SpM population were tightly link to the homologous genes of TaELF3 on 1D. It indicates the eps gene functions in heading and contributes to the frost impact in the location of 1D.
Heading related gene TaHD1 also designated as CONSTANS2 (CO2), is the homolog of CO in wheat, and the competitor of VRN2 [47,48]. TaHD1 are identified on 6A, 6B and 6D. In our study, one of the frost QTL together with the anthesis QTL were on the region of TaHD1-6D physical map location, implying that the TaHD1 gene contributes to the phenotype.
In summary, six DH populations planted at two different locations in 2018 were hit by a severe frost event at the critical reproductive phase. Evaluation of frost damage showed that the plant growth stage most susceptible to low temperature (<2°C) was during the young microspore stage (10-18 days before anthesis). Out of the 30 frost QTL detected, 18 major QTL were mapped onto 13 chromosomes. Most frost QTL overlapped or were closely linked to the QTL for anthesis and maturity Zadok stages as well as to anthesis-related genes. However, the frost tolerance contribution by these QTL clearly stems from the late flowering alleles, illustrating the frost escape mechanisms and indicating that they are not useful in wheat breeding for water-limited environments. The mutated Vrn-A1a, Vrn-D1a, Rht-B1b, and Rht-D1b alleles and a high-copy number Ppd-B1 allele contributed significantly to frost damage. Nevertheless, QTL or genetic factors outside the escape mechanisms were detected in the current study. Anthesis QTL were repeatedly detected in the proximity of the Vrn-B1a region and on chromosome 1B, whereas no frost QTL were detected on these two chromosomes. These striking results strongly imply that the Vrn-B1a and TaFT3-1B alleles on 5B and 1B should be utilised in breeding for frost tolerance, as the early-flowering phenotype associated with these two genes is frost tolerant. Meanwhile, the recessive non-dwarfing alleles Rht-B1a (rht1) and Rht-D1a (rht2) associated with normal plant height could also be used in breeding for reproductive frost tolerance without delaying the flowering time.

CrediT authorship contribution
Jingjuan Zhang, MD Shahidul Islam, Wujun Ma proposed the frost damage investigation. Jingjuan Zhang visually scored the frost damage, implemented the molecular work and field work, analyzed the data and wrote the paper. Ben Biddulph supervised frost damage phenotyping; Wujun Ma and Jorge E. Mayer contributed to the substantial paper writing and revisions. Junkang Rong surveyed the parental lines and provided Bethlehem and CSL seeds for DH population construction; Kefei Chen designed the field trials and statistically adjusted phenotype data. Zaid Alhabbar, Masood Anwar, and Rongchang Yang participated in the DH population construction; Jingjuan Zhang, Yun Zhao, Hang Liu, Maoyun She, Mirza Dowla, Sonia Afrin, and Nandita Roy constructed the genetic linkage map; Jingjuan Zhang, MD Shahidul Islam, Wujun Ma, Yun