Elucidation of the genetic architecture stabilizing the heading time in barley (Hordeum vulgare L.)
Maho Okuma, Kazusa Nishimura, Emiko Aoki, Masahiro Tsuchiya, Yuki Monden, Kenji Kato, Hidetaka Nishida
TL;DR
This study identifies genetic factors in barley that help maintain stable heading times despite climate changes, aiding in breeding early-maturing crops.
Contribution
The study discovers QHd.ouj-4H and HvPHYC as key genes influencing heading time stability in barley, modulated by HvCEN.
Findings
Three QTLs (QHd.ouj-2H, QHd.ouj-4H, QHd.ouj-5H) were identified for heading time in barley.
HvCEN and HvPHYC were linked to QHd.ouj-2H and QHd.ouj-5H, respectively, affecting heading time stability.
Interactions between QTLs show that HvCEN modulates the effects of QHd.ouj-4H and HvPHYC on stability.
Abstract
QHd.ouj-4H and HvPHYC were first found to affect the heading time stability under changing climate in barley, modulated by HvCEN, offering insights for stable early-maturing breeding. Heading time stability is essential for stable production of barley under the recent changes in climate. In this study, we sought to identify QTLs regulating heading time and stability, using RILs derived from a cross between two Japanese cultivars, Kashimamugi and Ishukushirazu, both of which are early-heading cultivars with spring growth habit, but differ in heading time stability. QTL analysis was performed by detecting genome-wide SNPs and InDels by MIG-seq, and detected three heading time QTLs (QHd.ouj-2H, QHd.ouj-4H, and QHd.ouj-5H). Among these, HvCEN and HvPHYC were considered to be the causative genes for QHd.ouj-2H and QHd.ouj-5H, respectively, while the gene underlying QHd.ouj-4H remains…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5- —JSPS KAKENHI
- —Okayama University
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsWheat and Barley Genetics and Pathology · Genetic Mapping and Diversity in Plants and Animals · Plant nutrient uptake and metabolism
Introduction
Barley (Hordeum vulgare L.) is the fourth most cultivated cereal crop in the world, and is used for feeding, brewing, and food. Global climate change has become a realistic threat in recent years, and is causing losses and instability of barley production in many parts of the world. For example, high temperatures at flowering time reduce pollen fertility, resulting in yield losses in wheat and barley (Jacott and Boden 2020). In addition, high temperatures at the ripening stage affect yield and damage grain quality (Mahalingam and Bregitzer 2019). High temperature stress during reproductive development can be avoided by adjusting flowering at the optimum time. Thus, heading time (HT) is a key trait for stable production under global climate change. To enhance barley breeding for stable production, it is essential to elucidate the molecular genetic mechanism for adjusting HT.
HT is a complex trait consisting of a vernalization requirement, photoperiod sensitivity, and earliness per se (Yasuda and Shimoyama 1965). Vernalization is a process by which exposure to prolonged periods of cold provides competence to a flower and protects the meristems from low temperatures. This trait is mainly controlled by three major genes, VRN-H1 (Yan et al. 2003), VRN-H2 (Yan et al. 2004), and VRN-H3 (Yan et al. 2006). VRN-H1 is located on the long arm of chromosome 5H and encodes an ortholog of the Arabidopsis APETALA1/FRUITFULL-like MADS-box transcription factor (Yan et al. 2003). VRN-H2 is a flowering repressor located on the long arm of chromosome 4H and encodes two tandem ZCCT-H genes with a CCT domain (Karsai et al. 2005) VRN-H3 is located on the short arm of chromosome 7H and encodes an ortholog of the Arabidopsis FLOWERING LOCUS T (Yan et al. 2006; Faure et al. 2007; Kikuchi et al. 2009). If at least one of the VRN genes is of the spring type (Vrn-H1, vrn-H2, or Vrn-H3), a floral transition occurs without exposure to prolonged cold periods, but if all VRN genes are of the winter type (vrn-H1, Vrn-H2, and vrn-H3), prolonged periods of cold are required. Barley cultivars of the former type are classified as spring type, and the latter as winter type. In barley and wheat, the vernalization requirement is important for adaptation to cold winter conditions, as indicated by the geographical distribution of spring and winter type accessions (Iwaki et al. 2001; Saisho et al. 2011). Photoperiod sensitivity is the response to suppress or accelerate flowering, depending on daylength. Barley is a long-day plant, which promotes floral transition under long-day conditions but delays it under short-day conditions. PPD-H1 and PPD-H2, located on the short arm of chromosome 2H and the long arm of chromosome 1H, respectively, have been identified as genes that control the photoperiod response in barley (Laurie et al. 1995; Turner et al. 2005). PPD-H1 encodes HvPRR37, the ortholog of Arabidopsis PSEUDO RESPONSE REGULATOR 7 (PRR7). Mutation at the CCT domain of PPD-H1 is associated with reduced sensitivity to photoperiod (Turner et al. 2005) and promotes flowering under long-day conditions as a dominant allele (Turner et al. 2005; Hemming et al. 2008; Sharma et al. 2020). PPD-H2 encodes HvFT3, a known flowering promoter under short-day conditions, but the details of its function in controlling HT are unclear (Laurie et al. 1995; Faure et al. 2007; Kikuchi et al. 2009; Casao et al. 2011). HvPHYC, which encodes the red and far-red photoreceptor located on the long arm of chromosome 5H, also affects HT (Nishida et al. 2013; Pankin et al. 2014). Two alleles of HvPHYC are known to affect HT: an early-heading allele (HvPHYC-e or haplotype 4) and a late-heading allele (HvPHYC-l or haplotype 1). In addition, haplotype 7 is a medium-heading allele (Nishida et al. unpublished). Of the three alleles, the late-heading allele is a wild type, and the others share an SNP with a nonsynonymous substitution in the GAF domain. The medium-heading allele has one additional SNP with a nonsynonymous substitution.
Earliness per se (EPS) genes control HT independent of vernalization and photoperiod. For example, HvELF3, located on the long arm of chromosome 1H, is the barley ortholog of the Arabidopsis circadian clock gene Early flowering 3 (ELF3), whose loss-of-function mutation causes early flowering regardless of daylength conditions (Faure et al. 2012; Zakhrabekova et al. 2012). Similarly, HvLUX1, located on the long arm of chromosome 3H, is the barley ortholog of the Arabidopsis circadian clock gene LUX ARRHYTHMO/Phytoclock1 (LUX/PCL1), and hvlux1 causes early flowering regardless of daylength conditions (Campoli et al. 2013). In addition, flowering repressor HvCEN, located in the centromere region of chromosome 2H, is an ortholog of Arabidopsis TERMINAL FLOWER 1 (TFL1) and Antirrhinum CENTRORADIALIS (Laurie et al. 1995; Comadran et al. 2012). Thirteen haplotypes of HvCEN have been reported, and one amino acid substitution (Pro135Ala) in this gene determines the growth habit (Comadran et al. 2012). In other words, Ala135 is associated with spring growth habit and Pro135 with winter growth habit. In recent years, the role of EPS genes in controlling HT has been revealed. Initially, EPS was defined as a trait independent of environmental factors, but it is now clear that some known EPS genes interact with environmental cues. For example, the expression of HvLUX1 was higher at 25 °C than at 15 °C, and the expression pattern also changed (Ford et al. 2016). Furthermore, the circadian clock genes and their response to temperature were mediated by HvELF3 (Ford et al. 2016).Recently, studies using various mapping populations (such as F2 and RILs) and diversity panels have revealed not only the solitary effect of these genes on heading, but also the interactions among them (Fernández-Calleja et al. 2021). It has been reported that VRN-H1, VRN-H2, and VRN-H3 form a feedback regulatory loop through epistatic interactions (Distelfeld et al. 2009). According to this model, high expression of VRN-H2 under long-day conditions during autumn suppresses the expression of VRN-H3, which upregulates VRN-H1 and promotes flowering. Subsequently, upregulation of VRN-H1 during winter suppresses VRN-H2, resulting in increased expression of VRN-H3 and flowering. An interaction among HvELF3, HvCEN, and VRN-H3 was also reported (Casas et al. 2021). When a barley cultivar harbors a ‘Logan’ type allele of HvELF3 (functional and early-heading allele) and a ‘Beka’ type allele of HvCEN (late-heading allele, Haplotype III), VRN-H3 does not affect HT. However, when a barley cultivar harbors other alleles of HvELF3 and HvCEN, VRN-H3 affects HT.
Based on the genetic studies of HT summarized above, the desirable alleles of HT-related genes can be easily selected by the PCR-based method. Spring type alleles of VRN genes and early-heading alleles of PPD genes have been selected for breeding early-maturing cultivars, which are effective in avoiding abiotic stresses such as high temperature, drought, or monsoon rains during the grain filling period. However, as a result, barley cultivars have lost the mechanisms to suppress ear development by sensing the environmental cues, and HT became unstable under the changing global climate. Genetic analysis of HT stability is of increasing importance, and several studies have been conducted using diversity panels of barley (Sato et al. 2020) and wheat (He et al. 2024). However, no genetic studies using segregating populations have been reported.
In this study, we sought to elucidate genetic variation in HT stability of barley over 12 consecutive years. HT data were analyzed to assess stability, and barley cultivars with stable and unstable HTs were identified. RILs (Recombinant Inbred Lines) derived from crosses of these cultivars were used to comprehensively identify loci controlling HT based on data from four growing seasons. This study aimed to assess the contribution of HT-related QTLs to the stability of HT across years, in addition to identifying QTLs associated with HT stability. To this end, we then evaluated the effect of these loci on HT stability, RILs data collected over four growing seasons (2015/2016, 2016/2017, 2019/2020, and 2020/2021). Through these analyses, we identified a novel QTL, QHd.ouj-4H, on chromosome 4H, that controls HT. We also found that HvPHYC and QHd.ouj-4H are involved in the regulation of HT stability, and that both interact epistatically with HvCEN, suggesting that the combination of these QTLs plays a key role in regulation HT stability.
Materials and methods
Plant materials and phenotyping of HT
The two early-heading cultivars with spring growth habit, Kashimamugi (KSM) and Ishukushirazu (ISH), were grown at the experimental field of Okayama University (34° 41′ 05.1" N, 133° 54′ 57.6" E) for 12 consecutive cropping seasons, from 2011/2012 to 2022/2023. They were directly sown in the field at a 10 cm spacing. Sowing was mostly done in November each year, but on December 25 in the 2019–2020 season to test the response to unusual growing conditions (Supplementary Fig S1). For evaluation and QTL analysis of HT stability, we developed a total of 191 RILs derived from a cross between KSM and ISH. The F5-generation RILs were sown on November 27, 2015, the F6 on November 25, 2016, the F7 on December 25, 2019, and the F8 on November 18, 2020. In each cropping season, HT, defined as the day when the tip of spike emerged from the flag leaf sheath, was recorded as days after 31 March (e.g., April 1 = day 1). The average HTs of three plants (F5 generation) or five plants (F6–F8 generations and parental cultivars) were used. QTL analysis of HT was conducted using 149 RILs selected randomly (Supplementary Table S1).
Evaluation of HT stability
Interannual variance of the mean HT was calculated for parental cultivars using 12 years’ data (2011/2012 to 2022/2023) and 4 years’ data (2015/2016, 2016/2017, 2019/2020, and 2020/2021), and for 140 RILs using 4 years’ data (2015/2016, 2016/2017, 2019/2020, and 2020/2021) from F5 to F8. HT stability, defined as the interannual variance of HT across years, was calculated for each cultivar and RIL, with larger values indicating instability and smaller values indicating stability. QTL analysis of HT stability was conducted using 149 RILs mentioned above.
Resequencing of parental cultivars and MIG-seq library construction for RILs
DNA was extracted from young leaves of KSM and ISH, the parental cultivars of RILs, by MagExtractor -Plant Genome- Nucleic Acid Purification Kit (Toyobo Co., Ltd., Shiga, Japan). DNBSEQ (MGI, Shenzhen, China) was then used to obtain genome sequences at about 30 × coverage of the barley genome size. In this study, to genotype RILs, we used multiplexed ISSR genotyping by sequencing (MIG-seq) (Nishimura et al. 2022), which is known to be effective in detecting polymorphisms in wheat with large genome sizes. First, leaf blades of 149 RILs (F7) were sampled, dried, and crushed using a multi-beads shocker (Yasui Kikai Co., Ltd., Osaka, Japan). DNA was then extracted by the SDS method (Edwards et al. 1991). The concentration of genomic DNA was measured using Nanodrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA) and then adjusted to 25 ng/µl for MIG-seq library construction. DNA of 149 RILs was used to construct the library, using a slightly modified version of the method described by Nishimura et al. (2024). Using KOD One PCR Master Mix (Toyobo Co., Ltd., Shiga, Japan) and a MIG-seq first PCR primer set (Supplementary Table S2), a first PCR was performed at 98 °C for 1 min, followed by 25 cycles at 98 °C for 10 s, 38 °C for 10 s, and 68 °C for 10 s, with a final extension at 68 °C for 5 min. The first PCR product, diluted to 1/50 concentration, indexing primer (Supplementary Table S3), and KOD One PCR Master Mix were used for the second PCR (indexing PCR) at 98 °C for 1 min, followed by 20 cycles at 98 °C for 5 s, 54 °C for 5 s, and 68 °C for 5 s, with a final extension at 68 °C for 5 min. The second PCR products were pooled, and purified with Ampure XP (Beckman Coulter, Inc., USA). The purified library was then mixed with KOD One PCR Master Mix, Illumina Primer P1 (5´-AATGATACGGGCGACCACCGA-3´), and Illumina Primer P2 (5´-CAAGCAGAAGACGGGCATACGA-3´) for reconditioning PCR. Reconditioning PCR conditions were as follows: 98 °C for 40 s, 54 °C for 15 s, 68 °C for 30 s, and final extension 72 °C for 10 min. Size selection for the MIG-seq library was then performed twice, using SPRIselect (Beckman Coulter, Inc., USA) at ratios of 1:0.75 and 1:0.56 (DNA:SPRIselect). After purification and size selection, the library was sequenced using Illumina NovaSeq X Plus.
QTL analysis for HT and HT stability in RILs
Short reads obtained by sequencing the MIG-seq library were filtered by fastp (Chen et al. 2018) with the following options: (−3 -q 20 -n 10 -f 17 -F 17). It was then mapped to the MorexV3_pseudomolecules_assembly by snap-aligner (Bolosky et al. 2021). Samtools (Danecek et al. 2021) was used to separate bam files by chromosome, based on the reference genome information. The samtools mpileup (Li et al. 2009) and bcftools call (Li 2011) were then used for variant calling. The vcf file was filtered by vcftools (Danecek et al. 2011) with the following options: (–max-missing 0.7 –minDP 5 –max-alleles 2 –min-alleles 2 –maf 0.1). These vcf files were used to construct a linkage map using Lep-Map3 (Rastas 2017). The HT-related genes that exhibited polymorphisms between the parents were genotyped in the RILs (F7) and integrated into the linkage map. Based on this linkage map, QTL analysis for HT of each generation from F5 to F8 was performed by R/qtl (Broman et al. 2003), using the simple interval mapping (SIM) and composite interval mapping (CIM) method. Similarly, we performed QTL analysis for HT stability. The logarithm the of odds (LOD) score threshold was determined based on 1,000 permutation tests. The loci where the LOD score exceeded the threshold were defined as the QTL candidate region.
Search for candidate genes of QHd.ouj-4H by SnpEff
The KSM and ISH reads by resequencing were first filtered by fastp with following options: (-z 9 –detect_adapter_for_pe -f 1 -t 1 −5 −3 -q 30 -l 50 -w 12). Their cleaned reads were then mapped to MorexV3_pseudomolecules_assembly by snap-aligner with the following options: (-t 4 -so -F s -F b -o -bam). From the bam files, only the QTL regions detected in this study were extracted, using samtools view. The samtools mpileup with options (-d 1,000,000 -q 20) and VarScan2 (Koboldt et al. 2012) with options (-Xmx16g –output-vcf –min-var-freq 0.01) were used for variant calling. For the vcf files, we added annotations using SnpEff (Cingolani et al. 2012b) and extracted only the mutations with a HIGH impact set by SnpSift (Cingolani et al. 2012a). The SnpSift filter option was: ("ANN[*].IMPACT has 'HIGH'"). The mutations that occurred in genes were searched on BARLEX (http://viewer.shigen.info/harunanijo/index.php) for detailed gene descriptions. A TBLASTN search was then performed on these genes against the Arabidopsis genome (TAIR10) on Ensembl Plants (https://plants.ensembl.org/index.html), and the gene with the largest alignment score on TBLASTN was designated as the barley homolog of the Arabidopsis gene. We then investigated whether Arabidopsis genes have been reported to affect flowering time, using ThaleMine v5.1.0 (https://bar.utoronto.ca/thalemine/begin.do). We investigated the expression levels of these genes based on Mascher et al. (2017) in BarleyExpDB (Li et al. 2023), and genes that were expressed in at least one developmental stages or tissues were considered as candidate genes.
Diagnostic PCR to genotype the major HT-related genes
The genotypes of VRN-H1, VRN-H2, VRN-H3, PPD-H1, and PPD-H2 were determined for KSM, ISH, and F5 RILs using allele-specific DNA markers reported previously (Nishida et al. 2013) (Supplementary Table S4). For the genotyping of HvCK2α, HvPHYC, and HvCEN, the primers were designed using the sequence reported previously (Comadran et al. 2012; Nishida et al. 2013). For QHd.ouj-4H, a novel HT-related QTL identified in this study, a PCR marker was designed based on InDels identified from the resequencing data of the parental lines and located in the genomic region flanking the significant SNPs detected by QTL analysis. The details of the primer sequence, PCR conditions, and restriction enzyme treatment are summarized in Supplementary Table S4.
QTL effect on HT and on HT stability
The diagnostic PCR in the F5 generation showed that 40 of 149 RILs used for QTL analysis were heterozygous at any of the three QTL loci or had recombination within the QHd.ouj-4H region. These 40 RILs were excluded from analysis of QTL effects, and an additional line was also excluded due to abnormally poor growth. Instead, additional 31 RILs homozygous at the three QTL loci with no recombination within the QHd.ouj-4H region were employed. The resulting 139 RILs were separated into 2 groups or 8 groups, based on QTL genotype, and HT and HT stability were compared between groups. HT differences were tested by Welch’s t test when comparing two groups, or Tukey–Kramer test for more groups. To assess HT stability, we transformed the interannual variance of HT, defined as an index of HT stability, after O'Brien (1981), so that it could be analyzed using analysis of variance (ANOVA) and Tukey–Kramer test.
Results
HT and HT stability of parental cultivars
The HT of ISH was always earlier than that of KSM in the 12 cropping seasons, irrespective of the sowing date (Supplementary Fig. S1). Correlation coefficients between sowing date and HT were insignificant, being r = 0.53 in KSM and r = 0.44 in ISH (Supplementary Fig. S2). The HT stabilities represented by interannual variance were 17.9 in KSM and 33.4 in ISH. ISH showed an early-heading phenotype and was more unstable in HT over the years.
Next, we genotyped the known HT-related genes for the two cultivars. KSM and ISH harbored homozygous recessive alleles of vrn-H1, vrn-H2, and vrn-H3, and homozygous dominant alleles of Ppd-H1 and Ppd-H2 (Supplementary Table S5). There was no polymorphism in HvCK2α between KSM and ISH. In contrast, the HvCEN and HvPHYC genotypes differed. KSM and ISH harbor early-heading and late-heading alleles at the HvCEN locus, respectively. The alleles of the HvPHYC locus were haplotype 7 in KSM and haplotype 4 in ISH, after Pankin et al. (2014). Hereafter, each allele of these two genes will be referred to as HvCEN-e, HvCEN-l, HvPHYC-e, and HvPHYC-l. In summary, it was revealed that the alleles of at least two known HT-related genes, HvPHYC and HvCEN, differed between the two cultivars.
QTL analysis of HT in RILs
QTL analysis was performed using 149 RILs. The average HT varied from −6.7 to 12.5 in 149 RILs of the F5 generation, −0.4 to 20.8 in the F6 generation, 0.4 to 25.4 in the F7 generation (late sowing in 2019/2020) (Supplementary Fig. S3a), and from −15.3 to 11.0 in the F8 generation (2020/2021) (Supplementary Fig. S3b). The average HTs of ISH and KSM were 10.4 and 15.4 in 2019–2020, and −7.2 and 1.8 in 2020–2021, respectively. In both generations, the frequency distributions of the HTs were continuous and showed clear transgressive segregation in the 149 RILs, suggesting the segregation of multiple genes.
To construct a linkage map, the pooled MIG-seq library of 149 RILs was sequenced, resulting in approximately 70 gigabases of raw read data volume. This data comprised 465,169,226 reads (an average of 3,121,941 reads per individual). After variant calling, 46,012,612 polymorphisms (42,463,288 SNPs and 3,549,324 InDels) were detected, and 45,954,293 polymorphisms were filtered by vcftools with the criteria described above. This filtered vcf file was used to construct a linkage map. The linkage map of RILs consisted of eight linkage groups. Chromosome 5H was separated into two linkage groups, 5H-1 and 5H-2. The total marker number of the linkage map was 6,865, with a total length of approximately 1,575 cM and an average distance between markers of approximately 0.2 cM (Supplementary Fig. S4 and Supplementary Table S6). The largest inter-marker distance was 16.6 cM, in linkage group 7H. The number of unique markers of the linkage map was 1,308. The order of the markers on the linkage map and the physical map corresponded well, and a suppression of recombination in the region near the centromere was also observed (Supplementary Fig. S5). These results indicate that the linkage map was successfully constructed by MIG-seq data, indicating that MIG-seq is effective for genome-wide polymorphism detection, not only in wheat but also in barley.
QTL analysis of HT detected three QTLs using SIM or CIM methods, designated QHd.ouj-2H, QHd.ouj-4H, and QHd.ouj-5H, in the F7 generation (Fig. 1a and Supplementary Table S7). QHd.ouj-2H (LOD = 4.4 in SIM method and LOD = 10.2 in CIM method) is located near the centromere of chromosome 2H, where HvCEN is located. In F7 generation, the peak SNP detected by CIM method mapped to the HvCEN. QHd.ouj-4H (LOD = 6.3 in SIM method) is located on the long arm of chromosome 4H. QHd.ouj-5H is located near the breakpoint of the linkage map of the long arm of chromosome 5H, and peaks of LOD scores were observed in both 5H-1 (LOD = 4.7 in SIM method and LOD = 8.2 in CIM method) and 5H-2 (LOD = 5.4 in SIM method and LOD = 5.6 in CIM method). The peak SNPs on 5H-1 in both F7 (SIM and CIM method) and F8 (SIM and CIM method) generations were located at the position of HvPHYC. Two genes, VRN-H1 and HvCK2α, were located in the gap region between 5H-1 and 5H-2. Two QTLs, QHd.ouj-2H and QHd.ouj-5H, were also detected in the F8 generation (Fig. 1b and Supplementary Table S7). The LOD score for QHd.ouj-2H was 12.8 in SIM method and in 13.4 in CIM method. For QHd.ouj-5H, the LOD scores on 5H-1 were 11.8 in SIM method and 12.7 in CIM method, and those on 5H-2 were 6.5 in SIM method and 6.8 in CIM method. Among these QTLs, QHd.ouj-2H and QHd.ouj-5H were commonly detected as significant QTLs in both generations. Tracing back across generations is influenced by heterozygosity, which leads to discrepancies between genotype and phenotype. As a result, the accuracy decreases; however, this approach, though primitive, was implemented. Nevertheless, QHd.ouj-2H and QHd.ouj-5H were detected in both the F5 and F6 generations, whereas QTL located on the short arm of chromosome 5H was detected only in the F5 generation (Supplementary Fig. S6).Fig. 1QTL likelihood curves of the LOD score of HT in F7 (a) and F8 (b) generations. Blue and red curves represent SIM and CIM method, respectively. The solid horizontal lines indicate the LOD significance thresholds for each method and generation: blue, 3.93 (F7) and 4.03 (F8); red, 6.59 (F7) and 6.32 (F8) for SIM and CIM, respectively
As mentioned above, parental cultivars of RILs are polymorphic at HvCEN and HvPHYC loci, which were located in the QTL regions of chromosomes 2H and 5H. Therefore, we assumed that HvCEN and HvPHYC were the causative genes for QHd.ouj-2H and QHd.ouj-5H, respectively. All subsequent analyses were performed on QTLs detected by SIM.
Evaluation of the effect of HT QTLs detected
HvCEN-e RILs headed significantly earlier by 3.4 (F7) or 3.5 days (F8) than HvCEN-l RILs (Supplementary Table S8). HvPHYC-e RILs headed significantly earlier by 5.3 (F7) or 6.2 days (F8) than HvPHYC-l RILs. The effects of these QTLs on HT were consistent with the results of the parental lines. Since HT-related genes are not known in the QHd.ouj-4H region, the QTL genotype was determined by the SNP with the highest LOD score in the F7 generation. KSM type homozygous RILs headed significantly earlier by 3.4 days in F7 generation than ISH type homozygous RILs (Supplementary Table S8), though the difference (1.7 days) was insignificant in F8 generation. From these results, the KSM type and ISH type alleles were designated QHd.ouj-4H-e and QHd.ouj-4H-l, respectively. A three-way ANOVA revealed no significant interaction between any of the QTLs in either generation (p > 0.4 in F7; p > 0.6 in F8).
HT stability in RILs
HT of KSM and ISK varied among 4 years, ranged from 1.8 to 15.4 and from −7.2 to 11.8, respectively. The HT stability represented by the variance of 4 year HT data was 47.1 and 62.8 in KSM and ISK, respectively, while it ranged from 19.8 to 90.0 among 191 RILs (Fig. 2). The continuous variation and transgressive segregation observed in RILs population indicated the involvement of multi-genes of which stable alleles dispersed between parents. Linear regression analysis after Finlay and Wilkinson (1963) clearly showed that HT of 4 generations of 191 RILs closely associated with overall average of HT in each generation (Supplementary Fig. S7a). The linear regression coefficient (β_1_) varied among RILs and ranged from 0.526 to 1.446 (mean = 1.000). The coefficient of determination ranged from 0.587 to 1.000 (mean = 0.946) among the RILs, and the correlation coefficient was insignificant in 121 RILs (df = 2, p < 0.01). Therefore, HT stability was assessed using variance rather than the regression coefficient. The correlation between variance and the regression coefficient was 0.967 (p < 0.01). Small regression coefficient indicated stable HT across the years. The linear regression coefficients were β_1_ < 0.783 (mean = 0.708) in the most stable 15 RILs (Supplementary Fig. S7b), while it was β_1_ > 1.155 (mean = 1.294) in the most unstable 15 RILs (Supplementary Fig. S7c). Thus, in this study, HT stability was evaluated by interannual variance instead of linear regression coefficient.Fig. 2. Frequency distribution of HT stability in 191 RILs. Blue and orange vertical lines indicate HT stabilities of ISH and KSM, respectively
Effect of three QTLs on HT stability
QTL analysis of HT stability detected no significant QTLs exceeding the threshold (Supplementary Fig. S8). Therefore, we investigated the effect of HT QTLs on HT stability. The average HT of the 139 RILs varied between cropping seasons, being 4.0 in 2015/2016 (F5), 9.7 in 2016/2017 (F6), 13.0 in 2019/2020 (F7), and −2.2 in 2020/2021 (F8). The interannual variance, reflecting the differing growing conditions, was 47.8. Parent–offspring correlation of HT was significant (p < 0.01), being r = 0.853 (F5–F6), r = 0.892 (F6–F7), and r = 0.910 (F7–F8). However, the correlation coefficient was not equal to 1, suggesting that responses to growing conditions and HT stability differed among RILs. An analysis of the relationship between parental HT and temperature across 11 growing seasons, excluding the 2019/2020 sown in late December, revealed that temperatures in mid-February exerted the strongest influence on HT. With respect to mid-February temperatures, an increase of 1 °C accelerated heading by 2.1 days in KSM and by 3.3 days in ISH (Supplementary Table S9 and Supplementary Fig. S9). Therefore, it was suggested that ISH is more sensitive to temperature fluctuations than KSM, resulting in a less stable heading time.
We first compared HT stability in 12 cropping seasons between parental cultivars. The HT stabilities were 17.9 in KSM and 33.4 in ISH, showing that ISH was more unstable (Supplementary Table S5). Similarly, in four cropping seasons (2015/2016, 2016/2017, 2019/2020, and 2020/2021), HT stability was 46.9 in KSM and 62.8 in ISH, confirming that ISH exhibited greater instability. To investigate the effect of heterozygosity at the three QTLs on HT, RILs of the F5 generation were genotyped for three QTL loci. In the F5 generation, the mean within-line variance of HT was 2.2 among lines heterozygous at any of the three QTLs, whereas it was 3.2 among lines homozygous at all three QTLs. Similarly, in the F6 generation, the mean within-line variance was 3.1 in heterozygous lines and 2.4 in homozygous lines; in the F7 generation, 3.2 and 2.2, respectively; and in the F8 generation, 4.1 and 2.3, respectively. As described above, the effect of heterozygosity on within-line variance of HT is not clear, and strong effect of ambient temperature was indicated. However, the F5 generation tended to show higher within-line variance compared with other generations, and within-line variance exceeded 10 in several RILs, presumably as a result of segregation of QTLs (Supplementary Fig. S10). Therefore, to eliminate the potential effects of heterozygosity on HT stability, only the RILs that were homozygous in the F5 generation were used for the subsequent analysis of HT stability. We then compared HT stability among RILs whose QTL genotype was the same as the parental cultivars, being 41.6 in the KSM type RILs (HvCEN-e QHd.ouj-4H-e HvPHYC-l) and 55.7 in the ISH type RILs (HvCEN-l QHd.ouj-4H-l HvPHYC-e) (Supplementary Table S10). From this result, it is possible that the HT QTLs detected in this study also affect HT stability.
Next, the effect of the single QTL on HT stability was evaluated. The average HT stabilities of RILs carrying HvCEN-e and HvCEN-l were identical, with both groups showing a mean value of 47.8 (Fig. 3a and Supplementary Table S11). In contrast, the average HT stability differed significantly depending on the QHd.ouj-4H genotype, being 45.1 and 50.3 in the RILs carrying QHd.ouj-4H-e and QHd.ouj-4H-l, respectively (p < 0.05; Fig. 3b and Supplementary Table S11). Similarly, the average HT stability differed significantly depending on the HvPHYC genotype, being 50.3 and 45.9 in the RILs with HvPHYC-e and HvPHYC-l, respectively (p < 0.05; Fig. 3c and Supplementary Table S11). Furthermore, significant interactions were detected for QHd.ouj-4H × HvCEN and HvPHYC × HvCEN (three-way ANOVA, p < 0.05). In the HvCEN-e background, QHd.ouj-4H had no significant effect on HT stability (Fig. 4 and Supplementary Table S11). Similarly, in the HvCEN*-l* background, HvPHYC had no significant effect on HT stability. These results suggest that the combination of QTLs is important for HT stability.Fig. 3. Allelic differences in HT stability at individual QTL loci (a–c) and in a combination of QTLs (d). (a–c) “*” and “n.s.” indicate significant (p < 0.05) and non-significant differences, respectively, according to the Tukey–Kramer test. “e” and “l” denote early and late genotypes, respectively. (d) Different letters indicate significant differences based on the Tukey–Kramer test at the 5% levelFig. 4Interactions between HvCEN and QHd.ouj-4H (a), and between HvCEN and HvPHYC (b), in HT stability
Finally, we analyzed the effects of interaction between three QTLs on HT stability, by comparing HT stability among RILs of eight QTL genotype combinations. The most significant difference was detected between HvCEN-l QHd.ouj-4H-e HvPHYC-e (39.2) and HvCEN-e QHd.ouj-4H-e HvPHYC-e (56.4) (Fig. 3d and Supplementary Table S10). These results suggest that the early alleles of HvCEN interact with QHd.ouj-4H to suppress HT instability caused by the QHd.ouj-4H-l allele. Likewise, the late alleles of HvCEN interact with HvPHYC to suppress HT instability associated with the HvPHYC-e allele.
Candidate genes of QHd.ouj-4H predicted by SnpEff
To investigate candidate genes of the novel QTL, QHd.ouj-4H, we used SnpEff to search for genes with a sequence polymorphism between KSM and ISH. We defined the QTL interval as the ± 5 cM region where the LOD value exceeded the threshold in the F7 generation. As a result of whole-genome resequencing of KSM and ISH, 93,676 polymorphisms (86,907 SNPs and 6,769 InDels) were detected in the QHd.ouj-4H interval of 8.1 Mb, and 38 polymorphisms (21 SNPs and 17 InDels) were considered HIGH annotation by SnpEff (Fig. 5). The 38 polymorphisms were located in 25 genes. We investigated the expression levels of these genes across 16 tissues of the Morex using the BarleyExpDB. Among them, 15 genes expressed in all developmental stages and tissues were identified as the most likely candidate genes for QHd.ouj-4H (Table 1 and Supplementary Table S12). Although these candidate genes have not been reported as HT-related genes, some are involved in development and growth in Arabidopsis. HORVU.MOREX.r3.4HG0414960 corresponded to Arabidopsis APYRASE 2 (AT5G18280), which encodes the enzyme with ATPase and ADPase activity.Fig. 5. Association between physical map and linkage map in the candidate region of QHd.ouj-4H. The red horizontal line indicates the threshold for QTL analysis, and the red vertical dotted line indicates genes annotated as HIGH by SnpEff (Cingolani et al. 2012b)Table 1. Candidate genes for QHd.ouj-4H identified by SnpEffGene IDPosition (bp)Variant annotationGene annotationArabidopsis gene IDHORVU.MOREX.r3.4HG0414960600,704,969splice_acceptor_variant&intron_variantApyraseAT5G18280**HORVU.MOREX.r3.4HG0415060600,879,254frameshift_variantrRNA N-glycosidaseNo hits foundHORVU.MOREX.r3.4HG0415100600,944,808frameshift_variantrRNA N-glycosidaseNo hits foundHORVU.MOREX.r3.4HG0415120601,039,430stop_gainedPol polyproteinNAHORVU.MOREX.r3.4HG0415250601,566,432frameshift_variant&start_lostPentatricopeptide repeat-containing proteinAT1G03100, AT1G04163**HORVU.MOREX.r3.4HG0415300601,627,867stop_gainedRNA-directed DNA polymerase-related family proteinNA601,628,086frameshift_variantHORVU.MOREX.r3.4HG0415360601,811,483frameshift_variant&stop_lost&splice_region_variantF-box family proteinAT5G02920601,811,523frameshift_variant601,811,532frameshift_variantHORVU.MOREX.r3.4HG0415720602,912,082stop_gainedEndoglucanaseAT1G64390602,912,240stop_gainedHORVU.MOREX.r3.4HG0415960603,674,507stop_gainedCCR4-NOT transcription complex subunit 1AT1G02080**HORVU.MOREX.r3.4HG0416570604,964,975start_lost&conservative_inframe_insertionMajor facilitator superfamily domain-containing proteinAT1G64650**HORVU.MOREX.r3.4HG0416670605,046,099frameshift_variantATPase, AAA family protein, expressedAT2G46620**HORVU.MOREX.r3.4HG0416940605,453,159splice_donor_variant&intron_variantChaperone DnaKAT5G09590605,454,521stop_gainedHORVU.MOREX.r3.4HG0417160605,875,982frameshift_variantFlavin-containing monooxygenaseAT5G61290**HORVU.MOREX.r3.4HG0417260606,085,065splice_acceptor_variant&intron_variantBeta-amylaseAT4G17090**HORVU.MOREX.r3.4HG0417290606,112,436frameshift_variant2-oxoglutarate (2OG) and Fe(II)-dependent oxygenase superfamily proteinAT4G10490**HORVU.MOREX.r3.4HG0417300606,121,464stop_gainedDihydroflavonol-4-reductaseAT1G09480NA indicates that no corresponding Arabidopsis gene was found in the TBLASTN results
Discussion
Although the importance of HT stability for stable production under the current climate change is widely recognized, genetic studies of this trait are fragmental in barley and wheat. By analysis of a diversity panel of barley, Sato et al. (2020) highlighted the contribution of vernalization requirement and PPD-H1 to HT stability. Using a diversity panel of wheat, He et al. (2024) revealed that the VRN and PPD genes affect HT stability under climate change scenarios based on Coupled Model Intercomparison Project Phase 5 (Lobell et al. 2015), through the phenological model. From a simulation-based analysis, they suggested that specific allelic combinations (Ppd-A1a Ppd-D1a Vrn-D1 or Ppd-A1b Ppd-D1b vrn-D1) caused an unstable phenotype of HT. Other allelic combinations (Ppd-A1a Ppd-D1a vrn-D1, Ppd-A1b Ppd-D1b Vrn-D1, or Ppd-A1b Ppd-D1a Vrn-D1) caused a stable phenotype of HT. They concluded that the combination of VRN and PPD genes, rather than individual VRN or PPD genes, plays a critical role in wheat flowering time stability. However, the materials used in the previous studies were not segregating populations, but diversity panels of barley (Sato et al. 2020) and wheat (He et al. 2024).
In this study, we first developed a segregating population (RILs) by crossing barley cultivars KSM and ISH, which were regarded as stable and unstable cultivars with HT stabilities of 17.9 and 33.4, respectively (Supplementary Table S5). Two parental cultivars commonly have a spring allele of vrn-H2 and early alleles of Ppd-H1 and Ppd-H2 (Supplementary Table S5). Therefore, RILs are not suitable for estimating the effect of VRN and PPD genes on HT stability, but are suitable for detecting QTLs related to HT stability in early-maturing cultivars with spring growth habit. By QTL analysis, using MIG-seq, three HT QTLs, QHd.ouj-2H, QHd.ouj-4H, and QHd.ouj-5H, were successfully identified, and their earliness effects were 3.4, 3.4, and 5.3 days, respectively (Fig. 1 and Supplementary Tables S7 and S8). HT QTLs with such a small effect could be used for fine tuning of HT. HvCEN and HvPHYC are considered as causal genes for QHd.ouj-2H and QHd.ouj-5H, respectively, since they were located in the QTL region and RILs with the respective early allele showed early heading (Supplementary Table S8), consistent with the results of Comadran et al. (2012), Nishida et al. (2013), and Pankin et al. (2014). No QTLs associated with HT stability were detected in our QTL analysis (Supplementary Fig. S8). This may be due to several factors. First, the 4-year dataset might not have been sufficient to capture the environmental variation affecting HT stability. Second, the population size of 149 RILs was relatively small for reliable QTL detection. Finally, since HT stability is a complex trait that may involve multiple loci and gene–gene interactions, it is likely difficult to detect QTLs associated with this trait using QTL analysis. Taken together, these limitations suggest that the stable detection of QTLs for HT stability will require evaluations across more diverse environments and a larger population size. These issues warrant further investigation in future studies. However, the significant involvement of HvPHYC and QHd.ouj-4H was indicated (Fig. 3 and Supplementary Table S11). Furthermore, interactions between HvCEN and QHd.ouj-4H, as well as between HvCEN and HvPHYC, were revealed (Fig. 4). HvCEN and HvPHYC have been reported as HT-related genes; however, their interaction contributing to HT stability was revealed for the first time in this study. The HvCEN-e allele suppressed the effect of QHd.ouj-4H-l on HT stability, whereas the HvCEN-l allele suppressed the effect of HvPHYC-e on HT stability (Fig. 4). Since early-maturing cultivars are cultivated to avoid abiotic stresses such as high temperature, drought, or monsoon rains during the grain filling period, these results facilitate the breeding of barley cultivars under the changing global climate.
Among the three QTLs, QHd.ouj-4H was detected on chromosome 4HL, where VRN-H2 for vernalization requirement (Yan et al. 2004; Karsai et al. 2005) and EPS-4L for earliness per se (Laurie et al. 1995) were located. Laurie et al. (1995) located EPS-4L near RFLP marker Xpsb37, which is closely linked with Ebmac0635 (Hv-Consensus2006-Marcel-4H). Based on the primer sequences, Ebmac0635 was located at Chr4H_580,928,328–580,928,507 of MorexV3_pseudomolecules_assembly. Although Morex lacks a VRN-H2 sequence, the physical position of VRN-H2 deletion was estimated near Chr4H_604,202,142–604,216,997 by physical mapping of the sequence of flanking marker SNF2P (Dubcovsky et al. 2005). Since QHd.ouj-4H was detected at Chr4H_619,974,683, EPS-4L and VRN-H2 would be excluded from the candidate of QHd.ouj-4H. In addition, a diagnostic PCR assay and short-read mapping against the Akashinriki genome (Jayakodi et al. 2020), which contains a VRN-H2 sequence, showed complete deletion of the VRN-H2 sequence in both KSM and ISH (Supplementary Table S5). Because of this deletion, both cultivars exhibit spring growth habit. Therefore, VRN-H2 cannot be the causal gene for QHd.ouj-4H.
From the results of SnpEff analysis, TBLASTN searches using whole-genome resequencing data, and expression profiles obtained from BarleyExpDB, 15 genes were considered candidate genes for QHd.ouj-4H (Table 1). Among them, HORVU.MOREX.r3.4HG0414960 is homologous to the Arabidopsis gene APYRASE 2 (APY2). In Arabidopsis, APY2:GUS expression is not detected in light-grown seedlings, suggesting that APY2 expression is regulated by light (Weeraratne et al. 2022). Moreover, polar auxin transport is significantly reduced in apy2 null mutants upon estradiol-induced APY1 RNA interference (Weeraratne et al. 2022). Liu et al. (2012) proposed that APY could alter the post-transcriptional activity of auxin transporters. Auxin has been shown to play a pivotal role not only in flower development, but also in the regulation of flowering time across diverse plant species. In Arabidopsis, it has been proposed that auxin homeostasis mediates the transition from floral meristem termination to gynoecium development (Yamaguchi et al. 2017). AUXIN RESPONSE FACTORS (ARFs) regulate auxin-mediated transcriptional activation and repression, and arf mutants exhibit delayed flowering (Okushima et al. 2005). In cotton, GhAUX22D interacts with GhFT, and auxin enhances the formation of the GhAUX22D–GhFT complex (Jian et al. 2025). Thus, HORVU.MOREX.r3.4HG0414960 might represent a particularly compelling candidate for the causal gene underlying QHd.ouj-4H. Identification of the causal gene for QHd.ouj-4H could deepen our understanding of HT stability and have practical implications for barley breeding. However, we cannot exclude the possibility that another gene located in the flanking region may be responsible. In addition, a chromosomal inversion was suggested within the candidate region (Fig. 5), and further fine-mapping will be essential. Based on this hypothesis, we propose a model in which QHd.ouj-4H contributes to HT stability. In several plant species, auxin transport velocity is influenced by temperature—it increases with rising temperature up to an upper limit (Morris 1979). The HvAPY2 allele in KSM might alter the post-transcriptional activity of auxin transporters, potentially rendering them temperature-insensitive or less responsive to temperature fluctuations. Consequently, auxin transport activity may remain nearly constant across temperature conditions, allowing KSM to maintain auxin transport velocity even at low temperatures, thereby achieving both earliness and stability. HvCEN, the causal gene for EPS-2S, affects the duration of the Z30–Z49 phase, which corresponds to the stem elongation stage (Castro et al. 2017). Wolbang et al. (2004) suggested that IAA plays an important role in GA biosynthesis and, consequently, in stem elongation. Taken together, these findings might suggest that the interaction between HvCEN and QHd.ouj-4H could affect stem elongation. Furthermore, HORVU.MOREX.r3.4HG0414960 exhibited higher expression in the vegetative tissue LEA (leaf tissue sampled from seedlings at 17 days after planting, dap), whereas its expression was reduced in the inflorescence tissues INF1 and INF2, which correspond to whole developing inflorescences sampled at 30 dap and 50 dap, respectively (Supplementary Table S12). This expression pattern supports the hypothesis that the interaction between HvCEN and QHd.ouj-4H contributes to the regulation of stem elongation. In Arabidopsis, ambient temperature regulates flowering through a pathway involving TFL1, and the delayed flowering under low temperature is partially alleviated in tfl1 mutants (Strasser et al. 2009). In addition, the expression of PHYC in bread wheat is affected by ambient temperature (Kiss et al. 2025). Integrating these findings, we propose that it is plausible that the interaction between HvCEN and HvPHYC identified in this study is modulated by ambient temperature. Further investigation will be necessary to clarify the molecular mechanisms by which temperature regulates this interaction.
Previous studies employed restriction enzyme-based methods such as GBS and RAD-seq for genome-wide genotyping of barley. The total number of markers were 1,391, 3,662, and 2,406 to 15,161 for GBS (Liu et al. 2014; Yao et al. 2018; Abed et al. 2021) and 12,998 and 1,894 for RAD-seq (Zhou et al. 2015; Wang et al. 2016), respectively. In this study, we employed MIG-seq (Nishimura et al. 2022) and constructed a linkage map with an average marker distance of 0.2 cM, using 6,865 markers (Supplementary Table S6). In this study, the linkage group of chromosome 5H was separated into two segments, 5H-1 and 5H-2. HvPHYC was located within the gap between 5H-1 and 5H-2 and was not included in the linkage map. This gap spanned approximately 16.0 Mb in physical position and 38.9 cM in genetic distance. Moreover, the read depth in this gap region was 18.9, which was lower than the average read depths in 5H-1 (25.2) and 5H-2 (31.5), resulting in the separation of the linkage maps. Although the number of polymorphic markers depends on the parents of the mapping population, and simple comparisons between studies are difficult, these results indicated that the marker density of this study was equivalent to that of restriction enzyme-based methods such as GBS and RAD-seq. In addition, MIG-seq is inexpensive, and it is also easier to construct libraries from low-quality DNA because it does not use restriction enzyme treatment in the library construction steps. Therefore, the use of MIG-seq will facilitate future genotyping of barley using NGS. It also may be used for linkage map construction and QTL analysis, and for efficient selection during backcrossing in barley, as suggested by Nishimura et al. (2024) in wheat.
To understand the environmental factors caused difference in HT and HT stability, the relationship between HT and average temperature was analyzed using 12 years data of KSM and ISH. As a result, average temperature of February, more specifically of mid-February (11th–20th), related closely with HT, indicating that HT variation could be explained by temperature variation in mid-February in barley cultivars insensitive to vernalization and short day-length. With excluding the data of 2019/2020 with extreme late sowing, correlation coefficient between HT and average temperature of mid-February was −0.768 in KSM (p < 0.01) and −0.860 in ISH (p < 0.01) (Supplementary Fig. S9), while correlation was insignificant in other 10 days during barley growing season (Supplementary Table S9). According to Yasuda (1984), vernalization requirement of fall sown barley was satisfied in the field by 15th December in spring type cultivars, and reached double-ridge stages by 16th January. Thereafter, the development became slow or rested, probably due to low temperature and/or short photoperiod. In middle of February, spike development and stem elongation started rapidly with an increase in temperature and photoperiod. KSM and ISH are spring type carrying vrn-H2 (Supplementary Table S5), and thus vernalization is not necessary for transition to reproductive growth. KSM, ISH, and the RILs were considered to reach early double-ridge stage by 5th January and to start rapid spike development in mid-February in the normal years. However, as shown in Supplementary Fig. S9, thermo-sensitivity seems different among KSM and ISH, causing the difference in the timing of the onset of rapid spike development and in HT. Furthermore, in Populus euphratica, PeAPY2 expression increased in response to cold treatment and peaked on day 7 of exposure (Deng et al. 2015). This cold-responsive expression pattern suggests that temperatures prevailing in mid-February may affect the expression of the gene underlying QHd.ouj-4H. Further study is required to clarify the relationship between HT and temperature and/or photoperiod in winter.
Supplementary Information
Below is the link to the electronic supplementary material.Supplementary file1 (PPTX 3054 KB)Supplementary file2 (XLSX 69 KB)
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Casao M, Karsai I, Igartua E, Gracia M, Veisz O, Casas A (2011) Adaptation of barley to mild winters: a role for PPDH 2. BMC Plant Biol 11:164. 10.1186/1471-2229-11-164