Identification of the vernalization gene VRN-B1 responsible for heading date variation by QTL mapping using a RIL population in wheat

Background Heading time is one of the most important agronomic traits in wheat, as it largely affects both adaptation to different agro-ecological conditions and yield potential. Identification of genes underlying the regulation of wheat heading and the development of diagnostic markers could facilitate our understanding of genetic control of this process. Results In this study, we developed 400 recombinant inbred lines (RILs) by crossing a γ-ray-induced early heading mutant (eh1) with the late heading cultivar, Lunxuan987. Bulked Segregant Analysis (BSA) of both RNA and DNA pools consisting of various RILs detected a quantitative trait loci (QTL) for heading date located on chromosomes 5B, and further genetic linkage analysis limited the QTL to a 3.31 cM region. We then identified a large deletion in the first intron of the vernalization gene VRN-B1 in eh1, and showed it was associated with the heading phenotype in the RIL population. However, it is not the mutation loci that resulted in early heading phonotype in the mutant compared to that of wildtype. RNA-seq analysis suggested that Vrn-B3 and several newly discovered genes, including beta-amylase 1 (BMY1) and anther-specific protein (RTS), were highly expressed in both the mutant and early heading pool with the dominant Vrn-B1 genotype compared to that of Lunxuan987 and late heading pool. Enrichment analysis of differentially expressed genes (DEGs) identified several key pathways previously reported to be associated with flowering, including fatty acid elongation, starch and sucrose metabolism, and flavonoid biosynthesis. Conclusion The development of new markers for Vrn-B1 in this study supplies an alternative solution for marker-assisted breeding to optimize heading time in wheat and the DEGs analysis provides basic information for VRN-B1 regulation study.


Background
Wheat is one of the most widely cultivated food crops worldwide, and successful adaptation of this important grain crop in different agro-ecological conditions is mainly determined by flowering and maturity time. Heading date is tightly connected with crop maturity and production. Identification and isolation of genes related to heading date variation will provide better understanding of the genetic pathways that control flowering time in plants and offer effective genetic resources for engineering high-yield varieties that can adapt to various conditions [1].
Plants have evolved multiple genetic pathways to regulate the flowering pathway, including vernalization requirements and photoperiod sensitivity, integrating both external and internal signals to adapt to changes in climatic conditions. Vernalization is the requirement of exposure to prolonged periods of low temperature that provides competence to flower in many plant species adapted to temperate climates [2]. In wheat, the vernalization-induced flowering pathway is mainly regulated by the VERNALIZATION gene series: VRN1 [3], VRN2 [4], VRN3 [5] and VRN-D4 [6]. Identification and characterization of these genes furthered our understanding of vernalization-mediated flowering regulatory networks in cereals. The MADS-box transcription factor VRN1 gene, was first identified in Triticum monococcum and functions in acceleration of flowering after vernalization [3]. VRN2 encodes a zinc-finger-CCT domain transcription factor and represses the flowering in cereals [4,7,8]. The CCT domain in VRN2 is a key element for vernalization in both wheat and barley, and mutation of the CCT domain eliminates a vernalization requirement [9,10]. VRN3 is homologous to the Arabidopsis photoperiod gene FLOW-ERING LOCUS T (FT) and its gene product functions as a mobile signaling protein, moving from leaves to the shoot apical meristem to accelerate flowering [5]. VRN-D4 originated from an insertion of a~290 kb region from chromosome arm 5AL into the proximal region and contains a copy of VRN-A1 [6].
VRN1 has been well reported in wheat and is orthologous to the Arabidopsis floral meristem-identity genes CAULIFLOWER (CAL), APETALA1 (AP1), and FURIT-FULL (FUL) [3,[11][12][13][14]. In hexaploid wheat, a homoeologous copy of the VRN1 gene is present in the proximal regions of chromosomes 5A, 5B, and 5D, and have been designated VRN-A1, VRN-B1 and VRN-D1, respectively. Allelic variation at the VRN1 locus is one of the main resources of genetic variation in vernalization requirements in wheat. A single dominant allele of one of the three homeoloci is sufficient to confer a spring growth habit [15][16][17]. A previous study showed that the dominant Vrn1 with spring growth habit contained large deletions in the first intron, and a 2.8-kb conserved segment within the deletion region of recessive vrn1 was important for vernalization requirement in wheat [16]. Additionally, the investigation of polymorphisms in the RNA immune precipitation fragment 3 (RIP3) of VRN-A1 first intron suggested that single nucleotide polymorphisms in RIP3 affected VRN-A1 transcript level, and are associated with differences in heading time of winter wheat cultivars from various geographic regions [17]. The differences between dominant and recessive alleles associated with the VRN-B1 locus lie near the first intron and mainly include a large deletion, SNPs, small deletions, and large insertion in the 5′-UTR region [16,[18][19][20][21][22].
The VRN1-VRN2-VRN3 regulatory feedback network integrates vernalization and photoperiod signals to precisely mediate the flowering time of wheat. In barley, high levels of VRN2 repress the expression of VRN3, inhibiting flowering. In wheat, the up-regulation of VRN1 transcriptional levels results in decreased expression levels of flowering repressor VRN2 after vernalization, which accelerates the induction of VRN3 expression levels, and thus induces flowering in the spring [23,24]. Epistatic analysis indicated that VRN2 interacts with VRN1 and shares a genetic pathway [25]. VRN3 up-regulates the expression of VRN1 through interactions with FDL2, a basic leucine-zipper (bZIP) transcription factor [26]. FDL2 is orthologous to the bZIP transcription factor FD in Arabidopsis that binds ACGT elements in the promoter of VRN1 [26,27]. The expression of VRN3 can also be affected by competition with CO and VRN2 to interact with the nuclear factor-Y (NF-Y) complex, which plays an important role in the integration of vernalization and photoperiod seasonal signals [28].
Many QTL have been identified for controlling heading time in wheat. Until now, more than 240 QTL related to heading time have been detected in wheat [29,30], and are important for marker-assisted selection in wheat breeding. However, most of genes underlying these loci remain unknown. In this study, we identified a QTL for heading time located on chromosome 5B by Bulked Segregant Analysis (BSA) and genetic linkage mapping, and showed that the mapped region contained the VRN-B1 gene that is involved in heading time variation in the RIL population but it is not the reason for early heading phenotype in the mutant (eh1) compared to that of wildtype (WT, Zhongyuan9). In order to identify genes regulated by VRN-B1, we identified genes that were differentially expressed between the early and late heading bulks, and investigated their associated pathways. This study provides important clues for marker assisted breeding for heading time in wheat and expands our knowledge on genetic regulation of the VRN-B1 gene.

Construction of the RIL population and heading time variation
The early heading mutant eh1 was obtained by mutagenesis of seeds of Zhongyuan9 wheat lines with 284 Gy γ-ray treatment. To map the QTL responsible for heading date (HD) variation, we developed 400 recombinant inbred lines (RILs) by crossing eh1 with Lunxuan987 (LX987), an elite cultivar with a minor defect of late heading for agricultural practice. The phenotype comparison between LX987 and eh1, WT and eh1 was shown in Fig. 1a. The HD of eh1 was 10-14 days earlier than WT over 3 years under field conditions (Fig. 1b, Table S1). The HD of LX987 is 4-8 days later than that of eh1 with significantly statistical difference (Fig. 1b, Table S1). Measurement of HD in RILs over a period of 3 years suggested that the maximum HD is 215 days and the minimum HD is 197 days (Tables S1, S2). Additionally, HD in the RIL population showed a continuous distribution (Fig. S1), indicating the flowering time in the RIL is a quantitative trait.
Mapping of an early heading gene using a BSR-seq approach Bulked Segregant RNA-Seq (BSR-Seq) is an efficient approach for gene mapping [31], especially in wheat plants that have large genomes. We performed BSR-Seq to map the QTL controlling HD in the RIL population. Young spikes of 21 early heading and 22 late heading lines, based on HD data collected over a 3 year period, were sampled when the early heading lines started to heading. The RNA of individual samples was isolated and then early or late heading RNA bulks were constructed. Meanwhile, young spikes of six individuals for two parent lines were also sampled for RNA-Seq analysis. In total, 357,941 SNPs were obtained from two parents and two bulks by RNA-Seq. A total of 160 SNPs with multiple alleles and 180,469 loci with reads < 4 were filtered out. Sixty-nine thousand two hundred twenty-two SNPs with the same genotype and 61,288 SNPs from bulks that differed from those in corresponding parents were also discarded. Finally, 46,802 reliable, high-quality SNPs distributed across all chromosomes were identified and used to perform further analysis.
The Δ (SNP-index) was calculated based on the reads and the genotypes of bulks and parents. The results showed that the region of 546.9 Mb -636.6 Mb on chromosome 5B with a SNP-index peak was the major QTL associated with HD in the RIL population (Fig. 2a). Calculation of the Euclidean Distance (ED) value was also used to identify loci controlling HD, and suggested that a total of 1113 genes were present in an 89.57 Mb region on chromosome 5B according to the gene information of Chinese Spring reference sequence (http:// www.wheatgenome.org/) (Fig. 2b). Both SNP-index and ED analyses detected a common interval on chromosome 5B, suggesting that an 89.57 Mb region of chromosome 5B harbors the QTL responsible for heading variation in the RILs.
The same region on chromosome 5B associated with HD was detected by BSA consisting of various RILs Based on 2 years phenotype data, three extremely early heading bulks consisting of 18-22 early heading lines for each bulk, and three extremely late heading bulks consisting of 16 late heading lines, respectively, were used for BSA analysis. A total of 6 bulks and 2 parents were genotyped by using 660 K SNP array. The SNPs showing different genotypes in a pair of early and late heading bulks, and the same genotypes between bulk and parent with the same phenotype, were selected and the SNP Values are means ± SD and different letters indicate the significant difference between the comparison groups at P < 0.05 numbers among 21 chromosomes were counted. The results showed that a large number of SNP were detected on chromosomes 2A and 5B between the early and late heading bulk 1 (Fig. S2A), on chromosomes 5B, 2B and 3A between the early and late heading bulk 2 ( Fig. S2B), on chromosomes 3B, 5B, 2B and 6A between the early and late heading bulk 3 (Fig. S2C), suggesting QTLs controlling HD exist on these chromosomes. The SNP frequencies distributed on chromosome 5B were analyzed. High peak regions from 540 to 600 Mb on chromosome 5B overlapped in the three pairs of bulks were observed (Fig. S2D), which is consistent with the BSR-Seq analysis.

Confirmation of the region by genetic mapping
To validate the BSA-based mapping region, 87 RNA-Seqderived SNPs on chromosome 5B were selected for development of KASP markers. Specific SNPs on the three wheat genomes were selected using PolyMarker [32], confirmed on parents and population lines, then 20 specific KASP markers were finally developed (Table S3). The genotypes of 400 lines in the RIL population were identified. We generated a linkage map spanning 112.83 cM on chromosome 5B, with an average genetic distance of 5.6 cM, using QTL IciMapping (Fig. 3). Based on the combined genetic and phenotypic data generated over 3 years, a stable QTL with a LOD score of more than 22.8 was detected in the region between the two markers, ch14 and ch8, with a genetic distance of 3.31 cM, corresponding to the physical position of 560.4-580.1 Mb (Fig. 3).
The Vrn-B1 is responsible for HD variation in the RILs The vernalization gene VRN-B1 located on chromosome 5B is a major locus associated with wheat heading and flowering [3]. In natural variants, the dominant Vrn-B1 allele with spring growth habit differs from the recessive allele vrn-B1 with winter growth habit due to a large deletion in the first intron of Vrn-B1 [16]. We found that the mapped region on chromosome 5B contained the VRN-B1 gene. To test whether the major QTL on chromosome 5B controlling HD is associated with the Vrn-B1 gene, we used the markers developed by previous study [16] but failed to get the PCR product in WT and eh1. Therefore, specific markers for Vrn-B1 were developed (Fig. S3A, Table S3). For the VRN-B1 gene markers, two pairs of genome-specific primers were designed to determine whether a large deletion in the first intron of Vrn-B1 was absent or present (Fig. S3A). Amplification products were detected both in eh1 and WT, but were not produced using DNA isolated from LX987 using the TavBI primer pair. However, a PCR product was amplified from LX987 DNA with the TavBII primer pair (Fig. 4a), indicating that the large deletion exists in   Values are means ± SD. Different letters indicate statistically significant difference among the comparison groups at P< 0.05 eh1 and WT but does not in LX987. The genotypes of VRN-B1 in the eh1, WT and LX987 suggested that the genotypic differences observed in VRN-B1 in the RIL population resulted from genetic background diversity between WT and LX987, but not from the γ-ray irradiation. We then identified the genotypes of VRN-B1 in the RIL population, and showed that 199 RILs with the Vrn-B1 allele displayed significantly (P < 0.001) earlier heading dates than 185 RILs that possessed the vrn-B1 allele. The Vrn-B1 lines headed approximately 2.7 days earlier than vrn-B1 lines on average (Fig. 4b). In addition, the heading time of the heterozygotes was 206 days, and the degree of dominance is 0.28, indicating that Vrn-B1 has a small dominant effect. Further sequencing the VRN-B1 in WT, eh1 and LX987 found that the A changed to G in the 3′ end and C changed to A in the middle of primer Intrl/B/F [16] in WT and eh1 (Fig. S3A, B). In addition to the large deletion in the first intron of WT and eh1, a 37 bp deletion located downstream of the large deletion was only detected in eh1 (Fig. S3A). We then developed specific markers for the 37 bp deletion (Fig. S3C) and identified the genotypes in the RIL population, and found that the 37 bp deletion co-segregated with the large deletion. Further analysis using an F 2 population containing 1060 individuals derived from cross of eh1 and WT indicates that the 37 bp deletion is not associated with HD variation (Fig. S4, P = 0.757).
To identify if natural variation of VRN-A1 existed in eh1 and LX987 affecting heading time, we used previously developed markers to distinguish dominant Vrn-A1 and recessive vrn-A1 [16]. The results showed that both eh1 and LX987 harbored the recessive vrn-A1 allele (Fig. S5A). Ppd-A1 is an important regulator of heading time that located on chromosome 2A. To determine the sequence variations of Ppd-A1 gene between eh1 and LX987, we sequenced the Ppd-A1 gene as well as its upand down-stream region. Sequence comparison showed no variations in the coding region of Ppd-A1 gene among eh1, WT and LX987 (Fig. S5C). A 131 bp deletion on the 3′ region of Ppd-A1 was detected in eh1 and this variation also existed in Chinese Spring (Fig. S5B, C). We also compared the heading time of lines carrying different variations at this locus in the RIL population and found that the 131 bp deletion did not affect heading time (Fig. S5D).
The difference of heading time between eh1, WT, and LX987 under vernalization or non-vernalization conditions To investigate whether the mutated early heading gene is related to vernalization, the heading time of eh1, WT and LX987 lines was compared starting with seeds that were or were not subjected to vernalization. When seeds were sown in autumn (vernalized during winter), during 3 years investigation the eh1 headed 10-14 days and 4-8 days earlier than WT and LX987 with significantly statistical difference, respectively [22]. When seeds were sown in spring without vernalization, LX987 failed to flower and eh1 headed approximately 13 days earlier than WT (Table 1, Fig. 5). A similar variation in heading date between eh1 and WT, with or without vernalization, suggested that the early heading phenotype observed for eh1 did not depend on vernalization.

Identification of differentially expressed genes (DEGs) by RNA-seq
The genotype of Vrn-B1 in RILs suggested that all the RILs used for the early heading pool harbored a dominant allele (Vrn-B1) and 21 of the 22 RILs used for the late heading pool were recessive (vrn-B1) ( Table S4). RNA-Seq data, to some extent, can be reliably used to investigate the transcriptomic changes due to the presence of different VRN-B1 alleles. Based on our criteria (fold change ≥2 and FDR < 0.05), a total of 9440 DEGs were identified between the two parents, in which 6078 were up-regulated and 3362 were down-regulated (Table S5). A total of 8571 DEGs, including 5860 up-regulated and 2711 downregulated genes, were identified between early and late heading pools (Table S5). Among all the DEGs, 4169 genes overlapped between the comparison groups of two parents and two bulks, including 3263 up-regulated and 906 down-regulated genes (Fig. 6a, Table S5). Among them, 103 up-regulated genes, including beta-amylase 1 (BMY1, TraesCS2B01G240100), 3-oxoacyl-[acyl-carrierprotein] synthase I (KAS12, TraesCS2B01G156800) and anther-specific protein (RTS, TraesCS3D01G436700), with a fold change > 100 and P < 0.0001 were identified (Table  S5). Additionally, we found the vernalization gene VRN-B3 (TraesCS7B01G013100) was up-regulated (P < 0.01) in both eh1 and the early heading bulk sample compared to LX987 and late heading bulk with a fold change of 4.63 and 3.36, respectively.

Expression comparison of selected DEGs and vernalization genes by RT-qPCR
To verify the reliability of gene expression levels generated by RNA-seq, nine genes including the Vrn-B3, were selected for validation using real-time quantitative PCR analysis. As shown in Fig. 7, the relative expression levels of selected genes by RT-qPCR analysis were similar with the expression trends observed in the RNA-seq data, with a high Pearson's correlation coefficient (R 2 = 0.79).
To compare the expression levels of vernalization genes between eh1 and LX987, RT-qPCR were used to quantify the relative expression of VRN-A1, VRN-B1, VRN-D1, and VRN2 genes. When eh1 started to heading, young spikes and leaves of eh1 and LX987 from 10 replicates were sampled, respectively. The expression levels of VRN-A1, VRN-B1 and VRN-D1 were significantly higher in the young spikes of eh1 compared to that of LX987 (Fig. S6). In the leaves, the expression amount of VRN-B1 and VRN-D1 were also significantly higher in eh1 than in LX987 (Fig. S6). In contrast, the expression levels of VRN2 gene in both young spikes and leaves of eh1 were significantly lower than that of LX987 (Fig. S6). To analyze the expression levels of these genes in eh1 and LX987 at the same developmental stage, young spikes and leaves of eh1 were also sampled 8 days earlier than that of LX987 when they showed similar developmental levels of spikes. A completely opposite results were obtained for the samples from the same developmental stage compared to that sampling at the same day (Fig. S7).

GO and KEGG pathway enrichment analysis
Gene ontology (GO) enrichment analysis was performed to investigate the functions of overlapped DEGs between two parents and two bulks. In total, 1537 GO categories were enriched and 42 of them were significantly enriched based on the corrected P-value < 0.05 (Table  S6). The top 20 significantly enriched GO terms suggested that biological processes including photosynthetic electron transport in photosystem I, reductive pentosephosphate cycle, carbohydrate metabolic process and fatty acid biosynthetic process, were affected by changes in heading time resulting from VRN-B1 variation (Table  S6). KEGG pathway enrichment of the overlapped DEGs indicated that significantly enriched pathways related to photosynthesis, including antenna proteins, fatty acid elongation, carbon fixation in photosynthetic organisms, starch and sucrose metabolism, carbon metabolism, and fatty acid metabolism, are involved in heading process (Fig. 6b, Table S7).

Discussion
Variability of HD is of great significance for the adaptation of wheat cultivars to local environments, and its optimization assists yield improvement during the wheat breeding process [33]. Induced mutation is an effective approach for creation of new germplasm [34]. The γ-ray induced early heading wheat mutant, eh1, used in this study provides novel genetic variation for HD. By using eh1 and the late heading cultivar LX987 as parents, a RIL population was constructed for BSA and genetic mapping analysis. By combining BSA and genetic mapping approaches, a stable QTL for HD with a LOD score exceeding 22.8 was detected, and flanking markers were located on a 20 Mb interval that was found to include the VRN-B1 locus [3]. Further analysis using two genome-specific primers suggested that the presence or absence of a large deletion in the first intron of VRN-B1 was associated with heading time variation in the RIL population. Consistent with these results, a previous study suggested that large deletions within the first intron of VRN1 resulted in a spring growth habit and affected flowering time [16]. Additionally, a 2.8 kb region in the first intron of vrn-B1 present in different recessive alleles showed high sequence conservation [16,35]. It was assumed that this critical region contained a binding site for a putative repressor that is down-regulated by vernalization; thus, the presence of large deletion in the first intron of Vrn-B1 ablates suppressor function on VRN-B1 [16,22,35]. Due to the different genotypes of VRN-B1 between the two parents of the RIL population in this study, it is reasonable to identify the major QTL (Figs. 2, 3). By using the two pairs of Vrn-B1 gene markers, we found that the lines used for construction of the early and late heading bulks for BSR-Seq were almost divided equally between Vrn-B1 and vrn-B1, respectively (Table S4), indicating that the results from BSR-Seq are reliable. The same genotypes of large deletion in the first intron of VRN-B1 between WT and eh1 (Fig. 4a), combined similar variations of HD between WT and eh1 with or without vernalization (Fig. 5, Table  1), indicated that the early heading phenotype in eh1 is not due to mutation of large deletion in VRN-B1. In addition, due to the deletion in the first intron of VRN-B1, the eh1 and WT was able to flower without vernalization (Table 1, Fig. 5). Since the variations in WT and eh1 at the position of primer Intrl/B/F [16] Fig. 6 Analysis of DEGs. a Venn diagrams illustrating the overlap in DEGs between the two heading bulks and parents. B1 and B2 indicate the early heading bulk and the late heading bulk, respectively, while P1 and P2 represent the early heading mutant and LX987, respectively. The number in each circle represents the total number of DEGs in each comparison group, and the number in the overlapping areas indicates the number of shared DEGs between two comparison groups. b KEGG pathways enriched in the DEGs overlapped in the B1_vs_B2 and P1_vs_P2 groups. The number of DEGs in each listed pathway is indicated by size the circle area, the rich factor reflects the degree of enriched DEGs in a given pathway and the circle color represents the range of the corrected P values leaded to failed amplification (Fig. S3A, B), the new markers developed in this study would be an alternative solution for marker-assisted breeding in wheat. The sequencing of VRN-B1 gene suggested that a 37 bp deletion located downstream of the large deletion was mutated in eh1 (Fig. S3A, C) and this 37 bp deletion cosegregated with the large deletion in the RIL population. The 37 bp deletion is similar to a 36 bp deletion region of previously reported Vrn-B1b, which was present in the spring wheat cultivar [18]. It has been reported that accessions containing the Vrn-B1b allele exhibited longer days to heading than those carrying the Vrn-B1a alleles [36]. In contrast, by phenotype analysis of F 2 individuals derived from cross of eh1 and WT, we found that the 37 bp deletion didn't result in HD variation (Fig.  S4). Since only 2 accessions carrying Vrn-B1b allele were used for comparison by Cho et al. [36], it is possible that the background differences of other loci affecting heading time among various accessions leaded to the detection of HD variation.
Due to the limited lines used for BSR-Seq, only one QTL located on chromosome 5B was detected by using this approach (Fig. 2). When we used different early heading or late heading lines to construct three pairs of bulks for BSA analysis, QTLs located on chromosomes 5B, 2A, 2B, 3A, 3B and 6A were detected ( Fig. S2A-C). The QTL on chromosome 5B detected by three pairs of bulks also included the VRN-B1 gene (Fig. S2D), which is in consistent with the results from BSR-Seq and genetic mapping (Figs. 2, 3). Previous studies also found QTLs associated with HD located on chromosomes 5B, 2A, 2B, 3A and 3B by genetic mapping [37][38][39] and the flowering regulator photoperiod genes Ppd-A1 and Ppd-B1 in wheat are located on chromosomes 2A and 2B, respectively [40]. Although a 131 bp deletion on the 3′ region of Ppd-A1 gene was detected in eh1 compared to LX987 (Fig. S5C), the locus detected on chromosome 2A was not Ppd-A1 gene due to no difference of heading date between RILs with and without this deletion (Fig.  S5D). Additionally, genome-wide association analysis also detected the QTL for HD on chromosome 6A [41,42]. It is possible that the genetic background diversity for HD between WT and LX987 leads to detection of several QTLs in the RIL population. Since different RILs contain distinct genes for HD, using different RILs for BSA analysis would be a feasible and comprehensive  Table S7. The RNA-seq values (blue) represent the ratio of expression fold change in eh1 relative to LX987. Data generated by qPCR (orange) are shown as means ± SD from three biological replicates and different letters indicate significant differences of gene expression levels between eh1 and LX987 by qPCR at P < 0.05 approach for QTL identification in such population. However, since natural variations for HD existed between eh1 and LX987, it is difficult to distinguish the mutation loci from these QTLs. This problem would be solved by further fine mapping or using the population derived from cross of mutant and WT.
Together, several studies have elucidated the genetic and phenotypic effects of different variations in the VRN1 locus over the past several decades. However, the impact of VRN1 deletions on global gene expression has not been widely reported [42]. Using BSR-Seq data, differentially expressed genes (DEGs) in young spikes were identified between two pools comprised of those with either the Vrn-B1 or vrn-B1 genes (Table S4). Although the DEGs detected in the RNAseq analysis might not be solely due to the large deletion in VRN-B1, the newly discovered genes and pathways by expression pattern analysis would be partly regulated by VRN1 or affected by different flowering time. One limitation is unreplicated data for RNA-seq and it would be partly offset by analysis of the overlapped DEGs between parents and bulks. It has been reported that the Vrn1 protein binds to the promoter of FT1/ Vrn3 by chromatin immunoprecipitation sequencing (ChIP-seq), and transcript levels of FT1/ Vrn3 were upregulated by high expression of Vrn1 in barley [43]. VRN1 also down-regulates the flowering repressor VRN2 gene [24]. Consistently, higher expression levels of VRN-A1, VRN-B1 and VRN-D1 in the spikes of eh1 compared to that of LX987 were detected while lower amounts of VRN2 were observed in eh1 (Fig. S6). Meanwhile, increased expression of Vrn-B3 was detected in eh1, or early heading bulks that had the Vrn-B1 genotype relative to those within the vrn-B1 group, suggesting Vrn-B1 may play a role in regulating Vrn-B3 in wheat. Genes that potentially involved in floral development including CONSTANS-like genes, MADS-box genes and MYB transcription factors were identified in our analysis (Table S5). These genes were also identified as VRN1-binding targets that potentially regulate flowering through RNA-Seq and ChIP-seq by Deng et al. [43]. Furthermore, extremely up-regulated genes with a fold change > 100, including Beta-amylase (BMY1) and anther-specific protein (RTS), may also be regulated by Vrn-B1 and play a key role in the relationship between heading and flowering. Additionally, functional enrichment analyses also indicated that there are several significantly enriched pathways among the DEGs, such as carbon metabolism, starch and sucrose metabolism, glycolysis, gluconeogenesis, and fatty acid biosynthesis and metabolism (Table S7). Deng et al. identified direct targets of VRN1 including genes involved in the biosynthesis or breakdown of jasmonic acid, abscisic acid and gibberellin [43]. Similarly, these pathways, except for gibberellin, were enriched in our GO analysis (Table S6), indicating these pathways are conserved in regulating flowering of cereals. Carbohydrates play an important role in the control of flowering transition in plants. The accumulation of sucrose, a major sugar in both leaf and apical exudates, in the meristem precedes the activation of energy-consuming processes such as mitotic activation [44,45]. The Arabidopsis late flowering mutant carbohydrate accumulation1 (cam1) contains higher starch levels than WT at the onset of flowering, but the late flowering phenotype is not caused by increased starch levels [46]. Interestingly, this study suggested that flowing time had an effect on expression patterns of the genes directly or indirectly participating in carbohydrate metabolism. Fatty acid metabolism is important for pollen development. The cell walls of pollen contain fatty substances produced in the tapetum of anthers [47]. In addition, the early flowering phenotype caused by overexpression of fatty acid amide hydrolase in Arabidopsis also supports the observation that fatty acids play an important role in floral induction [48]. Consequently, we observed that DEGs involved in fatty acid biosynthesis and metabolism were enriched in lines with flowering time variants. It also indicated that fatty acid metabolic pathways interact with the flower development pathway in wheat and may be affected by variations in VRN-B1 gene expression. The flavonoid biosynthesis pathway has been well documented and is known to be associated with flowering time variation in plants [49][50][51]. Flavonoids are a class of powerful antioxidant compounds whose synthesis is initiated by the phenylpropanoid pathway and subsequently routes into nine different metabolic branches, including those that produce flavones and flavonols [52]. Interestingly, our transcriptomic analysis indicated that pathways related to phenylpropanoid biosynthesis, phenylalanine metabolism, flavonoid biosynthesis, flavone biosynthesis, and flavonol biosynthesis were significantly enriched in our data set (Table S7) and further detailed analysis of DEGs for phenylpropanoid biosynthesis suggested that 98 upregulated-genes and 16 down-regulated genes participating in most of steps in the pathway were identified (Fig.  S8), indirectly suggesting that these pathways are involved in the early heading process or influenced by flowering time in wheat.

Conclusions
In the present study, a major QTL for heading time was identified in a wheat RIL population using a combination of BSA and genetic mapping. Development of molecular markers for Vrn-B1 suggested that a large deletion in the first intron was partially responsible for heading time variation in the RIL population. Additionally, our results also revealed several important genes and pathways associated with flowering. The molecular markers for VRN-B1 developed in this study will be deployed for marker assisted breeding programs to optimize heading time, and the candidate DEGs associated with Vrn-B1 gene regulation help us to better understand the molecular function of heading-related pathways in wheat.

Plant materials and phenotypic analysis
The early heading mutant (eh1) was obtained by mutagenesis of dry seeds of spring wheat line Zhongyuan9 with 284 Gy γ-ray irradiation. The eh1 mutant was crossed with a late heading Chinese elite cultivar, Lun-xuan987 (LX987), to produce F 6:7 recombinant inbred lines (RILs) by single-seed descent. Zhongyuan9 wheat lines and the eh1 were created in our lab, and the seeds of LX987 were kindly provided by Binghua Liu (Institute of Crop Sciences, Chinese Academy of Agricultural Sciences). The RIL population (400 lines) and parent lines were sown at the Zhongpuchang station of the Institute of Crop Sciences, Chinese Academy of Agricultural Sciences (Beijing, China) and grown under well-managed field conditions. Each line was planted with 15 plants in a row of 1 m. The heading time of each RIL was recorded from 3 years' experiment of 2016 to 2018 when two thirds of the plants in a line were headed and more than half of the spikes had emerged.

BSR-Seq analysis
When the early heading RILs were beginning to heading, young spikes from extremely early and late heading RILs, and two parents were flash frozen in liquid nitrogen and then stored at − 80°C. In total, 21 early heading and 22 late heading lines were used. RNA was extracted from these lines using TRNzol-A + Reagent (Tiangen Biotech, China). RNA samples were treated with DNase I (Takara) and cleaned using an RNA purification kit (Tiangen Biotech). An equal amount of RNA from each line was mixed to construct early and late heading bulks. The RNA quality was assessed using both a NanoDrop™ One (Thermo Scientific) and a Bioanalyzer 2100 system (Agilent Technologies). A total of 4 cDNA libraries, including two parents and two bulks, were successfully constructed using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB). RNA sequencing was carried out on an Illumina HiSeq platform according to the manufacturer's protocol and 150 bp paired-end reads were generated (Biomarker Technologies Corporation, Beijing). After removing low-quality reads, reads with adapters or those with more than 10% unidentified nucleotides from raw data, a total of 54.71 Gb of clean reads were obtained. The detailed information of number of reads per library, Q30 and GC content were shown in Table S8. The clean reads were then aligned to the reference genome Chinese Spring wheat v1.0 released by the International Wheat Genome Sequencing Consortium (IWGSC) (http://www.wheatgenome.org/) using STAR software with default parameters [53]. Duplicates were excluded based on their position in the reference genome using Picard (https://sourceforge.net/ projects/picard/). Local realignments, base recalibration, and SNP calling was conducted using GATK software [54]. High-quality SNPs were obtained according to the following filtering criteria: (i) loci were removed if multiple alleles exist; (ii) SNP loci with reads support < 4 were filtered out; (iii) SNPs with the same genotype in bulks were excluded; (iv) SNPs from late heading bulks that differed from the late heading parent were discarded. The sequencing datasets generated in this study are deposited in the National Centre for Biotechnology Information (NCBI) under the BioProject ID PRJNA517367 with the Sequence Read Achieve (SRA) accession SRP182626.

Association analysis of ED and SNP-index
The Euclidean Distance (ED) algorithm evaluates the correlation between the target region and the trait of interest based on the significant difference of markers from sequencing data in bulks [55]. Theoretically, markers close to the target region are different in two extreme phenotype bulks while others out of the region are consistent, and their ED value is close to 0. The higher the ED value is, the greater the difference between two bulks. The formula for ED is as follows: q X mut , X WT represents the frequency of X base in mutant and wild type bulk, respectively. In the analysis, the sequencing depth and ED value were calculated for each different SNP in the bulks. To eliminate the background noise, the fifth power of the original ED value was regarded as a correlation value and then the ED value was fitted using the SNPNUM algorithm. The SNPindex was defined as the ratio between mutant allele frequency and wild type allele frequency in bulks [56,57]. The index is equal to 1 when all SNP reads differ from the wild type allele, and is equal to 0 when all SNP reads are identical to the wild type allele. To exclude the false positives, SNPs located on the same chromosome were fitted according to their position in the genome. The regions above the association threshold were considered as the candidate interval related to the early heading phenotype. The number of genes in the candidate region was calculated according to the number of predicted genes in the corresponding region of reference genome.

Genomic DNA extraction
Genomic DNA was extracted from fresh leaf samples at the seedling stage as previously described [58]. Briefly, after grinding the frozen leaves into homogeneous powder by using a Vibration Mill Type MM301 (Restsch GmbH, Germany), 600 μL DNA extraction buffer was added and completely mixed.

Construction of a genetic linkage map
A linkage map of the candidate chromosome was constructed using QTL IciMapping version 4.1 [59,60]. Recombination frequency was converted into centimorgan (cM) distances with the Kosambi map function [61]. QTL mapping was conducted using the days to heading in three individual years with the inclusive composite interval mapping (ICIM) analysis. A value of phenotypic variance explained (PVE) by an individual QTL was determined using ICIM. Significant QTLs were identified at a logarithm of odds (LOD) threshold of 3.0.

Identification of VRN-B1
Two pairs of primers were designed to identify the presence or absence of a large deletion in VRN-B1. The forward and reverse primer pair termed TavBI were designed to flank the deletion region so that product could only be detected in the dominant Vrn-B1. For the primer pair termed TavBII, the forward primer was designed to anneal to the deletion region and the reverse primer was designed to anneal to a region of the genome outside of the deletion so that product can only be detected in the recessive vrn-B1. PCR reactions were performed using T 5 Super PCR mix (Tsingke, Beijing) with the following procedure: 98°C for 3 min, followed by 32 cycles of annealing (98°C for 10 s, 58°C for 10 s; 72°C for 1 min.). The ratio of dominant effect that represent the difference of heading time between heterozygote and the midpoint of two parents, and additive effect that represent the half value of the difference of two parents, was considered to be the degree of dominance. For sequencing Vrn-B1, 11 pairs of primers were designed to amplify the whole gene and then sequenced. The primers for sequencing of VRN-B1 and Ppd-A1 and marker primers for detection of a 37 bp deletion in VRN-B1 and a 131 bp deletion in Ppd-A1 are shown in Table S9.

Quantification of gene expression levels and DEGs analysis
Based on the BSR-Seq data, the number of clean reads mapped to each gene was counted using the Cuffquant and Cuffnorm modules in Cuffinks software. The expected gene expression levels from two parents and two bulks were calculated using Per Kilobase of transcript per Million fragments mapped (FPKM). The FPKM algorithm normalizes transcript lengths and the number of mapped reads, and is a commonly used method for estimating gene expression levels [62]. Differentially expressed genes (DEGs) analysis was performed using the EBSeq R package, which provides a statistical method based on an empirical Bayes hierarchical model. The resulting P-value was adjusted using the Benjamini-Hochberg approach in order to control the false discovery rate. Genes were considered to be differentially expressed between two groups if the adjusted P-values were < 0.05 (FDR < 0.05) by EBSeq and the Fold Change ≥2.0.

GO term and KEGG enrichment analysis
Gene Ontology (GO; http://www.geneontology.org) enrichment analysis of DEGs was performed using the GOseq R package (corrected P-value < 0.05), in which gene length bias was corrected. GO terms with corrected P-values less than 0.05 were considered significantly enriched. Based on the retrieving of KEGG pathways for DEGs (http://www.genome.jp/kegg/), the statistical significance of the enrichment of DEGs in each KEGG pathway was analyzed by using KOBAS software. KEGG mapper was used for analysis of the detailed DEGs in the most significantly enriched pathway (https://www. kegg.jp/kegg/mapper.html).

Quantitative real-time PCR analysis
The same tissues as used in RNAseq analysis (eh1 and LX987 spike samples) were used for quantitative real-time PCR verification. For quantifying the relative expression of vernalization genes, young spikes and leaves of eh1 and LX987 were sampled when eh1 started to heading. For analysis of the gene expression at the same developmental stage, young spikes and leaves of eh1 were also sampled 8 days earlier than that of LX987 when they showed similar developmental levels of spikes. RNA extraction and qRT-PCR were conducted as previously reported [63]. Generally, total RNA from spike or leaf sample was isolated using an RNeasy Plant Mini Kit (Qiagen). DNase I (Takara) and RNA purification kit (Tiangen) were used for elimination any DNA contamination in the extracted RNA. By using a TransScript First-Strand cDNA Synthesis SuperMix (TransGen) kit, the first strand cDNA was synthesized. The qRT-PCR was performed using the SsoFast EvaGreen Supermix Kit (Bio-Rad, USA) and a CFX 96 Real-TimeSystem (Bio-Rad, USA). This experiment was conducted with two technical repeats and three independent biological replicates. For verification of DEGs, genespecific primers were designed using Oligo 7 software. The details of primers were listed on Table S10. For expression analysis of VRN1 and VRN2 genes, primers described in previous studies [7,64] were used. ACTIN was used as an internal control to normalize the expression data. Relative expression levels were determined using the 2 -ΔΔCT method [65].