Identification of Quantitative Trait Loci Hotspots Affecting Agronomic Traits and High-Throughput Vegetation Indices in Rainfed Wheat

Understanding the genetic basis of agronomic traits is essential for wheat breeding programs to develop new cultivars with enhanced grain yield under climate change conditions. The use of high-throughput phenotyping (HTP) technologies for the assessment of agronomic performance through drought-adaptive traits opens new possibilities in plant breeding. HTP together with a genome-wide association study (GWAS) mapping approach can be a useful method to dissect the genetic control of complex traits in wheat to enhance grain yield under drought stress. This study aimed to identify molecular markers associated with agronomic and remotely sensed vegetation index (VI)-related traits under rainfed conditions in bread wheat and to use an in silico candidate gene (CG) approach to search for upregulated CGs under abiotic stress. The plant material consisted of 170 landraces and 184 modern cultivars from the Mediterranean basin. The collection was phenotyped for agronomic and VI traits derived from multispectral images over 3 and 2 years, respectively. The GWAS identified 2,579 marker-trait associations (MTAs). The quantitative trait loci (QTL) overview index statistic detected 11 QTL hotspots involving more than one trait in at least 2 years. A CG analysis detected 12 CGs upregulated under abiotic stress in six QTL hotspots and 46 downregulated CGs in 10 QTL hotspots. The current study highlights the utility of VI to identify chromosome regions that contribute to yield and drought tolerance under rainfed Mediterranean conditions.

Introduction

Wheat (Triticum aestivum L.) is the most common cultivated crop worldwide. It is grown on 216 million hectares of land, and its global production of 765 million tons of grain provides 19% of calories and 21% of the protein in the human diet (Faostat 2019, http://www.fao.org/faostat). To cover the expected food demand of a world population that will increase up to 60% by 2050, wheat production needs to be increased by 1.7% per year (Leegood et al., 2010). Achieving this objective will not be easy, considering the expected negative effects of climate change on wheat yield, particularly in areas, such as the Mediterranean basin, where a rise in temperature by 3–5°C and a decrease in the annual rainfall by 25–30% have been predicted (Giorgi and Lionello, 2008). An increasing frequency and severity of terminal drought stress will reduce grain weight, grain quality, and wheat yield (Araus et al., 2002; Slafer et al., 2005; Kulkarni et al., 2017). Therefore, there is a need to improve the identification of genotypes that are able to maintain acceptable levels of yield and yield stability in semiarid environments, which have been identified as the most sensitive to the effects of climate change (Rufo et al., 2021). The release of improved cultivars with enhanced drought adaptation will be critical for breeding programs focusing on wheat adaptability and stability under rainfed conditions (Graziani et al., 2014; Bhatta et al., 2018).

The recent progress in high-throughput phenotyping (HTP) based on the use of multispectral images acquired from unmanned aerial vehicles (UAVs) has increasingly improved the assessment of agronomic traits (Gracia-Romero et al., 2017; Xie and Yang, 2020; Gomez-Candon et al., 2021; Rufo et al., 2021) on large germplasm collections in a rapid, cost-effective, and high spatial resolution way (Duan et al., 2017), as it allows for the estimation of various plant traits using non-intrusive and non-destructive technology (White et al., 2012; Rufo et al., 2021). Remote sensing has attracted growing interest in breeding programs since it can deliver detailed information about biophysical crop traits in many situations to cope with the current phenotyping bottleneck (Araus and Cairns, 2014; Juliana et al., 2019; Bellvert et al., 2021). Some studies have demonstrated the use of vegetation indices (VI) to indirectly detect wheat plants under water stress due to a decrease in vegetative growth (Condorelli et al., 2018). Others have demonstrated the use of energy balance models to estimate the actual water status (Gomez-Candon et al., 2021). When VIs are derived from multispectral cameras, they are obtained from the combination of wavelengths located at the visible, red-edge and near-infrared (NIR) regions of the light spectrum (Kyratzis et al., 2017). These wavelengths allow for discerning differences in vegetative greenness, rate of senescence, photosynthetic efficiency, and stay-green duration (Stenberg et al., 2004; Babar et al., 2006; Lopes and Reynolds, 2012). It has been stated that anthesis (A) and milk grain are the most suitable growth stages for the assessment of agronomic traits on a plot-by-plot basis (Aparicio et al., 2002; Royo et al., 2003). The use of HTP as a suitable and accurate predictor of agronomic traits, such as phenology, grain filling duration, biomass, and yield will provide unique opportunities to increase the power of quantitative trait loci (QTL) discovery by increasing the number of genotypes included in the analysis (Juliana et al., 2019). This method will increase the frequency of rare alleles of potential interest to improve wheat adaptation to various environmental conditions.

The dissection of the genetic and molecular basis of complex traits, such as yield and drought stress tolerance through complementary approaches, such as QTL mapping and genome-wide association studies (GWAS) or association mapping (AM) is essential in breeding programs. GWAS is based on linkage disequilibrium (LD; Flint-Garcia et al., 2003), and it is a powerful approach that provides high mapping resolution due to the higher recombination events analyzed in comparison with biparental mapping (Soriano et al., 2017; Qaseem et al., 2019). The AM has been used to identify genomic regions related to drought and heat tolerance in durum and bread wheat (Maccaferri et al., 2016; Valluru et al., 2017). Several studies have been conducted to investigate the genetic basis of grain yield and yield-related traits in bread wheat under water-stress conditions using AM (Edae et al., 2014; Gizaw et al., 2018; Qaseem et al., 2019; Mérida-García et al., 2020). The release of genome sequences for emmer wheat (Avni et al., 2017), bread wheat (IWGSC, 2018), and durum wheat (Maccaferri et al., 2019) and the availability of open databases of RNA sequencing (RNA-seq) experiments (Ramírez-González et al., 2018) have made it possible to use a candidate gene (CG) approach to find targets within QTL intervals without performing new functional studies.

The aim of the current study is to identify molecular markers linked to important agronomic traits, vegetation indices (VIs), and plant features related to drought resistance assessed by HTP, to define the most important QTL hotspots for such traits and to perform in silico detection of the underlying CG in those genomic regions.

Materials and Methods

Plant Material and Field Trials

A germplasm collection of 354 bread wheat (Triticum aestivum L.) genotypes from the MED6WHEAT IRTA panel described by Rufo et al. (2019) was used in this study, of which 170 corresponded to landraces and 184 to modern varieties collected and adapted to 24 and 19 Mediterranean countries, respectively (Supplementary Table 1). The panel is structured into six genetic subpopulations (SPs) and 38 genotypes that remained admixed (Rufo et al., 2019). SP1: west Mediterranean landraces (43 accessions); SP2: north Mediterranean landraces (59 accessions); SP3: east Mediterranean landraces (42 accessions); SP4: France–Italy modern germplasm (82 accessions); SP5: Balkan modern varieties (24 accessions); and SP6: CIMMYT-ICARDA-derived varieties (62 accessions).

The field trials were conducted at Gimenells, Lleida (41°38' N and 0°22' E, 260 m.a.s.l), northeastern Spain, under rainfed conditions for three consecutive seasons: 2016, 2017, and 2018. Average minimum and maximum monthly temperatures and rainfall were calculated from daily data recorded for a weather station close to the experimental fields. Climatic data (rainfall and temperature) for the last 15 years corresponding to the weather station in Gimenells, Lleida (Spain) were downloaded from https://ruralcat.gencat.cat/web/guest/agrometeo.estacions. Experiments followed a non-replicated augmented design with two replicated checks (cv. “Anza” and “Soissons”) at a ratio of 1:4 between checks and tested genotypes and in 3.6 m² plots with eight rows spaced 0.15 m apart. The sowing density was adjusted to 250 germinable seeds m⁻², and the sowing dates were December 2, 2015; November 21, 2016; and November 15, 2017, whereas harvesting dates were July 7, 2016; July 5, 2017; and July 5, 2018. Weeds and diseases were controlled following standard practices at the site.

Agronomic Data

The following traits were measured across the 3 years (2016, 2017, and 2018) according to the protocol described by Rufo et al. (2021): grain yield (GY, t ha⁻¹), number of spikes per m² (NSm²), number of grains per m² (NGm²), thousand kernel weight (TKW, g), aboveground biomass at physiological maturity (biomass, t DM ha⁻¹), harvest index (HI), plant height (PH, cm) and early vigor (estimated as green area, GA). The HI was calculated as the ratio between grain and plant weights in a 1-m long row sample. The PH was measured at maturity for three main stems per plot and was measured from the soil to the top of the spike, excluding the awns. Early vigor was calculated by integrating the green area (GA) values obtained by ground-based red, green, and blue (RGB) images taken every 14 days as described by Casadesús and Villegas (2014) from emergence until the detection of the first node. Finally, days from sowing to anthesis (A) (DSA, GS65) and grain filling duration (GFD, GS87) were measured on each plot based on the growth stage (GS) scale of Zadoks et al. (1974). The GSs were achieved when at least 50% of the plants in each plot reached the stage.

Image Acquisition

Image acquisition was conducted with a multispectral camera (Parrot Sequoia) (Parrot, Paris, France) installed onboard an UAV (DJI S800 EVO hexacopter, Nanshan, China). Images were acquired during 2 years, on April 21, 2017 and May 19, 2017 and on April 17, 2018 and May 18, 2018. Flights were always conducted at ~12:00 solar time under sunny conditions and with a wind speed below 12 m/s. The UAV flew at a height of 40 m above ground level (agl) and with a flight plan of 80/60 frontal and side overlap. The multispectral camera has four spectral bands located at wavelengths of 550 ± 40 nm (green), 660 ± 40 nm (red), 735 ± 10 nm (red edge), and 790 ± 40 (near infrared). The camera yields a resolution of 1,280 × 960 pixels. All images were radiometrically corrected through an external incident light sensor that measured the irradiance levels of light at the same bands as the sensor, as well as with in situ spectral measurements in ground calibration targets (black, white, soil, and grass). Spectral measurements were conducted with a Jaz spectrometer (Ocean Optics, Inc., Dunedin, FL, USA). Jaz has a wavelength response from 200 to 1,100 nm and an optical resolution of 0.3–10.0 nm. The calibration of the spectrometer measurements was taken using a reference panel (white color Spectralon™). Geometrical correction was conducted by using ground control points (GCPs) and measuring the position in each with a handheld global positioning system (GPS) (Geo7x, Trimble GeoExplorer series, Sunnyvale, CA, USA). All images were mosaicked using AgisoftPhotoscan Professional version 1.6.2 (Agisoft LLC., St. Petersburg, Russia) software and geometrically and radiometrically corrected with QGIS 3.2.0 (USA, http://www.qgis.org). Then, six spectral vegetation indices (VIs) were carefully selected based on their significance in relation to certain plant physiology features in wheat (Table 1). In addition, the leaf area index (LAI) was measured using a portable ceptometer (AccuPAR model LP-80, decagon devices Inc., Pullman, WA, USA) from 13:00 to 15:00 (local time) on each image acquisition date in 64 different plots of each set of landrace and modern set. Then, the LAI was estimated for each plot in the whole collection through the MTVI2, following the methodology described by Rufo et al. (2021) and Gomez-Candon et al. (2021). All VIs were assessed in 2017 and 2018 through UAV multispectral images at two growth stages: (1) when all the plots reached anthesis (A) (VI_{_}A) and (2) postanthesis (PA) at the milk and dough developmental stages (VI_PA).

TABLE 1

Table 1. Spectral vegetation indices assessed in this study.

Genotyping

The panel was genotyped with 13,177 SNP markers using the Illumina Infinium 15K Wheat SNP Array at Trait Genetics GmbH (Gatersleben, Germany), and 11,196 markers were ordered according to the SNP map developed by Wang et al. (2014). To reduce the risk of errors in further analyses, markers and accessions were analyzed for the presence of duplicated patterns and missing values. After excluding markers with more than 25% missing values and with a minor allele frequency (MAF) lower than 5%, a total of 10,090 SNPs were used for mapping purposes.

Statistical Analyses

Phenotypic data were fitted to a linear mixed model considering the check cultivars as the fixed effect, and the row and column number and accessions as random in the model for each environment following the MIXED procedure of the SAS-STAT statistical package (SAS Institute Inc., Cary, NC, USA).

$\begin{array}{l}y=x\beta +z\gamma +\epsilon \end{array}$

where β is an unknown vector of fixed-effects parameters with known design matrix x, γ is an unknown vector of random-effects parameters with known design matrix z, and ε is an unknown random error vector whose elements are no longer required to be independent and homogeneous. Restricted maximum likelihood (REML) was used to estimate the variance components and to produce the best linear unbiased predictors (BLUPs) for agronomic traits and VIs (Supplementary Table 2).

To assess the differences between years and genetic subpopulations, one-way ANOVAs were conducted for the whole collection. The broad sensed heritability (H²) was estimated following Knapp et al. (1985).

$\begin{array}{l}{H}^2=\frac{{\sigma }_G^2}{{\sigma }_G^2+{\sigma }_E^2+{\sigma }_{GE}^2}\end{array}$

where ${\sigma }_G^2$ is the genotypic variance, ${\sigma }_E^2$ is the variance due to the environmental (year) effect, and ${\sigma }_{GE}^2$ is the is the variance for the interaction of genotype with environment.

Least squares means were calculated and compared using the Tukey's HSD test. Pearson's correlation coefficients were calculated among the evaluated traits. Mean phenotypic values across the 3 years were used to perform a hierarchical cluster analysis by the Ward method (Ward, 1963). Analyses of variance and mean differences were carried out using the JMP v14.2.0 statistical package (SAS Institute, Inc., Cary, NC, USA), considering a significance level of alpha = 0.05.

Marker Trait Associations

A GWAS with 10,090 SNP markers was conducted on the whole germplasm collection using Tassel 5.0 software (Bradbury et al., 2007) for all agronomic and VI traits per year and across the three growing seasons. A mixed linear model (MLM) was fitted using a principal component analysis (PCA) matrix with six principal components as the fixed effect and a kinship (k) matrix as the random effect (PCA + K model) at the optimum compression level based on the groups defined by the kinship matrix. Compression levels range from “no compression” (compression = 1) when each genotype belongs to its own group, to “maximum compression” (compression = n) when all genotypes belong to the same group. In addition, the (A) date was incorporated as a cofactor in the analysis, as reported in previous studies (Crowell et al., 2016; Condorelli et al., 2018; Soriano et al., 2021). Manhattan plots were generated using the R script, CMplot (https://github.com/YinLiLin/CMplot). A false discovery rate (FDR) threshold (Benjamini and Hochberg, 1995) was established at –log₁₀ p > 4.8 (p < 0.05), using 3,696 markers according to the results of the linkage disequilibrium (LD) decay (Rufo et al., 2019). Besides, a frequently used threshold was established at –log₁₀ P > 3, as previously reported in the literature (Wang et al., 2014, 2020; Condorelli et al., 2018; Mangini et al., 2018; Sukumaran et al., 2018). Confidence intervals (CIs) for MTAs were estimated for each chromosome according to the LD decay reported by Rufo et al. (2019) using the formula reported by Chardon et al. (2004).

$\begin{array}{l}{S}_i^2={\left(\frac{CI}{3.92}\right)}^2\end{array}$

where CI corresponded with the LD decay for each chromosome. To simplify the MTA information, the associations were grouped into QTL hotspots. To define a hotspot, the density of MTAs along the chromosome was calculated as the QTL overview index (Chardon et al., 2004) for each cM of the genetic map reported by Wang et al. (2014).

$\begin{array}{l}U=\frac{\frac{nbQTL}{nbE}}{Totallengthofmap}\end{array}$

where nbQTL is the number of QTLs and nbE is the total number of experiments.

Gene Annotation and in silico Gene Expression Analysis

Gene annotation for the target region of QTL hotspots was performed using the gene models for high-confidence genes reported for the wheat genome sequence (IWGSC, 2018), available at https://wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations. Physical distances were estimated using the genetic distances from the markers flanking the CIs of each QTL hotspot.

In silico expression analysis and the identification of upregulated gene models were carried out using the RNA-seq data available at http://www.wheat-expression.com/ (Ramírez-González et al., 2018) for the following studies: (1) drought and heat stress time-course in seedlings, (2) spikes with water stress, and (3) seedlings with polyethylene glycol (PEG) to simulate drought.

Gene Ontology (GO) data were retrieved from the high-confidence gene annotation at https://wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations.

Results

Environmental Conditions

The experimental site has a typical Mediterranean climate characterized by an irregular pattern of yearly rainfall distribution, low temperatures in winter that rise sharply in spring, and high temperatures continuing until the end of the crop cycle. Figure 1 represents a graphical summary of the rainfall and maximum and minimum temperatures during the crop cycle across the 3 years of field trials and the average of the last 15 years. Although precipitation values were representative of long-term data from the region for each growing season, the year 2017 was considered exceptionally dry due to the low rainfall received. The year 2018 was characterized as the wettest from December (sowing) to June (physiological maturity) with 269 mm of rainfall, whereas the first and second growing seasons with 207 and 105 mm, respectively, of rainfall were rather dry, suffering severe water scarcity during the grain filling period with only 5 mm of precipitation.

FIGURE 1

Figure 1. Monthly rainfall (mm) and minimum (Tmin) and maximum (Tmax) temperatures during the growth cycle of each growing season.

Phenotypic Analyses

A summary of the genetic variation is shown in Table 2, Supplementary Table 3 for agronomic and vegetation indices (VIs)-related traits. Agronomic traits showed coefficients of variation (CV) ranging from 36.6% for grain yield to 6.3% for days to anthesis (A). VIs showed higher CVs during postanthesis (PA) with values at A ranging from 17.4% for leaf area index (LAI) to 2.1% for normalized difference vegetation index (NDVI), and PA ranging from 54.0% for LAI to 15.3% for green normalized difference vegetation index (GNDVI). Agronomic traits showed higher values for heritability than VI for most of the traits. For the agronomic traits, heritability ranged from 0.9 for yield to 0.1 for green area (GA), whereas for VIs, heritability ranged between 0.41 and 0.05 for TCARI/OSAVI and RDVI, respectively at A, and between 0.46 and 0.03 for TCARI/OSAVI and MTVI, 2 respectively during PA.

TABLE 2

Table 2. Summary statistics of the agronomic traits, leaf area index (LAI), and VIs.

The results of ANOVAs for the agronomic traits measured during the three growing seasons are shown in Table 3. The percentage of variability explained by year was the highest for GA (81.6%) and GS65 (67.9%), while the sum of squares of subpopulation (SP) was the highest for yield (76.5%), NGm² (65.0%), harvest index (HI) (62.6%), plant height (PH) (61.4%), and GS87 (59.0%). Finally, the highest percentage explained by the interaction between year and SP was found fo rbiomass, NSm² and thousand kernel weight (TKW), reaching 86.2, 84.9, and 71.0%, respectively. Significant differences were found between SPs for all traits. The year and the year × SP interactions were also significant for all traits, except for HI.

TABLE 3

Table 3. Analysis of variance for grain yield, harvest index (HI), biomass, number of spikes per square meter (NSm²), number of grains per square meter (NGm²), thousand kernel weight (TKW), plant height (PH), green area (GA), number of days from sowing to anthesis (A) (GS65), and grain filling duration (GFD, GS87) for the three years of field trials.

Table 4 shows the results of the ANOVA for the VIs and LAI estimated through the MTVI2 at the A and PA stages during 2017 and 2018. Differences between SPs and between years, as well as the year × SP interaction, were statistically significant for all traits in both years. The sum of squares of the year accounted for 1.3% (NDVI) to 92.9% (RDVI) of the variation at A, whereas at PA, the percentages ranged from 10.0% (TCARI/OSAVI) to 92.3% (LAI). The percentages of the total variation explained by SP ranged from 2.3% (RDVI) to 11.0% (TCARI/OSAVI) at A, while they ranged from 1.0% (MTVI2) to 11.2% (TCARI/OSAVI) PA. Year was the most important criterion for explaining the variations in LAI, RDVI, MSAVI, MTVI2, and GNDVI in the two growth stages. SP explained the least percent of variation at both growth stages for all traits. The year x SP interaction accounted for 4.8% (RDVI) to 89.9% (NDVI) of the model variance at A, with the highest values for NDVI and TCARI/OSAVI. The variance explained by the year × SP interaction at PA ranged from 6.6% (LAI) to 78.8% (TCARI/OSAVI).

TABLE 4

Table 4. Analyses of variance for the LAI estimated through MTVI2 and all the VIs calculated at the anthesis (A) and postanthesis (PA) stages in 2017 and 2018.

The mean values of phenotypic traits for each year and SP are shown in Table 5. Yearly means showed that the highest yield was in 2016, a year in which the yield components NSm², NGm², and TKW reached intermediate values between those obtained in the two subsequent years. The shortest duration of the preanthesis period and the longest grain filling duration (GFD) were also observed in 2016. On the other hand, the lowest yield, NSm² and NGm², the heaviest grains and the shortest GFD were observed in 2017. The GA reached the highest value in 2016, which was characterized as the wettest year during the period from January–March, i.e., the stem elongation stage, when the trait was measured. In contrast, 2017 was the driest year in the same period, which showed the lowest value for GA. In 2017, PH showed maximum values but biomass showed the lowest values at maturity. Finally, in 2018, biomass, the number of spikes, and grains per unit area showed high values, and the cycle until A was the longest.

TABLE 5

Table 5. Mean values of grain yield, harvest index (HI), biomass, number of spikes per unit area (NSm²), number of grains per unit area (NGm²), thousand kernel weight (TKW), plant height (PH), green area (GA), number of days from sowing to anthesis (A) (GS65), and grain filling duration (GFD, GS87) in a set of 170 landraces and 184 modern cultivars of bread wheat for each growing season and genetic subpopulation.

Significant differences in agronomic traits between SPs highlighted the division of the whole set into landraces and modern cultivars (Table 5). Modern SPs (SP4, SP5, and SP6) showed higher values of grain yield and yield components, HI, and biomass than landrace SPs. The highest value for grain yield was observed for SP4, in agreement with its higher number of spikes and grains per unit area. The SP4 showed the lowest grain weight among modern SPs but was not significantly different from the heaviest grains observed in landraces (SP1 and SP2). As expected, landraces were taller than modern cultivars. The SP3 showed the lowest value for GA. For phenology, SP2 took the longest time to reach the A stage, whereas SP6 took the shortest time. In contrast, the GFD was the shortest for SP2 and the longest for SP6. Modern SPs showed a longer GFD than landraces.

The mean values of the VIs and LAI (estimated by MTVI2) at A and at PA for 2017 and 2018 and the different SPs are shown in Table 6. All traits had higher values at A, except for TCARI/OSAVI. For all traits, differences between years were statistically significant at the two stages. The LAI, RDVI, and MSAVI showed the highest mean values in 2018. The mean values for TCARI/OSAVI were the highest in 2017. The year 2017 showed the highest values of MTVI2 and GNDVI at A, but these VIs and NDVI were minimal at PA the same year. Due to saturation of the reflectance, NDVI became insensitive at high LAI values (LAI > 3) in both years at A. LAI, NDVI, RDVI, MSAVI, and MTVI2 significantly differed between landrace and modern cultivar SPs at A, with higher values being recorded in the landraces. However, no pattern was found for VI traits among SPs PA. SP2 and SP4 had higher mean values for all traits PA, with the exception of TCARI/OSAVI.

TABLE 6

Table 6. Mean values of the LAI estimated by MTVI2 and all the VIs at anthesis (A) and postanthesis (PA) stages in 2017 and 2018 as well as for each genetic subpopulation and the group of admixed genotypes from a set of 170 landraces and 184 modern bread wheat cultivars.

Correlation coefficients between agronomic traits, VIs, and LAI were calculated (Figure 2), showing highly significant coefficients among agronomic traits as yield with NGm² (r = 0.90) and PH (r = −0.71). Interestingly, when analyzed VIs-related traits against agronomic traits, highly significant coefficients (r > 0.61) were found between GS65 and RDVI, MTVI2, GNDVI, LAI, NDVI_PA, and MSAVI_PA, and between GA and MSAVI, GNDVI, LAI, MTVI2, and RDVI_A. In addition, NGm² and PH showed a moderate significant correlation with GNDVI_PA and NDVI_A, respectively (r = 0.46).

FIGURE 2

Figure 2. Pearson correlations between agronomic traits, vegetation indices (VIs), and LAI+. GS65, number of days from sowing to A. GFD, grain filling duration; HI, harvest index; NSm2, number of spikes per square meter; NGm2, number of grains per square meter; TKW, thousand kernel weight; LAI, leaf area index; PH, plant height; GA, green area; A, anthesis; PA, postanthesis. Significant correlations at P < 0.0001 were established for r > 0.45 and r < −0.45.

To quantify the relation between trait variation and population structure, multiple linear regressions were carried out between population structure (q_i) coefficients (Rufo et al., 2019) (Table 7) and phenotypic performance for landrace and modern sets separately and both sets were combined. The landrace R² values ranged from 0.10 for MSAVI_A to 0.39 for GA, while the modern R² values ranged from 0.10 for MTVI2_A to 0.64 for GNDVI_A. When the regressions were conducted on the combined data set, the R² values ranged from 0.11 for biomass to 0.60 for NGm². The traits yield and GNDVI_PA showed high R² values (>0.35) for each set separately and for the combined set. The highest R² values were found in modern set regressions for GNDVI_A, GNDVI_PA, NDVI_PA, and GS65. Among the components of yield, thousand kernel weight (TKW) showed the highest R² values in landrace set regressions, while in modern set regressions, NGm² showed the highest R² values.

TABLE 7

Table 7. Relationship between trait variation and population structure (q-values) for landrace and modern sets separately and the combined set.

The bidimensional clustering shown in Figure 3 represents the relationships among accessions and their mean phenotypic performances (3 years for agronomic traits and 2 years for VIs). The horizontal cluster grouped accessions according to their phenotypic similarity based on the traits in the vertical cluster. Horizontal clustering separated two main clusters: Cluster A was composed only of landraces and Cluster B included modern cultivars and two landraces: cv “TRI 11548” from Iraq and cv “1170” from Turkey. Cluster A was characterized by lower yield and yield components, except NSm², lower biomass, a shorter GFD but longer GS65, and taller plants than Cluster B, but Cluster A had higher values for VIs at A except for GNDVI_A.

FIGURE 3

Figure 3. Bidimensional clustering showing the phenotypic relationships between the 354 bread wheat genotypes based on the analyzed traits indicated in the vertical cluster at bottom. Red and green colors in the columns indicate high and low values, respectively. Dark, higher values; light, lower values; white, intermediate values.

Each of these two clusters was separated into two subclusters, A1 and A2 for landraces and B1 and B2 for modern cultivars. Subcluster A1 was represented mainly by south Mediterranean landraces (77%), including those from the east and west regions, whereas Subcluster A2 contained most of the north Mediterranean landraces (62%). East Mediterranean landraces were in a single cluster within A1, whereas west Mediterranean landraces were distributed in other clusters within A1. Differences among subclusters, A1 and A2 were due to higher NSm² and TCARIOSAVI_PA in A1 and longer cycles until anthesis in A2, along with higher values for GA and VIs at PA. Regarding modern cultivars, Subcluster B1 was composed mainly of genotypes carrying the CIMMYT/ICARDA genetic background (SP6) (62%) and included the two landraces allocated to Cluster B mentioned above. Moreover, Subcluster B2 included 91% of the cultivars from SP4 (French and Italian modern cultivars). Most of the modern Balkan cultivars (SP5) were grouped in Subcluster B1. Subcluster B1 was characterized by higher TCARIOSAVI_PA and GA, whereas Subcluster B2 was characterized by higher NSm², GNDVI_A and the rest of the VIs assessed in PA.

Marker-Trait Associations

A summary of the results of the genome-wide association study (GWAS) for all traits per year and for the mean values across years is reported in Figure 4. Due to the low number of marker trait associations (MTAs) showing significance above a false discovery rate (FDR) threshold at –log₁₀ P > 4.8, a common threshold of –log₁₀ P > 3, as reported in the literature (Wang et al., 2017, 2020; Condorelli et al., 2018; Mangini et al., 2018; Sukumaran et al., 2018), reported a total of 2,579 MTAs (Supplementary Table 4). Manhattan plots for each of the traits and year are represented in Supplementary Figure 1. The year 2017 presented the highest number of MTAs, 74% of the total number of MTAs, whereas 2018 and the mean across years presented the lowest number of MTAs, 3 and 4%, respectively (Figure 4A). During 2016, only MTAs related to agronomic traits were reported, accounting for 19% of the total MTAs across years, as no multispectral images were captured during that year.

FIGURE 4

Figure 4. Summary of MTAs. (A) Percentage of MTAs per year and trait. (B) Number of MTAs per chromosome. (C) Phenotypic variance explained (PVE). MTAs, marker–trait associations; HI, harvest index; LAI, leaf area index estimated by MTVI2; NSm², number of spikes per square meter; NGm², number of grains per square meter; TKW, thousand kernel weight (g); PH, plant height; GA, green area from emergence until the first node; A, anthesis stage; PA, postanthesis.

The number of MTAs per chromosome for all years and for the mean values across years ranged from 9 on chromosome 4D to 354 on chromosome 1B (Figure 4B). Genome B accounted for 48% of the total MTAs, followed by genomes A and D with 41 and 11%, respectively. The percentage of MTAs with a phenotypic variance explained (PVE) lower than 0.10 was 97.5%, which agreed with the highly quantitative nature of the analyzed traits (Figure 4C).

A total of 815 MTAs were identified for seven agronomic traits (Supplementary Table 5). Yield showed the highest number of MTAs (368), most of them (268) from 2016, whereas only one association was found for NSm² with the mean across years. MTAs for TKW were found mainly during 2016 (96%), and those for PH were found mainly during 2017 (68%).

A total of 1,764 MTAs over –log₁₀ P > 3 were identified for 15 VI traits (Supplementary Table 5). Among them, 1,718 were detected a tor before A of green area (GA), and only 46 MTAs were identified at PA. Ninety-six percent of the MTAs were identified during 2017, which was the year characterized by the lowest rainfall. TCARIOSAVI_A was the trait with the highest number of MTAs (1,243), followed by MTVI2_A with 350.

To identify the genomic regions mostly involved in trait variation, QTL hotspots were identified using the QTL overview index defined by Chardon et al. (2004) for each cM of the genetic map reported by Wang et al. (2014). Confidence intervals (CIs) were calculated using the linkage disequilibrium (LD) decay for each chromosome reported by Rufo et al. (2019).

A total of 209 peaks were identified using the mean of the overview index across the 21 chromosomes (0.7) as the threshold, whereas using a high threshold (3.5), a total of 41 peaks were detected (Figure 5). These 41 peaks were reduced to 28 QTL hotspots (Supplementary Table 6), 12 in genomes A and B and 4 in genome D. To simplify the search for candidate genes (CGs), quantitative trait loci (QTL) hotspots were excluded when (1) the centromere was included within the hotspot or the CI was higher than 35 Mb and (2) MTAs corresponded only to 1 year of field experiments. Eleven QTL hotspots grouping 295 MTAs remained for subsequent analysis (Table 8). As shown in Figure 5, hotspots defined by the QTL overview index correspond to genome regions with a higher number of MTAs.

FIGURE 5

Figure 5. QTL overview index. The index values are represented along chromosomes as a blue line. Yellow and red dashed lines represent the thresholds for average (0.7) and higher values (3.5), respectively. Green bars below the QTL overview index represent the number of significant MTAs per 10 cM (–log₁₀ P > 3). QTL, quantitative trait loci; MTAs, marker–trait associations.

TABLE 8

Table 8. QTL hotspots identified for agronomic and remotely sensed VI-related traits.

In silico Analysis of CGs

A search for CGs to study the relative gene expression levels under abiotic stress conditions and different tissues and developmental stages was performed within the QTL hotspot regions reported in Table 8 using the positions of flanking markers in the “Chinese Spring” reference genome (IWGSC, 2018) at https://wheat-urgi.versailles.inra.fr/Tools/JBrowse. A total of 1,342 gene models were detected, and to classify this information, Gene Ontology (GO) for 1,025 of the gene models (76%) was downloaded from https://wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations (Figure 6; Supplementary Table 7). Seven hundred ninety-one CGs were classified according to molecular function (MF), 183 according to biological process (BP), and 51 according to cellular component (CC). The most represented CGs according to molecular function were “protein binding” (31%), “protein kinase activity” (13%), and “nucleic acid binding” (11%). According to BP, 30% of the CGs were involved in “defense response” and 19% in “transport.” Finally, according to CG, 27% of the product of CGs were in the nucleus, 22% in the membrane, and 14% in the cytoplasm and cell wall.

FIGURE 6

Figure 6. Gene Ontology (GO) classification of gene models within QTL hotspots. (A) GO hierarchy. (B) Molecular function. (C) Biological process. (D) Cellular component.

Subsequently, a search for differentially expressed genes (DEGs) under three abiotic stress conditions as reported in http://www.wheat-expression.com was carried out. These conditions included (1) drought and heat stress time-course in seedlings, (2) spikes with water stress, and (3) seedlings treated with polyethylene glycol (PEG) to simulate drought, and DEGs were analyzed in four tissues (roots, shoots/leaves, spikes, and grains) during different developmental phases (seedling, vegetative, and reproductive). A total of 12 CGs that were upregulated under abiotic stress were found in six QTL hotspots and 46 were found downregulated in 10 QTL hotspots (Figure 7).

FIGURE 7

Figure 7. Upregulated and downregulated CGs under abiotic stress conditions in four tissues and three developmental phases. Values are based on log₂ tpm. CGs, candidate genes; tpm, transcripts per million.

Among the different upregulated DEGs, a defensin in hotspot QTL1A.1 showed the highest expression under abiotic stress conditions and was expressed in most of the tissues and all the developmental phases; it also showed the highest expression levels for each of the phases. All DEGs reported expression in the spikes with a range from 0.02 tpm for cytochrome b in QTL1B.2 to 3.48 tpm for defensin in QTL1A.1. Only four DEGs were expressed in the roots and five in the leaves/shoots and grain. Only zinc finger protein-like 1 in QTL2A.2 was expressed in all four plant tissues, showing the highest expression in roots. Regarding the developmental phase, no expression was reported in any stage for two DEGs, cytochrome b in QTL1B.2 and enoyl-[acyl-carrier-protein] reductase in QTL3D.1. The reproductive phase had the highest number of DEGs (nine out of 10 showing expression), whereas 6 were expressed in the seedlings and only four were expressed during the vegetative phase.

Among the downregulated DEGs, the hotspot QTL1B.2 showed the highest number of downregulated DEGS (16), whereas QTL2A.2 did not show any of them. Three DEGs showed expression in all tissues, whereas any of them were downregulated in all of tissues. About six DEGs were expressed only in roots under non-stressed conditions, three in leaves/shoots, seven in spikes, and nine in grains, whereas DEGs non-expressed under abiotic stress in only one tissue corresponded to one in roots and six in grains. According to the developmental phase, 12 DEGs were expressed in all of them, whereas any DEG was downregulated in all of them. One DEG was expressed only in the seedling stage at normal conditions, whereas two were expressed only in the vegetative and 22 in the reproductive stages.

Discussion

The current study was conducted under typical Mediterranean environmental conditions, with a pattern of increasing temperatures during the spring and an irregular distribution of rainfall across years. A genome-wide association study (GWAS) panel of 354 bread wheat genotypes, including Mediterranean landraces and modern cultivars, was grown for 3 years under these conditions in northeastern Spain. Given that the decrease in the genetic diversity of wheat occurred during the second half of the twentieth century, associated with the introduction of high-yielding semidwarf cultivars (Autrique et al., 1996), landraces are considered a natural reservoir of genetic variation within the species and an invaluable source of new alleles to widen the genetic variability in breeding populations, particularly for traits regulating adaptation to suboptimal environments (Lopes et al., 2015). Recent studies have demonstrated the scarce use of wheat landraces in breeding programs in the past, as suggested by the high genetic diversity and defined population structure among landrace and modern cultivar subpopulations (Soriano et al., 2016; Rufo et al., 2019).

Phenotypic Performance

The high heritability reported for the agronomic traits, reaching 0.9 for yield and 0.8 for HI, NGm² and PH indicated that genetic differentiation among landraces and modern cultivars played a predominant role in determining the variation for these traits.

The ANOVAs showed a large effect of subpopulation (SP) on the phenotypic expression of the agronomic traits, whereas year showed the largest effect for most of the VIs, followed by the year × SP interaction, with the SP effect being the lowest. The variability in agronomic traits was mostly caused by the different agronomic performances of wheat landraces and modern cultivars, as reported in previous studies (Soriano et al., 2016; Royo et al., 2020). On the other hand, the high year effect on VIs was likely due to the contrasting water availabilities during the 2 years in which images were acquired by unmanned aerial vehicles (UAVs) in the experimental fields. This was not an unexpected result given that the decrease in the rate of growth of wheat caused by drought stress results in a severe reduction in total aboveground biomass (Royo et al., 2004; Casadesús and Villegas, 2014).

Yearly variation in weather conditions, particularly water input, resulted in a yield range from 7.4 t/ha in 2017, the driest year, to 7.9 and 8.0 t/ha during the years with higher rainfall. Even with the low water input, the average experimental yields were higher than expected in a severe drought environment. The number of spikes and grains per unit area were the highest in 2018, the wettest year, but were the lowest in 2017. Grain weight showed the opposite pattern, suggesting that under drier and hotter conditions, cultivars filled their grains at a higher rate (1.29 g/day in 2017 and 1.06 g/day in 2016 and 2018), thus showing a shorter grain filling duration (GFD) in 2017. The high yields recorded, considering the rainfed conditions of the field trials, could be attributed to the high soil fertility (~3% of organic matter) and the superficial subsoil water layer at this site (Royo et al., 2021).

From a genetic viewpoint, a clear separation was observed between landraces and modern cultivars for most of the agronomic traits, which can be attributed to the improvement achieved by breeding. As expected, the yield was negatively correlated to plant height (PH) as reported previously by Royo et al. (2020). Among landraces, those from northern Mediterranean countries characterized by high rainfall and lower temperatures (Royo et al., 2014) showed higher yields due to an increase in the number of grains per unit area and grain weight. These genotypes showed longer cycles until A and a shorter grain filling duration, although this last trait was not statistically significant. Landraces from the eastern Mediterranean countries showed lower yields, a lower number of grains, and lighter grains but an increase in the number of spikes per unit area compared with landraces from the northern Mediterranean countries. Similar results for the east Mediterranean landraces were previously reported in durum wheat by Soriano et al. (2018) and Roselló et al. (2019b), suggesting an adaptation of landraces from this area to warmer environments, which has been associated with the allelic constitution of vernalization and photoperiod genes (Royo et al., 2020). The results of the current study are in agreement with the previous research reporting a tendency for wheat to increase the number of ear-bearing tillers as an adaptation strategy under heat stress (Hütsch et al., 2019) and to increase the number of spikes per unit area in genotypes adapted to dry and warm areas compared to genotypes adapted to wetter and colder areas (Royo et al., 2014, 2020). Among modern cultivars, significant differences were mainly found between SP4 (cultivars from France and Italy) and the other two SPs (Balkans and CIMMYT-ICARDA-derived germplasm). These results suggest that breeding in France and Italy was in the direction of increasing yield by increasing the number of spikes and grains per unit area, whereas the other SPs showed higher thousand kernel weight (TKW). In addition, the regression results of the modern set suggested a high impact of genetic population structure on the number of grains per unit area. Cultivars derived from CIMMYT and ICARDA germplasms reached A earlier, up to 8 days earlier compared with Balkan cultivars and 6 days earlier compared to French and Italian cultivars, which was in line with the high R² values obtained in the relation between the modern set structure and GS65. This earliness can help these cultivars from warmer regions avoid heat stress at the end of flowering.

All traits related to high-throughput phenotyping (HTP) showed significant differences between years before and after anthesis (A), showing higher values for most of the VIs in 2018 than in the previous year. These higher values agree with the rainfall recorded for both years, which was significantly lower in 2017 than in 2018, mostly during the grain filling period. Furthermore, the difference in the mean values between growth stages was much higher in 2017. This result could be explained by the water scarcity particularly affecting the PA stage, which results in an important loss of chlorophyll content during the grain filling period; therefore, VIs using bands mostly placed in the near-infrared (NIR) and green regions showed lower values (Adamsen et al., 1999). It was supported by the high and positive correlations values between GA and GNDVI and LAI at post A stage, which was the case in 2018. Even though water stress affects the growth of wheat, the effects are higher during the grain filling duration (GFD) (Moragues et al., 2006). Thus, the leaf area index (LAI) and green normalized difference vegetation index (GNDVI) values decreased at the end of the growing cycle due to a low chlorophyll content associated with senescence during the grain filling period (Rufo et al., 2021). In addition, Gitelson et al. (2002) reported that the sensitivity of the green band was higher than that of the red band when the vegetation fraction was more than 60%, so vegetation indices using green wavelengths perform better at high LAI values, which in wheat under Mediterranean conditions are the highest at booting (Aparicio et al., 2000; Royo et al., 2004; Kyratzis et al., 2017; Rufo et al., 2021). This agreed with the high and positive correlation values between GS65 and GNDVI and LAI, indicating that more days until A provides a high green LAI at postanthesis (PA) stage in wet years as 2018. TCARI/OSAVI had higher PA for both years. This agreed with the results of Zarco-Tejada et al. (2005), who reported that in advanced growth stages, chlorophyll indices, such as TCARI performed better due to being less sensitive to the loss of turgor and leaf drop. In fact, these authors also stated that the different patterns of the indices across growth stages suggested that chlorophyll-related indices are more suitable closer to harvest, while structural indices related to canopy light scattering and growth are better for early stages.

Differences in the mean values of SPs were found in the two growth stages, with the highest values mainly found at A based on the differences among years, thus highlighting the effect of PA senescence on the chlorophyll content. Landraces and modern cultivars showed significant differences in the LAI and structural vegetation indices (VIs) at A, and these values were higher in the landraces. As reported in previous studies in durum wheat (García Del Moral et al., 2005; Soriano et al., 2018), landraces are characterized by their tolerance to water scarcity and their superior water use efficiency before A compared to modern cultivars (Subira et al., 2015). The SPs showing the highest mean values for the LAI and VIs at PA were those including landraces from the north of the Mediterranean basin (SP2) and modern cultivars from France and Italy (SP4). Landraces from SP2 are better adapted to colder and wetter environments than landraces originating in the southern part of the Mediterranean basin. This adaptation pattern has been associated with the greatest early soil coverage and more aboveground biomass along the whole cycle length (Royo et al., 2014, 2021). For this reason, the canopy remains green much longer in landraces from the northern Mediterranean countries than in those from the southern Mediterranean countries (Royo et al., 2014). The same pattern was found in modern cultivars, with GNDVI values remaining higher than those of landraces after A and being significantly different among modern subpopulations. These results agreed with those from the relationship between structure and GNDVI_PA, where the modern set showed the highest R² values according to the differences found in GNDVI mean values among modern SPs. These results and the capacity to discern between landrace and modern SPs regarding the VI values at A proved the accuracy of HTP in characterizing populations. Several studies have stated the potential of remote sensing for assessing agronomic traits by screening hundreds of plots in a short period of time, minimizing replications (Araus et al., 2018; Gracia-Romero et al., 2019; Juliana et al., 2019). Furthermore, various authors have stressed the suitability of using VIs measured early in the season for grain yield forecasting (Aparicio et al., 2000).

Bidimensional clustering was helpful to jointly visualize the results obtained by Tukey's tests. Moreover, clustering of agronomic and HTP data revealed similarity with the separation obtained by Rufo et al. (2019) using SNP markers and SPs defined based on the structured collection. In both the cases, a clear differentiation among landraces and modern cultivars was observed, which resulted in separation into two main clusters. Within the landrace cluster, A separation was observed between landraces from the northern and southern Mediterranean countries, thus including landraces from SP2 in one cluster and those from SP1 and SP3 in the other cluster, with different groupings among them. Modern cultivars of SP6 (CIMMYT-ICARDA) clustered separately from the French and Italian cultivars (SP4), whereas modern cultivars from the Balkans grouped mostly with SP6. Although these two SPs were separated genetically, no significant differences were found for the agronomic traits, except for phenology, and regarding the VIs, no differences were found at A. Two landraces (TRI 11548 and 1170) were included within modern cultivars from CIMMYT-ICARDA and the Balkans. These two landraces were characterized by a longer GFD, higher HI, and lower number of spikes per unit area than the average for landraces. Landrace TRI 11548 from Iraq also showed higher yield and grain weight than other landraces, so it probably resulted from a selection made in an early landrace population.

Marker Trait Associations

Dissecting the genetic basis of complex traits in plant breeding is essential to tackle molecular-based approaches for crop improvement. Several efforts have been previously made to identify quantitative trait loci (QTLs) and marker–trait associations (MTAs) associated with traits of interest to carry out marker-assisted selection (MAS) approaches and the introgression of alleles of commercial interest in adapted phenotypes.

The highest number of MTAs related to agronomic traits was found in 2016, while 96% of MTAs related to VIs, GA, and the LAI were identified in 2017. It has been reported that under contrasting conditions, the G × E interaction could affect the identification of stable associations among different environments (Mwadzingeni et al., 2017), which could explain the difference in the number of significant associations among the 3 years of field trials. The highest number of associations for yield and TKW in 2016 could be due to the moderate amount of water input (rainfall) during the spring, together with the longest grain filling duration, as reported in previous studies where grain weight predominantly enhanced yield in wet environments (García Del Moral et al., 2003; Moragues et al., 2006; Royo et al., 2006). Moreover, Royo et al. (2020) found that genotypes with longer GFDs could have greater opportunities to increase grain weight in favorable growing seasons than in warmer and drier seasons. The elevated number of VI-related MTAs found in the driest year (2017) could be explained by the higher variability in traits related to leaf biochemical properties or canopy structural attributes within the set of genotypes grown in environments with water scarcity (Rufo et al., 2021). The highest number of MTAs was identified for PH in 2017, when the coefficients of variation (CV) was higher for this trait. Qaseem et al. (2019) suggested that taller genotypes under drought stress could increase yield accumulation and convert more assimilates into grain. Of the 1,764 MTAs detected for VIs, 1,718 were found at A, with 1,243 for TCARI/OSAVI. This result could be explained by the significant differences found between landraces and modern SPs at A for traits related to HTP. The highest variability was found when comparing SP4 with the rest of the SPs for TCARI/OSAVI, which could explain the elevated number of MTAs for this trait. The distribution of the MTAs across genomes agreed with the results of Rufo et al. (2020), with a similar number in the A and B genomes (41 and 48%, respectively) and the remaining 11% in the D genome. These results are consistent with those of previous studies (Chao et al., 2010; Wang et al., 2014; Gao et al., 2015), which attributed these values to the lower genetic diversity and higher LD found in the D genome of bread wheat compared with genomes A and B (Rufo et al., 2019).

QTL Hotspots

To reduce the complexity of the high number of identified MTAs, QTL hotspots were defined using the QTL overview index proposed by Chardon et al. (2004). Although this statistic was initially used for classical biparental QTL analysis, we adapted it to GWAS using the confidence intervals (CIs) of the MTAs as the distance of LD decay for each of the chromosomes. As reported in Figure 5, QTL hotspots defined by the high-value threshold of the overview index corresponded to genome regions with a higher MTA density, thus supporting the suitability of this approach in GWAS. To identify genome regions previously mapped in locations similar to our QTL hotspots and to detect new loci controlling agronomic traits and VIs, a comparison with previous GWAS studies and/or meta-QTL analysis reporting yield- and VI-related traits was conducted. Seven of the 11 QTL hotspots have been described previously in the literature. When compared with the meta-QTL analysis reported by Liu et al. (2020) in bread wheat, the QTL hotspots QTL1B.2 and QTL2D.1 were located at similar positions as MQTL1B.7 and MQTL1B.8 and MQTL2D.3 and MQTL2D.4, respectively, controlling grain yield, grain number, and TKW under drought and heat stress. QTL1B.2 was also in the homologous region of QTL IWB50693 in durum wheat controlling spike length (Anuarbek et al., 2020), QSN.caas-1BL controlling NSm² identified by Gao et al. (2015) in bread wheat and IWB3330 controlling the normalized chlorophyll pigment ratio index (NCPI) identified by Gizaw et al. (2018) in bread wheat. Gao et al. (2015) also found three QTLs for TKW, chlorophyll content, and NDVI located in a common region with the hotspot QTL5B.2. This QTL hotspot was also detected in a similar region as QTL yield/root_5B.1 controlling grain yield and shoot length identified by Rufo et al. (2020). QTL1B.2 and QTL5B.2 were found to have homology with several studies. This was an expected result, since they were the longest hotspots including the highest number of MTAs. The genomic regions for QTL hotspots QTL5B.4 and QTL7A.1 were also found in common with three QTLs identified by Anuarbek et al. (2020) controlling the number of fertile spikes and TKW in durum wheat under rainfed conditions. Hotspots QTL1A.1 and QTL2A.2 shared a common position with mtaq-1A.2 reported by Roselló et al. (2019a) in durum wheat and QTL yield/root_2A2 identified by Rufo et al. (2020) in bread wheat, respectively, controlling root-related traits and grain yield in bread wheat.

The detection of these regions in common with other studies opens the opportunity to a pyramid of different QTLs with pleiotropic effects in future breeding approaches. Moreover, the use of the reference genome sequence makes it possible to rapidly identify common molecular markers to be used in MAS.

Candidate Genes

Gene annotation from the “Chinese spring” reference genome sequence (IWGSC, 2018) allowed us to identify 1,342 gene models within the 11 QTL hotspots. Candidate gene (CG) mining was performed by searching for differentially expressed genes (DEGs) upregulated and downregulated under drought conditions in different tissues and developmental stages through in silico analysis at http://www.wheat-expression.com.

Four CGs that were upregulated under drought stress have been previously reported in the literature to be involved in stress resistance. Among them, in QTL hotspot 1A.1, a defensin protein (TraesCS1A01G013600) was found to show the highest expression under drought stress. According to Kumar et al. (2019), although defensins are mainly involved in antifungal responses, the defensin gene Ca-AFP from chickpea in transgenic Arabidopsis plants was overexpressed under drought stress and induced a higher germination rate, root length, and plant biomass. Two gene models enhancing drought and heat stress tolerance were found in QTL hotspot 2A.2: the gene model TraesCS2A01G550300 encoding a zinc finger protein, as reported by (Yoon et al., 2014) in poplar, and the gene model TraesCS2A01G552400 encoding a MYB transcription factor, which was described by Zhao et al. (2018). TaMYB31 from wheat is transcriptionally induced by drought stress in transgenic Arabidopsis plants (Zhao et al., 2018). Finally, in QTL hotspot 5B.2, a squamosa-binding protein was identified (TraesCS5B01G286000). These protein families have been found to be involved in several biological processes. Cao et al. (2019), in expression studies of the Squamosa binding protein from wheat TaSPL16, found that this gene was highly expressed in young panicles but expressed at low levels in seeds, in agreement with the expression profile of TraesCS5B01G286000 found in our study. The ectopic expression of TaSPL16 in Arabidopsis produced a delay in the emergence of vegetative leaves and early flowering and affected yield-related traits. Other gene models upregulated under drought stress, as reported in the RNA-seq analysis from Ramírez-González et al. (2018), such as NADH-quinone oxido reductase, cytochrome b, histone deacetylase 2, RNA-binding protein, enoyl-[acyl-carrier-protein] reductase, double-stranded RNA binding protein 3, and histone-lysine N-methyltransferase, have not been related to drought stress tolerance in the literature, and further experiments are required to assess their expression under drought stress conditions. On the other side, 46 gene models were shown to be downregulated under drought stress. However, the decrease in the expression level seems to be more associated with the breakdown of the physiological functions due to drought than a causal effect in response to the stress.

Conclusions

The use of local landraces in breeding programs is considered a valuable approach to broadening the genetic variability of crops lost during the breeding process and improving traits of commercial importance. The results reported in the present study evidenced the selection for grain yield, HI, NG m², and GNDVI_PA during the breeding process.Whereas differentiation among landraces were found for agronomic and VIs (grain yield, TKW, GNDVI_PA, and GA), in modern cultivars SPs differentiation were mainly due to GS65, GNDVI, and NDVI_PA. The use of a statistical approach as the QTL overview index for the definition of QTL hotspots resulted successful for the identification of consensus genome regions including most of the stable marker trait associations across years. The results of this study will be useful for our wheat breeding program by the selection of the appropriate genotypes carrying favorable alleles for the differential traits that will be useful for designing new crosses.

Using in silico approaches allowed gene mining in QTL hotspots, thus facilitating CG identification.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JS: conceptualization and supervision. RR, AL, ML, JB, and JS: methodology. RR, JB, and JS: formal analysis. RR, AL, and JS: data curation. RR: writing—original draft preparation. RR, ML, JB, and JS: writing—review and editing. ML and JS: project administration and funding acquisition. All authors have read and agreed to the published version of the manuscript.

Funding

This study was funded by the projects AGL2015-65351-R and PID2019-109089RB-C31 from the Spanish Ministry of Science and Innovation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Authors acknowledge Conxita Royo for his helpful revision of the manuscript and the contribution of the CERCA Program (Generalitat de Catalunya).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2021.735192/full#supplementary-material

Supplementary Figure 1. Manhattan plots for the analyzed traits during the 3 years and the mean across years.

Supplementary Table 1. Genotypes used in this study.

Supplementary Table 2. Phenotypic values of the analyzed traits for the 3 years of field trials.

Supplementary Table 3. Boxplots and summary statistics for the analyzed traits. Coefficients of variation for agronomic traits among the 3 years are included.

Supplementary Table 4. Total number of marker–trait associations (MTAs) per year and trait with –log₁₀ P > 3.

Supplementary Table 5. Marker trait associations (MTAs).

Supplementary Table 6. QTL hotspots.

Supplementary Table 7. Gene models within QTL hotspots.

References

Adamsen, F. J., Pinter, P. J., Barnes, E. M., LaMorte, R. L., Wall, G. W., Leavitt, S. W., et al. (1999). Measuring wheat senescence with a digital camera. Crop Sci. 39, 719–724. doi: 10.2135/cropsci1999.0011183X003900030019x

PubMed Abstract | CrossRef Full Text | Google Scholar

Anuarbek, S., Id, S. A., Pecchioni, N., Laid,ò, G., Maccaferri, M., Tuberosa, R., et al. (2020). Quantitative trait loci for agronomic traits in tetraploid wheat for enhancing grain yield in Kazakhstan environments. PLoS ONE. 15:e0234863. doi: 10.1371/journal.pone.0234863

PubMed Abstract | CrossRef Full Text | Google Scholar

Aparicio, N., Villegas, D., Araus, J. L., Casadesús, J., and Royo, C. (2002). Relationship between growth traits and spectral vegetation indices in durum wheat. Crop Sci. 42, 1547–1555. doi: 10.2135/cropsci2002.1547

CrossRef Full Text | Google Scholar

Aparicio, N., Villegas, D., Casadesus, J., Araus, J. L., and Royo, C. (2000). Spectral vegetation indices as non-destructive tools for determining durum wheat yield. Agron. J. 92, 83–91. doi: 10.2134/agronj2000.92183x

CrossRef Full Text | Google Scholar

Araus, J. L., and Cairns, J. E. (2014). Field high-throughput phenotyping: the new crop breeding frontier. Trends Plant Sci. 19, 52–61. doi: 10.1016/j.tplants.2013.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Araus, J. L., Kefauver, S. C., Zaman-Allah, M., Olsen, M. S., and Cairns, J. E. (2018). Translating high-throughput phenotyping into genetic gain. Trends Plant Sci. 23, 451–466. doi: 10.1016/j.tplants.2018.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Araus, J. L., Slafer, G. A., Reynolds, M. P., and Royo, C. (2002). Plant breeding and drought in C3 cereals: what should we breed for? Ann. Bot. 89, 925–940. doi: 10.1093/aob/mcf049

PubMed Abstract | CrossRef Full Text | Google Scholar

Autrique, E., Nachit, M. M., Monneveux, P., Tanksley, S. D., and Sorrells, M. E. (1996). Genetic diversity in durum wheat based on RFLPs, morphophysiological traits, and coefficient of parentage. Crop Sci. 36, 735–742. doi: 10.2135/cropsci1996.0011183X003600030036x

CrossRef Full Text | Google Scholar

Avni, R., Nave, M., Barad, O., Baruch, K., Twardziok, S. O., Gundlach, H., et al. (2017). Wild emmer genome architecture and diversity elucidate wheat evolution and domestication. Science 357, 93–97. doi: 10.1126/science.aan0032

PubMed Abstract | CrossRef Full Text | Google Scholar

Babar, M. A., Reynolds, M. P., van Ginkel, M., Klatt, A. R., Raun, W. R., and Stone, M. L. (2006). Spectral reflectance to estimate genetic variation for in-season biomass, leaf chlorophyll, and canopy temperature in wheat. Crop Sci. 46, 1046–1057. doi: 10.2135/cropsci2005.0211

CrossRef Full Text | Google Scholar

Bellvert, J., Nieto, H., Pelechá, A., Jofre-Cekalović, C., Zazurca, L., and Miarnau, X. (2021). Remote sensing energy balance model for the assessment of crop evapotranspiration and water status in an almond rootstock collection. Front. Plant Sci. 12:288. doi: 10.3389/fpls.2021.608967

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamini, Y., and Hochberg, Y. (1995). Controlling the false discovery rate: apractical and powerful approach to multiple testing. J. Royal Stat. Soc. 57, 289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

Bhatta, M., Morgounov, A., Belamkar, V., and Baenziger, P. (2018). Genome-wide association study reveals novel genomic regions for grain yield and yield-related traits in drought-stressed synthetic hexaploid wheat. Int. J. Mol. Sci. 19:3011. doi: 10.3390/ijms19103011

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, R., Guo, L., Ma, M., Zhang, W., Liu, X., and Zhao, H. (2019). Identification and functional characterization of Squamosa promoter binding protein-like gene Taspl16 in wheat (Triticum aestivum l.). Front. Plant Sci. 10:212. doi: 10.3389/fpls.2019.00212

PubMed Abstract | CrossRef Full Text | Google Scholar

Casadesús, J., and Villegas, D. (2014). Conventional digital cameras as a tool for assessing leaf area index and biomass for cereal breeding. J. Integr. Plant Biol. 56, 7–14. doi: 10.1111/jipb.12117

PubMed Abstract | CrossRef Full Text | Google Scholar

Chao, S., Dubcovsky, J., Dvorak, J., Luo, M.-C., Baenziger, S. P., Matnyazov, R., et al. (2010). Population- and genome-specific patterns of linkage disequilibrium and SNP variation in spring and winter wheat (Triticum aestivum L.). BMC Genomics 11:727. doi: 10.1186/1471-2164-11-727

PubMed Abstract | CrossRef Full Text | Google Scholar

Chardon, F., Virlon, B., Moreau, L., Falque, M., Joets, J., Decousset, L., et al. (2004). Genetic architecture of flowering time in maize as inferred from quantitative trait loci meta-analysis and synteny conservation with the rice genome. Genetics 168, 2169–2185. doi: 10.1534/genetics.104.032375

PubMed Abstract | CrossRef Full Text | Google Scholar

Condorelli, G. E., Maccaferri, M., Newcomb, M., Andrade-Sanchez, P., White, J. W., French, A. N., et al. (2018). Comparative aerial and ground based high throughput phenotyping for the genetic dissection of NDVI as a proxy for drought adaptive traits in durum wheat. Front. Plant Sci. 9:893. doi: 10.3389/fpls.2018.00893

PubMed Abstract | CrossRef Full Text | Google Scholar

Crowell, S., Korniliev, P., Falcão, A., Ismail, A., Gregorio, G., Mezey, J., et al. (2016). Genome-wide association and high-resolution phenotyping link Oryza sativa panicle traits to numerous trait-specific QTL clusters. Nat. Commun. 7, 1–14. doi: 10.1038/ncomms10527

PubMed Abstract | CrossRef Full Text | Google Scholar

Duan, T., Chapman, S. C., Guo, Y., and Zheng, B. (2017). Dynamic monitoring of NDVI in wheat agronomy and breeding trials using an unmanned aerial vehicle. F. Crop. Res. 210, 71–80. doi: 10.1016/j.fcr.2017.05.025

CrossRef Full Text | Google Scholar

Edae, E. A., Byrne, P. F., Haley, S. D., Lopes, M. S., and Reynolds, M. P. (2014). Genome-wide association mapping of yield and yield components of spring wheat under contrasting moisture regimes. Theor. Appl. Genet. 127, 791–807. doi: 10.1007/s00122-013-2257-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Flint-Garcia, S. A., Thornsberry, J. M., and Buckler, E. S. (2003). Structure of linkage disequilibrium in plants. Annu. Rev. Plant Biol. 54, 357–374. doi: 10.1146/annurev.arplant.54.031902.134907

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, F., Wen, W., Liu, J., Rasheed, A., Yin, G., Xia, X., 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. doi: 10.3389/fpls.2015.01099

PubMed Abstract | CrossRef Full Text | Google Scholar

García Del Moral, L. F., García Del Moral, M. B., Molina-Cano, J. L., and Slafer, G. A. (2003). Yield stability and development in two- and six-rowed winter barleys under Mediterranean conditions. F. Crop. Res. 81, 109–119. doi: 10.1016/S0378-4290(02)00215-0

CrossRef Full Text | Google Scholar

García Del Moral, L. F., Rharrabti, Y., Elhani, S., Martos, V., and Royo, C. (2005). Yield formation in Mediterranean durum wheats under two contrasting water regimes based on path-coefficient analysis. Euphytica 146, 203–212. doi: 10.1007/s10681-005-9006-2

CrossRef Full Text | Google Scholar

Giorgi, F., and Lionello, P. (2008). Climate change projections for the Mediterranean region. Glob. Planet. Change 63, 90–104. doi: 10.1016/j.gloplacha.2007.09.005

CrossRef Full Text | Google Scholar

Gitelson, A. A., Kaufman, Y. J., and Merzlyak, M. N. (1996). Use of a green channel in remote sensing of global vegetation from EOS- MODIS. Remote Sens. Environ. 58, 289–298. doi: 10.1016/S0034-4257(96)00072-7

CrossRef Full Text | Google Scholar

Gitelson, A. A., Kaufman, Y. J., Stark, R., and Rundquist, D. (2002). Novel algorithms for remote estimation of vegetation fraction. Remote Sens. Environ. 80, 76–87. doi: 10.1016/S0034-4257(01)00289-9

CrossRef Full Text | Google Scholar

Gizaw, S. A., Godoy, J. G. V., Pumphrey, M. O., and Carter, A. H. (2018). Spectral reflectance for indirect selection and genome-wide association analyses of grain yield and drought tolerance in North American Spring Wheat. Crop Sci. 58, 2289–2301. doi: 10.2135/cropsci2017.11.0690

CrossRef Full Text | Google Scholar

Gomez-Candon, D., Bellvert Rios, J., and Royo, C. (2021). Performance of the two-sourceenergy balance (TSEB) model as a tool for monitoring the response of durum wheat to drought by high-throughput field phenotyping. Front. Plant Sci. 12:658357. doi: 10.3389/fpls.2021.658357

PubMed Abstract | CrossRef Full Text | Google Scholar

Gracia-Romero, A., Kefauver, S. C., Fernandez-Gallego, J. A., Vergara-Díaz, O., Nieto-Taladriz, M. T., and Araus, J. L. (2019). UAV and ground image-based phenotyping: a proof of concept with durum wheat. RemoteSens 11:1244. doi: 10.3390/rs11101244

CrossRef Full Text | Google Scholar

Gracia-Romero, A., Kefauver, S. C., Vergara-Díaz, O., Zaman-Allah, M. A., Prasanna, B. M., Cairns, J. E., et al. (2017). Comparative performance of ground vs. aerially assessed RGB and multispectral indices for early-growth evaluation of maize performance under phosphorus fertilization. Front. Plant Sci. 8:2004. doi: 10.3389/fpls.2017.02004

PubMed Abstract | CrossRef Full Text | Google Scholar

Graziani, M., Maccaferri, M., Royo, C., Salvatorelli, F., and Tuberosa, R. (2014). QTL dissection of yield components and morpho-physiological traits in a durum wheat elite population tested in contrasting thermo-pluviometric conditions. Crop Pasture Sci. 65:80. doi: 10.1071/CP13349

CrossRef Full Text | Google Scholar

Haboudane, D., Miller, J. R., Pattey, E., Zarco-tejada, P. J., and Strachan, I. B. (2004). Hyperspectral vegetation indices and novel algorithms for predicting green LAI of crop canopies: modeling and validation in the context of precision agriculture. Remote Sens. 90, 337–352. doi: 10.1016/j.rse.2003.12.013

CrossRef Full Text | Google Scholar

Haboudane, D., Miller, J. R., Tremblay, N., Zarco-Tejada, P. J., and Dextraze, L. (2002). Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture. Remote Sens. Environ. 81, 416–426. doi: 10.1016/S0034-4257(02)00018-4

CrossRef Full Text | Google Scholar

Hütsch, B. W., Jahn, D., and Schubert, S. (2019). Grain yield of wheat (Triticum aestivum L.) under long-term heat stress is sink-limited with stronger inhibition of kernel setting than grain filling. J. Agron. CropSci. 205, 22–32. doi: 10.1111/jac.12298

CrossRef Full Text | Google Scholar

IWGSC (2018). Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science 361:eaar7191. doi: 10.1126/science.aar7191

PubMed Abstract | CrossRef Full Text | Google Scholar

Juliana, P., Montesinos-López, O. A., Crossa, J., Mondal, S., González Pérez, L., Poland, J., et al. (2019). Integrating genomic-enabled prediction and high-throughput phenotyping in breeding for climate-resilient bread wheat. Theor. Appl. Genet. 132, 177–194. doi: 10.1007/s00122-018-3206-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Knapp, S. J., Stroup, W. W., and Ross, W. M. (1985). Exact confidence intervals for heritability on a progeny mean basis. Crop Sci. 25, 192–194. doi: 10.2135/cropsci1985.0011183X002500010046x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kulkarni, M., Soolanayakanahally, R., Ogawa, S., and Uga, Y. (2017). Drought response in wheat : key genes and regulatory mechanisms controlling root system architecture and transpiration. Efficiency 5, 1–13. doi: 10.3389/fchem.2017.00106

PubMed Abstract | CrossRef Full Text | Google Scholar

Kumar, M., Yusuf, M. A., Yadav, P., Narayan, S., and Kumar, M. (2019). Overexpression of chickpea defensin gene confers tolerance to water-deficit stress in Arabidopsis thaliana. Front. Plant Sci. 10:290. doi: 10.3389/fpls.2019.00290

PubMed Abstract | CrossRef Full Text | Google Scholar

Kyratzis, A. C., Skarlatos, D. P., Menexes, G. C., Vamvakousis, V. F., and Katsiotis, A. (2017). Assessment of vegetation indices derived by UAV imagery for durum wheat phenotyping under a water limited and heat stressed Mediterranean environment. Front. Plant Sci. 8:1114. doi: 10.3389/fpls.2017.01114

PubMed Abstract | CrossRef Full Text | Google Scholar

Leegood, R. C., Evans, J. R., and Furbank, R. T. (2010). Food security requires genetic advances to increase farm yields. Nature 464:831. doi: 10.1038/464831d

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Mullan, D., Zhang, C., Zhao, S., Li, X., Zhang, A., et al. (2020). Major genomic regions responsible for wheat yield and its components as revealed by meta-QTL and genotype–phenotype association analyses. Planta 252:65. doi: 10.1007/s00425-020-03466-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, M. S., El-Basyoni, I., Baenziger, P. S., Singh, S., Royo, C., Ozbek, K., et al. (2015). Exploiting genetic diversity from landraces in wheat breeding for adaptation to climate change. J. Exp. Bot. 66, 3477–3486. doi: 10.1093/jxb/erv122

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, M. S., and Reynolds, M. P. (2012). Stay-green in spring wheat can be determined by spectral reflectance measurements (normalized difference vegetation index) independently from phenology. J. Exp. Bot. 63, 3789–3798. doi: 10.1093/jxb/ers071

PubMed Abstract | CrossRef Full Text | Google Scholar

Maccaferri, M., El-Feki, W., Nazemi, G., Salvi, S., Can,è, M. A., Colalongo, M. C., et al. (2016). Prioritizing quantitative trait loci for root system architecture in tetraploid wheat. J. Exp. Bot. 67, 1161–1178. doi: 10.1093/jxb/erw039

PubMed Abstract | CrossRef Full Text | Google Scholar

Maccaferri, M., Harris, N. S., Twardziok, S. O., Pasam, R. K., Gundlach, H., Spannagl, M., et al. (2019). Durum wheat genome highlights past domestication signatures and future improvement targets. Nat. Genet. 51, 885–895. doi: 10.1038/s41588-019-0381-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangini, G., Gadaleta, A., Colasuonno, P., Marcotuli, I., Signorile, A. M., Simeone, R., et al. (2018). Genetic dissection of the relationships between grain yield components by genome-wide association mapping in a collection of tetraploid wheats. PLoS ONE 13:e0190162. doi: 10.1371/journal.pone.0190162

PubMed Abstract | CrossRef Full Text | Google Scholar

Mérida-García, R., Bentley, A. R., Gálvez, S., Dorado, G., Solís, I., Ammar, K., et al. (2020). Mapping agronomic and quality traits in elite durum wheat lines under differing water regimes. Agronomy 10:144. doi: 10.3390/agronomy10010144

CrossRef Full Text | Google Scholar

Moragues, M., García Del Moral, L. F., Moralejo, M., and Royo, C. (2006). Yield formation strategies of durum wheat landraces with distinct pattern of dispersal within the Mediterranean basin: II. Biomass production and allocation. F. Crop. Res. 95, 182–193. doi: 10.1016/j.fcr.2005.02.008

CrossRef Full Text | Google Scholar

Mwadzingeni, L., Shimelis, H., Rees, D. J. G., and Tsilo, T. J. (2017). Genome-wide association analysis of agronomic traits in wheat under drought-stressed and non-stressed conditions. PLoS ONE 12:e0171692. doi: 10.1371/journal.pone.0171692

PubMed Abstract | CrossRef Full Text | Google Scholar

Qaseem, M. F., Qureshi, R., Shaheen, H., and Shafqat, N. (2019). Genome-wide association analyses for yield and yield-related traits in bread wheat (Triticum aestivum L.) under pre-anthesis combined heat and drought stress in field conditions. PLoS ONE 14:e0213407. doi: 10.1371/journal.pone.0213407

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, J., Chehbouni, A., Huete, A. R., Kerr, Y. H., and Sorooshian, S. (1994). A modified soil adjusted vegetation index. Remote Sens. Environ. 48, 119–126. doi: 10.1016/0034-4257(94)90134-1

CrossRef Full Text | Google Scholar

Ramírez-González, R. H., Borrill, P., Lang, D., Harrington, S. A., Brinton, J., Venturini, L., et al. (2018). The transcriptional landscape of polyploid wheat. Science 361:aar6089. doi: 10.1126/science.aar6089

PubMed Abstract | CrossRef Full Text | Google Scholar

Roselló, M., Royo, C., Sanchez-Garcia, M., and Soriano, J. M. (2019a). Genetic dissection of the seminal root system architecture in mediterranean durum wheat landraces by genome-wide association study. Agronomy 9:364. doi: 10.3390/agronomy9070364

CrossRef Full Text | Google Scholar

Roselló, M., Villegas, D., Álvaro, F., Soriano, J. M., Lopes, M. S., Nazco, R., et al. (2019b). Unravelling the relationship between adaptation pattern and yield formation strategies in Mediterranean durum wheat landraces. Eur. J. Agron. 107, 43–52. doi: 10.1016/j.eja.2019.04.003

CrossRef Full Text | Google Scholar

Roujean, J. L., and Breon, F. M. (1995). Estimating PAR absorbed by vegetation from bidirectional reflectance measurements. Remote Sens. Environ. 51, 375–384. doi: 10.1016/0034-4257(94)00114-3

CrossRef Full Text | Google Scholar

Rouse, J. W., Haas, R. H., Schell, J. A., Deering, D. W., and Harlan, J. C. (1974). Monitoring the Vernal Advancement and Retrogradation (Green Wave Effect) of Natural Vegetation. Greenbelt, MD: NASA/GSFC. Available online at: https://ntrs.nasa.gov/search.jsp?R=19750020419 (accessed June 16, 2020).

Google Scholar

Royo, C., Ammar, K., Villegas, D., and Soriano, J. M. (2021). Agronomic, physiological and genetic changes associated with evolution, migration and modern breeding in durum wheat. Front. Plant Sci. 12:674470. doi: 10.3389/fpls.2021.674470

PubMed Abstract | CrossRef Full Text | Google Scholar

Royo, C., Aparicio, N., Blanco, R., and Villegas, D. (2004). Leaf and green area development of durum wheat genotypes grown under Mediterranean conditions. Eur. J. Agron. 20, 419–430. doi: 10.1016/S1161-0301(03)00058-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Royo, C., Aparicio, N., Villegas, D., Casadesus, J., Monneveux, P., and Araus, J. L. (2003). Usefulness of spectral reflectance indices as durum wheat yield predictors under contrasting Mediterranean conditions. Int. J. Remote Sens. 24, 4403–4419. doi: 10.1080/0143116031000150059

CrossRef Full Text | Google Scholar

Royo, C., Dreisigacker, S., Ammar, K., and Villegas, D. (2020). Agronomic performance of durum wheat landraces and modern cultivars and its association with genotypic variation in vernalization response (Vrn-1) and photoperiod sensitivity (Ppd-1) genes. Eur. J. Agron. 120:126129. doi: 10.1016/j.eja.2020.126129

CrossRef Full Text | Google Scholar

Royo, C., Fanny, A. E., LvaroAe, A., Martos, V., Abdelhamid, A. E., Ae, R., et al. (2006). Genetic changes in durum wheat yield components and associated traits in Italian and Spanish varieties during the 20th century. Euphytica 9, 259–270. doi: 10.1007/s10681-006-9327-9

CrossRef Full Text | Google Scholar

Royo, C., Nazco, R., and Villegas, D. (2014). The climate of the zone of origin of Mediterranean durum wheat (Triticum durum Desf.) landraces affects their agronomic performance. Genet. Resour. Crop Evol. 61, 1345–1358. doi: 10.1007/s10722-014-0116-3

CrossRef Full Text | Google Scholar

Rufo, R., Alvaro, F., Royo, C., and Soriano, J. M. (2019). From landraces to improved cultivars: assessment of genetic diversity and population structure of Mediterranean wheat using SNP markers. PLoS ONE 14:219867. doi: 10.1371/journal.pone.0219867

PubMed Abstract | CrossRef Full Text | Google Scholar

Rufo, R., Salvi, S., Royo, C., and Soriano, J. (2020). Exploring the genetic architecture of root-related traits in mediterranean bread wheat landraces by genome-wide association analysis. Agronomy 10:613. doi: 10.3390/agronomy10050613

CrossRef Full Text | Google Scholar

Rufo, R., Soriano, J. M., Villegas, D., Royo, C., and Bellvert, J. (2021). Using unmanned aerial vehicle and ground-based RGB indices to assess agronomic performance of wheat landraces and cultivars in a mediterranean-type environment. Remote Sens. 13, 1–18. doi: 10.3390/rs13061187

CrossRef Full Text | Google Scholar

Slafer, G. A., Araus, J. L., Royo, C., and García Del Moral, L. F. (2005). Promising eco-physiological traits for genetic improvement of cereal yields in Mediterranean environments. Ann. Appl. Biol. 146, 61–70. doi: 10.1111/j.1744-7348.2005.04048.x

CrossRef Full Text | Google Scholar

Soriano, J. M., Malosetti, M., Rosell,ó, M., Sorrells, M. E., and Royo, C. (2017). Dissecting the old Mediterranean durum wheat genetic architecture for phenology, biomass and yield formation by association mapping and QTL meta-analysis. PLoS ONE 12, 1–19. doi: 10.1371/journal.pone.0178290

PubMed Abstract | CrossRef Full Text | Google Scholar

Soriano, J. M., Sansaloni, C., Ammar, K., and Royo, C. (2021). Labelling selective sweeps used in durum wheat breeding from a diverse and structured panel of landraces and cultivars. Biology 10, 1–20. doi: 10.3390/biology10040258

PubMed Abstract | CrossRef Full Text | Google Scholar

Soriano, J. M., Villegas, D., Aranzana, M. J., García Del Moral, L. F., and Royo, C. (2016). Genetic structure of modern durum wheat cultivars and mediterranean landraces matches with their agronomic performance. PLoS ONE 11:e0160983. doi: 10.1371/journal.pone.0160983

PubMed Abstract | CrossRef Full Text | Google Scholar

Soriano, J. M., Villegas, D., Sorrells, M. E., and Royo, C. (2018). Durum wheat landraces from east and west regions of the mediterranean basin are genetically distinct for yield components and phenology. Front. Plant Sci. 9, 1–9. doi: 10.3389/fpls.2018.00080

PubMed Abstract | CrossRef Full Text | Google Scholar

Stenberg, P., Rautiainen, M., Manninen, T., Voipio, P., and Smolander, H. (2004). Reduced simple ratio better than NDVI for estimating LAI in Finnish pine and spruce stands. Silva Fenn. 38, 3–14. doi: 10.14214/sf.431

CrossRef Full Text | Google Scholar

Subira, J., Álvaro, F., García del Moral, L. F., and Royo, C. (2015). Breeding effects on the cultivar × environment interaction of durum wheat yield. Eur. J. Agron. 68, 78–88. doi: 10.1016/j.eja.2015.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Sukumaran, S., Reynolds, M. P., and Sansaloni, C. (2018). Genome-wide association analyses identify QTL hotspots for yield and component traits in durum wheat grown under yield potential, drought, and heat stress environments. Front. Plant Sci. 9:81. doi: 10.3389/fpls.2018.00081

PubMed Abstract | CrossRef Full Text | Google Scholar

Valluru, R., Reynolds, M. P., Davies, W. J., and Sukumaran, S. (2017). Phenotypic and genome-wide association analysis of spike ethylene in diverse wheat genotypes under heat stress. New Phytol. 214, 271–283. doi: 10.1111/nph.14367

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S., Wong, D., Forrest, K., Allen, A., Chao, S., Huang, B. E., et al. (2014). Characterization of polyploid wheat genomic diversity using a high-density 90,000 single nucleotide polymorphism array. Plant Biotechnol. J. 12, 787–796. doi: 10.1111/pbi.12183

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, S.-X., Zhu, Y.-L., Zhang, D.-X., Shao, H., Liu, P., Hu, J.-B., et al. (2017). Genome-wide association study for grain yield and related traits in elite wheat varieties and advanced lines using SNP markers. PLoS ONE 12:e0188662. doi: 10.1371/journal.pone.0188662

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Guan, P., Xin, M., Wang, Y., Chen, X., Zhao, A., et al. (2020). Genome-wide association study identifies QTL for thousand grain weight in winter wheat under normal- and late-sown stressed environments. Theor. Appl. Genet. 134, 143–157. doi: 10.1007/s00122-020-03687-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Ward, J. H. (1963). Hierarchical grouping to optimize an objective function. J. Am. Stat. Assoc. 58, 236–244. doi: 10.1080/01621459.1963.10500845

CrossRef Full Text | Google Scholar

White, J., Andrade-Sanchez, P., and Gore, M. (2012). Field-based phenomics for plant genetics research. F. Crop. Res. 133, 101–113. doi: 10.1016/j.fcr.2012.04.003

CrossRef Full Text | Google Scholar

Xie, C., and Yang, C. (2020). A review on plant high-throughput phenotyping traits using UAV-based sensors. Comput. Electron. Agric. 178:105731. doi: 10.1016/j.compag.2020.105731

CrossRef Full Text | Google Scholar

Yoon, S.-K., Park, E.-J., Choi, Y.-I., Bae, E.-K., Kim, J.-H., Park, S.-Y., et al. (2014). Response to drought and salt stress in leaves of poplar (Populus alba × Populus glandulosa): expression profiling by oligonucleotide microarray analysis. Plant Physiol. Biochem. 84, 158–168. doi: 10.1016/j.plaphy.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Zadoks, J. C., Chang, T. T., and Konzak, C. F. (1974). A decimal code for the growth stages of cereals. Weed Res. 14, 415–421. doi: 10.1111/j.1365-3180.1974.tb01084.x

CrossRef Full Text | Google Scholar

Zarco-Tejada, P. J., Ustin, S. L., and Whiting, M. L. (2005). Temporal and spatial relationships between within-field yield variability in cotton and high-spatial hyperspectral remote sensing imagery. Agron. J. 97, 641–653. doi: 10.2134/agronj2003.0257

CrossRef Full Text | Google Scholar

Zhao, Y., Cheng, X., Liu, X., Wu, H., Bi, H., and Xu, H. (2018). The wheat MYB transcription factor TaMYB31 is involved in drought stress responses in arabidopsis. Front. Plant Sci. 9:1426. doi: 10.3389/fpls.2018.01426

PubMed Abstract | CrossRef Full Text | Google Scholar