Identification and validation of quantitative trait loci for fertile spikelet number per spike and grain number per fertile spikelet in bread wheat (Triticum aestivum L.)

A major and stable QTL for fertile spikelet number per spike and grain number per fertile spikelet identified in a 4.96-Mb interval on chromosome 2A was validated in different genetic backgrounds. Fertile spikelet number per spike (FSN) and grain number per fertile spikelet (GNFS) contribute greatly to wheat yield improvement. To detect quantitative trait loci (QTL) associated with FSN and GNFS, we used a recombinant inbred line population crossed by Zhongkemai 13F10 and Chuanmai 42 in eight environments. Two Genomic regions associated with FSN were detected on chromosomes 2A and 6A using bulked segregant exome sequencing analysis. After the genetic linkage maps were constructed, four QTL QFsn.cib-2A, QFsn.cib-6A, QGnfs.cib-2A and QGnfs.cib-6A were identified in three or more environments. Among them, two major QTL QFsn.cib-2A (LOD = 4.67–9.34, PVE = 6.66–13.05%) and QGnfs.cib-2A (LOD = 5.27–11.68, PVE = 7.95–16.71%) were detected in seven and six environments, respectively. They were co-located in the same region, namely QFsn/Gnfs.cib-2A. The developed linked Kompetitive Allele Specific PCR (KASP) markers further validated this QTL in a different genetic background. QFsn/Gnfs.cib-2A showed pleiotropic effects on grain number per spike (GNS) and spike compactness (SC), and had no effect on grain weight. Since QFsn/Gnfs.cib-2A might be a new locus, it and the developed KASP markers can be used in wheat breeding. According to haplotype analysis, QFsn/Gnfs.cib-2A was identified as a target of artificial selection during wheat improvement. Based on haplotype analysis, sequence differences, spatiotemporal expression patterns, and gene annotation, the potential candidate genes for QFsn/Gnfs.cib-2A were predicted. These results provide valuable information for fine mapping and cloning gene(s) underlying QFsn/Gnfs.cib-2A.


Introduction
Bread wheat (Triticum aestivum L.), one of the most widely adapted food crops in the world, provides more than 20% of the calories consumed by human beings (Gao 2021). By 2050, it is predicted that the global population will increase to nine billion, which requires an increase of at least 70% in total food production to meet mankind's demand (http:// www. fao. org). However, the understanding of genes and gene networks that determine grain yield is insufficient due to the low heritability of this trait. Grain yield is controlled by three major components: thousand-grain weight (TGW), grain number per spike (GNS), and spike number per unit area (Wang et al. 2018;Zhang et al. 2022). GNS can be further divided into fertile spikelet number per spike (FSN) and grain number per fertile spikelet (GNFS), which are the focus of the present study. Communicated by David D Fang. In view of the increased potential of FSN and GNFS to improve wheat yield, it is critical to identify, verify and pyramid more genes controlling FSN and GNFS from wheat germplasm. Up to now, several genes controlling or associating with fertile spikelet number or grain number per spikelet have been cloned. For example, WAPO1, an ortholog of rice gene ABERRANT PANICLE ORGANIZATION 1 (APO1), on chromosome arm 7AL could affect spikelet number by regulating inflorescence meristem (Kuzay et al. 2019;Voss-fels et al. 2019). GNI1, encoding a homeodomain leucine zipper class I (HD-Zip I) transcription factor, regulates floret fertility by inhibiting rachilla growth and development. The reduced-function allele of GNI-A1 increases the grain number per spikelet (Sakuma et al. 2019). TaMOC1, an ortholog of rice MOC1, involves in spikelet development and thus associates with spikelet number (Zhang et al. 2015). TaDEP1, an ortholog of rice DEP1, regulates spike length, grain density and reduces the spikelet number (Huang et al. 2009). Meanwhile, wheat FRIZZY PANICLE (Dobrovolskaya et al. 2015), PHOTOPER-IOD RESPONSES locus (Ppd-1) (Boden et al. 2015), bh t -A1 (Wolde and Schnurbusch 2019) and TEOSINTE BRANCHED 1 (TB1) (Dixon et al. 2018) could increase spikelet number. In addition, regulation of Q expression (Greenwood et al. 2017) and miR172 activity (Debernardi et al. 2017) increase grain number per spikelet.
Similar to other yield-related traits, fertile spikelet number and grain number per fertile spikelet in wheat are complex quantitative traits determined by the interaction of genetic and environmental factors. Thus, studies on FSN and GNFS have preliminary focused on quantitative trait loci (QTL) mapping. To date, a large number of FSN and GNFS QTL have been found on almost all wheat chromosomes using different genetic or natural populations (Li et al. 2007(Li et al. , 2021Cui et al. 2012Cui et al. , 2017Yu et al. 2014;Gao et al. 2015;Luo et al. 2016;Deng et al. 2017;Wang et al. 2018;Koppolu and Schnurbusch 2019;Ma et al. 2019;Farokhzadeh et al. 2020;Mo et al. 2021;Ding et al. 2022;Zhang et al. 2022). However, few major QTL have been detected and verified in different genetic backgrounds and multiple environments, which limits the elucidation and further utilization of them. Thus, it is essential to identify and validate new QTL/genes for FSN and GNFS.
In the present study, QTL for FSN and GNFS were identified using the bulked segregant exome sequencing (BSE-Seq) analysis and genetic mapping. Major QTL was further verified in different genetic backgrounds and the candidate gene(s) were predicted.

Plant materials and field management
Two populations generated using the single-seed descent method were employed in this study. They were: (1) a recombinant inbred line (RIL) population derived from a cross of Zhongkemai 13F10 (ZKM13F10)/Chuanmai 42 (CM42) or 13CM: 316 F 7 lines; (2) an F 2 population derived from a cross of Chuanmai 104/SH352 or CS352: 657 individuals. The first population, 13CM, was used for QTL mapping, and the other population, CS352, was used for validating major QTL detected in the mapping population.
CM42 (Syncd768/SW3243//Chuan6415) is a core cultivar widely used in Winter Wheat Breeding Program in China. It is characterized by desirable yield-related traits including high grain weight and high FSN due to its long spike. ZKM13F10 (ZKM138/PW18) is a stable breeding line selected by our lab characterized by higher grain number per spikelet. Chuanmai 104, derived from the cross of CM42/Chuannong 16 (CN16), inherits major elite traits (including high FSN) from CM42. SH352 is a wheat line with fewer spikelets and was used for validating the major QTL.
Page 3 of 13 69 SPSS 24.0 (IBM SPSS, Armonk, NY, USA) and Microsoft Excel 2020 were used to analyze the phenotypic variation, normal distribution and Pearson's correlation of wheat lines grown in different environments. The best linear unbiased estimation (BLUE) for target traits was calculated by QTL IciMapping v4.2 (http:// www. isbre eding. net/ softw are/). Broad sense heritability (H 2 ) was calculated according to the method described by Smith et al (1998). Student's t-test (P < 0.05) was employed to detect significant differences in a given trait between parents and different RIL groups using SPSS. By calculating Pearson's correlation coefficients, the correlations between FSN, GNFS, and other traits based on the BLUE datasets in the 13CM population were also analyzed. Based on the genotypes of the flanking markers, lines carrying different alleles at the major QTL from the mapping population were measured and compared for the other traits using the Student's t-test (P < 0.05) with SPSS.

BSE-Seq analysis
High-quality genomic DNA from parents and 13CM lines were extracted from tender leaves by cetyltrimethylammonium bromide (CTAB) method and checked by agarose gel electrophoresis. Based on the FSN values obtained in the six environments (the mean of eight random individual plants in each line at E1-E6), lines in each environment were rearranged from small to large. The lines within two tails (about 11%) simultaneously in at least four of the six environments were selected to bulk extreme pools. Two pools (FSN-H and FSN-L) were bulked using equal quantities (1 μg) of DNA from 30 individuals in the tails. Exon capture, sequencing and, analysis of these four libraries (ZMK13F10, CM42, FSN-H, and FSN-L) were performed by Bioacme Biotechnology Co., Ltd (Wuhan, China, http:// www. whbio acme. com).
The raw sequence data were deposited in the Genome Sequence Archive (Chen et al. 2021) in National Genomics Data Center (CNCB-NGDC Members and Partners 2022), China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA008821) that are publicly accessible at https:// ngdc. cncb. ac. cn/ gsa. After removing the reads containing sequencing adapters, low-quality bases, or undetected bases, clean data was used for subsequent analysis . It was aligned to the Chinese Spring (CS) reference genome (IWGSC RefSeq v1.0) using the alignment tool BWA (Li 2013). The single nucleotide polymorphism (SNP) and InDel sites were detected and extracted using BCFtools (Danecek et al. 2016), and then annotated using ANNOVAR, which mainly included different regions of the genome and different types of exon regions (Kai et al. 2010).
Two methods, Euclidean distance (ED) and SNP-index, were used to screen the SNP and InDel sites with significant differences between the progeny mixed pools in this study (Abe et al. 2012;Hill et al. 2013;Li 2011;Takagi et al. 2013).

Development of molecular marker
According to the BSE-Seq result, SNPs between ZKM13F10 and CM42 were screened. In addition, polymorphic SNPs simultaneously detected in the parents and the mixed pools within the associated regions (Chinese Spring IWGSC Ref-Seq v2.1) were converted into KASP markers. These markers were developed according to the previously described method by Ji et al (2022) for QTL confirmation and linkage analysis. To evaluate the effects of the major QTL on corresponding traits, tightly linked markers were used. To further validate the major QTL in a different genetic background, these markers were also used to trace the targeted QTL in the validation population.

QTL mapping
Markers that co-localized with others and had > 20% missing values were discarded. The linkage map of the target region in 13CM population was constructed using JoinMap v4.0, and the LOD score was between 5 and 10. Phenotypic values observed from 13CM population in eight environments and BLUE datasets were used. QTL mapping was conducted by inclusive composite interval mapping (ICIM-ADD) based on Bi-parental Populations (BIP) module with QTL IciMapping v4.2. The step speed of QTL detection was 1.0 cM, and the LOD score threshold was set to 2.5. In addition, the interaction of the QTL × environment (QE) was analyzed using QTL IciMapping v4.2 based on the Multi-Environment Trails (MET) module with pre-justed parameters: LOD = 2.5, PIN = 0.001, and step = 4 cM. QTL repeatedly detected in more than three environments (including the BLUE dataset) were considered as stable. The QTL with more than 10% of explained phenotypic variation were treated as major ones, and those with shared flanking markers or located within 1 cM region were denoted as a single one. QTL were named according to the rules of International Rules of Genetic Nomenclature (http:// wheat. pw. usda. gov/ ggpag es/ wgc/ 98/ Intro. htm), where 'cib' represents 'Chengdu Institute of Biology'.

Haplotype analysis
Furthermore, haplotypes at the crucial region of the major QTL were analyzed using the re-sequencing data of 145 landmark cultivars in China (http:// wheat. cau. edu. cn/ Wheat_ SnpHub_ Portal/). Based on this result, haplotype analysis was conducted using 117 wheat accessions collected by our lab (including 66 landraces and 51 cultivars in China). The 117 wheat accessions were planted in Shifang during the 2014-2015 wheat growing season. Each plot consists of two rows with a length of 1.5 m and a row spacing of 0.25 m, with 20 seeds in each row. Eight individuals from each plot were randomly selected to obtain the phenotypic data of FSN.

Prediction of candidate gene(s) for major QTL
According to the mapping result, the physical location of the flanking markers was obtained by blasting against (E-value of 1e-10) Chinese Spring reference genome (IWGSC Ref-Seq v2.1). The high-confidence genes between these markers were extracted from the JBrowse website (https:// urgi. versa illes. inra. fr/ jbrow seiwg sc). The annotations and functions of these genes were analyzed using UniProt (https:// www. unipr ot. org/). Gene expression data in different tissues was obtained from expVIP (http:// www. wheat-expre ssion. com). The expression patterns of candidate genes in various tissues were normalized using the ZeroToOne method and then presented in the HeatMap drawn by TBtools (Chen et al. 2020). The orthologous analysis of these genes was carried out on the EnsemblPlants website (https:// plants. ensem bl. org/ Triti cum_ aesti vum/ Info/ Index). Meanwhile, SNPs in this region were collected using BSE-Seq data.

Phenotypic analyses
There were significant differences (P < 0.05) in GNFS between ZKM13F10 and CM42 among the three environments and the BLUE dataset (Fig. S1, Table 1). CM42 had more fertile spikelet per spike than that of ZKM13F10. In the 13CM population, the FSN and GNFS showed wide and significant variations. Based on the BLUE datasets, the FSN ranged from 17.17 to 22.09, and the GNFS ranged from 2.44 to 3.53, respectively. The estimated H 2 for FSN and GNFS were 0.84 and 0.86, respectively, indicating that both traits were mainly controlled by genetic factors. In a variety of environments and BLUE datasets, FSN and GNFS were distributed continuously, indicating that they were typical quantitative traits controlled by multiple genes (Fig. 1). The frequency distributions of FSN and GNFS in different environments were indicative of polygenic inheritance. Significant Pearson's correlations ranging from 0.123 (P < 0.05) to 0.803 (P < 0.01), and 0.252 (P < 0.01) to 0.809 (P < 0.01) were measured for FSN and GNFS among multiple environments, respectively (Table S1).

Correlation analysis between FSN, GNFS, and yield-related traits
Using BLUE datasets of FSN, GNFS and, yield-related traits to evaluate their correlation (Fig. 2). FSN was significantly (P < 0.01) and negatively correlated with GL, GW, TGW, and PHT, and significantly and positively correlated with GNS, SC, and SL. A significant and positive correlation (P < 0.01) was detected between GNFS and GNS. Significant and negative correlations were observed between GNFS and GL, GL/GW, TGW, and PTN. Among them, the correlation coefficient between FSN and GNS was 0.52 (P < 0.01), and the correlation coefficient between GNFS and GNS was 0.83 (P < 0.01). In addition, a negative but not significant correlation (r = − 0.047) was observed between FSN and GNFS.

BSE-Seq analyses
BSE-Seq analyses of the four bulked pools were employed to detect genomic regions associated with FSN (Table S2), and then the results were compared with the CS reference genome. A total of 417,028,268 raw reads were generated, and after filtering the clean reads were 412,735,792 ranging from 77,283,512 to 130,971,296 for a single pool, indicating that the amount of data is suitable for the subsequent analysis (Table S3). At least 99.49% of the captured reads could be aligned to the CS reference genome, and the mean sequencing depth ranged from 20.03 × to 50.33 × . The coverage ≥ 5 × , representing the proportion of bases with sequencing depth no less than 5 × in the exon region, varied from 60.38% to 75.41% in each pool, suggesting that the quality of BSE-Seq was high and sufficient in the present study. After calling with BCFtools, a total of 5,969,324 SNPs were identified, with the SNPs ranging from 58,929 to 630,245 in each chromosome. The SNP and InDel density on chromosome 2A were 0.76 and 0.04 per 1 Kb, respectively (Table S4). According to the ED and SNP-index methods (threshold values were 0.1149 and 0.3591, respectively), the genomic regions for FSN were found on chromosomes 2A and 6A, and 2A, 2B, and 6A, respectively (Fig. S2, Table S5). The intersection of associated regions between the two methods was observed on chromosomes 2A and 6A with physical intervals of 651.81-665.84 Mb and 446.73-479.11 Mb, respectively.

Genetic map construction and QTL mapping
To confirm the preliminarily identified genomic regions using BSE-Seq, polymorphic SNPs were converted into KASP markers in expanded region (chr2A: 546.84-748.20 Mb; chr6A: 411.94-491.21 Mb) (Table S6 and Table S12). The phenotypic data evaluated in eight environments and their corresponding BLUE datasets were used for QTL mapping, while BLUE datasets were considered as an additional environment. For chromosome 2A, 15 KASP markers were developed to construct a genetic map spanning 42.8 cM in length. According to this map, a major additive QTL (namely QFsn.cib-2A) was detected in seven environments as well as the combined analysis (BLUE). It explained 6.66-13.05% of the phenotypic variance with LOD values ranging from 4.67 to 9.34. As expected, the favorable allele of QFsn.cib-2A was contributed by CM42. This QTL was located in a 4.6-cM interval between the markers KASP08 and KASP11.
For chromosome 6A, 10 KASP markers were developed to construct a genetic map spanning 7.8 cM in length. Based on this genetic map, an additive QTL (namely QFsn.cib-6A) was detected in six environments as well as the combined analysis (BLUE). It explained 3.84-9.78% of the phenotypic variance with LOD values ranging from 2.64 to 7.04. The favorable allele of QFsn.cib-6A was contributed by ZKM13F10.
Seven QTL, including the four QTL identified in the single-environment analysis, were detected in the QE interaction analysis. Except for QFsn.cib-2A, the PVE(A) of the remaining QTL were significantly larger than that of PVE(AE), indicating these QTL were environment stable (Table S7). No epistatic QTL was found in this study (data not show).

Effects of QFsn/Gnfs.cib-2A on yield-related traits
To detect the correlations between the major QTL QFsn/ Gnfs.cib-2A and yield-related traits, RILs of the 13CM were divided into two groups according to their marker profiles. Based on the BLUE datasets, the comparative analysis between the two groups showed that there were significant differences for GNS and SC, indicating that QFsn/Gnfs.cib-2A could be interrelated with spike morphology. No significant differences in GL, GW, GL/GW, PHT, PTN, SL and TGW were observed between the two groups (Fig. S4).

Validation of QFsn/Gnfs.cib-2A in a different genetic background
The validation population (CS352) was employed to evaluate the effects of QFsn/Gnfs.cib-2A in a different genetic background. CM104, a RIL line derived from CM42/CN16, inherited major elite traits from CM42 including a high fertile spikelet number. KASP marker (KASP09), tightly linked to QFsn/Gnfs.cib-2A, was used to genotype. As expected, polymorphism between CM104 and SH352 was detected with this marker. According to their marker profiles, lines of the CS352 were divided into three groups, lines with alleles to be homozygous for CM42; lines with alleles to be homozygous for non-CM42, and lines possessing heterozygous alleles (Fig. 5a). Significant differences (P < 0.01) on FSN and GNFS were detected between groups with different allele based on Student's t-test. The differences in FSN and GNFS between homozygous CM42 allele and non-CM42 allele ranged from 5.89-7.66% and 4.20-6.49%, respectively ( Fig. 5b and c).

Comparison of QFsn/Gnfs.cib-2A to those reported in previous studies
To identify whether QFsn/Gnfs.cib-2A overlaps with previously detected QTL, we compared their physical intervals on chromosome 2A based on the CS reference genome (Table S8)  and CM42 represent lines with alleles from ZKM13F10 and CM42, respectively; NS, *, **, and ***, not significant, P < 0.05, P < 0.01, and P < 0.001, respectively was detected in three environments and the mean data, indicating it is a stable QTL. However, its phenotypic variation is low (1.56-3.25%), suggesting it is a minor QTL (Hu et al. 2020). Meanwhile, QSpn.kibr-2AS is also a minor QTL (PVE = 5.00%), which was detected in an F 2:3 population (Katkout et al. 2014). Given that QFsn/Gnfs.cib-2A identified in this study was simultaneously responsible for FSN and GNFS, we suggest that QFsn/Gnfs.cib-2A may be a novel QTL.
As GNI-A1 on chromosome 2AL was reported to be associated with floret fertility, we thus integrated its linked marker (Xhuj013) into our linkage map by genotyping the 13CM lines (Sakuma et al. 2019). Remapping GNFS showed that QFsn/Gnfs.cib-2A was not linked to GNI-A1, indicating that GNI-A1 was not the candidate gene of QFsn/Gnfs. cib-2A.

Candidate gene analysis of QFsn/Gnfs.cib-2A
According to the CS reference genome (IWGSC RefSeq v2.1), there were 87 prediction genes within the interval of QFsn/Gnfs. , including 44 high-confidence genes (Fig. S5, Table S9). Analysis of the spatial expression patterns (Borrill et al. 2016) showed that there were 11 genes highly and four genes specifically expressed in spike, indicating that they might be involved in spike development (Fig. S5d). Furthermore, two genes, TraesCS2A03G0972300 and TraesCS2A03G0978700 were likely related to spike development (Table S9) based on gene annotation and homolog gene functions in rice. According to BSE-Seq data, a non-synonymous SNP (Leu to Ser) was identified in the coding region of TraesCS2A03G0972300. Twelve SNPs and two InDels were detected in the promoter region of TraesCS2A03G0978700 between ZKM13F10 and CM42 (Table S10).

Haplotype analysis
Haplotypes at the crucial region of QFsn/Gnfs.cib-2A were analyzed based on the re-sequencing data of 145 historically diverse elite wheat cultivars in China. As shown in Fig. S6, there were mainly six haplotypes within this region. In 117 wheat accessions, five SNPs that can effectively distinguish the six haplotypes were selected for haplotype analysis (Fig.  S7). As expected, six haplotypes, named haplotype-I, -II, -III, -IV -V, and -VI, were detected (Fig. S7, Table S11). According to the association analysis result, significant differences (P < 0.05) in FSN were detected between haplotype-III (contained CM42) and haplotype-I (contained ZKM13F10) in cultivars (Fig. 6). These results indicate that the candidate genes might be located within chr2A region: 661,956,566-663,563,612 (CS IWGSC RefSeq v1.1).

Relationships between FSN, GNFS and yield-related traits and pleiotropic effects of QFsn/Gnfs.cib-2A
In the present study, fertile spikelet number per spike and grain number per fertile spikelet were significantly and Red circle represents lines carrying the same allele as CM42 (Homozygous CM42, n = 228); blue circle represents lines carrying the same allele as ZKM13F10 (Homozygous ZKM13F10, n = 278); green circle represents heterozygous alleles (Heterozygous, n = 151). Effects of QFsn/Gnfs.cib-2A on FSN (b) and GNFS (c) in the CS352 population. ** and ***, P < 0.01 and P < 0.001, respectively positively correlated with GNS (Fig. 2). This finding was consistent with that of Ma et al (2019) and Sakuma et al (2019). FSN and GNFS are the two components of GNS in wheat. More fertile spikelets combined with more fertile florets can significantly increase the GNS (Gao et al. 2019). Significant and positive relationships between FSN and SL and SC were detected in the present study, which was consistent with the result of Li et al (2021) and Ma et al (2019). The increase in SL and SC provides more space for holding spikelet.
GNS and TGW are two essentials but usually counteracting determinants of grain yield (Quarrie et al. 2006;Zhai et al. 2017). As the two components of GNS, FSN and GNFS were significantly and negatively correlated with TGW in the present study (Fig. 2). This phenomenon may attribute to source limitations during grain filling or more grains with lower grain weight (Quintero et al. 2018).
CM42 is an elite variety characterized by high yield in the Yangzi River region of China. In the past 20 years, CM42 has been used as one of the parents to select more than 40 varieties. In this study, we found that QFsn/Gnfs.cib-2A has pleiotropic effects on GNS and SC, which may contribute to the high yield of CM42. Therefore, the fine mapping and cloning of this QTL will help to optimize the spike morphology and improve the grain yield.

QFsn/Gnfs.cib-2A is a target of artificial selection in wheat improvement
During wheat improvement, the traditional landraces were continually replaced by cultivars to adapt to various application scenarios. As a result, the genetic distribution was changed by artificial selection, especially for genes controlling agronomic traits (Hedden 2003). Recently, the whole-genome resequencing of wheat landraces and cultivars provides us with valuable information on this process (Guo et al. 2020a;Zhou et al. 2020;Hao et al. 2020;Chen et al. 2019). It is worth noting that haplotype frequency has changed significantly in some regions during the improvement process. These selected haplotypes associated with elite performance in agricultural traits rapidly raised their frequency in later-developed cultivars, and become dominant haplotypes. For instance, according to pedigree-specific haplotype analysis, Xiaoyan 6 carrying haplotype Hap-6A-1 was strongly favored in artificial selection because of its significant dwarfing and early maturing characteristics. Therefore, all new derivatives from Xiaoyan 6 maintained the haplotype Hap-6A-1 (Hao et al. 2020).
In the present study, for QFsn/Gnfs.cib-2A, only one accession was detected in both haplotype-I (1/66) and haplotype-III (1/66) in landrace. In cultivar, ten and twelve accessions were detected in haplotype-I (10/51) and haplotype-III (12/51), respectively. This result suggested that the two haplotypes (haplotype-I and -III) were enriched during artificial selection. As we know, the main objective of selection is focused on yield, and GNS is under a positive selection. These findings suggest that QFsn/Gnfs.cib-2A is a target of artificial selection during wheat improvement.

Potential candidate gene(s) for QFsn/Gnfs.cib-2A
Within the physical interval of QFsn/Gnfs.cib-2A, there were 44 high-confidence predicted genes in the CS genome (Table S9). Based on the haplotype analysis and spatiotemporal expression pattern, ortholog analysis, functional analysis, and sequence differences, TraesCS2A03G0978700 was predicted as the candidate gene. TraesCS2A03G0978700, an ortholog of rice OsMKKK10 and Arabidopsis MAPK KINASE KINASE4 (MAPKKK4), encodes a protein kinase (Table S9). OsMKKK10 can interact with and phosphorylate OsMKK4, while phosphorylated OsMKK4 in turn phosphorylates OsMAPK6. Furthermore, OsMKKK10, OsMKK4, and OsMAPK6 form an OsMKKK10-OsMKK4-OsMAPK6 signaling pathway and participate in the formation of rice panicle morphology (Guo et al. 2018(Guo et al. , 2020bXu et al. 2018). The mutation of MAPKKK4 in Arabidopsis thaliana participates in the biological process of inflorescence development, and stomatal complex morphogenesis (Li et al. 2019a;Wang et al. 2022;Guo et al. 2022;Liu et al. 2022). Multiple SNPs/InDels were detected in the upstream region of TraesCS2A03G0978700 (Table S10), which may change the expression level and affect the function of the gene (Li et al. 2019b). These results indicated that TraesC-S2A03G0978700 might be the candidate gene for QFsn/