Dissection of a grain yield QTL from wild emmer wheat reveals sub-intervals associated with culm length and kernel number

Genetic diversity in wheat has been depleted due to domestication and modern breeding. Wild relatives are a valuable source for improving drought tolerance in domesticated wheat. A QTL region on chromosome 2BS of wild emmer wheat (Triticum turgidum ssp. dicoccoides), conferring high grain yield under well-watered and water-limited conditions, was transferred to the elite durum wheat cultivar Uzan (T. turgidum ssp. durum) by a marker-assisted backcross breeding approach. The 2B introgression line turned out to be higher yielding but also exhibited negative traits that likely result from trans-, cis-, or linkage drag effects from the wild emmer parent. In this study, the respective 2BS QTL was subjected to fine-mapping, and a set of 17 homozygote recombinants were phenotyped at BC4F5 generation under water-limited and well-watered conditions at an experimental farm in Israel and at a high-throughput phenotyping platform (LemnaTec-129) in Germany. In general, both experimental setups allowed the identification of sub-QTL intervals related to culm length, kernel number, thousand kernel weight, and harvest index. Sub-QTLs for kernel number and harvest index were detected specifically under either drought stress or well-watered conditions, while QTLs for culm length and thousand-kernel weight were detected in both conditions. Although no direct QTL for grain yield was identified, plants with the sub-QTL for kernel number showed a higher grain yield than the recurrent durum cultivar Uzan under well-watered and mild drought stress conditions. We, therefore, suggest that this sub-QTL might be of interest for future breeding purposes.


Introduction
The meteorological concept of drought can be defined as "an extended period of time characterized by a deficiency in a region's water supply" (Below et al., 2007). Stress might be defined as an altered physiological condition that alters an equilibrium, leading to a strain, i.e., a biochemical and/or physical change, which can lead to injury, disease, or aberrant physiology of the plant (Gaspar et al., 2002). Plants can experience drought stress either when water supply becomes limited or when the transpiration rate becomes too high (Fahad et al., 2017). Drought stress effects on plants are, e.g., the reduction of water content, diminished leaf water potential and turgor loss, closure of stomata, decrease in cell enlargement/growth due to low turgor pressure, and arrest of photosynthesis, which can finally result in the disturbance of metabolism and death. Increased periods of drought were shown to have major negative effects on the yield (Lesk et al., 2016;Fahad et al., 2017). Potential yield losses differ according to the magnitude of stress and the developmental stage of the plant. Extreme yield losses were reported to be up to 92% (Farooq et al., 2014). Plant responses to drought stress can be divided into long-and short-term responses on the biochemical, molecular, and finally physiological levels. These responses lead to different strategies to deal with drought stress, i.e., escape, tolerance, avoidance, and recovery Fang and Xiong, 2015).
Wheat (Triticum spp.) is among the top five grown crops in the world (Oyewole 2016;FAO 2019) and provides about 20% of the total calories consumed by humans (Poole et al., 2021). Future plant production faces the challenge of feeding the growing global population, which is predicted to reach about 9.8 billion people by 2050 (UN 2017). This challenge is increased by climate change, because per degree-Celsius increase in the global mean temperature wheat yield will be decreased by 6.0% (Fahad et al., 2017;Zhao et al., 2017). Developing high-yielding cultivars under drought with yield stability between environments is, thus, of prime importance (Langridge and Fleury, 2011).
The genetic diversity of tetraploid and hexaploid wheat has been depleted by the limited number of allopolyploidization events, domestication, and modern breeding (Dempewolf et al., 2017). Using crop wild relatives (CWR) such as wild emmer wheat in pre-breeding turned out to be an efficient tool to exploit their genetic diversity (Xie and Nevo, 2008;Huang et al., 2016;Dempewolf et al., 2017). A significant increase in the use of CWR has been noted since 1980 (Dempewolf et al., 2017). However, since abiotic stress tolerance is complex, the introgression of quantitative trait loci (QTL) from CWR aiming to confer tolerance to abiotic stress into elite cultivars is difficult (Xie and Nevo, 2008;Dempewolf et al., 2017). Traits that are advantageous for the CWR's fitness can be detrimental for breeding purposes. An increased culm length, for instance, might allow the plant to compete for light (Lane et al., 2000) or to promote the distribution of the plants' own seeds in the ecosystem-even if this leads to a trade-off in grain yield (GY). Avoiding trans-, cis-, or linkage drag effects from the CWR is, therefore, one of the most difficult issues in pre-breeding.
In 2017, the near-isogenic line NIL-U-2B-1 carrying a wild emmer QTL region on chromosome 2BS was shown to produce more GY than its recurrent elite durum parent, under water-limited and well-watered conditions (Merchuk-Ovnat et al., 2016;Merchuk-Ovnat et al., 2017). It also turned out that negative traits such as an increased culm length or delayed heading were introduced, possibly due to trans-, cis-, or linkage drag effects from the CWR (Merchuck-Ovnat et al., 2016). A recent analysis of NIL-U-2B-1 revealed that more than two-thirds of chromosome 2B and additional fragments of chromosomes 2A, 3A, and 5A were also introgressed from the donor wild emmer parent (Deblieck et al., 2020). Consequently, the original RIL mapping population, previously used for mapping the QTL on chromosome 2B, was re-genotyped with the 15k iSelect chip (TraitGenetics) to narrow down the previously identified QTL regions by a higher marker density of~4,000 SNPs (Fatiukha et al., 2021), which is much higher than the density used for the original mapping (Peleg et al., 2008;Peleg et al., 2009). The existence of the corresponding 20 cM large QTL-region for GY on chromosome 2BS was again confirmed under both water-limited and well-watered conditions (Fatiukha et al., 2021). Therefore, this study aimed at the fine mapping of a 15.67 cM region of this 20 cM large QTL, to identify the smallest sub-QTL region affecting grain yield, and to diminish the observed trans-, cis-, or linkage drag effects from the CWR. For this purpose, we used the iSelect marker Tdurum_contig27976_414 and wsnp_Ex_c6537_1133876, flanking the corresponding QTL-interval (Fatiukha et al., 2021) and 10 additional markers within the QTL, to establish a population of sub-NILs, carrying segmental chromosomal substitutions of the target region.

Plant material
All plant materials were developed at the Hebrew University of Jerusalem (HUJI) in Rehovot (31°54′N, 34°47′E, 54 m above sea level). In past works, recombinant inbred lines (RILs), derived from a cross between durum wheat (cv. Langdon) and wild emmer wheat (acc# G18-16) were phenotyped under contrasting water availabilities and used for QTL mapping (Peleg et al., 2009;Fatiukha et al., 2021). Subsequently, the wild emmer QTL allele on chromosome 2BS conferring higher grain yield and harvest index (HI) was introgressed into the elite tetraploid durum wheat cultivar Uzan to develop the near-BC 3 F 3 isogenic line (NIL) NIL-U-2B-1 as previously described (Merchuk-Ovnat et al., 2016;Deblieck et al., 2020). In this study, BC 3 F 3 NIL-U-2B-1 (pollinator) was crossed with Uzan (Supplementary Figure S1) and BC 4 F 3 plants were screened with the respective flanking molecular markers to identify heterozygous recombinants. The heterozygous recombinant BC 4 F 3 plants were selfed and BC 4 F 4 single seed descendants were again genotyped to identify homozygous recombinant plants for seed multiplication. Finally, BC 4 F 5 seeds of 17 homozygous recombinant plants were used for phenotypic experiments in 2017, 2018, and 2019.

Development of molecular markers and genotyping by sequencing (GBS)
The 15k iSelect data of G18-16, Langdon, Uzan, and NIL-U-2B-1 were previously published (Merchuk-Ovnat et al., 2016;Soleimani et al., 2020). Tdurum_contig27976_414 (32.88 cM) and wsnp_Ex_c6537_1133876/IAAV980 (48.55 cM) which flank the QTL-interval on chromosome 2B (Fatiukha et al., 2021) were converted into competitive allele-specific PCR (KASP) markers (Wilkinson et al., 2016). In addition to this, a total of 10 other PCR-based molecular markers, including KASP, cleaved amplified polymorphic sides (CAPS), and simple sequence repeats (SSRs) were developed along the QTL-interval at a distance of~1 cM to genotype and identify heterozygous recombinant BC 4 F 3 sub-NILs. All molecular markers were tested and validated on the DNA of F 7 lines of the original mapping population (Peleg et al., 2008;Peleg et al., 2009;Fatiukha et al., 2021). Furthermore, genotyping by sequencing (GBS) was applied to NIL-U-2B-1, the parental lines, the wild emmer wheat acc. G18-16, the durum elite parent Uzan, and homozygous sub-NILs with recombination close to the marker Gene-1741_103, which was previously described to be the nearest marker to the target QTL interval (Fatiukha et al., 2021). This additional step served to identify new recombination events within the recombinant sub-NILs to further improve finemapping efficiency.
GBS libraries were prepared as described previously (Elshire et al., 2011) and sequenced (150 bp paired-end, Illumina MiSeq). The sequencing produced millions of reads. These were de-multiplexed according to the barcodes and the adapters/barcodes using the CASAVA pipeline 1.8 (Illumina, Inc.). Trim Galore software from Babraham Bioinformatics (2012) was used for adapter and quality trimming of the amplified genomic sequences. After this first filtering, the trimmed sequences were aligned to the draft genome sequence of the wild emmer acc. Zavitan (v1) (Avni et al., 2017) using BWM-MEM (version: 0.7.7.-r1140) (Li, 2013), and variant-calling was performed with the samtool and bcftools (version: 0.1.19-96) (Li et al., 2009). High-quality bi-allelic SNPs were filtered, and the imputation of missing SNPs was conducted with Beagle (Browning and Browing, 2016). Aligned sequencing reads were used for SNP detection after a quality check (Q score >20). Multi-allelic SNPs, SNPs with minor allele frequency (MAF) < 5%, missing values ≥5%, or heterozygosity ≥90% were further excluded and a highquality SNP genotyping dataset was compiled. SNPs which could clearly be assigned to a unique position on the physical genome of wild emmer were kept for analysis. GBS data were finally evaluated with the GenoTypeMapper (GTM) (Deblieck et al., 2020).

Plant growth conditions and experimental design
A total of five plants per genotype were phenotyped under contrasting water-limited (WL) and well-watered (WW) conditions from 2017-2019.
In 2017 and 2018, plants were grown at the HUJI experimental farm at Rehovot. Seedlings were first placed in moist germination paper at 4C°in a dark vernalization room for 2 weeks and then transplanted at the beginning of December into an insect-proof screen house rain-protected by a polyethylene top. The soil was a brown-red sandy loam (Rhodoxeralf) composed of 76% sand, 16% clay, and 8% silt (Merchuk-Ovnat et al., 2016). The different water regimes were simulated in a factorial (genotype x irrigation regime) split plot block design. Each block consisted of two main plots (for the two irrigation regimes), split into subplots for genotypes. Each subplot consisted of a single row with five plants, 10 cm apart (50-cm long plots). In each bed, two 40 cm spaced rows were planted with 100 cm between each pair of rows (Supplementary Figure S2). Seasonal rainfall was simulated by applying water once or twice a week from planting in December to heading in April/May, leading to a total seasonal water application of 350 mm (for WL conditions) and 650 mm (for WW conditions) in 2017. However, drought stress appeared to be mild in 2017. To increase drought stress, 384 mm and 201 mm of water were applied to the segmental sub-NILs in 2018. Plants that grew poorly in a certain area of the screen house were excluded from the experiment, keeping three replicates per genotype.
In July 2019, one plant per genotype was grown per pot on the high throughput phenotyping (HTP) facility (LemnaTec-129 Scanalyzer 3D) (http://www.lemnatec.com) of the Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) in Gatersleben in Germany (Supplementary Figure S3). To avoid the confounding effect of flowering time on drought evaluation, sub-NILs with a much later flowering time were not included in the HTP experiment. The controlled greenhouse is described by Neumann et al. (2015) along with the potting procedure. In short, pots (2 L) were filled in a standardized way by weighing in the same amount of soil as a standard substrate ("Klasmann Substrate no. 2" (http://www.klasmann-deilmann.com)) allowing a slow development of drought as described by Neumann et al. (2015). After sowing, 7 g of fertilizer was added to each pot (19% N, 9% P 2 O 5 , and 10% K 2 O) to supply the plants with nutrients throughout the life cycle. After 10 days, the plants were thinned out to one plant per pot. Directly after the daily image recording, each pot was weighed individually and watered to a previously defined target weight according to Dhanagond et al. (2019). Plants were grown on the platform from sowing until maturity and watered daily to the corresponding target weights.
The aim of the experiment was to mimic as good as possible the average temperatures and slowly progressing drought conditions of the screen house experiment conducted in Israel. Drought thresholds for severe drought stress were defined as 20% plant-available water (PAW) based on the results in barley (Dhanagond et al., 2019), while mild drought was defined as 30% PAW. No drought stress was applied until 30 days after sowing (30 DAS). The temperature during this predrought phase was set to 12°C at night and 16°C during the day. Supplementary greenhouse lights provided light for 15 h per day. To simulate drought stress, irrigation of plants intended for stress treatment was changed 31 DAS from 90% PAW to 30% PAW (mild drought threshold) and the temperature was raised to 20°C during the day and 16°C at night. At 62 DAS, the temperature was further increased to 24°C during the day and 20°C at night. From 64 DAS, watering was reduced to 20% PAW (severe drought level). This drought level and temperature regime lasted until maturity. At maturity, the plants were subjected to a detailed assessment of agronomic traits. Furthermore, the raw images were inspected, and the date of heading (BBCH 55) was determined for each plant.
Phenotypic data, statistical analysis, and QTL-detection Table 1 summarizes the type of traits that were recorded at maturity for each single plant (replicate) at the HUJI experimental farm and HTP platform (LemnaTec-129 Scanalyzer 3D) in Gatersleben. All statistical analyses were performed using R (version 3.4.1). The data of each trait, replication, year, and treatment were inspected separately. First, extreme outliers were filtered out, if they fell outside of an interval of plus or minus three times the standard deviation of the mean. Then, quantile-quantile (QQ) plots and density plots were used to evaluate whether the data were normally distributed. Data points which were not normally distributed at the QQ-plots residuals were removed manually. Subsequently, arithmetic means of the traits were calculated for each year/treatment and Shapiro-Wilk tests (Shapiro and Wilk, 1965) were carried out to confirm that data were normally distributed and suitable for further statistical analyses and tests. T-tests and ANOVA were applied to prove for each trait and/or year whether and to what extent they differ under both conditions. Descriptive statistics, i.e., density and correlation plots were calculated to analyze how the traits changed and correlate under the respective irrigation regimes. Finally, mean values of each trait, year, and treatment were used separately to calculate QTLs for each irrigation regime with the software MultiQTL2.6 (http://www. multiqtl.com). *15k iSelect markers that were used to develop PCR-based molecular markers are written in italics. Genotypic information that was obtained via GBS, is marked with "GBS" **. Physical positions were obtained from the wild emmer acc. Zavitan (v1) (Avni et al., 2017). Please note that no physical position was obtained for the marker Tdurum_contig27976_414. *** Genetic positions were obtained from the original mapping population (Fatiukha et al., 2021). **** Plants that showed recombination events between Kukri_c6227780 and Rac875_c2138_474 were subjected to GBS (please read material and methods section for more detailed information).

Marker development
A total of 10 PCR-based molecular markers were developed at a distance of 1-2 centimorgan (cM) along the QTL-interval of 15.67 cM (Supplementary Tables S1-3). Within the interval, six additional polymorphic regions were identified by applying GBS to a subset of plants that showed recombinations close to the marker Gene-1741_103, thus, in between Kukri_c6227780 and Rac875_c2138_474 ( Table 2). All markers share the same order in the genetic map and sequenced genome of the wild emmer acc. Zavitan (Table 2, Supplementary Table S1, Deblieck et al., 2020). For Tdurum_contig27976_414, no physical position could be annotated.

Phenotypic data of homozygous recombinants sub-NILs
A total of 600 BC 4 F 3 plants were derived from a cross of Uzan with NIL-U-2B-1 (Merchuk-Ovnat et al., 2016). Eighty-two heterozygous recombinant BC 4 F 3 plants were identified after genotyping with the flanking Tdurum_contig27976_414 (32.88 cM) and wsnp_Ex_c6537_1133876/IAAV980 (48.55) markers (Supplementary Table S1). Subsequently, a sample of six to eight descendent plants of each heterozygote recombinant BC 4 F 3 plant was screened and a total of 96 BC 4 F 4 homozygote recombinants were obtained. After further genotyping the BC 4 F 4 plants with the 10 markers (Supplementary Table S1,  Supplementary Table S2, Supplementary Table S3), 29 of them with representative recombinations along the QTL interval were selected and used for seed multiplication. Finally, 10 plants per genotype, i.e., 290 BC 4 F 5 plants, were subjected to phenotyping in 2017. Two very distinct groups of genotypes with different numbers of days from planting to heading (DPH) were observed (Figure 1). One group comprised 17 plants and required a mean time of 66.8 days for heading, whereas the other group comprised 12 plants and flowered 24.4 days later, thus, at 91.2 days. Grouping of these data was additionally confirmed with a Tukey's test (Tukey, 1949) and independent segregation of DPH was confirmed with a chi-square test (Pearson, 1900). A comparison of GBS data of these plants revealed that sub-NILs with delayed DPH-values, i.e., 1029B, 1115A, 1329A, and 1929C, share a single 669 kbp fragment on chromosome 2A ranging from 39980886 to 40649713 BP, which is absent in early flowering sub-NILs, e.g., 1004B, 1539A, 1121B, 1663G, 1766G, 1761D, 1735F, 1336E, 1145C, and Uzan (Supplementary Table S4). As indicated previously, this region harbors the photoperiod sensitive gene PpdA-1 from the wild emmer parent (Takenaka and Kawahara, 2012), ranging from 40487317-40489398 bp (Deblieck et al., 2020). Remarkably, this huge effect was not observed in the original mapping population (Fatiukha et al., 2021) although Langdon and Uzan share the same GBS-marker alleles along the PpdA-1 locus (Supplementary Table S4).
DPH correlates with grain yield (GY) and many important yield-related traits and might increase/decrease the respective Frontiers in Genetics frontiersin.org 06 QTL-LOD scores (Fatiukha et al., 2021). We, therefore, decided to exclude plants with increased DPH values from the phenotypic experiments in 2018 and 2019 ( Figure 1). In addition to this, two sub-NILs did not deliver enough phenotypic data in 2018. For these genotypes, less than three replicates per genotype/ treatment and year were available, which was considered unreliable and not representative, leading to a final set of 17 segmental sub-NILs. For each of these sub-NIL, 10 plants were subjected to phenotyping, i.e., phenotypic data were resolved for 170 plants in 2017, 2018, and 2019 (Table 2).
Almost all traits showed a normal distribution in both environments and all 3 years (Supplementary Table S5,  Figure S4). However, in a few cases such as the main spike length (MSpL), the number of spikes per plant (Sppp), and grain yield per spike (GYpSp), phenotypic data were not normally distributed, and therefore, (log-) transformed (Hackett, 1997). These data are clearly highlighted, e.g., in Supplementary Table S5. ANOVA revealed that the means of all traits, except for the number of main spike spikelets (MSpSp), were different between the 3 years (Supplementary Table S6).
While in 2017 the drought stress effect was mild but significant for most of the traits except for culm length (CL), thousand kernel weight of the main spike (MSpTKW), and   (Figure 3). In addition to this, HI correlates more negatively with CL and slightly more positively with CKN under WW conditions. These drought stress effects were observed independently in 2018 and 2019. Genotypic and phenotypic data of all years and both water regimes are given in Supplementary Table S7.

Sub-QTL-detection
Sub-QTLs under well-watered and water-limited conditions were calculated using data from 2017, 2018, and 2019, respectively. As mentioned previously, CL, MSpTKW, and HI did not show significant differences between drought and control conditions in 2017. Therefore, the data of these traits were not considered to calculate drought-specific QTLs. Table 3 shows some general information about the QTLs obtained, while Table 4 shows the logarithm of odds (LOD) values and dominance (d) effects of significant QTLs (p < 0.05) along the whole target interval. While QTLs for CL and TKW could be detected in both environments, QTLs for CKN or HI appear under drought stress or control conditions, respectively.
The QTLs for HI, CKN, and CL co-localize at the very distal parts of the QTL interval, leading to a biased set of segmental sub-NILs (i.e., northern or southern recombinants) which lack or harbor the respective part of the QTL (Table 2,  Table 4 and Figure 4). Sub-NILs from the upper part of the wild emmer QTL tend to have shorter culms, an increased number of kernels, and slightly higher GY than sub-Nils than the lower part of the QTL (Table2; Figure 4).
Furthermore, different QTLs related to traits of the main spike, e.g., for MSpSP, MSpTKW, and MSpL, were detected under WL and WW conditions between the markers Frontiers in Genetics frontiersin.org 08    LOD, logarithm of odds; d, sign of the QTL, substitution effect d represents the effect of G18-16 on the respective trait. It is the difference of the mean of homozygote sub-NILs, with the G18-16 and Uzan allele, respectively, CKN, calculated kernel number; GY, grain yield; HI, harvest index; MSPL, main spike length; TKW, thousand kernel weight; P, Probability (based on a permutation test with 1,000 repeats); PEV, Percentage of explained variance. No physical position was obtained for the marker Tdurum_contig27976_414. *** Genetic positions were obtained from the original mapping population (Fatiukha et al., 2021). Significant QTLs, with positive or negative effects on the trait and a LOD >2.0 are colored in yellow and blue.
Frontiers in Genetics frontiersin.org Gene_1741_1 and Tdurum_contig_68806_677 (Table 4). Sub-NILs with the wild emmer allele of the marker Tdurum_contig_68806_677 appear to have a smaller main ear with fewer main ear spikelets (MSpSp) than the elite parent Uzan ( Figure 5). Differentiating again between those sub-NILs which carry the upper part of the QTL interval and the G18-16 allele of the marker Tdurum_contig27976_414, but not the G18-16 allele of the flanking marker Gene-1741_103, revealed that these lines (i.e., 1663G, 1688C, 1488A, 1767E, and 1174B) again have a significantly higher HI, GY, and more CKN under WW conditions in all 3 years and under mild drought stress conditions in 2017 ( Figure 6 and Table 4, Table 2). The

FIGURE 4
Grain yield (GY), culm length (CL), calculated kernel number (CKN), and main spike length (MSpL) of segmental sub-NILs with G18-16 or Uzan alleles of the iSelect marker Tdurum_contig27976_414. A and B Alleles were derived from G18-16 or Uzan, respectively. Values obtained for the elite parent Uzan (U) were illustrated separately. Significant differences between groups of genotypes with the different alleles of Tdurum_contig27976_414 were calculated with a two-sided t-test. p-values ≤ 0.05 or 0.01 were marked with one or two stars, respectively.

Discussion and conclusion
In general, the following observations should be considered when comparing the results of this article to the previous studies by Peleg et al. (2009), Merchuk-Ovnat et al. (2016 and Fatiukha et al. (2021). First, it should be noted that the genetic background Main spike length (MSpL), main spike spikelets (MSpSp), and main spike seeds per spikelet (MSpSpSp) of segmental sub-NILs, carrying different alleles of the iSelect marker Tdurum_contig_68806_677. A and B Alleles were derived from G18-16 or Uzan, respectively. Values obtained for the elite parent Uzan (Uz) were illustrated separately. Significant differences between groups of genotypes with the different alleles of Tdurum_contig68806_677 were calculated with a two-sided t-test. p-values ≤ 0.05 or 0.01 were marked with one or two stars, respectively.
Frontiers in Genetics frontiersin.org of the NIL-U2B-1 comes from cv. Uzan and not from cv. Langdon, the original parent of the mapping population (Peleg et al., 2009;Fatiukha et al., 2021). Although the effect of this QTL in NIL-U-2B-1 was confirmed previously (Merchuk-Ovnat et al., 2016), the different genomic backgrounds of Uzan against which the effect of the QTL region was determined in this study might have an impact on the detected QTLs. Second, the drought stress conditions that were used in 2018 and 2019 were probably more severe than in previous studies. Peleg et al. (2009) Figure C shows the GY values of each of the sub-NIL that belongs to T2. Significsant differences between the T1, T2, and T3 groups of genotypes were calculated with a two-sided t-test. p-values ≤ 0.05 or 0.01 were marked with one or two stars, respectively.
Frontiers in Genetics frontiersin.org 15   Figure S4). This may be due to environmental factors, e.g., the size of the pots and light intensity. However, as mentioned in the results section, the drought stress effect in  Table S5). In contrast to the previous studies by Peleg et al. (2009) andFatiukha et al. (2020), no QTL for GY could be confirmed, but different segments of the introgressed wild emmer wheat region on chromosome 2B of NIL-U-2B-1 showed impacts on CL, CKN, and GY. The sub-NILs with the upper part of the emmer wheat QTL-interval had a shorter CL, more CKN, and a slightly higher GY, while sub-NILs with an introgression at the lower part of the QTL had longer stems, less CKN, and reduced GY ( Figure 4; Table 4, Table 5). Furthermore, the central part of the QTL interval has a negative impact on MSPL and the number of MSPSP. Removing sub-NILs with this MSPL-QTL from the set of sub-NILs that harbor the upper part of the QTL interval leads to five sub-NILs (i.e., sub-NIL 1663G, 1688C, 1488A, 1767E, and 1174B) with significantly (p < 0.001) higher HI, CKN, and GY than Uzan (Table2, Table 5, Figure 6). The increased number of kernels seems to have a decisive effect on the GY of these five lines. This effect was significant under mild drought stress in 2017 and controlled conditions in 2017 and 2019, but not under more severe drought stress conditions in 2018 and 2019. It can, therefore, be stated that the upper part of the introgressed fragment on chromosome 2B in the sub-NILs 1663G, 1688C, 1488A, 1767E, and 1174B might be of value for future breeding programs. The observed increased GY under mild drought stress conditions in 2017 is consistent with the results of previous studies (Peleg et al., 2009;Merchuk-Ovnat et al., 2016;Fatiukha et al., 2021). It might, therefore, be hypothesized that this region has a yield stabilizing effect under mild drought stress conditions. However, further trials with these sub-NILs would be useful to characterize up to which degree of drought stress the effect of this introgression remains advantageous.
The identification of candidate genes from the wild emmer parent that might be responsible for the increased kernel number remains difficult because a large fragment of about 50 million bp (Mbp) was introgressed upstream of Tdurum_contig27976_414 in NIL-U-2B-1 (Table 2, Deblieck et al., 2020). In addition to that, it is not clear which trait(s) are associated with the increased CKN of these plants. Plants with the upper part of the QTL region have an increased MSPL and an increased Sppp (T2 in Table 5). These differences are not significant (Table 5), and therefore, any further conclusions at this point of time remain speculative. However, it might be worth mentioning the presence of the Ppd-B1 gene from wild emmer on chromosome 2B in NIL-U-2B-1, which has previously been demonstrated to specifically influence the number of seeds per spikelet and not DPH in durum wheat (Arjona et al., 2018). Alternative QTLs in the same region, but close to Ppd-B1, have also been identified in wheat in previous studies (Gao et al., 2015;Shi et al., 2017).
Referring to the QTL interval that might be related to MSpL, MSpSp, or MSpSpSp, 111 high-confidence genes were annotated between Gene-1741_103 and Tdurum_contig_68806_677 (Table 2,  Table 4, Supplementary Table S9). One of these genes, the ethylene responsive factor (ERF) (TRIDC2BG016990), is a homolog of the well-analyzed Frizzy panicle (FZP) gene (LOC4344233) or Branched Silkless gene (BD1) in rice and maize, respectively (Colombo et al., 1998;Komatsu et al., 2001;Komatsu et al., 2003;Wang et al., 2020;Li et al., 2021). Wheat fzp lines were recently shown to have an increased number of spikelets, longer spike length, and reduced TKW (Li et al., 2021). Interestingly, these traits can be attributed to the elite parent allele (Table 4). It might, therefore, be hypothesized that a favorable FZP allele was selected in the course of evolution and introgressed into cultivars such as Uzan. Additional experiments, such as sequence and/ or expression analysis, are required to prove or reject this hypothesis.
Finally, a CL-QTL was detected at the lower part of the QTL interval (Table 4). In accordance with this result, Zanke et al. (2014) confirmed that the iSelect marker wsnp_Ex_c6537_11338763 and markers downstream of it are associated with an increased plant height. The authors mention the existence of a gene in wheat that is orthologous to the GIBBERELLIN INSENSITIVE DWARF1 (GID1)like receptor in rice (Zanke et al., 2014). This region was transferred from G18-16 into the NIL (Deblieck et al., 2020) and might have an impact on the detected CL-QTL. Remarkably, Merchuk-Ovnat et al. (2016) described that an alternative NIL, NIL-U-2B-3, carrying a smaller introgression of G18-16 on chromosome 2B than NIL-U-2B-1, downstream of Ku_c7740_879 (Supplementary Table S1, Merchuck-Ovnat et al., 2016) did not show such a strong increase in culm length relative to Uzan (Supplementary Table S5, Merchuck-Ovnat et al., 2016). Aligning Ku_c7740_879 to the reference genome of Zavitan (Avni et al., 2017) revealed that the marker is located at Frontiers in Genetics frontiersin.org