Identification and Validation of a Novel Locus Controlling Spikelet Number in Bread Wheat (Triticum aestivum L.)

Spikelet number is an important target trait for wheat yield improvement. Thus, the identification and verification of novel quantitative trait locus (QTL)/genes controlling spikelet number are essential for dissecting the underlying molecular mechanisms and hence for improving grain yield. In the present study, we constructed a high-density genetic map for the Kechengmai1/Chuanmai42 doubled haploid (DH) population using 13,068 single-nucleotide polymorphism (SNP) markers from the Wheat 55K SNP array. A comparison between the genetic and physical maps indicated high consistence of the marker orders. Based on this genetic map, a total of 27 QTLs associated with total spikelet number per spike (TSN) and fertile spikelet number per spike (FSN) were detected on chromosomes 1B, 1D, 2B, 2D, 3D, 4A, 4D, 5A, 5B, 5D, 6A, 6B, and 7D in five environments. Among them, five QTLs on chromosome 2D, 3D, 5A, and 7D were detected in multiple environments and combined QTL analysis, explaining the phenotypic variance ranging from 3.64% to 23.28%. Particularly, QTsn/Fsn.cib-3D for TSN and FSN [phenotypic variation explained (PVE) = 5.97–23.28%, limit of detection (LOD) = 3.73–18.51] is probably a novel locus and located in a 4.5-cM interval on chromosome arm 3DL flanking by the markers AX-110914105 and AX-109429351. This QTL was further validated in other two populations with different genetic backgrounds using the closely linked Kompetitive Allele-Specific PCR (KASP) marker KASP_AX-110914105. The results indicated that QTsn/Fsn.cib-3D significantly increased the TSN (5.56–7.96%) and FSN (5.13–9.35%), which were significantly correlated with grain number per spike (GNS). We also preliminary analyzed the candidate genes within this locus by sequence similarity, spatial expression patterns, and collinearity analysis. These results provide solid foundation for future fine mapping and cloning of QTsn/Fsn.cib-3D. The developed and validated KASP markers could be utilized in molecular breeding aiming to increase the grain yield in wheat.


INTRODUCTION
Bread wheat (Triticum aestivum L.) is one of the most important stable foods in the world, providing more than 20% of the calories and protein in our daily diet (Yao et al., 2019;Liu et al., 2020). It was reported that world population is expected to reach nine billion by 2050, which will require raising overall food production by at least 70% to fulfill future food demand (Yang et al., 2012). Conversely, arable land and water resources for agriculture continue to dwindle, and climate change also causes adverse effects on crop productivity (Evans et al., 2008;Butterworth et al., 2010;Ceccarelli et al., 2010). Thus, increasing wheat yield is an urgent task to fulfill global food and nutritional security.
Grain yield of wheat is a complex trait affected by a multitude of genetic and environmental factors and is usually constituted of three major components, i.e., thousand kernel weight (TKW), grain number per spike (GNS), and spike number per unit area (Wang et al., 2018;Kuzay et al., 2019). GNS potential is largely established by the spike architecture. Similar to other Triticeae species, wheat contains an unbranched inflorescence, and a number of spikelets are directly attached to the inflorescence axis in a distichous phyllotaxis, with a terminal spikelet at its apex (Koppolu and Schnurbusch, 2019). The spikelet is made of an indeterminate number of florets attached to a secondary axis (rachilla) (Hanif and Langer, 1972;Wolde et al., 2019). Therefore, the number of fertile spikelet number per spike (FSN) and fertile floret per spikelet ultimately determined GNS (Guo et al., 2017;Golan et al., 2019). Previous studies have showed that floret fertility is severely affected by environmental and genetic factors as an abortion process of floral structures existing continuously during the whole floral developmental process (Guo and Schnurbusch, 2015;Sakuma et al., 2019). By comparison, spikelet number is determined at an early stage of wheat reproductive development, exhibiting less environmentally sensitive with high heritabilities Gao et al., 2019). Thus, understanding the genetic factors underlying spikelet number would provide the prerequisite information necessary to improve wheat yield.
Similar to other yield-related traits, spikelet number is a complicated quantitative trait [quantitative trait locus (QTL)] controlled by multiple genes in wheat. Thus, determination of the number, chromosomal localization, and genetic effects of these polygenes is desirable for obtaining optimal genotype in breeding practice (Ma et al., 2007;Chen et al., 2020). However, due to the complexity of wheat genome, only a few genes associated with spikelet number or spike morphology have been cloned through homologous cloning approach. For instance, WAPO1, an ortholog of the rice gene APO1, on chromosome 7AL could regulate spikelet number by effect inflorescence development (Kuzay et al., 2019;Voss-fels et al., 2019). The domestication gene Q on chromosome 5A encodes a member of AP2 transcription factor family and can regulate rachis fragility, glume tenacity, head length, and spikelet density (Faris et al., 2003;Simons et al., 2006;Debernardi et al., 2017;Greenwood et al., 2017;Jiang et al., 2019). The photoperiod sensitivity gene Ppd-1 on chromosome 2D is associated with multiple pollinator traits, including spikelet number, spike length (SL), and number of days to heading (Okada et al., 2019); TaMOC1, an ortholog gene of rice MOC1, is significantly associated with spikelet number ; FZP gene can drive supernumerary spikelet (Dobrovolskaya et al., 2015); TaDEP1, an ortholog of rice DEP1, relates to increase SL and reduces the spikelet number (Huang et al., 2009).
In addition, numerous QTLs for spikelet number have been identified on almost all the 21 chromosomes in wheat using different genetic populations and natural varieties (Cui et al., 2012;Yu et al., 2014;Gao et al., 2015;Luo et al., 2016;Zhai et al., 2016;Deng et al., 2017;Liu et al., 2018b;Wang et al., 2018;Koppolu and Schnurbusch, 2019;Ma et al., 2019;Yao et al., 2019;Kuang et al., 2020). Nevertheless, only a few of QTLs have been genetically verified, offering a foundation for fine mapping and map-based cloning, which greatly restrict the dissection of the molecular mechanisms underlying spikelet number as well as improvement of spikelet number in wheat breeding. Thus, the identification and verification of novel QTL/genes for spikelet number are vital.
In this study, we constructed a high-density genetic map using the Wheat 55K single-nucleotide polymorphism (SNP) array and a doubled haploid (DH) population derived from a cross between two wheat elite varieties K1 and CM42. QTL mapping for two spikelet number-related traits, total spikelet number per spike (TSN), and FSN was performed, and subsequently, a major QTL on chromosome 3DL was further validated in different genetic backgrounds.

Plant Materials and Field Trials
Quantitative trait locus mapping was conducted using a population of 187 DH lines derived from a cross of K1 × CM42 (K1/CM42). CM42, kindly provided by Prof. Wuyun Yang (Sichuan Academy of Agricultural Science, China), is a backbone parent with excellent agronomic traits such as high TKW and high yield (Yang et al., 2009). K1, bred by our lab, is characterized by high spikelet number per spike and GNS. Two populations, including a recombinant inbred line (RIL) populations derived from a cross of K1 × Kechengmai4 (K1/K4, 70 F 7 RILs) and a F 2 populations derived from a cross of K1 × Yangfumai2 (K1/YFM2, 75 F 2 lines), were employed to further validate the major QTL.
The K1/CM42 lines along with parents were grown in five replicates following randomized complete-block design at the following two experimental sites: Shuangliu of Sichuan province (103 • 52 E, 30 • 34 N) and Shifang of Sichuan province (104 • 11 E, 31 • 6'N), China. The field trials were carried out in two crop seasons of 2017-18 and 2018-19 at Shuangliu, and three crop seasons of 2016-17, 2017-18, and 2018-19 at Shifang. Each plot consisted of five 200-cm rows with an inter-row spacing of 20 cm. The sowing density was 50 seeds per row. The K1/K4 lines were generated by single-seed descent and cropped at Shuangliu and Shifang in 2018-19 crop seasons. Each plot consisted of two 200cm rows with an inter-row spacing of 20 cm. The K1/YFM2 population was cropped at Shuangliu in 2019-20 crop seasons in 200-cm rows with 15-cm space between individuals.

Phenotypic Evaluation and Statistical Analysis
The TSN, FSN, plant height (PH), and SL of the K1/CM42 and K1/K4 populations were manually measured at maturity by using plants from middle row of each plot. For the K1/CM42 population, the primary spike of 10 random plants from each family in 2017-18 and 2018-19 crop seasons and three random plants from each family in 2016-17 crop seasons were selected. For the K1/K4 population, the primary spikes of five random plants from each family were selected. For the K1/YFM2 population, the TSN, FSN, PH, and SL of primary spike of each individual were measured. Subsequently, all selected spikes were harvested and manually threshed for evaluating GNS, TKW, grain width (GW), and grain length (GL) using SC-G software (WSeen, Hangzhou, China).
The frequency distribution of traits in each environment and their correlation were calculated using SPSS version 20.0 for Windows (IBM SPSS, Armonk, NY, United States). The best linear unbiased predictions (BLUPs), which were used for QTL detection, correlation analyses, and effect analyses, were calculated using the R package "lme4." Analysis of variance (ANOVA) was performed using the data from 2017-18 and 2018-19 crop seasons in the K1/CM42 population by SPSS. The broad-sense heritability (H 2 ) was calculated according to the method described by Smith et al. (1998) and Muqaddasi et al. (2019). Student's t test used to evaluate the significance of difference was performed by SPSS.

Genotyping
Genomic DNA of each line in the K1/CM42 population and their parents was extracted using the 2 × cetyl trimethylammonium bromide (CTAB) method and hybridized on the Wheat 55K SNP array by China Golden Marker (Beijing, China). DNA integrity was checked and confirmed on agarose gels, and DNA quantity was measured by spectrophotometry.

Genetic Map Construction and Quantitative Trait Locus Detection
IciMapping 4.1 (Meng et al., 2015) and JoinMap 4.1 were used for genetic construction and QTL detection in the present study. First, the function of "bin" in IciMapping 4.1 was used to remove redundant markers based on their segregation patterns in the mapping population with the parameters of "Missing Rates" and "Distortion Value" being set as 20% and 0.01, respectively. Then, the function of "Population" in JoinMap 4.1 was used to create groups with limit of detection (LOD) score values ranging from 2 to 10. Finally, the Kosambi mapping function was used to order the bin markers with the parameters being set as LOD ≥ 5 and round = 3 in JoinMap 4.1. To reduce the complexity of QTL mapping analyses, only one marker was selected as a delegate from each bin to construct the linkage map. QTL detection in each environment was performed by IciMapping 4.1 with the inclusive composite interval mapping (ICIM), and a test of 1,000 permutations was used to identify the LOD threshold that corresponded to a genome-wide false discovery rate of 5% (P < 0.05). The missing phenotype was deleted in QTL analysis.

Development of Kompetitive Allele-Specific PCR Markers
Probe DNA sequences for selected SNP were subjected to design Kompetitive Allele-Specific PCR (KASP) markers using the Triticeae Multi-omics Center 1 . Sequence probes for the FAM signal and the HEX signal were added to the primers of two parental genotypes, respectively. The KASP assays were performed in Bio-Rad CFX96 real-time PCR system with 10-µl reaction volumes consisting of 5 µl of 2 × master mix, 0.2 µl of primer mix, 3 µl of ddH 2 O, and 2 µl of DNA sample (50-150 ng/µl). Thermal cycling conditions were 94 • C for 15 min hot-start activation, followed by a touchdown phase of 10 cycles (94 • C for 20 s, touchdown at 61 • C initially and then decreased by 0.6 • C per cycle for 60 s) and finally 26 cycles of regular PCR (94 • C for 20 s, 55 • C for 60 s, and rest plates at 37 • C for 1 min). Further cycling and resting were performed if the clustering is not significant: 94 • C for 20 s, followed by 57 • C for 60 s (3-10 cycles per step). The information of DNA sequences flanking the SNP was kindly provided by Prof. Jizeng Jia's Lab (Institute of Crop Sciences, Chinese Academy of Agricultural Sciences).

Prediction of Candidate Gene
The physical positions of flanking markers were obtained by blasting against (E-value of 1e-10) the genome assembly of Chinese Spring (CS) 2 . Genes among the mapping interval were extracted from IWGSC RefSeq v1.1 annotation 3 . The annotations and functions of a given gene were analyzed using UniProt 4 . Gene collinearity among wheat, barley, rice, and maize was performed using the Triticeae-Gene Tribe 5 . The temporal and spatial gene expression patterns were extracted from Wheat Expression Browser 6 . To further explore the potential candidate genes, we designed the corresponding specific primers 1,7 to isolate and sequence these candidate genes (Tsingke Biotech Co., Ltd., China). The primer sequences were listed in Supplementary  Table 2. Twenty-five-microliter PCR reaction volumes were as follows: 12.5 µl of 2 × EasyTaq PCR SuperMix (TransGen Biotech Co., Ltd., China), 3 µl of DNA sample, 1 µl of primer mix, and 7.5 µl of ddH 2 O.

Phenotypic Variation
Significant differences (P < 0.001) on TSN and FSN between K1 and CM42 across five environments and their corresponding BLUP data were detected. K1 showed more TSN and FSN than did parent CM42 (Figure 1 and Table 1). In addition, we observed the development of spikelet, and we found that K1 showed more TSN at a very early stage of spike development than that of CM42, as well as at heading and maturity stages, indicating that spikelet number is determined at an early stage of wheat reproductive development (Figure 1). For the K1/CM42 lines, the TSN and FSN ranged from 16.3 to 26.8 and from 14.2 to 25.8, respectively, showing wide and significant variation and approximately normal distribution (Table 1 and Figure 2). In addition, the TSN and FSN values across all environments and BLUP data were significantly positively correlated (P < 0.001), showing high broad-sense heritabilities (0.85 for TSN and 0.82 for FSN), indicating that they were mainly under genetic control ( Table 1).
Phenotypic correlations between TSN, FSN, and other yield related traits are listed in Table 2. Significant and positive correlations (P < 0.001) between TSN and FSN and GNS, and significant and negative correlation (P < 0.01) between TSN and GL, and FSN and GL were observed in the K1/CM42 population. The highest positive correlation was detected between TSN and FSN (r = 0.96). No significant correlation was detected between TSN and TKW, TSN and PH, TSN and SL, FSN and TKW, FSN and PH, and FSN and SL were observed in the K1/CM42 population.

Linkage Map Construction
Among 53,063 SNPs from the Wheat 55K SNP array, 14,645 SNP markers (27.6%) were polymorphic between the two parental lines. Of them, 822 SNPs had more than 20% missing data points, 10 SNPs were not anchored on the physical linkage map, and 745 SNPs were not anchored on the genetic linkage map. Finally, the rest 13,068 SNP markers (call rate ≥ 80%) were used for linkage analysis and map construction. These markers were divided into 2,406 bins: 1,384 bins contain one SNP marker only, 1,022 bins contain multiple SNPs, and the largest bin contains 626 SNP markers on 4B. Linkage analysis showed that 2,406 bin markers were mapped on 27 linkage maps for the 21 chromosomes of wheat, in which two linkage groups were constructed for each of chromosomes 1A, 1D, 5A, 5B, 6A, and 6D. Based on the genetic information, all SNP markers including bin markers and redundant markers were integrated onto the genetic map with a total length of 3,091.39 cM and an average interval distance of 1.28 cM per bin and 0.24 cM per marker. All of the mapped markers were located to the A (33.75%), B (45.98%), and D (20.26%) genomes with a total length of 1,009.07, 895.54, and 1,186.79 cM, respectively. In addition, the lengths of constructed linkage maps range from 17.67 (1A-1) to 206.47 cM (5D); the markers on linkage maps range from 20 (1D-2) to 1228 (4B); the average interval distance between adjacent markers ranges from 0.08 (4B) to 4.37 cM (1D-2); and the average interval distance between adjacent bin markers ranges from 0.63 (2B) to 5.14 cM (1D-2) ( Table 3 and Supplementary Table 1). Based on the physical locations of these SNPs on the CS genome, the marker order was relatively consistent with that in the wheat genome assembly on most of the chromosomes (1A, 1B, 2A, 3B, 3D, 4A, 4B, 5A, 5B,5D, 6A, 6D, 7A, 7B, and 7D) (Supplementary Figure 1).

Quantitative Trait Locus Mapping
Phenotypic data of TSN and FSN evaluated in five environments and their corresponding BLUP values were used for QTL mapping, where the BLUP values were treated as an additional environment. A total of 27 QTLs for TSN and FSN were detected, which were located on chromosomes 1B, 1D, 2B, 2D, 3D, 4A, 4D, 5A, 5B, 5D, 6A, 6B, and 7D, respectively. Among them, five QTLs defined as stable were detected in multiple environments and combined QTL analysis, explaining the phenotypic variance ranging from 3.64% to 23.28%. The rest QTLs were detected in one environment, explaining 2.47-18.4% of the phenotypic variance (Table 4).
For TSN, 14 QTLs were detected on chromosomes 1B, 1D, 2B, 2D, 3D, 5A, 5B, 5D, 6A, 6B, and 7D. Of them, a major QTL QTsn.cib-3D was stably detected in four environments and combined QTL analysis, explaining 6.5% to 23.15% of the phenotypic variance with LOD values ranging from 6.1 to 18.51. Another stable QTL QTsn.cib-2D was detected in two environments and combined QTL analysis, explaining 3.64-10.23% of the phenotypic variance. The positive alleles of the two loci were all contributed by K1. For the remaining 12 QTLs, they were detected in one environment and explained 2.47-15.04% of the phenotypic variance. The positive alleles of eight loci were from K1, and four loci were from CM42 ( Table 4).
For FSN, 13 QTLs were identified on chromosomes 1B, 1D, 2D, 3D, 4A, 4D, 5A, 5B, 5D, and 7D. Of them, QFsn.cib-3D was a major QTL and stably detected in four environments and combined QTL analysis, explaining 5.97-23.28% of the phenotypic variance with LOD values ranging from 3.73 to 14.94. QFsn.cib-7D was detected in two environments and combined QTL analysis, explaining 7.51-12.58% of the phenotypic variance. The positive alleles of the two loci were contributed by K1. QFsn.cib-5A was detected in two environments with the phenotypic variance ranging from 5.86% to 18.4%, and the positive allele was contributed by CM42. The rest of the 10 QTLs were detected in single environment, accounting for 4.42-18.26% of the phenotypic variance. The positive alleles of seven and three loci were contributed by K1 and CM42, respectively (Table 4).

Validation of QTsn/Fsn.cib-3D in Different Genetic Backgrounds
Since the favorable allele of QTsn/Fsn.cib-3D was from K1, the other two populations (K1/K4 and K1/YFM2) containing a common parent K1 were used to validate the effects of QTsn/Fsn.cib-3D in different genetic backgrounds. The KASP marker KASP_AX-110914105 closely linked to QTsn/Fsn.cib-3D was used to assess the alleles from the common parent K1. For the K1/K4 population, 27 lines were found to be homozygous for the K1 allele, and the remaining 41 lines were found to be homozygous for the non-K1 allele (Figure 5). For the K1/YFM2 population, 17 plants were found to be homozygous for the K1 allele, 19 plants were found to be homozygous for the non-K1 allele, and the remaining 39 plants were found to be heterozygous allele ( Figure 5). Subsequently, Student's t-test was performed to analyze significant differences on TSN and FSN (P < 0.05) between classes with different allele. As expected, the lines with the homozygous K1 allele had significantly (P < 0.01) higher TSN (5.56-7.96%) and FSN (5.13-9.35%) than those carrying the non-K1 homozygous allele ( Figure 5). Moreover, a significant difference on TSN (6.62%) and FSN (7.71%) was also detected between lines with K1 homozygous allele and heterozygous allele (Figure 5).

Candidate Genes Analysis of QTsn/Fsn.cib-3D
Alignment of the flanking markers of QTsn/Fsn.cib-3D showed that it was corresponding to a physical interval of 549.49-555.2 Mb on chromosome arm 3DL (Figure 6 and Table 4). According to CS reference genome, there were 76 predicated genes in this interval (Figure 6). An analysis of the spatial expression patterns (Borrill et al., 2016;Ramírez-González et al., 2018) showed that 40 genes were expressed in various tissues, and several of them were abundantly expressed in spikes, indicating that they are supposed to be involved in spike development (Figure 6). Among these genes, TraesCS3D02G439000, TraesCS3D02G439200, TraesCS3D02G442000, TraesCS3D02G443900, and TraesCS3D02G445400 are likely associated with spike growth and development according to the gene annotation (Supplementary Table 3). To further analyze potential polymorphism of these candidate genes between K1 and CM42, and K4 and YFM2, we designed gene-specific primers to isolate and resequence the five genes (Supplementary Table 2). For TraesCS3D02G439000, TraesCS3D02G439200, and TraesCS3D02G442000, no sequence variation in coding sequences (cds) between K1 and the other three parents was detected (Supplementary Figures 4-6).
For TraesCS3D02G443900, three SNPs in the cds between K1 and the other three parents were detected, and two of them were non-synonymous SNPs (Figure 7 and Supplementary  Figure 7). Interestingly, we isolated and assembled the complete coding sequences of TraesCS3D02G445400 from K1 and YFM2 using these designed specific primers. However, only partial coding sequences could be obtained from CM42 and K4 (Supplementary Figure 8). We speculated that significant differences may exist in the sequences of TraesCS3D02G445400 between K1 and CM42 and K4. Based on these sequences, two synonymous SNPs in the coding sequence between K1 and the other three parents, and extra six synonymous SNPs and three non-synonymous SNPs between K1 and YFM2 were detected, respectively (Figure 7 and Supplementary Figure 8).

DISCUSSION
Grain number per spike as one of three major components determining wheat yield is ultimately determined by the number of FSN and fertile floret per spikelet. In the present study, by assessing two spikelet number-related traits TSN and FSN in a DH population among five environments, we found that TSN and FSN were significantly and positively correlated with GNS ( Table 2) with a high heritability (0.85 for TSN and 0.82 for FSN) ( Table 1). The results were consistent with previous knowledge (Che et al., 2018;Würschum et al., 2018;Zhang et al., 2019) that they are environmentally stable yield components and mainly determined by genetic factors. Therefore, increasing spikelet number could be considered as an effective strategy for increasing GNS and hence the grain yield.

Comparison of Constructed Genetic Linkage With Previous Studies
Quantitative trait locus analysis is a well-established and widely used tool for dissecting the genetic basis of complex traits, and a high-density genetic map plays a fundamental role in QTL analysis (Salvi and Tuberosa, 2005;Cooper et al., 2009;  Note. Add, additive effect (positive values indicate that alleles from CM42 are increasing the trait scores, and negative values indicate that alleles from K1 are increasing the trait scores); Chr., chromosome; Pos., position; PVE, phenotypic variation explained; LOD, limit of detection. Yamamoto et al., 2009). In the present study, by using the wheat 55K SNP array, a high-density genetic map containing 13,068 SNPs with a total length of 3,091.39 cM was constructed. The marker order was relatively consistent with that in the wheat genome assembly on most of the chromosomes (Supplementary Table 1 and Supplementary Figure 1). It is well known that the SNP markers in the Wheat 55K SNP array were well chosen from the 660K SNP array. Compared with previously reported genetic maps, we extracted 2,739 and 4,860 common mapped markers from the Wheat 55K and 660K SNP array, respectively (Supplementary Table 4) Liu et al., 2018a). The genetic positions of these markers were relatively consistent in the three genetic maps (Supplementary Figures 2, 3), indicating that the genetic map constructed in present study was accurate and credible and could be effectively used for QTL analysis.

Comparison of Stable Quantitative Trait Loci for Spikelet Number With Previous Studies
Over the past decades, QTLs that control spikelet number have been extensively studied and identified on almost every wheat chromosome. In this study, five stable QTLs for TSN and FSN were detected on chromosomes 2D, 3D, 5A, and 7D. Of them, QTsn/Fsn.cib-3D associated with TSN and FSN was detected as a major QTL for spikelet number. This QTL was mapped in a 4.455-cM interval and physically corresponding to 549.49 and 555.2 Mb on chromosome arm 3DL according to the physical location of flanking markers (Table 4 and Figure 3). As far as we know, several QTLs for spikelet number on chromosome 3D including QSns.sau-3D, QSsn.czm-3D, QSspn.cau-3D, and QSpn.ipk-3D had been reported in previous studies. QSsn.czm-3D and QSspn.cau-3D were located on the short arm of chromosome 3D flanked by 3DS_2577014_1698-IWA4559 and Xcau.3D-5-Xcau.3D-6, respectively (Cui et al., 2012;Chen et al., 2020).
Effects on Yield-Related Traits and Potential Use of QTsn/Fsn.cib-3D Similar to previous studies Cui et al., 2012), TSN is significantly and positively correlated with FSN.
Thus, QTLs, especially the major QTLs, for TSN and FSN were usually identified in same locus. For instance, seven stable QTLs for spikelet number were identified by Chen et al. (2020), and five of them were associated with TSN and FSN simultaneously. The loci of QTsn/Fsn.czm-1A and QTsn/Fsn.czm-1D were identified by Zhou et al. (2017), showing significant effects on TSN and FSN. A QTL cluster for FSN and TSN identified on chromosome 5D was located in the same region as that detected by Li et al. (2007) and Cui et al. (2012). Moreover, as one of the subcomponent of GNS, an increase of spikelet number could directly increase GNS (Hai et al., 2008). Therefore, it is not surprising that the major QTL QTsn/Fsn.cib-3D identified in present study showed significant and positive effects on TSN, FSN, and GNS (Figures 4, 5). Additionally, we further analyzed the effects of QTsn/Fsn.cib-3D on other yield-related traits including TKW, SL, PH, GW, and GL in the K1/CM42 populations. QTsn/Fsn.cib-3D showed moderate negative influence on GL, but no influence on TKW, PH, SL, and GW (Figure 4). These results suggested that QTsn/Fsn.cib-3D combined with developed KASP markers has potential application values in wheat breeding.
Potential Candidate Genes for QTsn/Fsn.cib-3D In the interval of QTsn/Fsn.cib-3D, there are 76 confidence genes in CS (Supplementary Table 3). Based on the expression analysis and gene annotation, five genes TraesCS3D02G439000, TraesCS3D02G439200, TraesCS3D02G442000, TraesCS3D02G443900, and TraesCS3D02G445400 are abundantly expressed in spikes, indicating they are likely associated with spike growth and development (Figure 6). Of them, TraesCS3D02G445400 encodes a glutathione S-transferase, which has been reported to be related to panicle and spikelet numbers, and plays an important role during the growth and development processes in rice (Kim et al., 2013;Hu et al., 2014). TraesCS3D02G442000 encodes an auxin response factor. Previous studies have revealed that auxin signal transduction acts critical roles in modulating various biological processes including growth, development, and stress resistance (Qiao et al., 2018;Brauer et al., 2019). TraesCS3D02G439200 is an ortholog of OsHAK5 of the rice and encodes a potassium transporter (Supplementary Figure 9). OsHAK5 plays a crucial role in maintaining rice architecture including PH, root, and tiller by altering cellular chemiosmotic gradients and regulating ATP-dependent polar auxin transport . TraesCS3D02G439000 encodes an F-box family protein, which is known to be involved in the vegetative and reproductive growth and development of many plants ( Van Den Burg et al., 2008;Ma et al., 2019). TraesCS3D02G443900 encodes a BTB domain-containing protein, which was reported to be involved in multiple functions including floral organ development and disease resistance in Arabidopsis, rice, maize, and tomato (Gingerich et al., 2005;Chen et al., 2014;Chern et al., 2014;Xu et al., 2016). To further explore the candidate gene, we isolated and assembled the coding sequences of the five genes from four parents. Sequence alignment showed that two non-synonymous SNPs and one synonymous SNP in coding sequence of TraesCS3D02G443900 were detected between K1 and CM42, and K4 and YFM2; two synonymous SNPs in the coding sequence of TraesCS3D02G445400 were found between K1 and CM42, and K4 and YFM2; and the extra six synonymous SNPs and three non-synonymous SNPs were found between K1 and YFM2, respectively (Figure 7). For the rest three genes TraesCS3D02G439000, TraesCS3D02G439200, and TraesCS3D02G442000, no sequence variation in coding sequences was simultaneously detected between K1 and the other three parents (Supplementary Figures 4, 5, 6). Thus, TraesCS3D02G443900 and TraesCS3D02G445400 could be considered as a focus for further work on fine mapping and gene cloning.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
TL undertook the field trials and subsequent analyses of all available data, including the phenotyping and population genotyping, and drafted this manuscript. GD assisted in the field trials. YT, YS, JW, JC, ZY, XQ, XP, HZ, and JL participated in phenotype measurement. MY discussed results and revised the manuscript. HL and YW designed the experiments, guided the entire study, participated in the data analysis, discussed the results, and revised the manuscript. All authors read and reviewed the manuscript.

ACKNOWLEDGMENTS
We thank Prof. Lifeng Gao of Institute of Crop Sciences, Chinese Academy of Agricultural Sciences, for providing us the sequence information of wheat 55K arrays for KASP marker development. The Triticeae Multi-omics Center (http: //202.194.139.32/) is also acknowledged for providing an integrated platform of tools and genomic data, bringing great convenience to our work.