Spikelet number of YC RIL population in different environments
The SNS and spike types of the YP and other parents were significantly different. Among the different growth stages, the spikelet number of CS ranged from 22 to 28, and the spikelet number of YP ranged from 45 to 85 (Figure 1). The mean number of SNS in YP was 70.5, while the mean number of SNS in CS was 26.8, the mean number of SNS in YP was significantly higher than that of CS.
The SNS trait in YC population was observed in different environments. The mean minimum number of SNS was 15.7, the maximum number was 94.8, the broad-sense heritability was 87.77% (Table 1). The correlation coefficients between the different environments were all significant, ranging from 0.655 (P<0.01) to 0.961 (P<0.01) (Table 2).
BSE-seq analysis
The exome capture sequencing of bulked segregation (BSE-Seq) analysis was performed to identify the genomic regions, and the results were compared with the CS reference genome v2.1. We use the ΔSNP-index algorithm to calculate the allele segregation of the SNPs and InDels between two extreme DNA pools. ΔSNP-index algorithm showed that abundantly candidate SNPs/InDels enriched on chromosome 2A、2D(Chr2A、Chr2D). Select ΔSNP-index greater than 95% as the threshold for screening (Shang et al. 2021). One 16.3 Mb region of Chr2A (63.8Mb-80.1Mb) and one 9.1 Mb region of Chr2D (68.1Mb-77.2Mb) exhibiting significant linkage disequilibrium were identified as the candidate region for SS trait (Figure 2). These SNPs/Indels in the candidate regions would help for further mapping the SS QTL.
Linkage map construction and QTLs identification
To confirm the preliminarily identified genomic regions responsible for SS trait, SNPs and InDels in the target regions were converted into KASP and InDel markers, and then construct the genetic map. In total, 13 KASP markers and one InDel marker were used for the construction of the genetic maps. Six KASP markers developed for Chr2A generated a linkage map spanning 18.1 cM, while seven KASP markers and one InDel marker developed for Chr2D generated a linkage map spanning 17.9 cM (Figure 3, Table S1). The phenotypic data of SS trait evaluated in the three environments were used for QTL mapping. Two stable QTLs named QSS.sicau-2A and QSS.sicau-2D were detected in three environments and BLUE. QSS.sicau-2A located between A2 and A17, explained 6.3%~15.6% of the phenotypic variance with the LOD values ranging from 3.05 to 7.94. QSS.sicau-2D located between D3 and D11, explained 23.4%~34.8% of the phenotypic variance with the LOD values ranging from 12.25 to 19.96 (Table 3). The favorable alleles of the two loci were all contributed by YP.
Effects of supernumber spikelet QTL in different genetic backgrounds
Three KASP markers (A4, D6, 2DYC) and one InDel marker were used to assess the effects of QTL on SNS in validation populations (A4 and InDel were used in YJ population, 2DYC and A4 were used in G2 population, D6 and A4 were used in NY2 population). In YJ population, the average number of spikelet in homozygous “AA, DD” genotypes was 40, that was significantly increased by 75.67% to 96.21% than that in homozygous “aa, dd”, “AA, dd” and “aa, DD” genotypes. Similarly, in G2 and NY2 populations, the average spikelet number in homozygous “AA, DD” genotypes were 36.25 and 36.47, respectively, that significantly increased by 54.1% ~ 61.8% and 73.5% ~ 79.21% than that in other three homozygous genotypes (Figure 4).
Candidate genes prediction
According to the CS genome (IWGSC Refseq v2.1), QSS.sicau-2A was located between 67.0 Mb and 78.3 Mb on Chr2A, with 145 high-confidence genes, including WFZP-A (TraesCS2A03G0239400). Sequence analysis of the WFZP-A, a 4 bp deletion of WFZP-A was identified resulting in a frame-shift in YP (Figure S1, Table S2), and the variation type is consistent with Zang734. These results suggest that WFZP-A is likely responsible for QSS.sicau-2A. QSS.sicau-2D was mapped between 70.4 Mb and 77.3 Mb on Chr2D, including 94 high-confidence genes. Eight spike-specific genes were detected by examining their expression profiles, that might be probably involved in spike growth and development (Table S3). Exome sequencing data of YP reveal one T for A substitution converts the codon (TAT) to a translational stop codon (TAA) in TraesCS2D03G0260700 (Figure S2, Figure 5). Additionally, one known branching gene WZFP-D (TraesCS2D03G0226400LC) is adjacent to QSS.sicau-2D, so we check the gene WZFP-D in YP and CS. However, no sequence variation between YP and CS was identified for WZFP-D, and no significant difference in gene expression (Figure 6, Table S2).
Effects of QSS.sicau-2A and QSS.sicau-2D on yield-related traits
Two QTLs, QSS.sicau-2A and QSS.sicau-2D, were identified in YC population. Their effects on the yield-related traits were analyzed by linking markers. The YC population was divided into four genotypes of “aa, dd” “AA, dd”, “aa, DD” and “AA, DD” by molecular markers. The TGW, GL and GW of homozygous “AA, DD” genotypes were significantly lower than that of other three genotypes, the GNS and GWP of homozygous “AA, DD” genotypes were significantly higher than that of other three genotypes (Figure 7a-7f). Additionally, the length and width of flag leaf were not affected by QSS.sicau-2A and QSS.sicau-2D (Figure 7g-7h).