QTL Mapping for Important Agronomic Traits Using a Wheat55K SNP Array-Based Genetic Map in Tetraploid Wheat

Wheat yield is highly correlated with plant height, heading date, spike characteristics, and kernel traits. In this study, we used the wheat55K single nucleotide polymorphism array to genotype a recombinant inbred line population of 165 lines constructed by crossing two tetraploid wheat materials, Icaro and Y4. A genetic linkage map with a total length of 6244.51 cM was constructed, covering 14 chromosomes of tetraploid wheat. QTLs for 12 important agronomic traits, including plant height (PH), heading date (HD), awn color (AC), spike-branching (SB), and related traits of spike and kernel, were mapped in multiple environments, while combined QTL-by-environment interactions and epistatic effects were analyzed for each trait. A total of 52 major or stable QTLs were identified, among which may be some novel loci controlling PH, SB, and kernel length-width ratio (LWR), etc., with LOD values ranging from 2.51 to 54.49, thereby explaining 2.40–66.27% of the phenotypic variation. Based on the ‘China Spring’ and durum wheat reference genome annotations, candidate genes were predicted for four stable QTLs, QPH.nwafu-2B.2 (165.67–166.99 cM), QAC.nwafu-3A.1 (419.89–420.52 cM), QAC.nwafu-4A.1 (424.31–447.4 cM), and QLWR.nwafu-7A.1 (166.66–175.46 cM). Thirty-one QTL clusters and 44 segregation distortion regions were also detected, and 38 and 18 major or stable QTLs were included in these clusters and segregation distortion regions, respectively. These results provide QTLs with breeding application potential in tetraploid wheat that broadens the genetic basis of important agronomic traits such as PH, HD, AC, SB, etc., and benefits wheat breeding.


Introduction
Wheat is one of the most important food crops. With the global population expected to be more than 9 billion by 2050, the demand for wheat is also expected to increase by 60% [1,2]. Improving food yields remains an urgent task to achieve global food and nutrition security.
In recent years, researchers studying crop breeding have focused on the ongoing exploitation of novel genetic resources to increase crop yield potential. During the "Green Revolution" around the 1960s, the discovery and application of several genes related to crop dwarfism and dry matter accumulation increased the harvest index (HI), which not only led to significant increases in grain yield per unit area but also prevented global famine [3][4][5]. For more than half a century, the production of crops has been increasing. However, the widespread use of a small number of parents in modern breeding programs and the large-scale promotion of only a few superior varieties have resulted in the loss of wheat genetic resources, an increasingly narrow genetic basis, and the production growth Seven traits (PH, HD, SL, TKW, LWR, AC, and SB) had relatively higher H 2 between 0.64 and 0.96, indicating that these traits were less influenced by environmental factors, while Area had the lowest H 2 at 0.20. The frequency distribution of HD and the six kernelrelated traits (TKW, KL, KW, LWR, Peri, and Area) as well as the spike-related traits KPS and SL were approximately normally distributed in most experimental environments and the BLUE dataset, and transgressive segregations (the absolute values of both skewness and kurtosis are close to 1) in the direction of both parents were observed in the population, which were typical of quantitative trait characteristics ( Figure S1). For PH, remarkable differences in height between population lines resulted in a bi-modal distribution, which allowed most lines to be classified as the dwarf and tall ( Figure S1), suggesting that the height of the population may be primarily influenced by one major gene [41]. AC and SB could be divided into two opposing phenotypes that exhibit partially similar features to qualitative traits (Figure 1b, Figure S1).

Correlation Analyses among Different Traits
Based on BLUE values, the phenotypic correlations between the 12 agronomic traits are presented in Table 2.   As a result, PH was significantly and negatively correlated (p < 0.01) with HD, while it has nonsignificant correlations with both spike and kernel traits. SL was significantly and positively correlated with other spike traits KPS and AC, and negatively correlated (p < 0.01) with kernel traits KL, LWR, and Peri. KPS was significantly and negatively correlated with the kernel-related traits TKW, KL, KW, Peri, and Area, and significantly positively correlated (p < 0.01) with SB. Except for no significant correlation with Area, LWR had a significant negative correlation with KW and TKW and a significant positive correlation with KL and Peri (p < 0.01). The other kernel traits (TKW, KL, KW, Peri, and Area) were significantly positively correlated with each other (p < 0.01). The correlation coefficient between KL and Peri reached 0.974. The correlation coefficients between Area and KL, KW, and Peri ranged from 0.854 to 0.943.

Genetic Map Construction
All 165 lines of the RIL population were genotyped using the Wheat55K SNP array, and after the rejection of missing values and random removal of redundant markers using the BIN function of IciMappingV4.1, 1356 BIN markers were finally screened from 53,063 markers to construct the genetic map. This genetic map covered all 14 chromosomes of tetraploid wheat, spanning 6244.51 cM in length with an average marker density of 4.61 cM/SNP. The A genome contained 602 SNPs and the B genome had 754 SNPs. The total genetic distance of the A and B genome maps was 3259.14 cM and 2985.37 cM, respectively, and the marker densities were 5.41 cM/SNPs and 3.96 cM/SNPs, respectively (Table S1).

Analysis of Segregation Distortions in the Genetic Map
It is well known that marker segregation distortion (SD) is a common phenomenon in genetic linkage maps construction and QTL mapping [42]. After excluding markers with severe distortion, a chi-square test was performed on the 1356 BIN markers involved in the construction of the genetic map, and if p < 0.05, the marker was considered to have undergone SD (significantly different from the expected 1:1 Mendelian separation ratio). The results show that a total of 837 (61.73%) BIN markers showed the distortion from Mendelian segregation. Among these distorted markers, 345 (57.31%) and 492 (65.25%) BIN markers showed SD in A and B genomes, respectively; the number of markers showing SD for each chromosome ranged from 19 (chr2A) to 90 (chr2B and 5B), and chromosomes 7A and 7B had the minimum (43.33%) and maximum (75%) marker distortion rate, respectively (Table S2).
In the RIL population genetic map, regions consisting of five or more consecutive genetic markers with the same direction of segregation were generally referred to as segregation distortion regions (SDRs). The analysis of the SDR is also an essential metric for evaluating the genetic map. While 44 SDRs were detected on 13 chromosomes except for 2A, chromosome 4A had the largest number of SDRs (6 SDRs) (Table S2). These SDRs comprised 470 SD markers with a total length of 2251.68 cM, accounting for 36.06% of the length of the genetic map, and most of them were biased towards the male parental (Y4) genotypes. No-tably, the four adjacent SDRs (AX-111495990~AX-110057517, AX-110935969~AX-110121036, AX-109346617~AX-110705049, and AX-111565010~AX-108891576) on chromosome 4A deviate in the same direction, and this distribution form was also present on other chromosomes. The region on chromosome 7A from AX-110526555 to AX-110484953 consisted of 18 distorted markers biased towards the male parental genotypes and was the longest SDR with a length of 274.5 cM. The shortest SDR was located on chromosome 6A (AX-108852957~AX-109515848), which was 3.16 cM in length and consisted of six distorted markers (Table S2).

QTL Mapping Analysis
The ICIM-BIP method was used to detect additive QTLs for the 12 target traits, and 152 additive QTLs were eventually detected (Table S3). The LOD values of the QTLs ranged from 2.51 to 54.49 and were distributed on 14 chromosomes, explaining 0.45-66.27% of the phenotypic variation (Table S3). Among the detected additive QTLs, the ones with LOD > 3 and explained phenotypic variation > 10% were defined as major QTLs. Additive QTLs that can be detected in two or more individual environments and BLUE dataset were also considered stable QTLs. There were 3, 3, 6, 4, 4, 7, 3, 5, 5, 3, 6, and 3 major or stable QTLs detected for HD, PH, SL, KPS, TKW, KL, KW, LWR, Peri, Area, AC, and SB, respectively (Table 3, Figure 2). These major or stable QTLs with LOD values ranging from 2.56 to 54.49 were distributed on 12 chromosomes and explained 2.40-66.27% of the phenotypic variation (Table 3).  Figure 2. Chromosome positions of the major or stable additive QTLs for the 12 traits described in Table 1. The underlined QTLs indicate that they are the stable QTLs identified in this study. Different colors are used to distinguish different traits on the same chromosome. Ten PH additive QTLs were mapped to chromosomes 2B (2), 3A (2), 4A (2), 6A (2), 6B (1), and 7A (1) ( Table S3). The LOD values of individual QTLs ranged from 3.04 to 19.53, and individual QTLs explained 3.78% to 23.02% of the PH phenotypic variation (Table S3). QPH.nwafu-6A.1 and QPH.nwafu-6A.2 had LOD values > 3, explained > 10% of the PH phenotypic variation, and were considered the major QTLs. QPH.nwafu-6A.1 and QPH.nwafu-2B.2 were detected in at least 3 years and were considered stable QTLs (Table 3). Icaro contributed to the effect of a decreased PH for the two stable QTLs.
Sixteen additive QTLs for KPS were mapped to chromosomes 1B (2) (Table 3). QKPS.nwafu-2B.3 was mapped in Y21 and BLUE dataset, with an average LOD value of 7.24, and the favorable allele was contributed by Icaro, explaining an average of 8.08% of the phenotypic variation ( Table 3).
The 19 additive QTLs of TKW were mapped to chromosomes 1A (1) (Table S3). A total of seven major QTLs were detected, with LOD values ranging from 4.43 to 10.70, explaining 6.61% to 20.41% of the phenotypic variation (Table 3).

QTL Cluster Analysis
In most crops, some QTLs associated with agronomic traits are distributed in clusters on chromosomes, manifesting as QTLs that control different traits located within the same marker interval or adjacent regions [11,43,44]. This distribution pattern was also observed in the tetraploid wheat RIL population via analysis of QTL mapping results in this study. A total of 31 QTL clusters were mapped to 12 chromosomes except for 1A and 5A (Table S4). These QTL clusters contained 97 QTLs, of which 38 were either major or stable QTLs, accounting for about 73% of the number of major and stable QTLs detected using the ICIM-BIP method. Twenty QTL clusters contained major or stable QTLs. The QTL cluster at 260-268 cM on chromosome 2B (C21) contained three QTLs (QTKW.nwafu-2B.2, QKPS.nwafu-2B.4, and QArea.nwafu-2B.3) associated with TKW, KPS, and Area, and three stable QTL QSB.nwafu-2B.1, QKPS.nwafu-2B.3 and QSB.nwafu-2B.2 were also located in this cluster (Table S4). Two stable QTL QSL.nwafu-6A and QPH.nwafu-6A.1 and one major QTL QHD.nwafu-6A.1 were genetically located close to each other and were also detected to be in the same cluster (C8) ( Table S4). The QTL clusters at 260-268 cM on chromosome 2B (C21) and 118-143 cM on chromosome 5B (C26) both contained the largest number of six QTLs controlling four traits including SB, KPS, TKW, and Area and five traits including Area, SB, LWR, TKW, and KPS, respectively (Table S4).
According to the findings, various yield traits such as KPS, SL, SB, and kernel-related traits were significantly correlated with each other (Table 2). Moreover, the positions of QTL for these traits tended to be close to each other (Table S3). Most QTL clusters (C2, C3, C6, C7, C10-17, C19-28, C30, and C31) also tended to contain multiple QTLs for yield traits (Table S4), which may be caused by the presence of pleiotropism at some loci, where a single gene directly or indirectly affects multiple traits in the spike and/or the kernel. The significant negative correlation between PH and HD in the correlation analysis based on the BLUE dataset (p < 0.01, Table 2) also led to the speculation that the genes or QTLs controlling PH and HD might be tightly linked in position or one gene affecting both PH and HD traits, thus showing a close genetic correlation; this speculation was further verified by the tightly linked distribution of QTLs for PH and HD on C5, C8, and C9 (Table  S4). Moreover, we have found the same phenomenon in many previous reports [45][46][47].

Combined QTL-Environment Interaction Analysis
The MET-Add method was used for cQTL analysis, and a total of 184 additive cQTLs for PH (23) (Table S5). Most of these cQTLs were minor QTLs. Moreover, 41 of the 52 major or stable QTLs detected using the ICIM-BIP method were similarly detected using the MET-Add method (Tables 3 and 4), showing that they might be stable in expression and less susceptible to environmental influences.
Among the cQTLs for PH, PVE (A) (phenotypic variation explained by additive and dominant effects) was 2.61-51.32%, while PVE (A by E) (additive and dominance by environment effects for corresponding QTLs) was 0.28-2.25% (Table S5), with a significant decrease in PVE (A by E) compared to PVE (A), demonstrating that PH was primarily influenced by genetic factors. This also coincided with the estimated Hˆ2 of 0.96 for PH based on the BLUE dataset (Table 1). In contrast, except for AC and SB, the other nine traits were relatively susceptible to environmental factors.

Epistatic Analysis
A total of 539 digenic epistatic QTLs (eQTLs) were detected on 14 chromosomes, explaining 0.31% to 7.31% of phenotypic variation (Figure 3, Table S6). For PH, HD, SL, KPS, TKW, KL, KW, LWR, Peri, Area, AC, and SB, 96, 21,3,41,34,41,39,60,38,64,98, and 4 digenic epistatic QTLs were detected, respectively (Table S6). Of these, 51.8% of eQTLs had negative epistatic effects (Table S6), with negative values indicating that the recombinant genotype outperformed the parental genotype in terms of epistatic impact. The epistasis of these six traits (PH, SL, KL, KW, LWR, and AC) were characterized by additive interactions between additive QTLs and random loci and between multiple random loci. Some of the major or stable QTLs detected using the ICIM-BIP method also showed epistatic effects (Table S6). These results showed that these traits were influenced by both epistatic and additive genetic effects. However, for HD, KPS, TKW, Peri, Area, and SB, numerous eQTLs were not detected in the BIP and MET-Add analyses, suggesting  The epistasis of these six traits (PH, SL, KL, KW, LWR, and AC) were characterized by additive interactions between additive QTLs and random loci and between multiple random loci. Some of the major or stable QTLs detected using the ICIM-BIP method also showed epistatic effects (Table S6). These results showed that these traits were influenced by both epistatic and additive genetic effects. However, for HD, KPS, TKW, Peri, Area, and SB, numerous eQTLs were not detected in the BIP and MET-Add analyses, suggesting that these random loci may not have a direct effect on the phenotype but can indirectly influence it through interactions between loci.

Effects of Stable Quantitative Trait Loci on Related Traits
There were 184 cQTLs for the 12 traits identified by combined QTL and environment interactive effect analysis (Table S5). Four of them were identical to the four stable QTLs (QPH.nwafu-2B.2, QAC.nwafu-3A.1, QAC.nwafu-4A.1, and QLWR.nwafu-7A.1), which were identified via single-environment analysis on chromosomes 2B, 3A, 4A, and 7A, and had relatively high PVE and narrow genetic localization intervals (Figure 2, Tables 3 and 4). To assess the effects of different haplotypes on the traits of interest, lines with the corresponding homozygous alleles in the population were obtained based on the genotyping results of the two flanking markers of these QTLs. Depending on whether the homozygous alleles of the QTL locus were from Icaro or Y4, the RILs were assigned into two groups accordingly: Icaro type and Y4 type. In addition, the corresponding phenotypes of the two groups were analyzed.
For PH, one stable QTL QPH.nwafu-2B.2 and one stable major QTL QPH.nwafu-6A.1 were identified. Moreover, significant differences in PH were observed between the two groups of lines corresponding to each locus (p < 0.01, Figure 4). The PH of Icaro type lines was significantly lower than those Y4 type lines, and the effect of the QPH.nwafu-6A.1 was greater than that of QPH.nwafu-2B.2 in reducing the magnitude of PH. To distinguish more carefully the effects of the two loci on PH as well as HD and TKW, the lines were further divided into four groups based on the different genotype combinations of the two loci ( Figure 5). Compared with lines carrying Y4 alleles at both QPH.nwafu-2B.2 and QPH.nwafu-6A.1, each locus carrying Icaro alleles was able to reduce PH by 3.97% to 9.76% and 18.98% to 24.92%, respectively, and clustering of the two further contributed to this effect, reaching 26.53% to 28.10% (Figure 5a). For HD, the two loci show different effects. For the QPH.nwafu-6A.1, the HD of lines carrying Icaro alleles were significantly delayed in all environments and BLUE dataset (Figure 5b). However, for the QPH.nwafu-2B.2, the HD of lines carrying Icaro alleles was slightly earlier; in addition, neither locus had a significant positive or negative effect on TKW (Figure 5c).
interactive effect analysis (Table S5). Four of them were identical to the four stable QTLs (QPH.nwafu-2B.2, QAC.nwafu-3A.1, QAC.nwafu-4A.1, and QLWR.nwafu-7A.1), which were identified via single-environment analysis on chromosomes 2B, 3A, 4A, and 7A, and had relatively high PVE and narrow genetic localization intervals (Figure 2, Tables 3 and 4). To assess the effects of different haplotypes on the traits of interest, lines with the corresponding homozygous alleles in the population were obtained based on the genotyping results of the two flanking markers of these QTLs. Depending on whether the homozygous alleles of the QTL locus were from Icaro or Y4, the RILs were assigned into two groups accordingly: Icaro type and Y4 type. In addition, the corresponding phenotypes of the two groups were analyzed.
For PH, one stable QTL QPH.nwafu-2B.2 and one stable major QTL QPH.nwafu-6A.1 were identified. Moreover, significant differences in PH were observed between the two groups of lines corresponding to each locus (p < 0.01, Figure 4). The PH of Icaro type lines was significantly lower than those Y4 type lines, and the effect of the QPH.nwafu-6A.1 was greater than that of QPH.nwafu-2B.2 in reducing the magnitude of PH. To distinguish more carefully the effects of the two loci on PH as well as HD and TKW, the lines were further divided into four groups based on the different genotype combinations of the two loci ( Figure 5). Compared with lines carrying Y4 alleles at both QPH.nwafu-2B.2 and QPH.nwafu-6A.1, each locus carrying Icaro alleles was able to reduce PH by 3.97% to 9.76% and 18.98% to 24.92%, respectively, and clustering of the two further contributed to this effect, reaching 26.53% to 28.10% (Figure 5a). For HD, the two loci show different effects. For the QPH.nwafu-6A.1, the HD of lines carrying Icaro alleles were significantly delayed in all environments and BLUE dataset (Figure 5b). However, for the QPH.nwafu-2B.2, the HD of lines carrying Icaro alleles was slightly earlier; in addition, neither locus had a significant positive or negative effect on TKW (Figure 5c).  As a stable locus, QLWR.nwafu-7A.1 was identified under the Y21, Y22, and BLUE datasets (Table 3). For the QLWR.nwafu-7A.1, the two groups of RILs also showed significant differences in LWR (p < 0.01, Figure S2). The LWR of Y4 type lines had a significantly larger LWR compared with those Icaro type lines; however, the TKW was not significantly affected. For the QAC.nwafu-3A.1 and QAC.nwafu-4A.1, lines carrying from Icaro alleles and thousand-kernel weight (TKW) (c). + and − represent lines with the alleles from the Icaro and Y4 of the target locus, respectively; * and ** Significance at p < 0.05 and p < 0.01, respectively; ns: Not significant at p > 0.05. As a stable locus, QLWR.nwafu-7A.1 was identified under the Y21, Y22, and BLUE datasets (Table 3). For the QLWR.nwafu-7A.1, the two groups of RILs also showed significant differences in LWR (p < 0.01, Figure S2). The LWR of Y4 type lines had a significantly larger LWR compared with those Icaro type lines; however, the TKW was not significantly affected. For the QAC.nwafu-3A.1 and QAC.nwafu-4A.1, lines carrying from Icaro alleles had a higher penetrance of black AC, and clustering of the two loci would further enhance the penetrance ( Figure S4).
Spike-branching phenotype on tetraploid wheat is an important additional spikelet trait associated with a significant increase in KPS; thus, we also performed a preliminary effect analysis on the two identified SB stable loci. For the QSB.nwafu-2B.1 and QSB.nwafu-2B.2, the lines were also further classified into four groups based on the different genotype combinations of the two loci. Compared with lines carrying Y4 alleles at both QSB.nwafu-2B.1 and QSB.nwafu-2B.2, each locus carrying Icaro alleles was able to increase KPS, and clustering of the two could further contribute to this effect ( Figure S3); however, the increase did not all reach statistical significance in some of the environments. In addition, TKW was not significantly influenced by different genotypes (Figure S3), indicating that these two QTLs could increase KPS without affecting kernel weight.

Analysis of the Relation of Mapped QTLs to Those Found in Previous Studies
We identified loci for some of the major agronomic traits, both those already reported and confirmed, and novel loci with promising applications. These loci can subsequently develop specific molecular markers capable of efficiently targeting genes for marker-assisted selection in wheat breeding. It should be noted that most studies to date have focused on common wheat, whereas tetraploid wheat may exhibit different compositions for these complex traits. Thus, we focused mainly on the previous reports on tetraploid wheat compared with the QTLs in Table 3 (Table 5).
Plant height is a quantitative trait controlled by multiple genes/QTLs. In the present study, we identified two stable QTLs, QPH.nwafu-6A.1 and QPH.nwafu-2B.2, that simultaneously affect HD and PH ( Figure 5). QPH.nwafu-6A.1 and QHD.nwafu-6A.1 were mapped to the same interval (197.82-201.86 cM on 6A), which overlapped with QSL.nwafu-6A (196.86-201.86 cM on 6A), indicating that these three traits might be contributed by a single gene/locus on 6A chromosome. However, the semi-dwarf durum cultivar Icaro has previously been reported by Tang [48] to contain Rht18 (412-445 Mb) within the physical interval corresponding to the QPH.nwafu-6A.1. Furthermore, Icaro type lines at the QPH.nwafu-6A.1 locus were approximately 27% shorter in height compared to the tall progeny ( Figure 5), which is essentially the same as the result reported by Tang et al. [49] that Rht18 reduced PH by 26%. These results all seem to indicate that QPH.nwafu-6A.1 is the Rht18 that has been mapped and has a significant impact on the PH phenotype. Rht18 was localized in the region from 412 to 445 Mb on chromosome 6A, which contains the GA2oxA9 encoding a GA 2-oxidase, and overexpression of GA2oxA9 reduced bioactive GA content and PH [50,51]. Both QPH.nwafu-2B.2 and another semi-

Analysis of the Relation of Mapped QTLs to Those Found in Previous Studies
We identified loci for some of the major agronomic traits, both those already reported and confirmed, and novel loci with promising applications. These loci can subsequently develop specific molecular markers capable of efficiently targeting genes for marker-assisted selection in wheat breeding. It should be noted that most studies to date have focused on common wheat, whereas tetraploid wheat may exhibit different compositions for these complex traits. Thus, we focused mainly on the previous reports on tetraploid wheat compared with the QTLs in Table 3 (Table 5).
Plant height is a quantitative trait controlled by multiple genes/QTLs. In the present study, we identified two stable QTLs, QPH.nwafu-6A.1 and QPH.nwafu-2B.2, that simultaneously affect HD and PH ( Figure 5). QPH.nwafu-6A.1 and QHD.nwafu-6A.1 were mapped to the same interval (197.82-201.86 cM on 6A), which overlapped with QSL.nwafu-6A (196.86-201.86 cM on 6A), indicating that these three traits might be contributed by a single gene/locus on 6A chromosome. However, the semi-dwarf durum cultivar Icaro has previously been reported by Tang [48] to contain Rht18 (412-445 Mb) within the physical interval corresponding to the QPH.nwafu-6A.1. Furthermore, Icaro type lines at the QPH.nwafu-6A.1 locus were approximately 27% shorter in height compared to the tall progeny ( Figure 5), which is essentially the same as the result reported by Tang et al. [49] that Rht18 reduced PH by 26%. These results all seem to indicate that QPH.nwafu-6A.1 is the Rht18 that has been mapped and has a significant impact on the PH phenotype. Rht18 was localized in the region from 412 to 445 Mb on chromosome 6A, which contains the GA2oxA9 encoding a GA 2-oxidase, and overexpression of GA2oxA9 reduced bioactive GA content and PH [50,51]. Both QPH.nwafu-2B.2 and another semi-dwarfing gene Rht4, which is tightly linked to the microsatellite marker WMC317, are located on 2BL. However, BLAST analyses revealed that the positions of the two loci are not identical (Figure 6a). In addition, Icaro type alleles at the QPH.nwafu-2B.2 locus had a tendency to accelerate HD and reduce PH, which is similar to the function of the photoperiod response locus Ppd-D1a [52]. However, it also differs from the position of the homologous locus of Ppd-D1a on chromosome 2B (Ppd-B1) [53]. Therefore, QPH.nwafu-2B.2 may be a new locus for PH. Both the shortened heading date and the reduced plant height facilitate the transport and distribution of assimilates produced by the plant to the seeds, spikes and other "sink", and the increase in assimilates admitted to the "sink" can help improve crop yields under a certain number of "sources"; therefore, QPH.nwafu-2B.2 is a highly promising locus for application, and further research is necessary.
The stable major QTL QSL.nwafu-3A detected under the Y21, Y22, and BLUE dataset shared the same position as the locus reported by Giraldo et al. [57]. QKPS.nwafu-4B.3 was consistent with the locus reported by Kidane et al. [58]. QKPS.nwafu-2B.3 and QKPS.nwafu-5B.3 had extensive overlap with the two loci reported by Mangini et al. [39], respectively. The TKW-related physical intervals reported by Peng et al. [54], Graziani et al. [59], and Roncallo et al. [45] contained QTKW.nwafu-1B.1. The physical interval of QTKW.nwafu-2B.3 contained multiple previously reported TKW loci [39,60,61]. The physical location of QTKW.nwafu-5B.2 had a large overlap with the locus of TKW reported by Roncallo et al. [45] and Peng et al. [54]. They may be influenced or controlled by the same gene. Ma et al. [62] identified TaGS5 on 3B, but this gene's location was different from that of QTKW.nwafu-3B.2, suggesting that there may be novel genes associated with TKW in the QTL. In addition, compared with the physical map of QTLs for other yield-related traits such as KL, KW, LWR, and kernel size (KS) [37,63], the physical regions of QTLs in this study did not overlap with these loci. Spike-branching (SB) phenotype can alter the capacity of the "sink" by changing spike traits (e.g., spikelet number per spike and KPS), which has the potential to increase grain yield. In this study, two stable loci of SB (QSB.nwafu-2B.1 and QSB.nwafu-2B.2) located on chromosome 2B were identified, explaining 10.00-11.20% and 8.97-23.16% of the phenotypic variation, respectively. Both loci had favorable alleles for SB traits contributed by Y4 (Table 3). Previous studies have shown that an amino acid substitution in the AP2/ERF domain of TtBH-A1, representing a mutant allele at the WFZP-A locus, causes spike branching in tetraploid T. turgidum wheat [64]. The expression of the trait is influenced by various environmental factors, in addition to the fact that wfzp-A/TtBH-A1 is the main genetic factor determining the spike branching of turgidum type in tetraploid wheat, and other genetic factors also play a role in regulating turgidum-type SB [65][66][67]. Wolde et al. [26] reported the new modifier QTL QSS.ipk-2BS for spike-branching on 2BS, which was highly associated with coding sequence variation of the homologous allele of TtBH-B1 (bh t -B1). However, BLAST analysis showed that QSB.nwafu-2B.1 and QSB.nwafu-2B.2 were located on the long arm of chromosome 2B, which was far from QSS.ipk-2BS. This suggests that these two QTLs are potentially new SB loci.

Segregation Distortion Regions in the Genetic Map
It has generally been assumed that many factors from gametogenesis to zygote formation, and then to the post-zygote stage can lead to segregation distortion [68]. Marker segregation distortion is a common biological phenomenon [69]. However, these markers are either simply discarded or incorporated unnoticed in most studies and applications [70]. Hackett and Broadfoot [71] concluded from simulation analysis that the presence of segregation distortion had little effect on the accuracy of marker order or genetic map length. Zuo et al. [70], based on a study of 519 soybean RIL populations from orthogonal and reciprocal crosses between LSZZH and NN493-1, found that the inclusion of distorted markers had a positive effect on increasing genome coverage and improving the concordance of linkage maps with genome. Therefore, based on the results of previous studies, we believed that when the exact chromosome to which an abnormal marker belongs is known and anchored during the creation of a genetic map, the appropriate inclusion of slightly distorted markers not only reduces the effect of marker positional disorder but also increases the marker density of the genetic map.
Forty-four SDRs were detected with a length of 2251.68 cM, accounting for 36.06% of the overall length of the genetic map, distributed over all chromosomes except 2A (Table S2). Eighteen major or stable QTLs were located in the SDRs, accounting for 34.62% of the total number of major and stable QTLs detected. A few major or stable QTLs for all traits except AC, SB, and KW were also distributed in SDRs (Table 3, Table S2). Of the 837 distorted markers, 555 (66.3%) and 282 (33.7%) were biased toward the male and female parental genotypes, respectively. In addition, 31 of the 44 SDRs were also biased toward the male parental genotype (Table S2). The distribution and biased direction of distorted markers were more inclined to the B genome and male parental genotypes, respectively. The occurrence of segregation distortion in hybrid progenies makes the allele frequencies from both parents change in the progeny population. The changing allele frequencies indicate the evolution of the population, which has important evolutionary significance [72,73]. Thus, changes in allele frequencies due to selection and other factors, the impact of such changes on estimating genetic linkage distances between loci, and how to efficiently select parental material and loci based on the direction of distortion bias should also be critical issues to consider in future breeding practices. Fifteen genes were annotated in the chromosome 2B region between 663.7 and 664.2 MB, some of which may be associated with PH (Figure 6a, Table S7); for example, TraesCS2B01G468100 encodes the WD40 protein, which is involved in several cellular processes including signal transduction, transcriptional modulation, and histone modification, and was considered a key regulator of plant developmental processes [74]. Moreover, five genes, such as TraesCS2B01G467800, TraesCS2B01G467700, and its homologous gene TRITD2Bv1G217740, all encode MYB transcription factors that play significant roles in plant developmental processes and stress responses. AtMYB68, a root growth-specific regulator in Arabidopsis thaliana, influences overall plant development in unfavorable situations (such as elevated temperatures) [75]; during the seedling stage, AtMYB38 and AtMYB18/LAF1 regulate the hypocotyl response to blue and far-red light, respectively [76,77].
A total of 183 genes were annotated in the 652.3-661.2 MB region of 7A, several of which may be associated with LWR traits (Figure 6d, Table S7). For example, TraesCS7A01G457500 and the homologous TRITD7Av1G244030 encode an F-box protein that is involved in many pathways of plant nutrition and reproduction [37,81]. TraesCS7A01G461700 and its homologous TRITD7Av1G245930 encode an Auxin response factor. ARF family TFs play a key role in transmitting auxin signals to alter plant growth and development [82]. In rice, increased auxin content induces auxin response factor-mediated activation of NO 3 − transporter and N-metabolism genes, resulting in improved N-use efficiency and grain yield [83].

Plant Materials and Field Trials
The mapping population was from the 165 F 6 generation RILs derived from a cross between the durum wheat variety named Icaro and a variety of T. turgidum ssp. turgidum known as Y4 that was selected and bred in our laboratory. Icaro has an earlier maturity, black awn, and disease resistance, and Y4 has larger and fuller seeds, higher PH, longer spike with spike-branching characteristics, white awn, and the stalks are so tough that they do not fall over easily.
RILs and their parents were grown in one row for each line at Yangling (108 • 07 E, 34 • 30 N) during the 2018-2019 (Y19), 2019-2020 (Y20), 2020-2021 (Y21), and 2021-2022 (Y22) growing seasons. Before sowing, tilling and basal fertilizer were applied to the experimental field; conventional weeding and winter irrigation were carried out during plant growth. Three replications were used in the randomized block design of the experiment. Each line was planted with 15 seeds evenly spaced along a 1.5 m row, with 0.25 m between each row.

Phenotypic Evaluation and Statistical Analysis
The following 12 important agronomic traits were investigated at different growth stages. Among them, HD was investigated at the early heading stage, and PH, SB, and AC were investigated at the maturing stage; the rest were measured after harvest.
HD was recorded when half of the plant spikes of each line emerged from the flag leaf sheath and converted to a representation of the number of days difference between it and the earliest HD in the population. Five representative plants of each line were randomly selected for manual measurements of PH, SL, and KPS. PH was measured from the ground to the top of the main spike (excluding awns). For SL, the length from the base of the rachis to the tip of the terminal spikelet was measured, excluding the awns. KPS was obtained from five spikes that had previously been measured for SL, individually threshed per spike and counted for kernel number. The other kernel-related traits, including TKW, KL, KW, LWR, Peri, and Area, were measured in the laboratory using an SC-G automatic grain analyzer (Wanshen Testing Technology Co., Ltd., Hangzhou, China). The AC and SB phenotypes for each line were typed according to a rank of 0 and 2, where 0 was given to a spike-branching or black-awn phenotype and 2 to a non-branching or white-awn phenotype.
The PH and AC data were collected from Y19 to Y22, and the HD, SL, KPS, TKW, KL, KW, LWR, Peri, Area, and SB data were collected from Y20 to Y22. For each trait, best linear unbiased estimation (BLUE) was calculated across environments using R package "lme4" and used for the estimation of genotype fixed effects and the calculation of broad-sense heritability (H 2 ). Based on the BLUE dataset, descriptive statistics, analysis of variance (ANOVA), and Pearson correlation analysis were performed for the target traits using SPSS 23.0 (https://www.spss.com, accessed on 2 September 2022). Population phenotype frequency distributions and the boxplot were drawn using Origin 2021 (https://www.originlab.com, accessed on 20 October 2022).

Genotyping and Genetic Map Construction
The genomic DNA (gDNA) of each RIL and parent lines was extracted from fresh leaves following the modified cetyl trimethylammonium bromide (CTAB) extraction method [84]. After DNA integrity, concentration and purity were confirmed, and the gDNA was hybridized to the Wheat55K SNP array containing 53,063 SNP markers. Genotyping was performed by Beijing CapitalBio Technology Co., Ltd. (https://www.capitalbiotech. com/, accessed on 7 September 2022).
The polymorphic SNPs with heterozygous genotypes between parents or allele frequencies <0.3 or >0.7 among the RILs were removed. The IciMappingV4.1 (https:// isbreeding.caas.cn/, accessed on 10 September 2022) BIN module was used to remove redundant SNPs and SNPs with missing rate >10%, distortion p value >0.01, and obtained BIN markers. The BIN markers were then used for map construction. Briefly, the BIN markers were firstly anchored to the corresponding chromosomes by combining data on the corresponding physical positions of the markers in the Chinese Spring (CS) reference genome (IWGSC RefSeq v1.0, https://urgi.versailles.inra.fr/blast_iwgsc/blast.php, accessed on 2 October 2022). Ordering and Rippling of markers were then performed using the nnTwoOpt algorithm and the SARF criterion, respectively. Recombination frequencies were converted into genetic distance (cM, centimorgan) using the Kosambi mapping function. The genetic linkage map was visualized using Mapchart V2.3.

Quantitative Trait Loci Mapping
QTL analysis was undertaken for the 12 agronomic traits using the inclusive composite interval mapping (ICIM) method in the biparental population BIP module of IciMap-pingV4.1 with the following parameters: Step = 1 cM, PIN = 0.0001, and the LOD threshold for a significant QTL was set to 2.5 [85,86]. For each trait, QTLs less than 1 cM apart or sharing a common flanking marker were regarded as a single QTL. The multi-environment trials (MET-Add) module of IcimappingV4.1 was used to identify cQTL (combined quantitative trait loci) and additive-by-environment (A by E) interactive effects [87]. Epistatic analysis was performed using the IciMappingV4.1 MET-EPI module, and eQTL (epistatic quantitative trait loci) detection parameters were set as follows: LOD = 5, step = 5 cM, PIN = 0.0001. QTLs were named based on the International Rules of Genetic Nomenclature, where "nwafu" represents Northwest Agriculture and Forestry University.
Flanking molecular marker sequences of the major or stable QTLs were compared with the CS reference genome (IWGSC RefSeq v1.0) and the durum wheat genome reference sequence (Durum Wheat Svevo Refseq Rel. 1.0) on the Triticeae Multiomics Center website (http://202.194.139.32/, accessed on 10 October 2022) to obtain the corresponding physical interval locations. Moreover, Uniprot was used (https://www.uniprot.org/, accessed on 10 October 2022) for further annotation and functional analysis of the candidate genes in the QTL intervals.

Data Availability Statement:
The data presented in this study are available in this article and the Supplementary Materials.

Conflicts of Interest:
The authors declare no conflict of interest.