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

Abstract

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 BC₄F₅ 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 (Chen et al., 2015; 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.

Materials and methods

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₃F₃ isogenic line (NIL) NIL-U-2B-1 as previously described (Merchuk-Ovnat et al., 2016; Deblieck et al., 2020). In this study, BC₃F₃ NIL-U-2B-1 (pollinator) was crossed with Uzan (Supplementary Figure S1) and BC₄F₃ plants were screened with the respective flanking molecular markers to identify heterozygous recombinants. The heterozygous recombinant BC₄F₃ plants were selfed and BC₄F₄ single seed descendants were again genotyped to identify homozygous recombinant plants for seed multiplication. Finally, BC₄F₅ 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₄F₃ sub-NILs. All molecular markers were tested and validated on the DNA of F₇ 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 fine-mapping 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 high-quality 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₂O₅, and 10% K₂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 pre-drought 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).

Results

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.

TABLE 2

Genotypic data of the recombinant inbred lines.

*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).

Phenotypic data of homozygous recombinants sub-NILs

A total of 600 BC₄F₃ plants were derived from a cross of Uzan with NIL-U-2B-1 (Merchuk-Ovnat et al., 2016). Eighty-two heterozygous recombinant BC₄F₃ 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₄F₃ plant was screened and a total of 96 BC₄F₄ homozygote recombinants were obtained. After further genotyping the BC₄F₄ 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₄F₅ 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).

FIGURE 1

DPH correlates with grain yield (GY) and many important yield-related traits and might increase/decrease the respective 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 2 and Supplementary 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).

FIGURE 2

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 harvest index (HI), the effect of drought stress was very clear for all traits in 2018 and 2019 (Figure 2, Supplementary Figure S4, Supplementary Table S5).

In general, plants under drought stress have less GY, less MSpTKW, a reduced calculated kernel number (CKN), and shorter culms than under WW conditions (Figure 2 and Supplementary Figure S4, Supplementary Table S5, Supplementary Table S6). Traits such as CKN, MSpSp, MSpL, and the number of spikes per plant (Sppp) correlate most positively with GY under WW and WL conditions (Figure 3). In addition to this, HI correlates more negatively with CL and slightly more positively with CKN under WW conditions.

FIGURE 3

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).

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 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 average GY of these lines under WW conditions was 12.05 g (2017), 9.48 g (2018), and 12.07 g (2019), whereas Uzan yielded 10.60 g, 7.92 g, or 6.92 g (Table 5). The average GY of these lines under WW conditions in 2017, 2018, and 2019 was thus 12, 16, and 42% higher than of Uzan, respectively (Table 5).

FIGURE 5

FIGURE 6

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 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) used 350 mm and Merchuk-Ovnat et al. (2016) used 360–290 mm of seasonal water application to mimic drought stress. In this study, 350 mm of seasonal water application was used only once in 2017 resulting in mild drought stress symptoms. About 201 mm of seasonal water application was used in 2018, and in 2019, water application was reduced from 64 DAS to 20% PAW to simulate severe drought stress conditions (Dhanagond et al., 2019). Average GY was significantly reduced under drought stress to 10.39 g, 2.75 g, and 4.33 g per plant in 2017, 2018, and 2019, respectively. Lower GY in 2018 and 2019 is very likely a consequence of increased drought stress. Furthermore, ANOVA analysis revealed that the means of all traits were different between all 3 years (Supplementary Table S6). This is probably due to different reasons. On the one hand, the temperatures in the screen houses in 2018 were very high, i.e., above 45°C (Supplementary Table S8). These high temperatures suppressed tillering and induced early heading so that the plants were shorter, and yield was lower than in 2017. On the other hand, the experimental setup in 2019 differed significantly from 2017 to 2018. Plants in the pot experiments at the HTP in 2019 appeared to be smaller, to have reduced biomass, and less TKW (Supplementary Table S6, Supplementary 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 2017 was mild but significant for most of the traits, except for CL, MSpTKW, and HI, and very clear for all traits in 2018 and 2019 (Figure 2, Supplementary Figure S4, Supplementary Table S5).

In contrast to the previous studies by Peleg et al. (2009) and Fatiukha 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 164,96 Mbp on chromosome 2BS. We, therefore, hypothesize that the CL-QTL has a size of approximately 26,87 Mbp and ranges from 138.09 (Table 2) to 164,96 Mbp (Merchuck-Ovat et al. .,2016) on chromosome 2BS. This interval contains about 201 annotated high-confidence genes and was transferred into NIL-U-2B-1, but not into NIL-U-2B-3 (Merchuck-Ovnat et al., 2016; Avni et al., 2017). Aligning the rice GID1 gene (LOC4338764) against chromosome 2B of the reference genome of Zavitan (Avni et al., 2017), revealed a homolog candidate gene (TRIDC2BG021600) at position 146.52 Mbp within the interval. However, further fine-mapping of the introgressed G18-16 fragment that flanks wsnp_Ex_c6537_11338763 is required to further narrow down this QTL and to confirm or confute this candidate gene. In addition to this, it should also be taken into consideration that the reduced height gene 4 (Rht4) was described to be located on the long arm of chromosome 2B (Zanke et al., 2014; Wu et al., 2021). Aligning the primer sequences of the most significantly linked SSR marker WMC317 to the reference genome of Zavitan (Avni et al., 2017) revealed that it is located at 762 Mbp and was not transferred from G18-16 into NIL-U-2B-1 (Deblieck et al., 2020). Therefore, the observed differences in plant height cannot be attributed to Rht4.

The main tasks of this study were to identify the smallest sub-QTL region affecting grain yield and to diminish the effects of additional introgressions from the wild emmer parent in NIL-U-2B-1. Both goals were achieved. On the one hand, the delayed DPH and increased CL of NIL-U-2B-1 could be attributed to trans- and linkage drag effects of G18-16 introgressions on chromosome 2A or 2B, respectively. On the other hand, a GY and kernel number increasing G18-16 introgression was identified upstream of Tdurum_contig27976_414. The effect of this on the 50 million Mbp region was confirmed under controlled conditions in all years and under mild drought stress in 2017 (Table 5; Figure 6). Unfortunately, it flanks the QTL interval that was subjected to fine mapping and could not be further narrowed down. However, since the size of the region is much smaller than the original introgression in NIL-U-2B-1 (>400 Mbp) (Deblieck et al., 2020) and since the five selected sub-NILs with this fragment (sub-NILs 1663G, 1688C, 1488A, 1767E, and 1174B) do not show delayed DPH or increased CL, these genotypes might already be crossed with new elite cultivars to establish new NILs. In this respect, the effect of wild emmer introgression might be analyzed in different genomic backgrounds, environments, and years.

In addition to this, it might be useful to further examine the locus in the background of the cv. Uzan. A new fine-mapping approach might, therefore, be accomplished by again backcrossing 1663G, 1688C, 1488A, 1767E, or 1174B with the recurrent parent. Since the region upstream of Tdurum_contig27976_414 is located almost at the end of the chromosome (Table 2; Deblieck et al., 2020), recombinations should appear frequently in this region and fine mapping should be feasible. Finally, repeated trials with different levels of drought stress, Uzan, and the five sub-NILs might reveal to what extent of drought stress the actual yield-increasing effect of this wild emmer fragment persists.

References

• 1

  ArjonaJ. M.RoyoC.DreisigackerS.AmmarK.VillegasD. (2018). Effect of ppd-A1 and Ppd-B1 allelic variants on grain number and thousand kernel weight of durum wheat and their impact on final grain yield. Front. Plant Sci.9, 888. 10.3389/fpls.2018.00888

  • CrossRef
  • Google Scholar

• 2

  AvniR.NaveM.BaradO.BaruchK.TwardziokS. O.GundlachH.et al (2017). Wild emmer genome architecture and diversity elucidate wheat evolution and domestication. Science357, 93–97. 10.1126/science.aan0032

  • CrossRef
  • Google Scholar

• 3

  BelowR.Grover-KopecE.DilleyM. (2007). Documenting drought-related disasters. J. Environ. Dev.16, 328–344. 10.1177/1070496507306222

  • CrossRef
  • Google Scholar

• 4

  BrowningB. L.BrowningS. R. (2016). Genotype imputation with millions of reference samples. Am. J. Hum. Genet.98, 116–126. 10.1016/j.ajhg.2015.11.020

  • CrossRef
  • Google Scholar

• 5

  ChenD.WangS.CaoB.CaoD.LengG.LiH.et al (2015). Genotypic variation in growth and physiological response to drought stress and Re-watering reveals the critical role of recovery in drought adaptation in maize seedlings. Front. Plant Sci.6, 1241. 10.3389/fpls.2015.01241

  • CrossRef
  • Google Scholar

• 6

  ColomboL.MarzianiG.MasieroS.WittichP. E.SchmidtR. J.GorlaM. S.et al (1998). BRANCHED SILKLESSmediates the transition from spikelet to floral meristem duringZea maysear development. Plant J.16, 355–363. 10.1046/j.1365-313x.1998.00300.x

  • CrossRef
  • Google Scholar

• 7

  DeblieckM.FatiukhaA.GrundmanN.Merchuk-OvnatL.SarangaY.KrugmanT.et al (2020). GenoTypeMapper: Graphical genotyping on genetic and sequence-based maps. Plant Methods16, 123. 10.1186/s13007-020-00665-7

  • CrossRef
  • Google Scholar

• 8

  DempewolfH.BauteG.AndersonJ.KilianB.SmithC.GuarinoL. (2017). Past and future use of wild relatives in crop breeding. Crop Sci.57, 1070–1082. 10.2135/cropsci2016.10.0885

  • CrossRef
  • Google Scholar

• 9

  DhanagondS.LiuG.ZhaoY.ChenD.GriecoM.ReifJ.et al (2019). Non-invasive phenotyping reveals genomic regions involved in pre-anthesis drought tolerance and recovery in spring barley. Front. Plant Sci.10, 1307. 10.3389/fpls.2019.01307

  • CrossRef
  • Google Scholar

• 10

  ElshireR. J.GlaubitzJ. C.SunQ.PolandJ. A.KawamotoK.BucklerE. S.et al (2011). A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS One6, e19379. 10.1371/journal.pone.0019379

  • CrossRef
  • Google Scholar

• 11

  FahadS.BajwaA. A.NazirU.AnjumS. A.FarooqA.ZohaibA.et al (2017). Crop production under drought and heat stress: Plant responses and management options. Front. Plant Sci.8, 1147. 10.3389/fpls.2017.01147

  • CrossRef
  • Google Scholar

• 12

  FangYujieXiongLizhong (2015). General mechanisms of drought response and their application in drought resistance improvement in plants. Cell. Mol. Life Sci.72 (4), 673–689. 10.1007/s00018-014-1767-0

  • CrossRef
  • Google Scholar

• 13

  FarooqM.HussainM.SiddiqueK. H. M. (2014). Drought stress in wheat during flowering and grain-filling periods. Crit. Rev. Plant Sci.33, 331–349. 10.1080/07352689.2014.875291

  • CrossRef
  • Google Scholar

• 14

  FAO (2021). World Food and Agriculture-Statistical Yearbook 2021.

  • Google Scholar

• 15

  FatiukhaA.DeblieckM.KlymiukV.Merchuk-OvnatL.PelegZ.OrdonF.et al (2021). Genomic architecture of phenotypic plasticity in response to water stress in tetraploid wheat. Int. J. Mol. Sci.22, 1723. 10.3390/ijms22041723

  • CrossRef
  • Google Scholar

• 16

  GaoF.WenW.LiuJ.RasheedA.YinG.XiaX.et al (2015). Genome-wide linkage mapping of QTL for yield components, plant height and yield-related physiological traits in the Chinese wheat cross zhou 8425B/Chinese spring. Front. Plant Sci.6, 1099. 10.3389/fpls.2015.01099

  • CrossRef
  • Google Scholar

• 17

  GasparT.FranckT.BisbisB.KeversC.JouveL.HausmanJ. F.et al (2002). Plant Growth Regul.37, 263–285. 10.1023/A:1020835304842

  • CrossRef
  • Google Scholar

• 18

  HackettC. A. (1997). Model diagnostics for fitting QTL models to trait and marker data by interval mapping. Heredity79, 319–328. 10.1038/hdy.1997.160

  • CrossRef
  • Google Scholar

• 19

  HuangL.RaatsD.SelaH.KlymiukV.LidzbarskyG.FengL.et al (2016). Evolution and adaptation of wild emmer wheat populations to biotic and abiotic stresses. Annu. Rev. Phytopathol.54, 279–301. 10.1146/annurev-phyto-080614-120254

  • CrossRef
  • Google Scholar

• 20

  KomatsuM.ChujoA.NagatoY.ShimamotoK.KyozukaJ. (2003). FRIZZY PANICLE is required to prevent the formation of axillary meristems and to establish floral meristem identity in rice spikelets. Development130, 3841–3850. 10.1242/dev.00564

  • CrossRef
  • Google Scholar

• 21

  KomatsuM.MaekawaM.ShimamotoK.KyozukaJ. (2001). The LAX1 and FRIZZY PANICLE 2 genes determine the inflorescence architecture of rice by controlling rachis-branch and spikelet development. Dev. Biol.231, 364–373. 10.1006/dbio.2000.9988

  • CrossRef
  • Google Scholar

• 22

  LaneD. R.CoffinD. P.LauenrothW. K. (2000). Changes in grassland canopy structure across a precipitation gradient. J. Veg. Sci.11, 359–368. 10.2307/3236628

  • CrossRef
  • Google Scholar

• 23

  LangridgeP.FleuryD. (2011). Making the most of 'omics' for crop breeding. Trends Biotechnol.29, 33–40. 10.1016/j.tibtech.2010.09.006

  • CrossRef
  • Google Scholar

• 24

  LeskC.RowhaniP.RamankuttyN. (2016). Influence of extreme weather disasters on global crop production. Nature529, 84–87. 10.1038/nature16467

  • CrossRef
  • Google Scholar

• 25

  LiH. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. Available at: https://www.semanticscholar.org/paper/Aligning-sequence-reads%2C-clone-sequences-and-with-Li/74574ee09030e8aadb48fa349eb9b054e2f95ceb.

  • Google Scholar

• 26

  LiH.HandsakerB.WysokerA.FennellT.RuanJ.HomerN.et al (2009). The sequence alignment/map format and SAMtools. Bioinformatics25, 2078–2079. 10.1093/bioinformatics/btp352

  • CrossRef
  • Google Scholar

• 27

  LiY.LiL.ZhaoM.GuoL.GuoX.ZhaoD.et al (2021). Wheat FRIZZY PANICLE activates VERNALIZATION1-A and HOMEOBOX4-A to regulate spike development in wheat. Plant Biotechnol. J.19, 1141–1154. 10.1111/pbi.13535

  • CrossRef
  • Google Scholar

• 28

  Merchuk-OvnatL.BarakV.FahimaT.OrdonF.LidzbarskyG. A.KrugmanT.et al (2016). Ancestral QTL alleles from wild emmer wheat improve drought resistance and productivity in modern wheat cultivars. Front. Plant Sci.7, 452. 10.3389/fpls.2016.00452

  • CrossRef
  • Google Scholar

• 29

  Merchuk-OvnatL.FahimaT.EphrathJ. E.KrugmanT.SarangaY. (2017). Ancestral QTL alleles from wild emmer wheat enhance root development under drought in modern wheat. Front. Plant Sci.8, 703. 10.3389/fpls.2017.00703

  • CrossRef
  • Google Scholar

• 30

  NeumannK.KlukasC.FriedelS.RischbeckP.ChenD.EntzianA.et al (2015). Dissecting spatiotemporal biomass accumulation in barley under different water regimes using high-throughput image analysis. Plant Cell Environ.38 (10), 1980–1996. 10.1111/pce.12516

  • CrossRef
  • Google Scholar

• 31

  OyewoleC. I. (2016). The wheat crop. Anyigab: Kogi State University.

  • Google Scholar

• 32

  PearsonK. (1900). X. On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling. Lond. Edinb. Dublin Philosophical Mag. J. Sci.50, 157–175. 10.1080/14786440009463897

  • CrossRef
  • Google Scholar

• 33

  PelegZ.FahimaT.KrugmanT.AbboS.YakirD.KorolA. B.et al (2009). Genomic dissection of drought resistance in durum wheat x wild emmer wheat recombinant inbreed line population. Plant Cell Environ.32, 758–779. 10.1111/j.1365-3040.2009.01956.x

  • CrossRef
  • Google Scholar

• 34

  PelegZ.SarangaY.SuprunovaT.RoninY.RöderM. S.KilianA.et al (2008). High-density genetic map of durum wheat x wild emmer wheat based on SSR and DArT markers. Theor. Appl. Genet.117, 103–115. 10.1007/s00122-008-0756-9

  • CrossRef
  • Google Scholar

• 35

  PooleN.DonovanJ.ErensteinO. (2021). Viewpoint: Agri-nutrition research: Revisiting the contribution of maize and wheat to human nutrition and health. Food Policy100 (2021), 101976. 10.1016/j.foodpol.2020.101976

  • CrossRef
  • Google Scholar

• 36

  ShapiroS. S.WilkM. B. (1965). An analysis of variance test for normality (complete samples). Biometrika52, 591. 10.2307/2333709

  • CrossRef
  • Google Scholar

• 37

  ShiW.HaoC.ZhangY.ChengJ.ZhangZ.LiuJ.et al (2017). A combined association mapping and linkage analysis of kernel number per spike in common wheat (Triticum aestivum L.). Front. Plant Sci.8, 1412. 10.3389/fpls.2017.01412

  • CrossRef
  • Google Scholar

• 38

  SoleimaniB.LehnertH.KeilwagenJ.PlieskeJ.OrdonF.Naseri RadS.et al (2020). Comparison between core set selection methods using different Illumina marker platforms: A case study of assessment of diversity in wheat. Front. Plant Sci.11, 1040. 10.3389/fpls.2020.01040

  • CrossRef
  • Google Scholar

• 39

  TakenakaS.KawaharaT. (2012). Evolution and dispersal of emmer wheat (Triticum sp.) from novel haplotypes of Ppd-1 (photoperiod response) genes and their surrounding DNA sequences. Theor. Appl. Genet.125, 999–1014. 10.1007/s00122-012-1890-y

  • CrossRef
  • Google Scholar

• 40

  TukeyJ. W. (1949). Comparing individual means in the analysis of variance. Biometrics5, 99–114. 10.2307/3001913

  • CrossRef
  • Google Scholar

• 41

  United Nations, Department of Economic and Social Affairs, Population Division (2017). World Population Prospects: The 2017 Revision, Key Findings and Advance Tables. Working Paper No. ESA/P/WP/248.

  • Google Scholar

• 42

  WangS.-S.ChungC.-L.ChenK.-Y.ChenR.-K. (2020). A novel variation in the FRIZZLE PANICLE (FZP) gene promoter improves grain number and yield in rice. Genetics215, 243–252. 10.1534/genetics.119.302862

  • CrossRef
  • Google Scholar

• 43

  WilkinsonP. A.WinfieldM. O.BarkerG. L. A.TyrrellS.BianX.AllenA. M.et al (2016). CerealsDB 3.0: Expansion of resources and data integration. BMC Bioinforma.17, 256. 10.1186/s12859-016-1139-x

  • CrossRef
  • Google Scholar

• 44

  WuQ.ChenY.XieJ.DongL.WangZ.LuP.et al (2021). A 36 Mb terminal deletion of chromosome 2BL is responsible for a wheat semi-dwarf mutation. Crop J.9, 873–881. 10.1016/j.cj.2020.06.015

  • CrossRef
  • Google Scholar

• 45

  XieW.NevoE. (2008). Wild emmer: Genetic resources, gene mapping and potential for wheat improvement. Euphytica164, 603–614. 10.1007/s10681-008-9703-8

  • CrossRef
  • Google Scholar

• 46

  ZankeC. D.LingJ.PlieskeJ.KollersS.EbmeyerE.KorzunV.et al (2014). Whole genome association mapping of plant height in winter wheat (Triticum aestivum L.). PLoS One9, e113287. 10.1371/journal.pone.0113287

  • CrossRef
  • Google Scholar

• 47

  ZhaoC.LiuB.PiaoS.WangX.LobellD. B.HuangY.et al (2017). Temperature increase reduces global yields of major crops in four independent estimates. Proc. Natl. Acad. Sci. U. S. A.114, 9326–9331. 10.1073/pnas.1701762114

  • CrossRef
  • Google Scholar