Global QTL Analysis Identifies Genomic Regions on Chromosomes 4A and 4B Harboring Stable Loci for Yield-Related Traits Across Different Environments in Wheat (Triticum aestivum L.)

Major advances in wheat production are needed to address global food insecurity under future climate conditions, such as high temperatures. The grain yield of bread wheat (Triticum aestivum L.) is a quantitatively inherited complex trait that is strongly influenced by interacting genetic and environmental factors. Here, we conducted global QTL analysis for five yield-related traits, including spike yield, yield components and plant height (PH), in the Nongda3338/Jingdong6 doubled haploid (DH) population using a high-density SNP and SSR-based genetic map. A total of 12 major genomic regions with stable QTL controlling yield-related traits were detected on chromosomes 1B, 2A, 2B, 2D, 3A, 4A, 4B, 4D, 5A, 6A, and 7A across 12 different field trials with timely sown (normal) and late sown (heat stress) conditions. Co-location of yield components revealed significant tradeoffs between thousand grain weight (TGW) and grain number per spike (GNS) on chromosome 4A. Dissection of a “QTL-hotspot” region for grain weight on chromosome 4B was helpful in marker-assisted selection (MAS) breeding. Moreover, this study identified a novel QTL for heat susceptibility index of thousand grain weight (HSITGW) on chromosome 4BL that explains approximately 10% of phenotypic variation. QPh.cau-4B.2, QPh.cau-4D.1 and QPh.cau-2D.3 were coincident with the dwarfing genes Rht1, Rht2, and Rht8, and haplotype analysis revealed their pleiotropic architecture with yield components. Overall, our findings will be useful for elucidating the genetic architecture of yield-related traits and developing new wheat varieties with high and stable yield.


INTRODUCTION
Common wheat (Triticum aestivum L.) is one of the most widely adapted food crops worldwide, providing approximately 30% of global grain production and 20% of the calories consumed by humans (FAO, 2017). The development of high-yield varieties is one of the important targets of modern wheat breeding programs worldwide because of the ever-growing global population and limited land for agricultural expansion (Lobell et al., 2011;Ray et al., 2012). Therefore, the identification, understanding and incorporation of QTL/genes that beneficially influence yield can facilitate the genetic improvement of varieties with high yield.
Grain yield in wheat is a complex quantitative trait that is strongly influenced by interacting genetic and environmental factors and can usually be broken down into three components: spikes per plant (SPP), grain number per spike (GNS), and thousand grain weight (TGW) (Quarrie et al., 2006;Gao et al., 2015). These yield components are sequentially fixed, influencing each other during the growth cycle, and are affected by other traits, such as plant height (PH), crop phenology, and biomass. Also they vary in terms of heritability (Del Moral et al., 2003;Ogbonnaya et al., 2017). SPP reflects wheat tiller capacity, which is one of the key characteristics affecting the yield potential of cereal crops (Naruoka et al., 2011). To date, wheat tillering QTL have been identified on chromosomes 1A, 1B, 1D, 2A, 2B, 2D, 3A, 3B, 4D, 5A, 5D, 6A, 6D, and 7A (Spielmeyer and Richards, 2004;Kuraparthy et al., 2007;Liu G. et al., 2014;Wang et al., 2016). Wheat spike is an important reproductive organ and is positively correlated with grain yield. Previous studies have identified more than 100 QTL for GNS distributed on all 21 chromosomes in wheat through linkage analysis and association analysis; these QTL were primarily located on chromosomes 1A, 1B, 1D, 2A, 2D, 3B, 3D, 4A, 5A, 6A, 7A, and 7D (Gao et al., 2015;Zhai et al., 2016;Cui et al., 2017;Shi et al., 2017;Zhou et al., 2017). Meanwhile, several genes related to GNS have been cloned and characterized using a homology-based approach (Zhang et al., 2011(Zhang et al., , 2015Zheng et al., 2014). Grain weight is another essential yield component that is more stably inherited than final yield and is part of domestication syndrome in cereal crops (Quarrie et al., 2005;Meyer and Purugganan, 2013). Recently, several major QTL for grain size have been cloned and characterized in rice, providing important information about the molecular basis of grain weight in crop plants (Wang et al., 2008;Huang et al., 2013;Zuo and Li, 2014). Briefly, these cloned rice genes for grain size are mainly involved in multiple signaling pathways, including ubiquitination-mediated proteasomal degradation, phytohormones, and G protein signaling, to regulate cell division and cell expansion (Zuo and Li, 2014;Li and Yang, 2017). In wheat, a wealth of QTL for grain weight have been identified to date on almost all wheat chromosomes based on linkage mapping in bi-parental genetic populations and genome-wide association study (GWAS), and a number of genes associated with grain weight have been isolated using comparative genetic analysis (Quarrie et al., 2005;Zhang et al., 2010;Liu G. et al., 2014;Gao et al., 2015;Zanke et al., 2015;Li and Yang, 2017;Nadolska-Orczyk et al., 2017;Sajjad et al., 2017;Shao et al., 2017;Zhai et al., 2017). Nonetheless, to the best of our knowledge, there are no reports of cloned major QTL for grain weight through map-based cloning approach in wheat, and the molecular roles of QTL/genes in the regulation of grain weight are still largely unknown.
Wheat is a typical cool season crop, and increasing temperature is a major limitation to further improve wheat yield potential (Wardlaw et al., 1989;Acuna-Galindo et al., 2015;Ni et al., 2017). It was estimated that each temperature increase of one degree Celsius reduces grain weight per spike (GWS) by 3-4% (Wardlaw et al., 1989) and global wheat production by 6% (Asseng et al., 2015). Correspondingly, high temperatures can reduce GNS at the pre-flowering and flowering stages, and reduce grain weight at early grain filling stage (Sharma et al., 2016;Shirdelmoghanloo et al., 2016;Valluru et al., 2017). Therefore, the identification of stable and robust QTL for yield-related traits under varying severities of heat-prone environments is useful for maintaining wheat adaptability and production stability against a backdrop of fluctuating climate change patterns (Ogbonnaya et al., 2017). Moreover, heattolerant QTL mapping and the identification of traits associated with heat tolerance of yield components for common QTL are pre-requisites for developing molecular markers suitable for heat tolerance breeding (Shirdelmoghanloo et al., 2016). Despite its importance, only a few QTL mapping studies for the heat tolerance of yield components in wheat have been reported (Mason et al., 2010;Paliwal et al., 2012;Shirdelmoghanloo et al., 2016;Ni et al., 2017).
A breakthrough in wheat production during modern breeding was the utilization of reduced height (Rht) loci, which led to the Green Revolution in the late twentieth century (Hedden, 2003). To date, more than 20 dwarfing genes have been reported in wheat, and Rht1 (Rht-B1b), Rht2 (Rht-D1b), and Rht8 are currently the three most commonly adopted dwarfing genes worldwide (Ellis et al., 2002;Zhang et al., 2006;Tian et al., 2017). The Rht1 and Rht2 belong to a group of genes known as gibberellic acid (GA) insensitive dwarfing genes and are located on chromosomes 4BS and 4DS, respectively (Peng et al., 1999;Achard et al., 2006). The two Rht-1 homoeoloci, Rht-B1 and Rht-D1, exert a pleiotropic effect on grain number, grain weight and yield in addition to reducing height (Zhang et al., 2013;Würschum et al., 2017). GA-responsive Rht8, located on chromosome 2DS, is another extensively used dwarfing gene, and it has been implemented in different environments because it has no effect on grain yield (Zhai et al., 2016;Tian et al., 2017). Furthermore, grain yield in wheat largely depends on plant architecture, particularly plant height (PH); thus, genetic loci associated with PH and yield components that are obtained by QTL mapping can provide a clear understanding of genetic relationships.
Recently, we developed a DH population derived from an elite cross of Nongda3338 (ND3338), and Jingdong6 (JD6) that exhibit contrasting phenotypes in PH, SPP, GNS, TGW, and GWS. Hence, the objectives of this study were to (i) evaluate the phenotypic performance of yield-related traits across different field trials with normal and late sowing mediated heat stress conditions; (ii) identify genomic regions with stable and robust QTL associated with yield-related traits; (iii) detect QTL controlling heat susceptibility index of thousand grain weight (HSITGW) for two contrasting treatments; and (iv) provide diagnostic markers to be deployed in marker-assisted selection (MAS) breeding for high-yield and heat-tolerant wheat varieties.

Plant Materials and Field Experiments
A DH population consisting of 203 individuals was developed through in vitro anther culture (De Buyser and Henry, 1980) of the F 1 hybrids from a cross between two Chinese elite winter wheat varieties, ND3338 and JD6. Briefly, the female parent ND3338 is a "core parental" breeding line for the North China Winter Wheat Breeding Program with high general combining ability developed by China Agricultural University, while JD6 is a variety released by Beijing Academy of Agricultural and Forestry Sciences. The results of functional molecular markers (Ellis et al., 2002) indicated that ND3338 has mutant (dwarf) Rht-B1b and Rht-D1b, while JD6 possesses wild-type Rht-B1a and Rht-D1a (Kabir et al., 2015). Additionally, the diagnostic molecular marker Xgwm261 for Rht8 was mapped on the short arm of chromosome 2D in our DH population with the 192-bp allele (Rht8c) from parent JD6.
The field experiments included two parts. First, the plants were grown at a conventional sowing time during the autumn across four different geographical locations in northern China: Beijing, Linfen, Shijiazhuang, and Urumqi. Dry-hot wind, defined as strong wind with high temperature and low humidity, often occurs during the grain filling period in the higher latitude areas of northern China (Teixeira et al., 2013;Liu B. et al., 2014;Ni et al., 2017). Second, two additional experiments were conducted under a timely sown condition and a latesown condition to expose plants to higher temperatures (heat stress), particularly during grain filling period, at two locations in China: Linfen and Sanyuan. Detailed environment characteristics are provided in Table 1. The method of late-sown trials was as described by Cheng et al. (2015). Meteorological data for the experiment sites are presented in Table S1. In each field environment, the 203 DH lines and their parents were planted in randomized complete blocks with three replicates. Each plot contained two rows that were 2 m long and 30 cm apart with a sowing rate at 30 seeds in each row. All fields were wellwatered by both rainfall and broad irrigation. Other management procedures of field trials followed local standard practices.

Phenotypic Evaluation and Statistical Analysis
Once the plants reached physiological maturity, 10 representative plants per genotype from each replication were used for phenotypic evaluation. PH (cm) was measured from the soil surface to the tip of the spike, excluding awns; SPP was also measured before harvesting. GNS, GWS, and TGW were measured after the seeds had naturally dried following harvest. Based on the TGW of two sowing dates, the HSI for each individual line was calculated using the formula by Fisher and Maurer (Fischer and Maurer, 1978): where Y heatstress and Y control are the TGW means for each genotype under heat-stressed and controlled conditions, respectively, and X heatstress and X control are the TGW means for all lines under heat-stressed and controlled conditions, respectively.
Basic statistical analyses, phenotypic correlation and Shapiro-Wilk tests for departure from normality were performed by SPSS software version 20.0 (SPSS, Chicago, USA). The adjusted mean (Best Linear Unbiased Prediction, BLUP) values across multiple environments were calculated using SAS v9.1.3 (SAS Institute Inc., North Carolina, USA) with the PROC MIXED procedure. Broad sense heritability (h 2 ) on a family basis was calculated with the PROC GLM procedure in SAS according to the following formula: h 2 = σ g 2 /(σ g 2 +σ ge 2 /n+σ 2 /nr), where σ g 2 is the genotypic effect, σ ge 2 is the genotype by environmental effect, σ 2 is the residual error, n is the number of environments and r is the number of replicates (Liu G. et al., 2014;Zhai et al., 2016). Linear regression analysis was conducted based on the additive effects in Microsoft Excel 2016. The biplot of principal component analysis (PCA) and additive main effects and multiplicative interactions (AMMI) analysis were performed in R software (v. 3.4.2) (R Core Team, 2013; Dixit et al., 2014).

Genotyping and Construction of Genetic Map
Young leaf tissues of the parents and DH lines at the seedling stage were used for total genomic DNA extraction with the modified CTAB method . Approximately 2,800 SSR and STS markers were used to detect polymorphisms between the two parents. Primer sequences for most SSR markers are publicly available at http://wheat.pw.usda.gov/GG2/index. shtml. The PCR system, DNA amplification condition and product fragment detection were determined as described by Liu G. et al. (2014). In addition, the two parents and 203 DH lines were also genotyped with the Illumina 90K iSelect wheat SNP assay  at the Genome Center at the University of California, Davis. SNP clustering and genotype calling were performed using GenomeStudio version 2011.1 software (Cavanagh et al., 2013). The integrated genetic map based on SSR and STS markers, SNP makers and gene specific markers were generated using the programs RECORD 2.0 (Van Os et al., 2005) and JoinMap 4.0 (Van Ooijen, 2006). SNP markers with large numbers of missing values (>20%) were discarded. The map construction procedure is based on the methodology described by Zhai et al. (2016). Genetic linkage maps were compared with consensus 90K SNP maps  to orient each linkage group with respect to the short (S) and long (L) chromosome arms, further checking the accuracy of the marker order. Additionally, genetic maps and QTL graphs were drawn using MapChart 2.2 software (Voorrips, 2002).

QTL Detection
Mean data of each trait for individual environments and the adjusted mean (BLUP) values across multiple environments were used for QTL analysis with Windows QTL Cartographer software version 2.5  through composite interval mapping (CIM). In QTL Cartographer, the parameters were as follows: model 6 (standard model), forward and backward regression, five control markers (co-factors), window size of 10 cM, and walk speed of 1 cM. An empirical genome-wide LOD threshold to identify significant QTL were calculated using 1,000 permutations for P ≤ 0.05. Confidence intervals were estimated based on positions ± 2 LOD (from the peak) method using QTL Cartographer. QTL with overlapping confidence intervals or QTL located within 10 cM region were considered equivalent. Only QTL that were significant at a LOD value ≥ 2.5 were accepted in this study. QTL names were denoted according to the International Rules of Genetic Nomenclature (http://wheat. pw.usda.gov/ggpages/wgc/98/Intro.htm). QTL for traits that colocalized within the same genomic region were assigned a common QTL name.

Bioinformatics Analysis
The SNP flanking sequences mapped in the integrated genetic map were aligned with respect to the newly released bread wheat Chinese Spring Reference sequence by International Wheat Genome Sequencing Consortium (IWGSC) (http:// www.wheatgenome.org/, IWGSC RefSeq v1.0) and the coding sequences (CDS) with high-confidence genes to obtain physical positions and candidate genes. In addition, the annotation of high-confidence genes was obtained from wheat genome database websites (https://wheat-urgi. versailles.inra.fr/Seq-Repository/Annotations and http://plants. ensembl.org/Triticum_aestivum/Info/Index). Synonymous or nonsynonymous SNPs were annotated, and synteny analyses with rice genomes were performed as described by Wang et al. (2014).

High-Density Integrated Genetic Linkage Map Construction
A coarse-scale linkage map was initially constructed over the whole wheat genome with 475 SSR and STS markers using 203 DH lines. To construct linkage map for the unlinked regions in the initial map, 81,587 SNPs from the wheat 90K SNP array were used for genotyping the ND3338/JD6 DH population , and 10,409 (12.76%) SNP markers that showed polymorphism between the parents were used for linkage analysis. Forty-six of these SNP markers were discarded because they had more than 20% missing data, and 176 SNP markers were not anchored on the linkage map. Finally, a high-density integrated genetic linkage map based on 475 SSR and STS markers, 10,187 transcript-derived SNP markers and two gene functional markers (Rht1 and Rht2) was constructed for QTL mapping (Tables S3, S4). The 10,664 markers represented a total of 2,017 unique loci (18.91%) distributed among 26 linkage groups (LGs) representing the 21 hexaploid wheat chromosomes (Tables S3-S5). All linkage maps covered 3,391.20 cM in length with an average density of 1.68 cM/locus, and the map length of the three A, B, and D genomes was divided approximately equally (Table S3). However, the distribution of markers in the genomes was not uniform with about five times as many polymorphic markers mapping to the A and B genomes than to the D genome, namely, 4,692, 4,530, and 965 markers for the A, B, and D genomes, respectively (Table S3). In total, 958 (47.50%) of the 2,017 unique loci in the ND3338/JD6 integrated linkage map had segregation ratios that deviated significantly (chi-square ≤ 0.05) from the expected 1:1 ratio (Table S5). Of these loci, 527 (55.01%) favored the male parent (JD6), and 431 (44.99%) loci favored the female parent (ND3338).

Phenotypic Performance Across Multi-Environments
Two parents and 203 DH lines were trialed in a multienvironment design, including different locations and growing seasons, to identify stable QTL for yield-related traits across different field trials with timely sown (normal) and late sown (heat stress) conditions ( Table 1). The DH population means and ranges of five yield-related traits (PH, SPP, GNS, TGW, and GWS) across the eight shared environments are shown in Table  S2. JD6 had higher PH, TGW, and GWS and lower SPP and GNS than ND3338 (Table 2, Figure S1). PH, SPP, TGW, and GWS displayed obvious deviations from normality in the population, whereas GNS exhibited normal distributions with the BLUP values of eight shared environments ( Table 2). All traits had h 2 over 0.80, and the highest h 2 was observed for PH, which reached 0.99, suggesting that PH was controlled by major effect genes in the DH population ( Table 3). Pearson's correlation coefficients of the five traits were estimated based on the BLUP values of eight shared environments, which showed that SPP was significantly and negatively corrected with PH, GNS, TGW, and GWS ( Table 4). TGW displayed strongly positive correlations with PH and GWS and a significantly negative correlation with GNS ( Figure 1, Table 4). ANOVA with the eight shared environments for five yield-related traits revealed that there were significant variations from the environments and DH lines by environment interactions (Table 3). PCA biplot for environmental variability showed the differences among the eight shared experimental sites (Figure 1). For contrasting field trials, means of TGW in the DH population showed significant reduction under the late-sown heat-stressed condition compared to the timely sown condition, which was consistent with meteorological data ( Figure S2, Table  S1). Meanwhile, the HSI of TGW was calculated to assess the heat tolerance performance of the DH population, which suggested that a part of the progeny had transgressive phenotypes and the DH lines showed a continuous normal distribution based on the Shapiro-Wilk test ( Figure S2).

QTL Mapping Analysis
In this study, a total of 226 QTL controlling five yield-related traits and the heat susceptibility index (HSI) were detected across 12 different environments using CIM (Tables S6-S8). QTL that were repeatedly detected in ≥ 3 individual environments and in the BLUP analysis were considered to be stable. According to this criterion, 50 stable QTL for PH, SPP, GNS, TGW, and GWS were identified; and of these, 39 stable QTL were mapped     Table 1.
within 12 genomic regions with corresponding physical intervals of Chinese Spring RefSeq v1.0 sequence on chromosomes 1B, 2A, 2B, 2D, 3A, 4A, 4B, 4D, 5A, 6A, and 7A ( Table 5, Table S6). The other 168 putative QTL for five yield-related traits are listed in Table S7, and 8 QTL for the HSI of TGW were identified in two conditions (Table S8). Detailed parameters of QTL detected for each trait and environment are as follows.

Grain Number Per Spike (GNS)
For GNS, a total of 34 QTL were identified across eight shared environments with a LOD score range of 2.51-13.48, explaining 3.13-21.36% of the variance (Tables S6, S7

Thousand Grain Weight (TGW)
For TGW, a total of 69 QTL were identified in 12 environments that include eight shared environments, two timely sown controlled environments and two late-sown heat-stressed environments with a LOD score range of 2.51-32.07, explaining 1.96-48.84% of the variance (Tables S6, S7). Of these, 13 stable QTL for TGW were mapped on chromosomes 2A, 2D, 4A, 4B, 5A, 6A, and 7A, and the favorable alleles were contributed by both parents, ND3338 and JD6 (Tables 5, 7, Table  S6), which was consistent with regression analysis (Figure 3).

Heat Susceptibility Index of TGW (HSITGW)
The analysis for the HSI was conducted to detect QTL directly related to stress performance traits. Using composite interval mapping, we detected 3 and 6 QTL for the HSITGW in E11  and E12 locations, respectively (Table S8). Individual QTL explained from 4.84% to 14.29% of the phenotypic variance for the HSITGW, and both parents contributed favorable alleles. QHsitgw.cau-4B.2 was repeatedly detected in both locations, which explained 9.23-9.92% of the total variation and co-localized with a minor QTL (QTgw.cau-4B.7) for TGW (Figure 2, Tables 7, 8, Tables S7, S8). The parent JD6 contributed the positive allele, enhancing heat tolerance. In addition, the localization of QHsitgw.cau-2D on chromosome 2DS coincided with QTgw.cau-2D.2, explaining the highest observed variation (14.29%) for HSITGW in the E12 location (Tables S7, S8).

Phenotypic Variations in Response to Environments
Phenotypic characterization of yield-related traits over multienvironment field trials is essential for assessing trait stability across environments and contributes to accurate identification of stable genomic regions (Ogbonnaya et al., 2017). With rising global temperature and climate change, heat stress is becoming an increasingly severe constraint on wheat production in many parts of the world (Paliwal et al., 2012;Ni et al., 2017). In the present study, the two parents, ND3338 and JD6, are elite winter wheat varieties from the Northern Winter Wheat Zone in China, and Beijing and Shijiazhuang sites are optimum growing zones. The Sanyuan site is located in the Yellow and Huai River Valleys Facultative Wheat Zone. The Linfen and Urumqi sites were hotspot regions of hot and dry winds during the grain filling period (Liu B. et al., 2014;Yang et al., 2017). Moreover, AMMI biplot based on the population means of five yield-related traits across eight shared sites depicted the environmental difference and genotype by environment interaction (Figure 1), indicating that the identified QTL for yield-related traits showed stability in a wide range of environments. Two additional experiments under high environmental temperatures due to a delayed planting date were conducted to further evaluate the performance of the parents and the DH population under heat stress. Consequently,
late-sown environments in Linfen and Sanyuan with a maximum temperature >35 • C during the grain filling stage resulted in 20.09-36.36% reduction of grain weight compared to the control, which is consistent with the results of AMMI analysis (Figure 5).
In this study, three major stable QTL (QPh.cau-4B.2, QPh.cau-4D.1, and QPh.cau-2D.3) for PH were detected on chromosomes 4BS, 4DS, and 2DS at positions within the three dwarfing genes Rht-B1, Rht-D1, and Rht8, respectively, which agrees with results from previous studies (Peng et al., 1999;Zhang et al., 2006;Gao et al., 2015). A stable locus (QPh.cau-5A.2) that had an effect on PH was found on chromosome 5A at 86.70-95.00 cM; similarly, the dwarfing gene Rht9 was previously reported to be located on chromosome 5AL (Wu et al., 2010). Furthermore, a minor but stable PH QTL, QPh.cau-6A.2, was detected across all environments in this study. In parallel, Tian et al. (2017) also found a major quantitative trait locus (QPH.caas-6A) in a similar region on chromosome 6A, which was designated Rht24. After that, Würschum et al. (2017) assessed the relevance of Rht24 using an association mapping approach based on a large panel of 1,110 winter wheat cultivars, suggesting that Rht24 was an important Rht gene of commercial relevance in worldwide wheat breeding. To the best of our knowledge, the stable QTL QPh.cau-1B.2, QPh.cau-2D.2, QPh.cau-3A.2, QPh.cau-3A.3, QPh.cau-4D.2, QPh.cau-6A.3, and QPh.cau-7B were likely novel QTL for PH owing to the high-density integrated genetic linkage map and special performance of the genetic background between ND3338 and JD6. These novel loci displayed relatively smaller additional effects compared with Rht1, Rht2, and Rht8; hence, near isogenic lines are ideally required in future work to determine their effects on PH and other agronomic traits dependent on the environment.
For GWS and yield components, we detected a QTL cluster for TGW and GWS located in the interval 112.90-133.80 cM on chromosome 2AL and a QTL cluster for GNS and GWS located in the interval 111.60-131.90 cM on chromosome 2B; similar genomic regions were previously reported to harbor QTL for GNS, TGW, and GWS using the different recombinant inbred line (RIL) populations (Cui et al., 2014;Liu G. et al., 2014;Gao et al., 2015.) The genomic regions with the highest numbers of stable QTL for yield components were on group-4 chromosomes near the Green Revolution genes Rht1 and Rht2, which agrees with results from earlier reports (Zhang et al., 2013;Liu G. et al., 2014). Furthermore, the stable QTL for GNS and TGW were found on chromosome 5A; and previous studies also identified stable loci for GNS and TGW at the similar interval based on linkage analysis and GWAS approach (Cui et al., 2014;Gao et al., 2015;Wu et al., 2015;Ogbonnaya et al., 2017). Likewise, another QTL-rich region associated with GNS, GWS, and TGW was identified on chromosome 7A in this study; consistently, Quarrie et al. (2006) reported the co-localized QTL cluster for yield and yield-related traits on chromosomes 7A in the Chinese Spring/SQ1 DH population under stressed and non-stressed conditions, some of which were also validated using near isogenic lines.
In this study, we identified 8 QTL for the HSI of TGW. Among these QTL, QHsitgw.cau-6D was identified on chromosome 6DL in the Sanyuan location, which is consistent with the study of Mason et al. (2010), which identified QTL for the HSI of kernel weight in a similar position using a RIL population derived from a cross between the heat-tolerant cultivar "Halberd" and heat-sensitive cultivar "Cutter." Notably, the locus QHsitgw.cau-4B.2 was detected in the two heat-stressed trials and was located on chromosome 4BL in the interval 115. 80-123.10 cM, spanning 14.28 Mb (4B: 620058668-4B: 634340181) in physical position (Tables S4, S8). Moreover, IAAV5323 is a SNP in the confidence interval within the protein coding region of the annotated gene TraesCS4B01G340600 encoding Zinc finger protein CONSTANS (Table S10), which is known to be involved in the regulation of multiple processes, including flowering time, hormone metabolism, and biotic and abiotic stresses (Griffiths et al., 2003;Weng et al., 2014;Kiss et al., 2017). Thus, we postulate that TraesCS4B01G340600 is a possible candidate gene for the QHsitgw.cau-4B.2 QTL. Moreover, to the best of our knowledge, QHsitgw.cau-4B.2 is a novel stable QTL for the HSI of TGW, and diagnostic molecular markers can be developed and deployed within breeding programs. Additionally, the SNP maker IACX4386 for QHsitgw.cau-1A corresponded to the gene TraesCS1A01G285000 encoding a 70-kDa heat shock protein; this result merits further investigation (Table S10).

QTL/Genes Controlling PH Showed Pleiotropic Effects on Yield Components
PH is an important agronomic trait in wheat, and wheat yield increases during the Green Revolution were achieved through the introduction of reduced height (Rht) dwarfing genes (Hedden, 2003;Zhang et al., 2006). QTL mapping in previous studies confirmed that PH is a complex trait controlled by the few major Rht loci and by minor QTL (Wu et al., 2010;Griffiths et al., 2012;Würschum et al., 2015Würschum et al., , 2017Tian et al., 2017). In this study, three major stable QTL (QPh.cau-4B.2, QPh.cau-4D.1, and QPh.cau-2D.3) for PH corresponded to Rht1, Rht2, and Rht8. Haplotype analysis revealed that Rht loci exhibited pleiotropic architecture, which affects not only PH but also yield component traits ( Table 6). Significant differences were detected between Haplotype 1 and Haplotype 3 for PH, GNS, and TGW with the BLUP values. The dwarf gene Haplotype 2 and Haplotype 3 showed pleiotropic effects on PH, GNS, and GWS compared with Haplotype 1 (Table 6). Furthermore, Haplotype 8 exhibited the lowest PH, TGW, and GWS but the most SPP, indicating that shorter plants lead to the reduction of total biomass yield. Notably, for GNS, Haplotype 8 displayed no significant difference from Haplotypes 1, 2, 4, and 5 (Table 6). Additionally, previous research found that the two Rht-1 semi-dwarfing genes only improved yields under optimal conditions, whereas tall isogenic lines without Rht-B1b or Rht-D1b yielded more than their Rht-1b carrying counterparts in biotic and abiotic stresses environments (Zhang et al., 2013;Würschum et al., 2017). Similarly, there is a significant difference in thermotolerance between Haplotype 1 and Haplotype 3 in our study (Table 6). Therefore, the choice of the Rht loci best suited for achieving the desired plant stature must account for pleiotropic effects on yield components in different geographic and climatic regions.
Co-location of Yield Components Revealed Significant Tradeoffs Between TGW and GNS on Chromosome 4A TGW and GNS are two crucial but counteracting determinants of wheat grain yield (Quarrie et al., 2006;Schulthess et al., 2017;Zhai et al., 2017). In the present study, we identified two stable QTL on chromosome 4A controlling TGW (QTgw.cau-4A.3) and GNS (QGns.cau-4A.4) with superior alleles coming from the opposite parent (Table 5), which exhibited a strong pleiotropic trade-off between TGW and GNS, consistent with phenotypic correlation analysis (Table 4). For TGW, QTL in the same region of chromosome 4AL were reported by Cui et al. (2016) and Gao et al. (2015); however, to the best of our knowledge, this is the first report about stable locus QTgw.cau-4A.3 for TGW under heat-stressed conditions. For GNS, interestingly, Cui et al. (2017) also identified a major stable QTL for kernel number per spike (KNPS) in 10 environments using a Wheat660K SNP array-derived high-density genetic map, and the overlapping confidence intervals spanned 3.23 Mb (4A: 680398739-4A: 683638403) in physical position, which is consistent with our results (4A: 632020468-4A: 684998591). Moreover, a QTL cluster for TKW and KNPS positioned in the interval 136.80-157.30 cM on chromosome 4AL corresponded to a similar physical interval (4A: 632864778-4A: 688093018) detected by Gao et al. (2015) using the Zhou 8425B/Chinese Spring RIL population, and the positive alleles increasing TGW and GNS were all contributed by Zhou 8425B, showing no negative pleiotropic trade-off. Hence, these coincidences confirmed the authenticity of this locus, which should be subjected to fine mapping and map-based cloning in the future. Only these techniques can distinguish between whether these loci have a pleiotropic effect or are closely linked. Additionally, we performed corresponding gene annotations and synteny analyses with rice genomes with the Chinese Spring IWGSC RefSeq v1.0 based on the SNP marker flanking sequences in the confidence interval (Table S9).

Dissection of a "QTL-Hotspot" Region on Chromosome 4B Was Useful in MAS for Grain Weight
Generally, at least three QTL for yield-related traits were identified on every chromosome in our study with chromosome 4B having the most QTL (Tables S6, S7). Thus, the shared genomic region with a pleiotropic effect or tightly linked loci affecting two or more traits on chromosome 4B is referred to as a "QTL-hotspot" region. Based on the positions of stable QTL for TGW, the "QTL-hotspot" region on chromosome 4B was artificially divided into three parts: region 4B.1, region 4B.2, and region 4B.3 (Figure 2). For region 4B.1, the stable QTL controlling TGW (QTgw.cau-4B.1 and QTgw.cau-4B.2) and SPP (QSpp.cau-4B.3 and QSpp.cau-4B.4) with favored alleles coming from the opposite parent were detected, exhibiting strong TGW-SPP tradeoffs, and this result was consistent with the phenotypic negative correlation. Similarly, Chen et al. (2014) identified a stable QTL for TGW on chromosome 4BS with the same SNP makers using the Shannong 01-35/gaocheng9411 RIL population, and it was considered to be a novel QTL for TGW in their results. For region 4B.2, the stable QTL (QTgw.cau-4B.3) was co-localized with one QTL for PH (QPh.cau-4B.2), and increasing alleles all came from JD6. Meanwhile, it was corroborated as a pleiotropic consequence of the dwarfing gene Rht1 in the present study, consistent with the results previously published in wheat (Liu G. et al., 2014;Schulthess et al., 2017). For region 4B.3, we found that the stable QTL (QTgw.cau-4B.4) was also constitutively expressed in two heat-stressed environments, and more importantly, it displayed no negative pleiotropic association with the other two yield component traits (GNS and SPP). Although previous research also reported the detection of QTL for grain weight or grain dimensions on chromosome 4B through linkage analysis and genome-wide association analysis (Quarrie et al., 2005;Liu G. et al., 2014;Huang et al., 2015;Wu et al., 2015;Zanke et al., 2015;Cui et al., 2016;Kumar et al., 2016), to the best of our knowledge, this is the first report on the dissection of the QTL cluster on chromosome 4B for TGW. We speculate that the "QTL-hotspot" region on chromosome 4B may have experienced artificial selection during past breeding practice and played an important role in grain yield determination and improvement, especially on grain weight. Thus, dissection of a "QTL-hotspot" region on chromosome 4B was useful for future map-based cloning and MAS-based QTL pyramiding for grain weight.

QTL Combinations for TGW and HSITGW Explored the Value of Breeding
To explore the effects of different QTL combinations on grain weight, the favored or unfavored allele effects from ND3338 and JD6 were simulated. The patterns of the relationships were similar, where TGW increased gradually along with increases in favored alleles from ND3338 or JD6 (Figure 3). Outstandingly, favored alleles from JD6 showed higher R-squared values in linear regression analysis (R 2 = 0.53) than the values for favored alleles from ND3338 (R 2 = 0.12), consistent with the total phenotypic variation explained by QTL analysis (Tables S6, S7). Moreover, AMMI analysis was conducted to rank genotypes based on TGW across 12 environments, and the top 10 high-TGW genotypes and top 10 low-TGW genotypes were identified to determine the genetic composition of stable QTL for TGW. Among the 13 stable QTL for TGW, top 10 high-TGW genotypes possess at least 7 favorable alleles from ND3338 or JD6, whereas top 10 low-TGW genotypes have no more than 5 favorable alleles (Table 7). Moreover, top 10 high-TGW genotypes all have the major robust QTL (QTgw.cau-4B.4) that harbors favored alleles from JD6 (Table 7), whereas the stable QTL (QTgw.cau-4B.1 and QTgw.cau-4B.3) that harbor unfavorable alleles from ND3338 were all observed in top 10 low-TGW genotypes ( Table 8). Additionally, 9 out of the top 10 high-TGW genotypes for QHsitgw.cau-4B.2 carried the favored alleles that increase thermotolerance, and only three out of the top 10 low-TGW genotypes carried the favored alleles (Tables 7, 8). Furthermore, the top 10 high-TGW genotypes exhibited relatively more stable performance across all seasons compared to the top low-TGW genotypes based on the principal component values (PC1 and PC2). Taken together, our results provide a basis for developing new heat-tolerant wheat varieties with high yield stability via a molecular design breeding strategy to fix good additive alleles.

CONCLUSION
Overall, this study identified a total of 226 QTL controlling five yield-related traits (i.e., PH, SPP, GNS, TGW, and GWS) and the heat susceptibility index (HSI) in Nongda3338/Jingdong6 DH population across 12 different field trials with normal and late sowing heat stress conditions. Of these, 39 stable QTL for PH, SPP, GNS, TGW, and GWS were mapped within 12 genomic regions with corresponding physical intervals of Chinese Spring RefSeq v1.0 sequence on chromosomes 1B, 2A, 2B, 2D, 3A, 4A, 4B, 4D, 5A, 6A, and 7A. Three QTL QPh.cau-4B.2, QPh.cau-4D.1 and QPh.cau-2D.3 corresponded to dwarfing genes Rht1, Rht2, and Rht8, which had the pleiotropic effect on yield component traits. A QTL-hotspot region on chromosome 4B for grain weight and a novel QTL for HSITGW on chromosome 4BL were detected. These results will contribute to our understanding of the genetic basis of yield-related traits, and could be used to improve grain yield and in wheat through MAS breeding after validation.

AUTHOR CONTRIBUTIONS
HP: conceived the project; PG, LL, and LJ: carried out experiments; MK, JZ, TL, and YZ: participated in field trials; PG: analyzed experimental results; PG and HP: wrote the manuscript; MX, ZH, YY, ZN, and QS: helped to revise the manuscript. All authors have read and approved the final manuscript.