Identification of a major and stable QTL on chromosome 5A confers spike length in wheat (Triticum aestivum L.)

Spike length (SL) is the key determinant of plant architecture and yield potential. In this study, 193 recombinant inbred lines (RILs) derived from a cross between 13F10 and Chuanmai 42 (CM42) were evaluated for spike length in six environments. Sixty RILs consisting of 30 high and 30 low SLs were genotyped using the bulked segregant analysis exome sequencing (BSE-Seq) analysis for preliminary quantitative trait locus (QTL) mapping. A 6.69 Mb (518.43–525.12 Mb) region on chromosome 5AL was found to have a significant effect on the SL trait. Fifteen competitive allele-specific PCR (KASP) markers were successfully converted from the single nucleotide polymorphisms (SNPs) in the SL target region. Combined with four novel simple sequence repeat (SSR) markers, a genetic linkage map spanning 21.159 cM was constructed. The mapping result confirmed the identity of a major and stable QTL named QSl.cib-5A in the targeted region that explained 7.88–26.60% of the phenotypic variation in SL. QSl.cib-5A was narrowed to a region of 4.84 cM interval corresponding to a 4.67 Mb (516.60–521.27 Mb) physical region in the Chinese Spring RefSeq v2.0 containing 17 high-confidence genes with 25 transcripts. In addition, this QTL exhibited pleiotropic effects on spikelet density (SD), with the phenotypic variances proportion ranging from 11.34 to 19.92%. This study provides a foundational step for cloning the QSl.cib-5A, which is involved in the regulation of spike morphology in common wheat.


Introduction
Common wheat (Triticum aestivum L.) is one of the most important food crops widely cultivated, and estimates suggest that the wheat production needs to increase 70% to fulfill future food demand. Therefore, the major objective of recent wheat breeding programs has been to improve grain yield (Cao et al. 2019;Su et al. 2016). Grain yield is determined by the number of spikes, grain number of per spike (GNS), and thousand grain weight (TGW). Among these factors, GNS and TGW are closely related Abstract Spike length (SL) is the key determinant of plant architecture and yield potential. In this study, 193 recombinant inbred lines (RILs) derived from a cross between 13F10 and Chuanmai 42 (CM42) were evaluated for spike length in six environments. Sixty RILs consisting of 30 high and 30 low SLs were genotyped using the bulked segregant analysis exome sequencing (BSE-Seq) analysis for preliminary quantitative trait locus (QTL) mapping. A 6.69  Mb) region on chromosome 5AL was found to have a significant effect on the SL trait. Fifteen competitive allele-specific PCR (KASP) markers were successfully converted from the single nucleotide polymorphisms (SNPs) in the SL target region. Combined with four novel simple sequence repeat (SSR) markers, a genetic linkage map spanning to spike-related traits, such as spike length (SL), spikelet density (SD), grain length (GL), and grain width (GW) (Cao et al. 2019). SL is usually positive for spikelet number per spike (SNS), resulting in more GNS. A lower SD usually leads to reduced SNS, decreasing yield production, while a higher SD increases the incidence of Fusarium head blight and pre-harvest sprouting (Huang et al. 2018). Therefore, understanding the genetic basis of desirable traits is important for breeding cultivars that can adapt to diverse environments and cropping systems (Yu et al. 2014a).
Previous results indicated that SL and SD are usually controlled by multiple genes (Jin et al. 2020;Cui et al. 2012;Yu et al. 2014b). Q, C and S are three major genes known to affect the spike morphological traits of wheat. The Q gene, on chromosome 5A, is widely studied with respect to cultivation because it confers the grain free-threshing ability. In addition, Q exerts a pleiotropic influence on many other traits, such as plant height, spike fragility, spike emergence time, spike length and spikelet compactness (SC) (Xie et al. 2018a). Due to the effect of the dominant allele at the compact (C) locus, which is positioned in the interval of Xgwm484-Xgwm539 on chromosome 2DL, the spikes of club wheat are more compact than those of common wheat and have pleiotropic effects on grain size, shape, and number (Johnson et al. 2008). S, on chromosome 3D, can determine whether a spike has round seeds and glumes (Faris et al. 2014). Moreover, spike development is affected by genes associated with heading date. Three categories of genes, including vernalization response genes (Vrn), photoperiod-response genes (Ppd), and earliness per-se genes (Eps), control heading date and are associated with spike development (Cao et al. 2019). Vrn genes have a profound influence on flowering time. The vernalization response is mainly controlled by the homologous genes of Vrn-A1 (Vrn-1), Vrn-B1 (Vrn-2) and Vrn-D1 (Vrn-3), which are on the 5A, 5B, and 5D chromosomes of wheat, respectively. However, Vrn-1 has the greatest effect on wheat vernalization. VRN1 encodes a MADS box transcription factor that binds the promoter of FLOWERING LOCUS T-like 1 to trigger flowering of cereal crops. RNA-seq revealed that VRN1 also targets genes that play central roles in spike architecture and hormone metabolism. For example, the elevated expression of VRN1 in cereals might reduce height and spike length (Deng et al. 2015). Wheat Eps are completely linked to spike development. The Eps-A1 gene affects the heading time, spike development and spikelet number of diploid wheat (Faricelli et al. 2010). By affecting the production time of terminal spikelets, Ppd genes play an important role in the photoperiod response and determine the heading date, thus influencing the spikelet number .
Furthermore, QTLs for SL and SD were identified on almost all of the wheat chromosomes of different genetic populations and natural varieties examined (Jin et al. 2020;Cui et al. 2012;Heidari et al. 2011;Luo et al. 2016;Yu et al. 2014b;Liu et al. 2018a). For example, nine QTLs for spike length on chromosomes 2D, 3A, 4A, 5B, 5D, and 6B were identified by using F 2 and F 2:3 populations (Cao et al. 2019). Liu et al. (2020) reported several stably expressed QTLs for spikelet density of common wheat in multiple environments. A major SL QTL associated with flanking markers Xcfd53 or Xgwm261 was validated in different wheat varieties and populations (Deng et al. 2017;Zhai et al. 2016;Wu et al. 2014). Reports have indicated that Rht8 has not only a significant effect on plant height but also has a regulatory effect on SL, SD, SNS, and TGW (Cui et al. 2012). Heidari et al. (2011) reported that Rht-B1 colocalized with QTLs for spike compactness and putative effects on other spike characteristics. Moreover, Zhai et al. (2016) identified two other QTL clusters on chromosomes 5A and 7B, which are both associated with QTLs for spike length and spike compactness.
Bulked segregant analysis (BSA) with pooled DNA samples obtained from individuals with extreme phenotypes in the RIL population provides a highefficiency approach to the identification of polymorphic markers associated with target traits (Zhai et al. 2016;Zou et al. 2016;Win et al. 2017;Ramirez-Gonzalez et al. 2015). Due to the large size of the wheat genome (~ 17 Gb) and highly repetitive content (IWGSC 2014;Gupta et al. 2008), whole-genome DNA resequencing for BSA is still costly (Martinez et al. 2020). Exons constitute approximately 0.95% of the wheat genome and code primary protein components that play essential roles in plant growth and development (Dong et al. 2020). Therefore, exome capture can effectively capture the protein-coding regions of the wheat genome to discover SNPs (Allen et al. 2013;Winfield et al. 2012) and catalog induced mutations (Henry et al. 2014;King et al. 2015). A recent report confirmed that exome capture resulted in a median mutation site coverage of 21 times and a 40-to 60-fold enrichment (Uauy et al. 2017). A pilot study assessing more than 10 million mutations in the protein-coding regions of 2735 mutant lines of tetraploid and hexaploid wheat revealed deleterious alleles for ∼ 90% of the captured wheat genes and demonstrated that exome capture followed by sequencing was a viable strategy to identify mutations in wheat (Krasileva et al. 2017;King et al. 2015;He et al. 2019). Therefore, BSE-Seq with wheat is likely a superior cost-efficient mapping method for identifying quantitative traits that exhibit clearly different phenotypes.
Synthetic hexaploid wheat (SHW), a carrier of various elite genes, is an important genetic resource in the improvement of common wheat (Yu et al. 2014a). Chuanmai 42 (CM42), an excellent cultivar, was developed by crossing synthetic hexaploid wheat Syn769 with Sichuan commercial wheat cultivars. CM42 has excellent characteristics, such as high yield, wide adaptability, and rust resistance (Yang et al. 2009). Hence, it has been used as a core parent in breeding programs. To identify the genetic loci controlling spike length in CM42, we performed QTL mapping in an RIL population derived from CM42 (longer SL) and 13F10. The objective was to gain new insights into the genetic basis of SL under different environments in wheat.

Materials and methods
Plant materials and field trials CM42 is a wheat cultivar with a high spike length and low spike density. 13F10 is a high-yield line with high SD selected by our laboratory (Supplemental Fig. 1). A population of 193 RIL lines, generated by the single seed descent method from a cross of 13F10 × CM42 (13CM), was used in this study. The parents and RILs of the 13CM population were grown in six different environments for phenotype evaluation: Shuangliu in (2017SHL, 2018SHL, and 2019SHL), Shifang in 2016. A randomized complete block design was used for RIL populations in all environments. Each plot comprised two rows with a 1.0-m row length, 0.2-m row spacing, and 10 seeds per row. Field management was performed according to customary practices for wheat production.
Phenotype evaluation of spike-related traits and statistical analyses The RILs along with their parents were measured for six traits including: SL, SD, SNS, GL, GW, and TGW. At maturity, six representative plants from each genotype were used for phenotypic evaluation. All data were collected from the main tillers. SL was the length of the main spike of an individual plant and measured from the base of the rachis to the tip of the terminal spikelet (excluding awns). SNS was determined by counting the number of fertile spikelets of the main spike. SD was calculated by dividing the SL by the SNS. GL and GW were determined from the weight of more than 200 random kernels with two technical repeats using SC-G software (Hangzhou Wanshen Detection Technology Co., Ltd., Hangzhou, China). TGW was calculated as 10-folds of the weight of 100 seeds with electronic balance, and three replicates for each line were set.
The phenotypic variation, normal distribution, and Pearson's correlation analyses of the wheat lines grown in different environments were performed using SPSS version 20.0 (IBM SPSS) and Microsoft Excel 2016. The best linear unbiased estimates (BLUEs) for target traits were calculated using QTL IciMapping V4.1 (http:// www. isbre eding. net/ softw are/). The estimated heritability (H 2 ) of traits was calculated according to the method described by Smith et al. (1998). Student's t test (P < 0.05) was performed by SPSS to detect significant differences in the parent and the RIL populations. Correlations between SL and SNS, GNA, SD, GL, GW, and TGW in the 13CM population based on the BLUE datasets were also analyzed by calculating Pearson's correlation coefficients.

BSE-Seq analysis
High-quality genomic DNA from the parents and the 13CM RIL individuals was extracted from tender leaves using the cetyltrimethylammonium bromide (CTAB) method followed by RNase-A digestion. The isolated DNA was quantified using an agarose gel. Based on the phenotypic investigation of RIL population, two DNA pools were constructed by mixing equal amount of DNA from 30 long-spike individuals (named SL-H) and 30 short-spike individuals (named SL-D), respectively. The selected proportion comprised 11% at each tail for extreme phenotypes, which was presumed to provide a 95% probability of detecting QTLs with large effects (Sun et al. 2010). Exome capture sequencing and analysis of these four libraries (13F10, CM42, SL-H, and SL-D) were performed by Bioacme Biotechnology Co., Ltd. (Wuhan, China, http:// www. whbio acme. com).
Based on the sequencing data, two extensive analysis methods were employed in this study. Euclidean distance (ED), with the advantage of linearity, is a better calculation method to measure the separation of alleles (Hill et al. 2013). Another widely used approach for analyzing BSA-Seq data is the SNPindex method (Li et al. 2020). When performing an ED analysis, we calculated the ED value of each site, took the 4th power of the original ED as the correlation value to eliminate background noise, and then used LOESS (locally weighted scatterplot smoothing) to fit the ED value. The SNP-index was used to screen the sites with significant differences between the mixed pools (Abe et al. 2012), and then, the absolute value of Δ(SNP-index) was used for local polynomial regression. ANNOVAR software was used for deep annotation of multiple databases (GO, KEGG) on the coding genes in the candidate interval ).

Development of molecular markers
According to BSE-Seq analysis, the variants sites (including SNPs and InDels) between 13F10 and CM42 were screened. Furthermore, only polymorphic SNPs between the two bulked pools (SL-H and SL-D) associated with the target QTLs in the candidate region (Chinese Spring reference genome IWGSC RefSeq v2.0) were converted into KASP markers for QTL confirmation and linkage analyses. SNP clustering and genotype calling of these KASP markers were performed by China Golden Marker Biotechnology Co., Ltd. (Beijing, China, http:// www. cgmb. com. cn). Finally, 15 pairs of KASP markers were successfully developed (Supplemental Table 5).
In addition, we developed SSR markers with parental polymorphisms by referring to the Chinese Spring reference genome IWGSC RefSeq v2.0 (http:// www. wheat genome. org). SSR primers were designed using the online primer design pipeline PolyMarker (http:// polym arker. tgac. ac. uk/), and they are listed in Supplemental Table 6. Polymerase chain reaction (PCR) amplification was performed in a 20-μL reaction volume containing 0.5 μL of template DNA (200 ng/μL), 1 μL of primer mixture (10 μM), 8.5 μL of RNasefree water, and 10 μL of 2 × Taq Master Mix. The PCR conditions were as follows: a denaturation step (94 °C for 4 min), followed by 28 cycles of 94 °C for 30 s, 51 °C-65 °C (primer-dependent) for 30 s and 72 °C for 30 s, with a final incubation at 72 °C for 5 min. The PCR products were separated by 8% polyacrylamide gel electrophoresis, and then, the gel was subjected to silver staining and the results statistically analyzed. In addition, based on the SSR marker alleles in the parental lines and RILs, the genotypes were categorized into two groups: genotypes with homozygous alleles from 13F10 (designated as "AA") and genotypes with homozygous alleles from CM42 (designated as "aa"). The average value of SL of the RILs was used to measure the QTL effect.
Furthermore, we used a KASP marker to identify the Vrn gene. The KASP marker (exon 7_C_T_Vrn-A1: A013151) was provided by China Golden Marker Biotechnology Co., Ltd. (Beijing, China, http:// www. cgmb. com. cn). The T3142C substitution mutation, located in the miR172 target site, differentiates the q allele from the Q allele. Genome-specific primers PJG14 and PJG18 have been successfully used to amplify a fragment spanning T3142C (Wolde et al. 2019).

Genetic linkage map construction and QTL mapping
Markers that were co-located with other markers and had > 20% missing values were discarded. The linkage map of the target QTL region in the 13CM RIL population was constructed by JoinMap 4.0, with LOD scores ranging from 5 to 10. The observed phenotypic values from the 13CM population in the six environments as well as BLUE data were used for QTL mapping analyses. IciMapping 4.1 was used to detect QTLs based on the Biparental Populations (BIP) module with inclusive composite interval mapping (ICIM). The walking speed for the QTLs was 1.0 cM, and the LOD score threshold was set at 2.5. Joint analysis for SL in multiple environments was performed with the MET module of QTL IciMapping 4.1. MapChart 2.2 was used to demonstrate the location of the QTLs on chromosome 5A. QTLs were named according to the rules of International Rules of Genetic Nomenclature.
Identification of candidate genes in QSl.cib-5A The physical positions and sequences of flanking markers of the QTL were obtained by blasting against the Chinese Spring genome (IWGSC RefSeq v 2.0). Genes predicted with high confidence in this candidate region were extracted from the JBrowse website (https:// urgi. versa illes. inra. fr/ jbrow seiwg sc). Candidate gene analysis for QSl.cib-5A was performed as indicated above using the closely linked markers on chromosome 5A. The gene expression data in different tissues were obtained from expVIP (http:// www. wheat-expre ssion. com).

Phenotypic analysis
The measurements of SL and SL-related traits in the RIL population as well as in the two parental lines are listed in Table 1. As expected, the SL of CM42 was much higher than that of 13F10. In addition, we also found that SD, GW, and TGW differed between 13F10 and CM42; the values of these traits were significantly lower in CM42 than in 13F10 (Table 1, Supplemental Fig. 1). In the RIL population, the continuous distributions of the six traits, together with the skewness and kurtosis for all traits in six environments of less than 1.0 in absolute value, indicated that these traits followed a normal distribution (Table 1, Supplemental Fig. 2). The heritability of the six traits ranged from 75.55 to 93.85%. Genetic factors played key roles in the phenotypic variation of the corresponding six traits in the mapping population. Significant Pearson's correlation coefficients ranged from 0.40 to 0.74 (P < 0.01) for the SL trait and from 0.53 to 0.79 (P < 0.01) for the SD trait in the six different environments (Supplemental Table 1). These results suggest that these traits are typically quantitatively controlled and suitable for QTL mapping.
The genetic correlation coefficients between the traits studied using the BLUE datasets are presented in Table 2. Correlation coefficients of the six traits ranged from − 0.843 to 0.898. Higher positive correlation coefficients were observed between SL and GL, SNS and SD, GL and GW, TGW and GL, TGW and GW, and higher negative correlation coefficients were found between SL and SD, SNS and GL, SNS and GW, SNS and TGW, SD and GL, SD and TGW. Among these traits, the significant negative correlation was observed between SL and SD.

BSE-Seq analysis
Exome capture and high-throughput sequencing of the two bulked pools (SL-H, SL-D) and the parents (13F10, CM42) were performed, and the results were Table 1 Parental and population means, ranges, and broad sense heritability of SL, SNS, SD, GL, GW, and TGW based on the BLUE dataset a Traits are spike length (SL), spikelet number per spike (SNS), spike density (SD), grain length (GL), grain width (GW), and thousand grain weight (TGW); b 13CM population represent '13F10/CM42'; H 2 : broad-sense heritability; ** Significant at P < 0.01, * Significant at P < 0.05  Table 2). These results indicated that these high-quality data can be used for subsequent SNP/InDel detection and analysis. After SNP/InDel variant calling with BCFtools, a total of 2,641,275 and 3,253,239 SNPs were obtained from the two mixed pools, of which 409,424 SNPs were synonymous (Supplemental Table 3). Furthermore, 316,148 InDels from parents and 352,585 InDels from the mixed pools were identified during the analysis, and 52,435 InDels were nonsynonymous among the four samples (Supplemental Fig. 3). As two common methods for BSA analysis, the Euclidean distance (ED) algorithm and SNP/InDel index algorithm were used to determine the candidate region for spike length. With the SNP-index algorithm, a total of 883 SNPs (approximately 55.67% of all SNPs) were identified on chromosome 5A between two extreme bulks under the significance threshold (LOESS fit of Δ(SNP-index) > 0.7416), of which 765 SNPs were concentrated in the 510-540 Mb chromosomal region (Fig. 1a, b). Based on this result, only one QTL was identified, and it mapped to chromosome 5A (Fig. 1c). Additionally, this QTL for SL was detected by ED analysis (Supplemental Fig. 4). Taking overlapping regions into account, the main peak of the different SNPs was found to be located in a physical region of 6.69 Mb (from 518.43 Mb to 525.12 Mb) on chromosome 5A (Supplemental Table 4). SNPs and InDels identified between 13F10 and CM42 were then screened against the SL-H and SL-D bulked exome sequences. Finally, a total of 88 high-quality SNPs and 7 InDels were found to span the candidate interval.

Linkage map construction and QTL analysis
To confirm the preliminarily identified QTL based on the extreme mixed pool of SL phenotypes, polymorphic SNPs in the target region of the parents were converted to KASP markers for QTL analysis. Furthermore, SSR markers were selected with parental polymorphisms in the candidate region by referring to the Chinese Spring genome IWGSC RefSeq v2.0. The polymorphic KASP markers and SSR markers in the whole 13CM population were genotyped. In total, 15 KASP markers (Supplemental Table 5) and 4 SSR markers (Supplemental Table 6) were successfully used for genetic map construction by genotyping 186 RILs. The resulting linkage map represented a segment of chromosome 5AL spanning 21.16 cM in length (Fig. 2a). Based on this genetic linkage map, a major additive QTL (named QSl.cib-5A) was mapped between markers SSR1017544 and SSR1018075 with logarithm of the odds (LOD) values from 3.38 to 12.68. QSl.cib-5A was located in a 4.84 cM genetic interval corresponding to 4.67 Mb (516.60-521.27 Mb) physical interval with the longspike allele from CM42 explaining 7.88-26.60% of the phenotypic variance in the six tested environments, and the BLUE data were calculated (Supplemental Table 7). No epistatic QTL was identified in this study. In addition, a QTL involved in SL was detected by multi-environmental analysis. The QTL was in the same marker interval as that detected in individual environment analysis, indicating that the loci interact with the environment (Supplemental  Table 8). However, the phenotypic variance between environments in the joint analysis was smaller than the phenotypic variance of the additive effect in a single environment, suggesting that the contribution of the additive effect to the phenotypic variation played a main role. Thus, the QTL was stably expressed and mainly controlled by genetic factors. SL was strongly and negatively correlated with SD (r = − 0.843, P < 0.01). QTL mapping for SD was conducted using the linkage map to determine whether the QTL region had significant effects on SD. A major QTL explaining 11.35-19.93% of the phenotypic variation was identified, with 13F10 contributing the favorable allele.
To verify the mapping region, the phenotype and location of recombination events were compared to the reference interval. By comparing the SL phenotype between plants with the 13F10 genotype and plants with the CM42 genotype, we suggest that QSl.cib-5A was closely linked to the SSR marker SSR1017719 and KASP marker A014150 (Fig. 3, Supplemental Table 7). To verify the influence of QSl.cib-5A, the markers SSR1017719 and A014150 were tested in the 13CM population. Finally, 193 RILs were categorized into two groups according to the genotypes of the two flanking markers. Student's t test showed that the SL of the lines carrying the favorable alleles from 'CM42' was significantly higher than that from '13F10' in all environments as well as in the BLUE dataset. In contrast, the lines with the favored allele (QSd.cib-5A) from '13F10' showed significantly increased SD in the six different environments (Supplemental Fig. 5). However, no significant effects were detected in the SNS, TGW, GL. or GW traits which contributed by QSl.cib-5A (Supplemental Fig. 6).

Analysis of candidate genes in QSl.cib-5A
The candidate genes in QSl.cib-5A were mapped to a 4.67-Mb physical interval, which contained 17 highconfidence genes with 25 transcripts. We analyzed sequence variation in the region between 13F10 and CM42 using BSE-seq data. A total of 17 SNPs and InDels in each parent homozygous for different alleles were identified. Analysis of the SNP/InDel effects suggested that two SNPs at the exonic level were nonsynonymous mutations (Supplemental Tables 10 and  11). One SNP in TraesCS5A03G0749800 located at the 686th position caused a change in the amino acid threonine in 13F10 to serine in CM42 and was predicted to have a large effect on gene function.

Discussion
Synthetic hexaploid wheat (SHW) lines carrying elite genes that are lost during wheat polyploidization show outstanding traits in term of disease resistance, abiotic stress tolerance, suitable quality, and anti-sprouting (Yu et al. 2014a). CM42, showing large kernels, great spike length, and resistance to stripe rust, was the first widely used variety derived from SHW in China wheat breeding programs (Yang et al. 2009). SL and SD are the most critical components of spike morphology and determine plant architecture and yield potential. Therefore, mapping genes/ QTLs controlling SL/SD in CM42 may provide new insight into spike architecture development.
Selective genotyping of individuals associated with the high and low phenotype distribution tails provides an effective approach for QTL mapping (Yang et al. 2020). In this study, a BSE-Seq approach was carried out to mine the QTLs associated with SL. As an effective method, BSE-Seq has been successfully used for mapping and gene identification in wheat. However, the traits identified by BSE-Seq to date are mainly controlled by a single gene, and few results have been reported for complex traits in RIL populations (Martinez et al. 2020;Dong et al. 2020;Breseghello and Sorrells 2006;Harrington et al. 2019). To determine the usefulness of BSE-Seq data in complex trait gene mapping, SL throughout an RIL population was analyzed. The result was a main peak   (Fig. 1), which confirmed the effectiveness of mixing pools constructed with representative samples of extreme phenotypes in a RIL population (Martinez et al. 2020;Dong et al. 2020). BSE-Seq led to a reduction in the data analyzed and in the cost related to sequencing a large genome, especially in wheat. As a result, BSA experiments using exome sequencing data led to the identification of major QTL loci associated with complex traits in wheat (Mo et al. 2018).
In the present study, QTLs for the SL trait were mapped to chromosome 5AL. Previously identified QTLs associated with the SL trait on chromosome 5AL are summarized in Supplemental Table 9. The IWGSC (2018) Chinese Spring reference sequence was used as a common coordinating system for comparison of QTLs identified in different studies. According to the physical position of the flanking markers, several QTLs were previously reported to be on chromosome 5AL, which is close to or overlaps the interval identified in our study (Luo et al. 2016;Zhai et al. 2016;Fan et al. 2019;Li et al. 2018). For example, Zhai et al. (2016) used the Y8679/J411 RIL population to identify a major and stable QTL cluster for SL and SD, in which the physical interval is located in a region of chromosome 5AL corresponding to about 307-439 Mb in the IWGSC reference sequences, partially overlaps with our reference interval. R5A, a major QTL qSl-5A.3 and its colocalized major QTL for SC (qSc-5A.2), was previously mapped to the interval of AX-110071854-AX-111139819 (corresponding to an interval of 479.12-541.91 Mb), which also covers the region identified in our current study (Fan et al. 2019). However, these aforementioned QTLs span a large physical interval as determined on the basis of the reference genome sequence. In addition, two QTL clusters (wPt-9094-wPt-9513 and Xgpw4457-Xgpw2273) on chromosome 5A were identified in an RIL population derived from a cross between Tibetan semiwild wheat Q1028 and common wheat Zhengmai 9023 (Luo et al. 2016). The QTL in Xgpw4457-Xgpw2273 was confirmed to be the same as the Q gene (Xu et al. 2018). However, the first QTL cluster, located in wPt-9094-wPt-9513, which corresponding to about 537 Mb in the IWGSC RefSeq v2.0, nearby our QTL QSl.cib-5A. Due to the position of the vernalization locus Vrn-A1, the author of the previous study suggested that this QTL may be the Vrn-A1 allele. In this study, KASP marker (Vrn-A1) was used to genotype the parents and RILs. The results indicated that there was no significant phenotypic difference in SL and SD was detected between the two alleles although Vrn-A1 had polymorphic in the parents and RILs. In addition, this KASP marker of Vrn-A1 was significantly correlated with GW in our population (Supplemental Fig. 8). These observations suggest that Vrn-A1 is unlikely to be the candidate gene for this QTL. The Q gene (651.81-651.82 Mb) is located on chromosome 5AL and encodes a member of the AP2 transcription factor, and microRNA target site variation increases Q allele expression (Liu et al. 2018b). Comparisons of our genetic map with those previous reports revealed that the QTL identified in our study was located far away from the Q gene (~ 130.54 Mb). Moreover, the Q gene of the two parents was sequenced, and the results showed that this gene in CM42 was the same as that in 13F10, with both carrying the Q locus (Supplemental Fig. 9). As a result, QSl.cib-5A is most likely a novel locus that warrants further investigation. Gene annotation of the corresponding 4.67 Mb genomic region in Chinese Spring RefSeq v2.1 revealed 17 high-confidence genes with 25 transcripts. According to gene annotation and analysis of the effect on SNPs in the mapped region (Supplemental Table 10), a nonsynonymous mutation in TraesCS5A03G0749800.1 encoding an AEE5 protein might have a great impact on gene function. In plants, the superfamily of carboxyl-CoA ligases and related proteins, collectively called acyl activating enzymes (AAEs), has evolved to include enzymes for many pathways of primary and secondary metabolism and for the conjugation of hormones to amino acids (Shockey and Browse 2011). Moreover, based on the expression profiles of expVIP, six genes were abundantly expressed in spike tissue during spike development (Fig. 2b). Among these candidates, TraesC-S5A03G0746000.1 is a transcription factor in the MYC2 family, which is one of the largest transcription factor families involved in a wide and diverse array of biological processes. A series of reports indicated that MYC2 participates in the regulation of plant growth and development, including morphogenesis, flowering time and fruit development (Wang et al. 2017;Kazan and Manners 2013;Wei and Chen 2018;Bai et al. 2019). TraesCS5A03G0748800.1, a cell division cycle 5-like protein (CDC5), is a DNA-binding protein of the MYB3R-and R2R3-type. It plays an important role in cell division and cell expansion, which are related to the loosening and rebuilding of cell walls (Xie et al. 2018b). TraesCS5A03G0748900.1 is a serine carboxypeptidase-like (SCP) gene that may control the width and filling of rice grains by regulating plant cell development (Li et al. 2011). Further work on these genes is necessary to identify the genes involved in the regulation of SL and SD.