The fine-mapping population was derived from a cross between two winter wheat lines: the maternal parent DM and the paternal parent Shi 4185. RIL population from this cross was previously established (Li et al., 2018; Xu et al., 2019). The paternal line Shi 4185, a cultivar released in Hebei Province characterized by smaller spikes, higher tillering capacity, and lower plant biomass; and the maternal line DM, which exhibits larger spikes, broader flag leaves, greater leaf area, and reduced tillering ability (Li et al., 2018; Xu et al., 2019). For this study, we utilized F₂ seeds from the same cross, obtained from the germplasm bank of the National Wheat Improvement Center in June 2022.

Field trials were conducted over three consecutive growing seasons (2022–2025) at two experimental locations: Boxin Agricultural Technology Co., Ltd. (Zibo, Shandong Province) and the Jiaozhou Demonstration Park (Qingdao, Shandong Province), China. In the 2022–2023 season, the population was planted in a double-row plot design (row length: 1.5 m) with two replicate blocks. Individual F₂ plants were screened to identify those heterozygous at the target QSL.caas-6BL.1 region. These RHLs were self-pollinated to generate F₃ seeds. In the subsequent 2023–2024 and 2024–2025 seasons, progeny rows from the selected RHLs (F₄ and F₅ generations) were planted in a three-row plot design with two replicate blocks at each location. SL was measured as the length of the main spike from the base of the rachis to the tip of the terminal spikelet, excluding awns. For each plant, three primary spikes were measured, and the average value was recorded as the representative SL. Plants were individually tagged, and phenotypic data were collected prior to tissue sampling for DNA extraction. Genomic DNA was extracted from leaf tissue using a high-salt and low-pH precipitation method. All field management followed local standard agronomic practices.

Marker development was based on the reference genome of Shi 4185 (Jiao et al., 2025). Whole-genome resequencing data of DM (unpublished data) were aligned to this reference and using Geneious Prime software (v2023.0.4), a 21 Mb target interval (700–721 Mb) on chromosome 6B (IWGSC RefSeq v2.1) was defined based on the initial mapping results (approximately 705 Mb). Polymorphic InDels and single-nucleotide polymorphisms (SNPs) were identified within this region. These variants were then used to develop InDel, CAPS, and KASP markers. Initial comparative genomic analysis between the parental lines enabled the development of seven polymorphic markers (3 InDel, 1 CAPS, 3 KASP) within the target interval. These markers were used to genotype an F₄ population of 492 individuals to identify plants harboring recombination events within the region. Two representative recombinant lines (L1–L2) were selected, as their recombination breakpoints collectively spanned the entire initial 21 Mb target interval. These plants were self-pollinated to generate five distinct F₄:₅ families. To enhance mapping resolution, a subsequent round of refined genomic comparison was performed, leading to the development of four additional high-quality polymorphic markers (1 InDel, 2 CAPS, 1 KASP) within critical regions indicated by preliminary analysis. This expanded set of 11 markers was deployed for high-density genotyping of approximately 230 individuals derived from the five F₄:₅ families (L3–L7). By analyzing the segregation patterns within each family, the recombination breakpoints in their respective F₄:₅ progenies were precisely determined and mapped to intervals between flanking markers, enabling a significant reduction of the target QTL to a narrower physical interval.

Based on phenotypic evaluation, the parental lines DM and Shi 4185 showed significant differences in young spike length at the stamen and pistil differentiation stage. To investigate the molecular basis, we selected a near-isogenic line (NIL-133-9) homozygous for the target QSL.caas-6BL.1 region from the cross population. Young spikes were collected from NIL plants homozygous for either the DM or Shi 4185 allele at the same developmental stage. For each genotype, spikes were collected from three different, independently grown plants, constituting three biological replicates. Samples were immediately frozen in liquid nitrogen and stored at -80 °C. Total RNA was extracted using TRIzol reagent. RNA integrity was verified with an Agilent 2100 Bioanalyzer. Libraries were constructed and sequenced on an Illumina platform with a PE150 strategy. Raw reads were quality-trimmed and aligned to the Chinese Spring reference genome (IWGSC RefSeq v2.1). Gene expression levels were quantified as TPM. Differential expression analysis was performed using thresholds of |log₂(fold change)| > 1 and adjusted P-value < 0.05. Total RNA for RNA sequencing was reverse transcribed using the Evo M-MLV Reverse Transcriptase Pre-mix Kit (AG). Gene-specific primers were designed via Primer-BLAST. qRT-PCR was performed on the QuantStudio 5 Real-Time Fluorescent Quantitative PCR System (Applied Biosystems) using ChamQ Blue Universal SYBR qPCR Master Mix (Vazyme). The wheat EF1α gene served as the internal control. Relative expression levels were calculated using the 2^(-ΔΔCt) method. Each sample underwent triplicate technical replicates.

Phenotypic data comparisons between genotypes were performed using Student’s t-test. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of differentially expressed genes (DEGs) were conducted using TBtools, with significance set at an adjusted P-value < 0.05. Homologous genes in rice (Oryza sativa Japonica Group) and maize (Zea mays) were identified using the Ensembl Plants database and cross-referenced with functional annotations from literature. All statistical analyses and visualizations were performed in Microsoft Excel, R (v4.3.2) or TBtools (v2.363; Chen et al., 2023).

To validate the major SL QTL QSL.caas-6BL.1, phenotypic characterization of the two parental lines (DM; Shi 4185) and the Residual heterozygous lines (RHLs) population (n=275) was conducted across two distinct environments (Zibo, Shandong Province in 2024; Qingdao, Shandong Province in 2025). The parental lines exhibited significant and consistent differences in SL (P < 0.05) across all tested environments. DM had an average SL of 7.09 ± 0.74 cm in 2024 and 6.60 ± 1.07 cm in 2025, while Shi 4185 had longer spikes with an average of 8.16 ± 0.96 cm in 2024 and 8.51 ± 0.82 cm in 2025 (Figure 1a; Supplementary Table 1). This consistent phenotypic difference confirms that Shi 4185 contributes the allele for increased spike length. In the RIL population, SL displayed continuous variation over a wide range: specifically, SL ranged from 3.80–10.50 cm in 2024 and 4.50–11.50 cm in 2025. The frequency distributions of SL in the population approximated normal distributions across both environments (Figure 1b; Supplementary Table 1). These results confirm that wheat SL is a quantitative trait controlled by multiple genes and validate the suitability of this population for subsequent QTL mapping of SL.

To finely map the major SL-regulating QTL QSL.caas-6BL.1, we identified kompetitive allelespecific PCR (KASP) markers on chromosome 6BL between DM and Shi 4185. Seven molecular markers (INDEL691, CAPS699, INDEL702, INDEL704, KASP704, KASP705, KASP711) were mapped into the genetic region of QSL.caas-6BL.1 (Supplementary Figure 1; Supplementary Table 2). Using these markers, we screened recombinants from the F₄ population (492 individuals) of the cross DM/Shi 4185. This screening identified 180 recombinant individuals carrying crossovers within the target interval (700–721 Mb) (Supplementary Figure 1, Supplementary Table 2). Within these, two key recombinant individuals (L1–L2) were selected for further analysis—their recombination breakpoints collectively spanned the entire 700–721 Mb region, and their progeny (F₅ lines) derived via self-pollinated were used for subsequent high-resolution mapping.

In the subsequent growing season, five F₅ families (L3–L7) were genotyped to precisely determine the recombination breakpoints of their corresponding F₄ parent individuals. To refine genotyping resolution, four additional markers (KASP702, CAPS703, INDEL708, CAPS710) were developed for critical regions identified in the initial mapping, establishing a high-density genotyping panel of 11 markers (Supplementary Figure 1; Supplementary Table 2). By analyzing segregation patterns within each family, the high-density genotyping data enabled precise mapping of recombination breakpoints as intervals between adjacent markers in these nine critical F₄:₅ recombinant lines (L1–L7). Figure 2 illustrates the distribution of these recombination breakpoints, which are positioned between markers and partition the target interval into multiple sub-intervals. Through phenotypic evaluation and genotype-phenotype co-segregation analysis of the F₄:₅ families, the causal genetic interval responsible for spike length variation was ultimately delineated to a 3.3 Mb physical region (708.2–711.5 Mb) (Figure 2).

The QSL.caas-6BL.1 locus was initially identified within a genomic region on chromosome 6BL known to harbor a cluster of QTLs for multiple yield-related traits. In the primary QTL analysis of the same DM × Shi 4185 population, Li et al. (2018) reported that the interval containing QSL.caas-6BL.1 co-localized with stable QTLs for thousand-kernel weight (QTKW.caas-6BL), spike dry weight (QSDW.caas-6BL), heading date (QHD.caas-6BL), and plant height (QPH.caas-6BL). This preliminary evidence suggested that the genetic region responsible for spike length might also influence other key components of grain yield, providing the impetus for the present fine-mapping study.

Following the fine-mapping of the major SL-regulating QTL QSL.caas-6BL.1 to a 3.3 Mb physical interval (708.2–711.5 Mb), we performed an integrated analysis of this genomic region. Based on the IWGSC RefSeq v2.1 genome annotation, 25 high-confidence genes within this interval were identified (Table 1). To explore genetic variations underlying SL differences between parents, we conducted comparative genomic analysis using genome assemblies of DM (unpublished data) and Shi 4185 (Jiao et al., 2025). This analysis revealed sequence variations in the open reading frames (ORFs) of 18 genes between the two parents, consisting of missense mutations derived from SNPs in 10 genes and frameshift or disruptive mutations caused by insertions/deletions (InDels) in 14 genes (detailed in Supplementary Table 3). To prioritize the most promising candidate genes for SL regulation, we conducted a systematic screening of sequence variations between DM and Shi 4185. This screening identified three genes—TraesCS6B03G1207300, TraesCS6B03G1211400, and TraesCS6B03G1211800—as strong candidates potentially directly involved in SL development based on the presence of protein-altering mutations. The remaining 22 genes were excluded from further consideration, as they either lacked functionally meaningful sequence variations or had no documented association with known spike development regulatory pathways.

To further precisely identify candidate genes regulating SL within the 3.3 Mb QSL.caas-6BL.1 interval, we performed RNA-seq analysis on young spikes of DM and Shi 4185 near-isogenic lines at the stamen and pistil differentiation stage. Three biological replicates were included for each genotype to ensure data reliability. Principal component analysis (PCA) of the transcriptomic data showed that replicates from the same parent clustered tightly, while clear separation was observed between the two parental genotypes (Supplementary Figure 2). This result indicated that genotype was the primary driver of transcriptomic variation. Additionally, a kinship matrix heatmap constructed from genome-wide expression data confirmed high intra-group reproducibility, further validating the rigor of the experimental design and the quality of the RNA-seq data (Figure 3a).

Subsequent differential expression analysis identified 5,046 differentially expressed genes (DEGs) (|log₂FoldChange| > 1 and P-adjust < 0.05) when mapped to the Chinese Spring reference genome (IWGSC RefSeq v2.1) (Supplementary Table 4). We then focused specifically on the 25 high-confidence genes within the 3.3 Mb fine-mapped QSL.caas-6BL.1 interval. Among these, 23 genes showed no significant differential expression between DM and Shi 4185 NILs (Figure 3b; Supplementary Tables 3–5). Of the remaining two DEGs (TraesCS6B03G1210400 and TraesCS6B03G1211800), although TraesCS6B03G1210400 shows differential expression, comparative genomic analysis (as shown in Supplementary Table 3) revealed no non-synonymous SNPs or insertion-deletion (InDel) mutations within its coding sequence. Therefore, TraesCS6B03G1211800 was the gene that exhibited both significant differential expression between parents and coding sequence variations (Figure 3b; Supplementary Tables 3–5), highlighting it as a high-priority candidate gene.

To characterize their functional roles, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on 5,046 DEGs (Figures 3c, d; Supplementary Tables 6, 7). For GO enrichment, the top 10 enriched terms in each of the three core categories—cellular component, biological process, and molecular function—were visualized. Among these, the molecular function category, “catalytic activity (GO:0003824)” was the most prominently represented term, reflecting the involvement of metabolic enzyme activities in spike development. In the cellular component category, terms related to photosynthetic structures—including “thylakoid part (GO:0044436)”, “thylakoid (GO:0009579)”, “photosynthetic membrane (GO:0034357)”, and “photosystem (GO:0009521)”—were significantly enriched, consistent with the critical role of spike photosynthesis in grain development. In the biological processes category, enriched terms were mainly associated with regulatory processes, such as “regulation of DNA-templated transcription (GO:0006355)” and “regulation of cellular biosynthetic process (GO:0031323)”, indicating that transcriptional regulation is a key mechanism underlying SL variation (Figure 3c; Supplementary Table 6). KEGG pathway analysis further revealed the top 20 significantly enriched pathways. These DEGs were mainly associated with metabolic pathways, including carbohydrate metabolism, biosynthesis of secondary metabolites, and amino acid metabolism (Figure 3d; Supplementary Table 7). These results align with the dynamic physiological and biochemical changes occurring during young spike development, providing a functional framework for understanding the molecular mechanisms of SL regulation mediated by QSL.caas-6BL.1.

To independently validate the reliability of our transcriptome data, we performed quantitative real-time PCR (qRT-PCR) analysis on six selected genes using the same RNA samples employed for RNA-seq. The gene set included two DEGs from the fine-mapped QSL.caas-6BL.1 interval (TraesCS6B03G1211800 as the prime candidate and TraesCS6B03G1210400) and four additional DEGs randomly picked from the genome-wide RNA-seq dataset for expanded validation. The qRT-PCR results showed a highly significant positive correlation with the RNA-seq data (Supplementary Figure 3; Supplementary Table 2).

To robustly prioritize the causal gene from the 25 candidates and mitigate the risk of false negatives from expression screening alone, we employed a multi-layered evidence integration strategy. Homology and functional domain analysis indicated several genes with plausible roles in development (Supplementary Figure 4). Notably, TraesCS6B03G1211800 encodes a NAC domain transcription factor, a family key to meristem regulation. Its rice ortholog (Os07g0225300) belongs to a clade containing known panicle architecture regulators (Anjum and Maiti, 2024). Furthermore, protein sequence and domain comparison with a curated set of cloned spike development genes in wheat, rice, and maize (Supplementary Table 8) revealed that the NAC domain of TraesCS6B03G1211800 shares high similarity with those of established regulators such as OsNAC121 (Os10g0571600) (Anjum and Maiti, 2024).

Comparative analysis with previously reported spike-related QTLs/genes confirmed the functional relevance of this genomic region on chromosome 6BL (Simons et al., 2006; Dixon et al., 2018; Du et al., 2020; Kuzay et al., 2022; Liu et al., 2022; Anjum and Maiti, 2024; Lin et al., 2024; Agata, 2025). This interval overlaps with a yield-QTL cluster on 6BL, which contains TraesCS6B03G1211100 (KAT-2B)—a validated grain weight regulator—further supporting the functional importance of this chromosomal segment. However, genetic co-segregation analysis in our key recombinant families showed that the spike length phenotype perfectly tracked the haplotype of TraesCS6B03G1211800, but not that of TraesCS6B03G1211100, genetically separating the SL QTL from the linked grain weight effect. Conclusive integration of molecular evidence pinpointed the final candidate. TraesCS6B03G1211800 was the only gene within the fine-mapped interval that satisfied all three decisive criteria: 1) significant differential expression in young spikes, 2) presence of a protein-altering missense mutation, and 3) genotype-phenotype co-segregation. Therefore, through systematic integration of cross-species homology, genetic dissection, and transcriptomic-genomic variation, we conclusively identify TraesCS6B03G1211800 as the prime candidate gene underlying the QSL.caas-6BL.1 QTL for spike length.