Generation of RNA-seq data
In this study, we constructed 24 RNA-seq libraries of the 8 treatments, WHd-A2.0, MHd-A2.0, WHd-L2.0, MHd-L2.0, WHd-A3.5, MHd-A3.5, WHd-L3.5 and MHd-L3.5, and each treatment included 3 biological repeats. After high-throughput sequencing, each sample contained 40 – 59.4 million reads, the data size ranged from 6 – 8.9 Gb, the Q30 value exceeded 87%, the GC content distribution was 50 – 56%, and sequencing alignment analysis results showed that 89.35 – 97.79% of the reads could be mapped to IWGSC RefSeq v1.1 (Table S1). Principal component analysis (PCA) of the eight raw sequencing datasets indicated that the samples could be clustered into four groups according to their genotype, which showed good repeatability between samples that could be used for subsequent analysis (Fig. S1).
Identification of differentially expressed genes
To identify genes that differed significantly between the early and late heading genotypes during the development of young panicles, we analyzed the differentially expressed genes (DEGs) with strict quality control. The DEGs of L2.0, L3.5, A2.0, and A3.5 were compared between the genotypes with an early heading time (WHd) and a late heading time (MHd). Finally, we identified 18,325 unique DEGs (Fig. 1a, Table S2), which were mainly distributed at the end of the chromosome arm (Fig. S2). Specifically, there were 12,941 DEGs in the apical meristem at A2.0, 9,990 of which were upregulated and 2,951 of which were downregulated (Fig. S3a); 1,557 DEGs at A3.5, 1,406 of which were upregulated and 151 of which were downregulated (Fig. S3b), 6,390 DEGs at L2.0, 2,735 of which were upregulated and 3,655 of which were downregulated (Fig. S3c), 1,040 DEGs at L3.5, 764 of which were upregulated and 276 of which were downregulated (Fig. S3d). Flowering time genes Vrn1-5A , Vrn3-7B , Ppd-1D [9, 10], and WSOC1  were also found to express differentially. We also screened the common differentially expressed genes at different developmental stages in the same tissue and found that 1,084 genes were common DEGs at A2.0 and A3.5 (Fig. 1b), and 215 genes were common DEGs at L2.0 and L3.5 (Fig. 1c). We speculated that these genes may play important roles in spikelet development at the W2.0 and W3.5 stages. Moreover, nine genes were selected from the above DEGs to assess the expression using qRT-PCR (Fig. 2a and 2b). The results showed that the expression patterns were consistent with the transcriptome analysis (Fig. 2c and 2d).
GO enrichment analysis of differentially expressed genes
To explore the regulatory pathways of DEGs, GO enrichment analyses were conducted between the early heading date and the late heading date wheat accessions. We finally identified 468 unique GO terms (Table S3, Fig. 3) that were involved in biological process (BP), molecular function (MF), and cellular component (CC) (Fig. 3a; Fig. S4; Fig. S5; Fig. S6). There were 222, 92, 208, and 127 significant GO terms in A2.0, A3.5, L2.0, and L3.5, respectively (Table S3). Interestingly, we discovered several flowering time-related GO terms, including GO:0005975 (carbohydrate metabolism) , GO:0005991 (trehalose metabolic process) , and GO:0019684 (photosynthesis, light reaction) . Furthermore, we found 5 common GO terms, e.g., GO:0003700 (transcription factor activity, sequence-specific DNA binding), GO:0001071 (nucleic acid binding transcription factor activity), GO:0003824 (catalytic activity), GO:0043565 (sequence-specific DNA binding), and GO:0016491 (oxidoreductase activity), which were generally expressed in all comparative groups of DEGs (Fig. 3b).
MapMan metabolic pathway analysis of differentially expressed genes
The metabolic pathways of four comparative groups of DEGs were visualized by MapMan software (Fig. 4a; Fig. S7; Fig. S8; Fig. S9), and a total of 932 unique metabolic pathways were ultimately identified (Table S4, Fig. 4b). At A2.0, the DEGs mainly functioned in DNA.synthesis/chromatin structure (bin:28.1), signalling.receptor kinases (bin:30.2), signalling.light (bin:30.11), and hormone metabolism.jasmonate (bin:17.7). At A3.5, the DEGs function in the RNA.regulation of transcription.NAC domain transcription factor family (bin:27.3.27), PS.lightreaction.photosystem II (bin:1.1.1), hormone metabolism.ethylene.signal transduction (bin:17.5.2), protein.synthesis.initiation (bin:29.2.3). At L2.0, the metabolic pathway was mainly involved in secondary metabolism.phenylpropanoids.lignin biosynthesis (bin:16.2.1), and signalling (bin:30), RNA.regulation of transcription (bin:27.3.99). At L3.5, the DEGs were involved in DNA.synthesis/chromatin structure.histone (bin:28.1.3), transport (bin:34), metal handling.binding, chelation and storage (bin:15.2), RNA.regulation of transcription.TCP transcription factor family (bin:27.3.29). Moreover, 214 common metabolic pathways were screened, including bin:1.1.1 (light reaction), bin:11.8 (lipid metabolism), and bin:22.214.171.124 (amino acid metabolism) (Fig. 4b).
Transcription factor classification and identification
To identify the transcription factors (TFs) involved in the heading time development process, the DEGs were subjected to transcription factor analysis using iTAK software. In our study, according to the criteria of the plant transcription factor database (http://planttfdb.gao-lab.org/), a total of 1,225 transcription factors were classified into 45 transcription factor families (Table S5). A large number of DEGs were bHLH (129, 10.51%), WRKY (109, 8.88%), NAC (91, 7.41%), AP2/ERF-ERF (88, 7.17%), and MYB (64, 5.21%) transcription factors (Table S5). Among them, we also identified several transcription factor families involved in the plant flowering process, including three LFY transcription factors , eighteen SBP transcription factors, thirty-six MADS-MIKC transcription factors and ten MADS-M-type transcription factors (Table S5).
Gene expression level analysis showed that transcription factors were expressed at the critical period that determines flowering time. For example, we found that the three LFY transcription factor genes (TraesCS2A02G443100, TraesCS2B02G464200, TraesCS2D02G442200) were highly expressed in the apical meristem, especially in WHd-A3.5 and WHd-A2.0 (Fig. 5a). Most of the SBP transcription factor genes were highly expressed in WHd-A2.0 and WHd-A3.5 (Fig. 5b). In the MADS-box gene family, we discovered that the TraesCS5A02G391700 (Vrn1-5A) gene always showed a high expression level in both wild and mutant wheat, TraesCS3B02G612600 was highly expressed in leaf tissue, and two genes (TraesCS3D02G284200, TraesCS3A02G284400) were highly expressed in the apical meristem (Fig. 5c).
Construction of the flowering gene regulatory network
In our study, a total of 18,352 DEGs were used to construct the co-expression network. We calculated the average gene connectivity under different soft-thresholding powers and found that when β = 14, the co-expression network had scale-free characteristics (Fig. S10a; Fig. S10b; Fig. S10c; Fig. S10d). Finally, a total of 17 co-expression modules were obtained by the “dynamic tree cut” method. Each branch in the cluster tree represents a gene set, and different modules were distinguished by different colors (Fig. S10e). The correlation value between co-expression modules and samples was also calculated (Fig. S11). The reported flowering time genes Vrn1-5A (TraesCS5A02G391700), Vrn3-7B (TraesCS7B02G013100), Ppd-1D (TraesCS2D02G079600), and WSOC1 (TraesCS4D02G341700) were selected to build a gene co-expression network, and genes connected to the reported flowering time genes were extracted from the co-expression modules. We finally detected 16, 336, 446 and 124 DEGs with biological connections to Vrn1-5A (Fig. 6a), Vrn3-7B (Fig. S12a), Ppd-1D (Fig. S13a), and WSOC1 (Fig. S14a), respectively. The completed gene list related to the reported flowering time genes is summarized in Table S6.
Analysis of gene expression patterns co-expressed with Vrn1-5A, Vrn3-7B, Ppd-1D, and WSOC1
To further study the gene expression pattern of the differentially expressed genes co-expressed with the wheat heading date genes (Vrn1-5A, Vrn3-7B, Ppd-1D, and WSOC1), we found that the 16 genes that were co-expressed with Vrn1-5A could be divided into three patterns according to their gene expression level. First, TraesCS2D02G181400 is a MADS-MIKC transcription factor that was highly expressed in WHd-A3.5, MHd-A3.5, MHd-L3.5, and WHd-A2.0 but expressed at low levels in WHd-L2.0, MHd-L2.0, MHd-A2.0, and WHd-L3.5 (Fig. 6b). TraesCS2D02G181400 was upregulated between WHd-A2.0 and MHd-A2.0 and downregulated between WHd-L3.5 and MHd-L3.5 (Fig. 6c). Additionally, the other five genes (TraesCS1D02G339300, TraesCS2A02G562100, TraesCS3A02G007200, TraesCS3B02G361000 and TraesCS4D02G337000) were expressed at low levels at all stages in both the early and late heading date accessions (Fig. 6b). In addition, the remaining ten differentially expressed genes (TraesCS3A02G231100, TraesCS3B02G152300, TraesCS3B02G213900, TraesCS3B02G260400, TraesCS5D02G386500, TraesCS6A02G146500, TraesCS7A02G439300, TraesCS7A02G489000, TraesCS7D02G429000 and TraesCS7D02G475300) showed higher expression levels in WHd-L2.0, WHd-L3.5, and MHd-L2.0 but presented lower expression in WHd-A2.0, WHd-A3.5, MHd-A2.0, MHd-A3.5, and MHd-L3.5 (Fig. 6b). The genes highly related to Vrn3-7B were mainly expressed in the leaf tissues (Fig. S12b). Genes associated with Ppd-1D tended to be expressed in the leaf (Fig. S13b). Some genes that were co-expressed with WSOC1 were preferentially expressed in leaf tissues (Fig. S14b). Furthermore, by gene function annotation, we found some vital candidate flowering time genes in the co-expression network, such as TraesCS1A02G220300, TraesCS2D02G181400, TraesCS3A02G143100 (Fig. 2), and TraesCS5A02G166100 (Table 1).