“SY95-71” (with a pedigree of Fan6//Fan6/Eronga83) was developed by the Triticeae Research Institute of Sichuan Agricultural University in the 1990s, and “CH7034” (with a pedigree of Jing411/Xiaoyan7430//Zhong8601) was developed by the College of Agronomy of Shanxi Agricultural University. The two advanced breeding lines have obvious differences in spike morphology. A population of 184 F_2:10 recombinant inbred lines (RILs) derived from the cross SY95-71 × CH7034 was developed to map genes for important traits.

The 184 RILs of the SY95-71 × CH7034 population were tested in six locations/years, including Chengdu (104°13′E, 30°47′N) in 2014 and 2015, Yuncheng (110°53′E, 35°00′N) in 2016, and Linfen (111°31′E, 36°54′N) in 2015 and 2016. The field experiments were arranged in a randomized complete block design with three replications in all of the trials. Each line was grown in a single 2-m row with 0.25 m between rows. Twenty seeds were planted in each row with plant spacing of 0.1 m. Nitrogen and superphosphate fertilizers in the six environmental trials were applied at sowing at rates of 120 and 130 kg ha^–1, respectively. Field management was carried out with standard agronomic practices according to wheat production throughout the growing season to obtain even crop stands (Ma et al., 2019). The number of spikelet nodes per spike was characterized in the main spikes of 15 randomly selected plants from each line at physiological maturity and was calculated as the mean. In addition, three plants from each RIL were planted in the greenhouse of Oklahoma State University in 2019 to determine whether the supernumerary spikelet phenotype was associated with the QTLs for the SNS trait. To avoid any confusion, the number of spikelet nodes per spike in a regular pattern is referred to as SNS, and the total number of spikelets per spike in a pattern with supernumerary spikelets is referred to as SSS.

The 184 RILs of the SY95-71/CH7034 population and their parents were genotyped using 17K SNP chips with 17,526 markers and 35K DArT SNP chips containing 36,420 markers at Diversity Arrays Technology Pty Ltd. (Canberra, Australia¹). The single-nucleotide polymorphism (SNP) marker sequences were used to BLAST the Chinese Spring sequence [International Wheat Genome Consortium (IWGSC) RefSeq v1.1] to determine the physical positions of each marker. The SNP markers that had less than 20% missing values and contained different recombinant information were selected in the final mapping.

The SNP markers were used to construct linkage groups using JoinMap 4.0 (Van-Ooijen, 2006). The genetic distance between markers was estimated using the Kosambi mapping function, and a QTL for SNSs was claimed when the value of the logarithm of the odds (LOD) threshold for significance was above 2.5 in the interval mapping program in MapQTL 6.0 (Van-Ooijen, 2009). A marker with the maximum LOD value was used to determine the chromosomal position of the QTL peak.

Those genes under or close to the peak of the identified QTLs were considered candidates causing the QTLs. Specific primers were designed to isolate and sequence the candidate genes from each of the two parents to determine the allelic variation.

A diagnostic molecular marker for allelic variation in WFZP-A1 was developed. The primers WFZP-Indel-adF2 (5′-GCCAGCCAACCTCACTTCACTTC-3′) and WFZP-Indel-aR2 (5′-TGGCCGAGGACGCGGCGT-3′) were used to amplify the 449-bp CH7034 allele and the 463-bp SY95-71 allele with the following PCR amplification conditions: 94°C for 3 min, 40 cycles of 94°C for 30 s, 64°C for 30 s, and then 72°C for 40 s, with a final extension of 72°C for 10 min. The PCR products was differentiated on a 2.5% agarose gel at 180 V for 1.5 h.

WAPO-A1 was isolated using the forward primer WAPO-A1F2 (5′-AGGAGATCAATGTTTTTAGGACAATAGGA-3′) and the reverse primer WAPO-A1R6 (5′-GAGCAGGGTGTAGCAACTAGT-3′) with the following PCR amplification conditions: 94°C for 3 min, 40 cycles of 94°C for 30 s, 58°C for 30 s, and then 72°C for 2 min, with a final extension of 72°C for 10 min. A 2-bp deletion was found to distinguish between the two alleles after the PCR products were digested with restriction enzyme HypCH4V. The primers WAPOA1-Indel-F1 (5′-CACAAACATATTGAATGAAAGTTCC-3′) and WAPOA1-Indel-R3 (5′-CATTATATCTTAAACTATATGATGTCC-3′) were used to amplify a fragment containing the 2-bp deletion with the following PCR amplification conditions: 94°C for 3 min, 38 cycles of 94°C for 30 s, 58°C for 30 s, and then 72°C for 30 s, with a final extension of 72°C for 10 min. The digested PCR products were 42 and 70 bp for the CH7034 allele but 110 bp for the SY95-71 allele. The PCR products were separated and scored on an 8% non-denaturing polyacrylamide gel.

WAPO-D1 was isolated using the forward primer WAPO-D1F1 (5′-CAGACGCAGAAAGGAAAACCCTAAGGCAACTT-3′) and the reverse primer WAPO-D1R2 (5′-TGCTAAGAACAGAACAACAAATATACTACCC-3′) with the following PCR amplification conditions: 94°C for 3 min, 40 cycles of 94°C for 30 s, 58°C for 30 s, and then 72°C for 3 min, with a final extension of 72°C for 10 min.

Quantitative real-time PCR (qRT–PCR) was used to determine the transcriptional levels of WAPO-A1 and WAPO-D1 in SY95-71 and CH7034 cells with SYBR Green PCR Master Mix, and actin was used as an endogenous control. Total RNA was extracted from the leaves using TRIzol^® reagent (Invitrogen, Carlsbad, CA, United States). cDNA was synthesized from 1 μg of RNA treated with deoxyribonuclease I using a SuperScript™ II Reverse Transcriptase kit and with an oligo (dT)₂₀ primer (Invitrogen). Primers for the gene transcripts were UFO-A-RT-F2 (5′-CTCACTCACTCTCACTCCACG-3′) and UFO-A-RT-R2 (5′-GGTGGTGAGGCAGTAGGTTC-3′) (Kuzay et al., 2019) for WAPO-A1 and WAPO1-D1-RT-F (5′-CTCACTCTCCCTCCCACCACA-3′) and WAPO1-D1-RT-R (5′-GGTGGTGAGGCAGTAGGTTC-3′) for WAPO-D1. qRT–PCR was performed on a CFX96™ Real-time PCR System (Bio–Rad Laboratories, Hercules, CA, United States) using iQTM SYBR^® Green Supermix (Bio–Rad Laboratories, Hercules, CA, United States), with actin used as an endogenous control. Gene transcript levels were calculated by the 2^–ΔΔCT method, where CT is the threshold cycle.

No supernumerary spikelet on a spikelet node was observed in either of the two parental lines SY95-71 and CH7034 when tested in different years/locations (Figure 1A). However, supernumerary spikelets were observed in some RILs derived from the two parental lines (Figure 1B), indicating that the supernumerary spikelets were generated due to recombination of genes from the two parental lines. The SNS trait in the SY95-71/CH7034 RIL population showed clear segregation, as illustrated in Figure 1C.

Statistical analysis showed that among the 184 RILs the SNS association under the field conditions was significant (T > | t₀.₀₉|, Supplementary Table 1), indicating that major gene(s) might be responsible for SNS segregation in this population, regardless of the difference in changes in climate or environments in the field.

A total of 17,512 SNP markers was eventually generated from the 17K SNP chips, and a total of 36,420 SNP markers was generated from the 35K DArT SNP chips for the 184 RILs of the SY95-71/CH7034 population. According to their sequences, these SNP markers were assembled into 21 chromosomes, forming genetic maps for the RIL population. The total length of the 21 chromosomes containing the 2,347 SNP markers was 4,192 cM, with a marker density of 1.78 cM per marker. Detailed information on the sequence and genetic distances of the SNP markers on the three SNS QTLs mapped in this study is provided in Supplementary Table 2. The SX plus 8-digit SNP codes served as the reference number for each SNP from 17K SNP chips, and the SX-D plus 8-digit SNP codes served as the reference number for each SNP from the 35K DArT chips.

A total of 1,738 SNP markers, including 797 from 17K SNP chips and 941 from 35K DArT SNP chips, were assembled into chromosome 2A, which was confirmed by the IWGSC RefSeq v1.1 sequences. On the basis of whole-chromosome QTL scanning using interval mapping (IM) analysis, a QTL, QSns.sxau-2A, was found on this chromosome (Figure 2A). The LOD value at the peak position of QSns.sxau-2A ranged from 4.89 to 11.13 in six experimental locations/years, and this QTL alone explained 11.65 or 24.20% of the total phenotypic variation. The 48 markers on chromosome 2A were selected as a linkage group to indicate the genetic distance and physical location of the QSns.sxau-2A locus (Figure 2B).

The WFZP-A1 gene at position 66,848 kb was close to marker SX-D1098510, which was under the peak of QSns.sxau-2A (Figure 2C); hence, WFZP-A1 was considered a candidate gene causing QSns.sxau-2A. Sequencing results revealed that the SY95-71 allele showed the same sequence as the WFZP-A1 gene of Chinese Spring, which is referred to as the WFZP-A1a allele, whereas the CH7034 allele had a 14-bp deletion in its coding region, which is referred to as the WFZP-A1b allele (Figure 2D). The 14-bp insertion in the WFZP-A1a allele could cause it to gain its function at the protein level, or the 14-bp deletion in the WFZP-A1b allele could cause it to lose its function at the protein level (Figure 2E). No association was found between the WFZP-A1 marker and the SSS trait in this study.

On average, over the six experimental locations/years, the CH7034 WFZP-A1b allele (22.6 ± 0.1, n = 593) set 1.7 more spikelet nodes per spike than the SY95-71 WFZP-A1a allele (20.9 ± 0.1, n = 450) at the QSns.sxau-2A locus (p = 8.68E-27) (Figure 2F). Based on the 14-bp indel, a diagnostic marker for WFZP-A1 was developed (Figure 2G).

A QTL analysis revealed that SNS was associated with 65 markers on chromosome 7A (QSns.sxau-7A). The 1,730 SNP markers including 844 from 17K SNP chips and 886 from 35K DArT SNP chips were assembled into chromosome 7A, but only 43 markers on chromosome arm 7AL were selected as a linkage group to indicate the genetic distance and physical location of the QSns.sxau-7A locus (Figure 3A). The LOD values at the peak position of QSns.sxau-7A ranged from 3.20 to 11.00 in five experimental locations/years, and this QTL alone explained 9.17–28.18% of the total phenotypic variation.

The WAPO-A1 gene (TraesCS7A02G481600) was considered a candidate gene causing QSns.sxau-7A, as it was close to marker SX-D1102483 that was under the peak of QSns.sxau-7A (Figures 3B,C). Primers that were reported to specifically amplify WAPO-A1 (Kuzay et al., 2019) were used to isolate this gene from two parental lines. No difference was found between the CH7034 and WAPO-A1b alleles in the sequenced 2,439-bp region, including the 938-bp promoter region, 1,457 bp from the start codon to the stop codon, and 42 bp at the 3′ end of the gene. Compared with CH7034, SY95-71 had a 2-bp deletion at positions –891 to –892 bp, which was only the difference between the two parental lines (Figure 3D). Compared with the four reported alleles (WAPO-A1a through -A1d), the WAPO-A1 allele in SY95-71 was novel and was thus named WAPO-A1e (Figure 3D).

The WAPO-A1 transcript levels in 1-cm young spikes were compared between the WAPO-A1b and WAPO-A1e alleles. The CH7034 WAPO-A1b allele showed 2.6 times higher transcript levels than the WAPO-A1e allele (p = 0.0001) (Figure 3E). This result suggests that either the WAPOA1b promoter might increase its transcript levels due to the absence of the 2-bp deletion or the WAPO-A1e promoter reduces its transcript levels due to the presence of the 2-bp deletion.

The CH7034 WAPO-A1b allele (22.7 ± 0.1, n = 554) showed 1.5 more spikelet nodes per spike than the SY95-71 WAPO-A1e allele (21.2 ± 0.1, n = 435) at the QSns.sxau-7A locus (p = 1.54E-19), based on the averaged SNS from the six experimental locations/years (Figure 3F). The results showed an association of the SNS phenotype with the 2-bp deletion genotype. A diagnostic marker was developed to distinguish WAPO-A1e that had the 2-bp deletion from WAPO-A1b that did not have the 2-bp deletion (Figure 3G).

The third SNS QTL was located on chromosome 7D (QSns.sxau-7D), where only SNP markers, including 8 from 17K SNP chips and 24 from 35K DArT SNP chips, were assembled into the whole chromosome, spanning 22 cM in genetic distance (Figure 4A). QSns.sxau-7D explained up to 17.92% of the total phenotypic variation in five of the six experimental locations/years, and the effect was not significant in one experimental location/year. The marker SX-D1124962 closest to the peak of QSns.sxau-7D was associated with the SNS trait, the SY95-71 allele (22.1 ± 0.1, n = 561) versus the CH7034 allele (21.8 ± 0.1, n = 536) (p = 0.091) (Figures 4B–D).

The complete gene of homoeologous WAPO-D1 on chromosome 7D (TraesCS7D02G468700) was isolated and sequenced, and the sequenced region, including 564 bp upstream from the start codon, 1,488 bp from the start codon to the stop codon, and 196 bp at the 3′ end of the gene, showed no difference between the SY95-71 and CH7034 alleles. Both SY95-71 and CH7034 had the same WAPO-D1 sequence as Chinese Spring. WAPO-D1 showed no significant difference in its transcript level between the two alleles (Figure 4E). The SNP at SX-D1124962 may be the closest marker to the gene causing QSns.sxau-7D.

Further analyses were performed on each of the three SNS QTLs as independent genetic factors. Eight genotypes generated from these three genes/QTLs, each with two alleles, were analyzed for their SNS phenotypes.

These eight genotypes were classified for the three QTLs in the order of WFZP-A1, WAPO-A1, and QSns.sxau-7D. The SY95-71 allele was designated A, and the CH7034 allele was designated B. The contribution of each genotype to the phenotype was assessed by comparing the SNS mean, given the lack of epistatic interactions. The genotype expected to generate the highest SNS should be the recombinant (B_B_A), containing all the WFZP-A1b and WAPO-A1b alleles from CH7034 and QSns.sxau-7Da from SY95-71, whereas the genotype expected to generate the lowest SNS should be the recombinant (A_A_B), containing all the WFZP-A1a and WAPO-A1e alleles from SY95-71 and QSns.sxau-7Da from CH7034. The incorporated B_B_A genotype generated 23.4 ± 0.2 (n = 149), which was 3.4 more spikelets than the incorporated A_A_B genotype 20 ± 0.2 (n = 124). The remaining six genotypes showed intermediate SNS scores (Figure 5A).

The recently released genome sequences of 10 wheat cultivars provided haplotype information for the frequencies of favorable alleles, WFZP-A1 and WAPO-A1, from CH7034. Only cultivar Norin61 had the same WFZP-A1b allele with the 14-bp deletion as CH7034. All other cultivars showed 100% identity to SY95-71 without the 14-bp deletion, except Spelt, which had a SNP in the WFZP-A gene region.

CH7034 showed 100% identity in the complete 2,439-bp sequence of the WAPO-A1b allele as three wheat cultivars: Mace, Norin61, and Stanley. None of these ten cultivars showed 100% identity to the complete 2,437 bp sequence of the SY95-71 WAPO1-A1e allele, but the 2-bp deletion in WAPO1-A1e was found to exist in six wheat cultivars: Arina, Jaggar, Julius, Lancer, Landmark, and SY_Mattis. These cultivars were found to have the WAPO1-A1a allele, but the sequence of the WAPO1-A1a allele did not include the 2-bp deletion region identified in a previous study (Figure 3D; Kuzay et al., 2019; Voss-Fels et al., 2019). The last cultivar, Spelt, had no 2-bp deletion, but had a WAPO1-A1 that was different from WAPO1-A1a through WAPO1-A1e.

Correlation analysis examined the correlation of the total SNS with fertile spikelets per spike and sterile spikelets per spike over the six experimental locations/years. The total SNS was positively correlated with fertile spikelets per spike (R² = 0.922, n = 184) (Figure 5B). The results indicated that fertile spikelets were increased through spikelet development in wheat.