Root and canopy traits and adaptability genes explain drought tolerance responses in winter wheat

Bread wheat (Triticum aestivum L) is one of the three main staple crops worldwide contributing 20% calories in the human diet. Drought stress is the main factor limiting yields and threatening food security, with climate change resulting in more frequent and intense drought. Developing drought-tolerant wheat cultivars is a promising way forward. The use of holistic approaches that include high-throughput phenotyping and genetic markers in selection could help in accelerating genetic gains. Fifty advanced breeding lines were selected from the CIMMYT Turkey winter wheat breeding program and studied under irrigated and semiarid conditions in two years. High-throughput phenotyping was done for wheat crown root traits and canopy senescence dynamics using vegetation indices (green area using RGB images and Normalized Difference Vegetation Index using spectral reflectance). In addition, genotyping by KASP markers for adaptability genes was done. Overall, under semiarid conditions yield reduced by 3.09 t ha-1 (-46.8%) compared to irrigated conditions. Genotypes responded differently under drought stress and genotypes 39 (VORONA/HD24-12//GUN/7/VEE#8//…/8/ALTAY), 18 (BiII98) and 29 (NIKIFOR//KROSHKA) were the most drought tolerant. Root traits including shallow nodal root angle under irrigated conditions and root number per shoot under semiarid conditions were correlated with increased grain yield. RGB based vegetation index measuring canopy green area at anthesis was better correlated with GY than NDVI was with GY under drought. The markers for five established functional genes (PRR73.A1 –flowering time, TEF-7A –grain size and weight, TaCwi.4A - yield under drought, Dreb1- drought tolerance, and ISBW11.GY.QTL.CANDIDATE- grain yield) were associated with different drought-tolerance traits in this experiment. We conclude that–genotypes 39, 18 and 29 could be used for drought tolerance breeding. The trait combinations of canopy green area at anthesis, and root number per shoot along with key drought adaptability makers (TaCwi.4A and Dreb1) could be used in screening drought tolerance wheat breeding lines.


Introduction
Wheat (Triticum aestivum L.) is one of the most important food crops contributing around 20% of calories in the human diet worldwide. However, climate change has resulted in more frequent and intense periods of drought which affect wheat production [1]. Worldwide drought is the most important factor affecting wheat yields [2] and model-based predictions indicate that there will be 9-12% wheat yield reduction with climate change in 21 st century without considering the benefits of CO 2 fertilization and adaptations [3]. Developing new drought-tolerant varieties is therefore important to achieve food security in the context of climate change. Identifying drought-tolerance traits for deployment in breeding is therefore crucial. A deeper root system and cooler canopy temperature are two important traits reported to be correlated with for drought tolerance in wheat [4][5][6]. As roots are difficult to study under field conditions, canopy temperature has been applied as an indirect way of assessing the root depth as higher water uptake and leaf transpiration is related to a cooler canopy. Canopy staygreen characters have also shown promise to help select drought-tolerant genotypes as they may confer extended photosynthesis, nutrient and water uptake under stress conditions [7,8].
Wheat root systems consist of seminal roots (up to 6) and crown roots (around 10-15) per plant the latter emerging from basal node of main shoots and tillers [9]. These two root systems function together to acquire water and nutrients from the soil [10][11][12]. Distribution of root length density (root length per unit soil volume; RLD) with depth is an important trait affecting water capture in wheat crops [7,10,13]. In synthetic wheat-derived lines, yield increase under water stress conditions was correlated with increase in root dry weight at depth in NW Mexico [14]. Higher allocation of plant assimilates to deeper roots has been correlated with a cooler canopy and increase in overall grain yield under drought conditions in wheat synthetic-derived material [5]. Narrower root angle in the top soil (steeper roots) has correlated with higher root density in deeper soil in wheat genotypes in Australia [14][15][16].
The high-throughput phenotyping of root system architecture under the field conditions presents a bottleneck in breeding for drought tolerance in wheat [17]. Previously, the soil-core break method [18] and 'shovelomics' [19] have been used for high-throughput field phenotyping in cereals. The core-break method involves the counting of roots visible on the cross-section to estimate root length density [20]. Shovelomics, focuses on crown root phenotyping, and involves the excavation of roots in the topsoil and measuring root traits manually or through image analysis. Results of direct measurements and visual scoring in maize showed correlations with root depth for crown root number and angle [19]. Shovelomics methods have quantified genetic variation in crown root angle and root length in maize [19,21,22], barley [23] and durum wheat [24]. We used a high-throughput shovelomics technique for phenotyping root crown architecture of the whole root crown in wheat.
The use of Normalized Difference Vegetation Index (NDVI) spectral reflectance index to study canopy growth and senescence dynamics is well established, but has some limitations [25,26]. Especially at high values of leaf area index (LAI), NDVI tends to saturate and does not show strong linear association with yield components [27,28]. Also, in order to obtain accurate results with NDVI, bright conditions with direct sunlight are required while taking the measurements (in the case of passive sensors). Senescence is a genetically programmed and environmentally influenced process [29] and the stay-green phenotype has been shown to be improving yields under drought [8,30]. NDVI has been used to measure stay green trait in wheat under drought [8]. More recently, RGB image-based vegetation indices have proved to be better correlated with grain yield than NDVI under similar circumstances and also are time-saving [28].
Marker-assisted selection is a very important component of molecular breeding to develop resilient cultivars by selecting and accumulating favorable alleles. In bread wheat several genes underpinning drought tolerance have been identified and molecular markers have been developed to select favorable alleles [31,32]. However, the distribution and association of alleles for functional genes like Dreb1, PRR-73, TaCwi-A1 and TEF-7 associated with drought tolerance is largely unknown in wheat cultivars from the most parts of the world. Transcription factors like DREB control the expression of several functional genes responsible for plant tolerance to drought, and have been proposed for use in plant improvement for drought tolerance [33]. Similarly, TaTEF-7A is transcript elongation factor gene responsible for grain number per spike [34], thousand grain weight and chlorophyll content at grain-filling stage under drought stress [35]. The gene TaPRR73 was found to be regulation of flowering date and can be used in breeding to develop cultivars adaptable for different geographical areas [36]. TaCwi-A1 gene produce cell wall invertase enzyme mainly responsible for sink tissue development and carbon allocation and also showed affecting grain weight and potential role in drought tolerance in wheat [37,38].
The present study reports associations between nodal root traits measured using the wheat shovelomics techniques along with vegetation indices in a set of 50 CIMMYT Turkey winter wheat cultivars and advanced lines. The experiment was conducted under irrigated (IR) and semiarid (SA) field conditions in Turkey in two cropping seasons. The germplasm was also screened for allelic variation of genes previously related to drought adaptability genes using breeder friendly KASP markers. Our aim was to quantify genetic variation in grain yield responses to drought and its physiological basis and to identify drought adaptability genes for potential use in marker-assisted selection for drought tolerance in wheat.

Experimental design and plot management
Two field experiments were sown on 15 Nov 2017 and 8 Nov 2018 at Bahri Dagdas International Agricultural Research Institute, Konya in 2017-18 and 2018-19 wheat growing season. Before sowing the experimental field was fallow. The soil type was sandy clay. Experiments were conducted in a randomized block, split-plot design, in which two irrigation treatments (IR: drip-irrigated and SA: semiarid/rain-fed) were randomized on main-plots, and 50 CIM-MYT winter wheat cultivars and advanced lines including 4 check cultivars were randomized on sub-plots in two replicates. The 50 winter wheat genotypes represent modern germplasm developed by the CIMMYT-International Winter Wheat Improvement Program and obtained from cooperators in Eastern Europe. These lines were selected from 100 genotypes tested in previous years and represented three groups: high yielding under irrigation; high yielding under drought and balanced performance under both environments. The check cultivars used were Gerek (widely grown drought resistant check), Katea (widely grown irrigated check), Konya (high yield potential irrigated check) and Nacibey (check for supplementary irrigation) (S1 Table). Plots were 7.0 m × 1.2 m with 6 rows 20 cm apart and 450 seeds were sown per square meter. Fertilizers applied were 100 kg ha −1 of phosphorus (P) and 140 kg ha −1 of nitrogen (as ammonium nitrate) per hectare at the time of planting, and an additional 50 kg ha −1 of N at tillering (GS35). Under the irrigated treatment drip-irrigation was given as 50 mm application each time. Irrigation was given twice during the crop growth season at the tillering and flowering stage. The irrigation was applied according to past experience at the field site so that under irrigated conditions there would not be significant water stress.

Crop measurements: Grain yield and yield components
In 2018, a 1.5 m row bulk sample was hand-harvested by cutting at ground level at physiological maturity (GS89) [39]. The fertile shoots (those with an ear) were counted and 5 primary (large ear and stem) fertile shoots were selected for dry matter partitioning analysis. In 2019, 10-20 shoots were selected for dry matter partitioning analysis from the sample taken for root traits analysis (explained in next section). All selected shoots were separated into ears and straw. Dry weight of the ears and the straw was recorded after drying at 80˚C for 48 h. The ears of the bulk sample were then hand threshed and grain weighed. All grains were counted by a Contador seed counter (Pfeuffer, Germany) and thousand-grain weight (TGW) was calculated. From these data the grain DM per fertile shoot, harvest index (HI; grain DM / aboveground DM), fruiting efficiency (grain weight per ear dry weight) and above ground dry matter (AGDM; GY/HI) were calculated. The grain yield was calculated by weighing grain from rest of the plot which was machine-harvested (adjusted to 85% dry weight).

Shovelomics root crown trait measurements
The methodology for root excavation was the same in both the years whereas root traits were assessed in different ways for each years. Root crowns were excavated from all sub-plots during late-grain filling. A spade of 25 cm width and 30 cm depth was inserted to 20 cm depth on either side of plants keeping the blade parallel to the row. A single sample was taken per plot. The soil was placed into a 10 L bucket filled with water overnight. The next day the root crowns were sprayed with low pressure water from a hose to remove remaining soil. Three plants per sample were selected for scanning or image analysis. In 2018, root images were acquired and analyzed using WinRHIZO Regular V. 2009c scanner and software (Regent Instruments Inc., Canada). The traits measured were root surface area (cm 2 ), root diameter (mm), and root volume (cm 3 ) per plant. In 2019, images of the roots were taken with a RGB camera (Sony a 6000). A single image per sample was taken with auto setting. Roots were placed on a black background to maximum contrast and sample ID label and reference scale (white square of 2 cm x 1 cm) was placed on one side of the roots as shown in S1 Fig. Images were analyzed using a modified method from York and Lynch [40]. A project for the ObjectJ plugin (https://sils.fnwi.uva.nl/bcb/objectj) for ImageJ [41] was created to allow the angles, numbers, stem diameter and roots diameter to be annotated and measured from the image of the plant-root samples (S1 Fig). The pixel dimensions were converted to physical units using measurements of the known-sized scale in every image. The traits measured were root number per shoot and plant, root diameter (mm), and root angle (˚). A polyline as showed in S1 Fig was used to measure the crown lengths of the outermost roots. The seminal root length, and the angle was measured for the outermost crown roots at approximately 5 cm depth by measuring the width then later calculating angle using trigonometry and the actual depth measurement to where width was measured (S1 Fig). For nodal root number, each nodal root axis was manually annotated, and the count recorded in an output file. The image analysis gave values for the number of pixels corresponding to root diameter, length and numbers. Using the 2 cm x 1 cm reference square, these pixel values were then converted to the relevant units for each root measurement in Excel and the angles were calculated as detailed in York and Lynch [40]. The trait definitions for the shovelomics root system architecture traits measured in 2019 are given in S3 Table.

NDVI and RGB based vegetation indexes and canopy temperature
In 2018, Normalized Difference Vegetation Index (NDVI) was measured using the handheld active sensor Trimble GreenSeeker spectroradiometer (Trimble Navigation Ltd, USA) to assess the canopy green area starting from booting (GS41). Modified Gompertz curves (Eq 1) were fitted to the NDVI values against thermal time (base temp. 0˚C after anthesis, GS61, [39]). The Gompertz (T) parameter was fitted as the thermal time for the NDVI to decrease to 37% of the NDVI value at GS61. The thermal time (t, measured in˚CD) when NDVI values were 90% and 10% of the value at GS61 was taken as the onset of senescence (SenSt) and end of senescence (SenEnd), respectively, and the duration from 90% NDVI to 10% NDVI remaining was considered as the senescence duration (SenDu).
where t is thermal time (base temp. 0˚C), D is the duration of senescence (SenDu), T is the timing of the inflection point at 37% NDVI value remaining from initial point at GS61, and K is the maximum NDVI at GS61. The senescence parameters were estimated for each sub-plot and then subjected to ANOVA. In 2019, RGB images and NDVI (GreenSeeker Trimble Navigation Ltd, USA) were taken every two weeks from tillering (GS35) to crop maturity (GS89). RGB image-based vegetation index-green area per meter square (GA m -2 ) was calculated using Eqs 2 to 6: whereas GSD is ground sampling distance (centimeters/pixel), SW is camera sensor width (mm), H is camera height from top of the canopy (m), FL is focal length of camera (mm), DW is width of single image footprint on the ground (m), DH is height of single image footprint on the ground (m), IW is image width (pixels), IH is image height (pixels), A is ground area in image (m 2 ) and GA is index value output after image analysis using BreedPix. BreedPix is open source software [42], implemented as part of the open-source CerealScanner plugin developed for ImageJ software [28,41]. Green Area per meter square at anthesis (GA An) and 2 weeks after anthesis (GA 2W) values are used in this paper. Canopy temperature was measured at anthesis using handheld infrared temperature meter (SBRMART GM320).

Genotyping
DNA was extracted from all genotypes using a modified CTAB method [43]. Allele-specific KASP markers for five different loci were used. The primer sequences and amplification conditions of each gene are described in [38]. The detailed genotyping procedures have been described in previous studies [38,44]. Briefly, two allele-specific primers carrying a standard FAM tail (5 0 -GAAGGTGACCAAGTTCATGCT-3 0 ) and HEX tail (5 0 -GAAGGTCGGAGTCAA CGGATT-3 0 ), with targeting SNP at the 3 0 end, and a common reverse primer were synthesized. The primer mixture included 46 μl ddH 2 O, 30 μl common primer (100 μM) and 12 μl of each tailed primer (100 μM). Assays were tested in 384-well format and set up as 5 μl reaction [2.2 μl DNA (10-20 ng/μl), 2.5 μl of 2XKASP master mixture and 0.056 μl primer mixture]. PCR cycling was performed using the following protocol: hot start at 95˚C for 15 min, followed by ten touchdown cycles (95˚C for 20 s; touchdown 65˚C and decreasing by -1˚C per cycle for 25 s) further followed by 30 cycles of amplification (95˚C for 10 s; 57˚C for 60 s). The extension step was not required because amplicon size is less than 120 bp. The plate was read in the Bio-Tek H1 system and data analysis was performed manually using Klustercaller software (version 2.22.0.5; LGC Hoddesdon, United Kingdom).

Marker-traits association analysis
In this paper we presented results for five key KASP markers out of 150 for their association with phenotypes (Table 4). KASP markers for which one of the alleles was represented at relatively higher frequency (>80%) than other alleles were not considered for MTA. These markers are regularly used in marker-assisted selection in CIMMYT's wheat breeding program. As some of the traits that we measured differed between years, the marker-traits associations (MTAs) were identified separately for each year. MTA analysis was done in R using a linear model (lm) function (Eq 7) to test the significant effect of the marker on traits by comparing the mean. BLUEs (best linear unbiased estimators) were calculated using Meta-R for a randomized block design for all the traits for individual years. BLUE values were used to test the significant effect of the marker on traits using following liner model.
Y is phenotyping value, μ is mean of the population, M is mean effect of j th marker, G k (M j ) genotype within marker variance (error variance).

Statistics
In both years, GenStat 19th edition (VSN International, Hemel Hempstead, UK) was used for carrying out analysis of variance (ANOVA) of traits applying a split-plot design with replications and genotypes regarded as random effects and the least significant difference (LSD) test was used to compare the means between specific treatments. A cross-year ANOVA was applied to analyze irrigation treatments and genotypes effects across years and the interaction with year, assuming irrigation treatments and genotypes were fixed effects and replicates and year were random effects. Pearson's correlation coefficient (r) and the linear regressions coefficient (R 2 ) were calculated to quantify associations between traits for mean values in individual years and cross year means using GenStat 19 th editions. Principal component analyses were done to produce biplots using R software package "factoextra." Forward stepwise multi-linear regression was applied to 50 genotypes for each treatment and year separately with GY as the dependent variables and root surface area (RoSuAr), root diameter (RoDiM), root volume (RoVol), NDVI at anthesis (NDVI), NDVI senescence start (SenSt), NDVI senescence duration (SenDu) as independent variables in 2018 and root angle (RoAng), root diameter (RoDiM), root dry weight per shoot (RoDrWtSh), canopy green area per meter square at anthesis (GA An) and NDVI at anthesis (NDVI) in 2019 using GenStat 19th Editions (VSN International 2017). The R 2 statistic values are presented calculated as: 100 × (1 -(residual mean square/total mean square)).

Drought effects on plant growth
Averaging across the 50 genotypes, the drought/semiarid (SA) conditions reduced the grain yield (GY) compared to irrigated (IR) conditions by 2.67 t ha -1 (-50.1%) in 2018 and 3.51 t ha -1 (-44.6%) in 2019 (P < 0.001; Table 1) with an average reduction over two years of 3.09 t ha -1 (-46.8%, P = 0.01, Fig 1). The cross-year ANOVA showed a significant Year x Genotype (Y x G) interaction (P<0.001, Table 1). Relative loss in GY under SA conditions ranged from -36.1% (genotype code no. 32) compared to -58.5% (genotype 9). The three-way interaction of Y x T x G (Year x Treatment x Genotype) was not significant.
Regression analysis showed a positive correlation amongst the genotypes for GY between IR and SA conditions (R 2 = 0.27, P<0.001, Fig 1a). Nevertheless, some genotypes marked differences in percentage reduction in GY; for example: genotype 33 with a 49.3% reduction compared to genotypes 16 and 25 with 34.8 and 40.7% reduction, respectively, under SA conditions (S2 Table). Overall, genotypes 39, 18 and 29 showed relative less yield reduction under SA conditions indicating ability to tolerate the drought whereas genotype 40, 35 and 29 were the most susceptible to drought. GY also showed positive correlation with AGDM (IR: R 2 = 0.21, P<0.001 and SA: R 2 = 0.22., P<0.001) and NDVI at anthesis (IR: R 2 = 0.18, P = 0.01 and SA: R 2 = 0.23., P<0.001) under both IR and SA conditions (Fig 2a and 2c).
For harvest traits, averaging across years, overall the AGDM was the component affected most by the semiarid conditions reducing from 19.7 to 12.1 t ha -1 (-38.6%, P = 0.01, Table 1); whereas TGW reduced from 36.6 to 32.1 (-12.2%, P = 0.03, Table 1). For AGDM, T x G interaction was not significant whereas for TGW it was (P = 0.03). Reduction in grains per ear under SA conditions ranged from -0.5% (genotype 16) to -49.0% (genotype 9). Variation in response to drought for TGW was from -0.9 to -25.7%. Heading date (HD) was advanced by two days in SA conditions compared to IR (P = 0.004). There was a negative association between GY and HD amongst cultivars under SA conditions (R 2 = 0.18, P = 0.002), but no association under IR conditions (Fig 2c). Overall drought reduced plant height by 31.6% but there was no correlation between GY and PH amongst cultivars under either SA or IR conditions. Harvest index showed a positive correlation with GY under IR (R 2 = 0.39, P<0.001) and SA (R 2 = 0.28, P<0.001) conditions (Fig 2b).

Root system traits and correlations with yield and yield components
In 2018, root traits were not significantly affected by the irrigation treatment. However, differences between the genotypes were observed in all the root traits (P<0.05,   Table 4).
In 2019, root diameter (RoDiM), root number per plant (RoNoPl) and root:shoot ratio (Ro: Sh ratio) showed significant differences between genotypes and treatments (Table 3). Overall phenotypic variation amongst genotypes for root angle was from 46.7˚to 68.0˚with mean of 56.6˚and 46.1˚to 63.8˚with mean of 56.3˚under IR and SA conditions, respectively. Root number per shoot (RoNoSh) showed a positive correlation with GY and HI (r = 0.32 and  (Table 4). AGDM showed a negative correlation with root diameter (r = -0.29) and root dry weight per plant (r = -0.49) under IR conditions. Wider root angle was also correlated with higher GY and AGDM under IR conditions (r = 0.29 and r = 0.33, respectively; Table 5). Narrower root angle was correlated with more roots per plant under SA conditions, whereas under IR conditions these correlations were not significant. Root dry weight per plant also showed a positive correlation with root diameter under IR conditions. There was a strong positive correlation between root dry weight and root number per plant under both SA and IR conditions.
Traits correlations were also explored in a subset of 30 genotypes with median flowering dates and which differed in flowering date by only 1 day in each of the two seasons. The correlation matrices for this group of genotypes are shown in S4 and S5 Tables. However, there was no meaningful change in the correlations between the group of genotypes with narrow flowering date (S4 and S5 Tables) as compared to the correlations for the full set of genotypes.

Canopy senescence and temperature traits
Overall, in 2018 NDVI at anthesis (NDVI, GS61) ranged from 0.56 to 0.71 and 0.40 to 0.61 under IR and SA conditions, respectively. ANOVA shows that there was a significant difference between genotypes and irrigation treatments along with significant G x T interaction for NDVI at anthesis, senescence start (SenSt) and senescence duration (SenDu) ( Table 2). There was a positive correlation between NDVI at anthesis and GY under both, IR (r = 0.34) and SA (r = 0.34) conditions (Table 4). Under IR conditions, NDVI at anthesis and senescence start (SenSt) showed positive correlations with GY (r = 0.34 and r = 0.38, respectively) whereas senescence duration (SenDu) showed negative correlation with GY (r = -0.27) ( Table 4). In 2019, NDVI at anthesis ranged from 0.64 to 0.79 with mean of 0.72 and 0.46 to 0.66 with mean of 0.57 under IR and SA conditions, respectively. There was significant difference between genotypes and treatments for canopy green area per meter square at anthesis and NDVI at anthesis. G x T interaction was significant only for NDVI at anthesis (Table 3). In terms of correlation between GY and vegetation indexes, canopy green area per meter square at anthesis showed stronger correlation with GY than NDVI at anthesis and these correlations were stronger under SA (r = 0.56 and r = 0.55, respectively) than IR (r = 0.36 and r = 0.26, respectively) conditions (Table 5, Fig 3). Under SA conditions both canopy green area at anthesis and NDVI at anthesis showed a negative correlation (r = -0.39 and r = -0.31, respectively) with canopy temperature (Table 5; Fig 3). Canopy green area and NDVI after two and three weeks of anthesis was not correlated with GY.

Correlation between different yield components under IR and SA condition
Under IR conditions, the first principal component (PC1) explained 25.3% variation and the group of traits explaining this variation were fruiting efficiency, spikelet's per ear, days to heading with positive effect whereas TGW showed a negative effect. The second principal component (PC2) explained 20.8% variation and the group of traits explaining this variation included grains per ear, harvest index, grain yield, and root diameter with positive effects whereas above ground dry matter showed a negative effect.

PLOS ONE
Drought tolerance responses in winter wheat Under SA conditions, PC1 explained 25.6% variation and traits correlated with positive effect-were grains per ear, harvest index, grain yield, and NDVI whereas days to heading showed a negative effect. PC2 explained 21.0% of variation and traits showing correlation with positive effects were-thousand grain weight, plant height, root diameter, whereas spikelets per ear showed a negative effect (Fig 4).

Table 5. Correlation matrix showing correlation coefficient (r) values for grain yield (GY), above ground dry matter (AGDM), harvest index (HI), thousand grain weight (TGW), days to heading (DH), canopy green area per meter square at anthesis (GA An) and after 2 weeks of anthesis (GA 2W), NDVI at anthesis (NDVI), root angle (RoAng), root diameter (RoDiM), root dry weight per plant (RoDrWtPl), root number per shoot (RoNoSht), and canopy temperature at anthesis (CT). Below diagonal IR
Overall the bi-plots showed genotype code numbers (S1 Table) 11, 28, 7, and 29 under IR and 16, 39, 7 and 18 under SA conditions performed relatively better for GY than other genotypes (Fig 4). There was no grouping of genotypes that indicated strong genotype outliers for correlations with PC1 or PC2 under both treatments.

Stepwise regression analysis
A forward stepwise multiple-linear regression analysis was done to investigate which traits amongst root and shoot traits contributing the most towards improvement in the GY in the SA treatment (S6 Table). The stepwise regression analysis identified NDVI at anthesis (NDVI), root surface area (RoSuAr), and NDVI senescence duration (SenDu) in 2018 and green area at anthesis (GA An), root numbers per shoot (RoNoSh), and root angle (RoAng) in 2019 were the most important traits under drought conditions. These traits explained 17 and 41% of phenotypic variation contributing to GY in 2018 and 2019 respectively and the regression model was not improved by the addition of any further traits.

Marker-traits associations
Marker, allele/haplotype effect, mean of traits for each allele and probability of a significant difference in the mean under IR and SA conditions are presented in Table 6 for 2018 and 2019. In 2018, results showed that marker TaCwi-4A had a significant influence on GY under SA. This marker was responsible for increase in GY by 17.1% under SA in the presence of the Hap-4A-C allele which was associated with higher yield under drought conditions [37] and present in 59% of genotypes. Marker Dreb1 (Hap-I) was responsible for increase root surface area (RoSuAr) and root volume (RoVol) only under SA conditions. Dreb1 allele TaDreb-B1a increased these traits by 9.4% and 13.3%, respectively, under SA conditions. Dreb1 (Hap-I) was also correlated with extended SenEnd under SA conditions. The PRR73-A1 gene had a significant influence on GY under SA but not under IR conditions. There was an increase of

Discussion
In this discussion we address the potential application of the physiological traits and molecular markers examined in the field experiments in trait-based breeding for drought tolerance in wheat. We consider first the effect of phenology and canopy traits on the yield responses of genotypes to drought, then the effects of root traits and lastly the association between the molecular markers and responses to drought.

Grain yield responses to drought and association with phenology
In our experiments, comparing the two irrigation treatments, there was a moderately to severe drought with overall 3.06 t ha -1 (-47%) yield reduction under SA conditions. This is representative of Mediterranean drought effects in semi-arid conditions reported for wheat with reductions of yield typically ca. 30-50% [45]. The cultivars responded differently to the drought stress as indicated by the significant irrigation x genotype interaction. Higher yield under IR conditions was correlated with greater yield loss under SA conditions amongst the cultivars. From the physiological standpoint, it is not surprising that absolute reduction in yield for a given reduction in water resource is strongly influenced by yield potential. This is because higher yield potential genotypes will tend to use more water during the season and have higher biomass that low yield potential genotypes under optimal conditions [46][47][48].
Drought had only a small effect on heading date (GS59) advancing on average by one day in 2018 and three days in 2019; indicating genotypes responded similarly to drought stress (Fig 2). In both SA and IR conditions, HD ranged by 10 days amongst cultivars. Correlations between heading date and grain yield were negative under SA conditions in both years, and there was also a negative correlation under IR conditions in 2019 although it was less strong than under SA conditions. Bi-plots for the cross-year means also confirm these effects (Fig 4). This correlation among genotypes did not change even when considering just the median 30 genotype as a narrower flowering date group with a difference of 1 day in their flowering dates (S4 and S5 Tables). Early flowering has been correlated with drought escape in wheat in environments subjected to severe early season drought stress, e.g., in northern Mexico [46]. Similarly, Worland et. al. [49] reported increased yield for Ppd-D1a early-flowering NILs by ca. 5%  (14), TEF.7A-lower TGW (39) and high TGW (11), TaCwi.4A-high yield in drought (29) and low yield in drought (20), Dreb1-drought tolerant (20) and drought susceptible (29), ISBW11.GY.QTL.CANDIDATE (GY.QTL)-higher yield (27) and lower yield (23).
https://doi.org/10.1371/journal.pone.0242472.t006 compared to Ppd-D1b controls in dry years. Each day's advancement in HD raised yield by 0.11 t ha -1 . In the present study soil depth was more than 2 m at the field site with a very low organic matter. Presumably a shorter pre-anthesis phase reduced water uptake during this phase, so that season-long water uptake was redistributed more favorably with regard to the post-anthesis period; therefore water uptake during grain filling was increased with earlier flowering. However, in one year a similar negative correlation between HD and yield was recorded under irrigated conditions. This indicated that the negative association may have been partly correlated with advanced flowering leading to cooler prevailing temperatures during grain filling and therefore more calendar days for grain filling [50]. Brdar et. al. [51] showed that during grain filling, an increase in 1˚C mean daily temperature higher than optimum can be responsible for decrease in ca. 2.8 mg of grain weight. Observations regarding flowering time and drought resistance are very much dependent on the exact timing of drought stress and we recognize that the present experiments would need to be repeated over more years before we could conclude with certainty that later heading date overall has a negative effect on yield losses under droughts in Turkey. For example, it may be that there is a trade-off between early flowering and the development of a smaller root system, as suggested by Foulkes et. al. [52].

Correlations between canopy senescence traits and responses to drought
In the present study, greater yield amongst cultivars was correlated with higher green area per meter square and NDVI at around heading or anthesis under both drought and irrigated conditions. This was also confirmed by the stepwise multi-linear regression analysis which identified green area per meter square and NDVI as amongst the most important traits contributing to GY under drought conditions. Both of these are high-throughput measurements. This likely reflected a correlation between NDVI and biomass at anthesis and hence grains m -2 under both treatments. Genetic variation in GY in wheat has previously been correlated with green canopy area duration under drought in wheat [47,[53][54][55][56], and sorghum [57]. The role of senescence dynamics-start, end and rate of senescence-is important as they relate to grain filling duration and post-anthesis water and N uptake under abiotic stress [56,58]. Under irrigated conditions, GY was correlated positively with onset of senescence (SenSt) in 2018. Genotypes having delayed onset of senescence may be able to accumulate more plant nutrients and carbohydrates during grain filing duration resulting in higher yield. Our results suggested that there was source limitation if grain growth even under irrigated conditions. This could have been due to some heat stress incurred in the experiments combined with the higher grain number in irrigated conditions leading to source limitation during the later stages of grain filling [59].
One of the objectives of this experiment was to compare the NDVI and RGB-based vegetation indexes as methods for measuring canopy green area. We found that the RGB-based vegetation indexes showed better correlation with GY than NDVI measured with the handheld Trimble GreenSeeker and in addition it increased throughput. Better correlation of grain yield with RGB-based vegetation indexes than NDVI was also reported in durum wheat by Kefauver et. al. [27] and Fernandez-Gallego et. al. [28]. The role of a cooler canopy at anthesis correlated with a deeper root system for drought tolerance is previously reported [5] and was also indicated in this experiment with a negative correlation between GA at anthesis, NDVI at anthesis and GY with canopy temperature (CT) under SA conditions in 2019.

Effect of root traits on responses to drought
Shovelomics is a high-throughput phenotyping method for field-grown crops and has been used to quantify genetic variation in root traits in maize [19,21,22], legumes [60] and barley [23]. Maccaferri et. al. [24] carried out field shovelomics for durum wheat recombinant inbred lines (RIL) for crown root length, number and angle and reported QTL. In our study we applied a shovelomics methodology for bread wheat to quantify variation in nodal root angle, length, roots plant -1 and roots shoot -1 and correlation with GY. In the present study, the range in nodal root angle of 45-65˚in the semiarid treatment was similar to that of 42.3-69.2r eported by Maccaferri et. al. [24] for the Colosseo × Lloyd durum wheat mapping population assessed at anthesis in the field under optimum agronomic conditions in Italy. Values under irrigation of 45.0-70.0˚were also similar to those reported by Maccaferri et. al. [24] under irrigation. The correlation of shallower root angle with higher GY under irrigated conditions in 2019 was possibly correlated with increased root density at shallower depth and increased recovery of fertilizer N uptake which is predominately distributed in the top 30 cm of the soil. Shallow roots may also have increased rate of uptake of irrigation water leading to more transpiration and cooler canopies, so avoiding heat stress. Generally shallower root angle is correlated with shallower root distribution in soil and narrow root angle correlated with relatively deeper root distribution in soil [16]. In durum wheat, Hassouni et. al. [61] reported 20 to 40% yield advantage under irrigated conditions with the shallow root types compared to deeper root types.
Fan et. al. [62] reported 46.2% and 68.3% of wheat roots were distributed in the upper 15 and 30 cm, respectively. In contrast, a steeper angle would be expected to correlate with relatively deeper roots and greater yield under SA conditions, as has been reported in wheat in Australia [14][15][16]. In our experiments soil depth was more than 2 m. Previous studies in maize found steeper root angle related to increased rooting depth under low nitrogen field environments in the USA and South Africa [19]. However, in our results we did not see a significant correlation between root angle and GY under drought. Nevertheless, under SA conditions in 2019 there was a positive correlation between nodal root number per shoot and GY and HI, but no correlation under IR conditions. These results suggest that the wheat ideotype for drought tolerance may be a plant with relatively few tillers but a high number of nodal roots per shoot correlated with a longer residence time per root on average during the season and increased rooting depth. The multi-linear regression model also predicted the best combination of traits for drought-tolerance breeding to include root surface area, root number per shoot and root angle. A field study in Pennsylvania in maize found that reduced nodal roots per plant led to increased root length at depth and 57% higher grain yield under water-stressed conditions [63]. Tillering influences carbon partitioning, and there is some evidence that reduced tillering increases rooting depth in wheat [64,65] and rice [66].
In 2018 terminal stress was severe as rainfall was lower towards the end of crop season than in 2019 (Table 7). In 2018 under SA conditions genotypes which had less root surface area and root volume per plant produced higher yield and biomass. It is feasible that less surface area and/or volume of surface crown roots may have been correlated with more roots distributed relatively deeper in 2018 under SA conditions. Overall under drought conditions fine roots with smaller root diameter may have an advantage in exploring an increase area of soil relative to the energy invested to grow them [67]. In 2019, when terminal drought was less severe, this correlation with root surface area and volume was not observed.
Shovelomics is becoming an increasingly popular method for the high-throughput phenotyping of roots in field-grown crops. The shovelomics method we have applied for phenotyping nodal root traits in winter wheat in the present study were shown to be a valuable technique. In addition, the negative correlation of canopy temperature with GY indicates that water uptake at depth was contributing to yield increase. As the shovelomics methodology measures root traits in top 20 cm of soil profile, this technique requires further validation that it is a reliable indicator for roots at depth. Shovelomics is a relatively high-throughput method, making it possible in the present study to sample, wash and measure root crowns from 200 plots in two person-days. The shovelomics phenotyping platform is significantly faster than field soil coring, which would take approximately one person-month to sample, wash and extract roots, and image the samples for 200 plots in the present study. This high-throughput shovelomics platform can potentially be used to phenotype large populations to identify QTL, search for candidate genes and develop molecular marker for marker-assisted selection [68].
There are examples of deploying QTL for root depth in other cereal species. In rice, the Dro1 gene related to steeper crown root angles and deeper rooting was identified by measuring nodal root traits in a high-throughput controlled environment study [69] and has been used to produce drought tolerant NILs which have been phenotyped in field conditions [70].

Association between molecular markers and responses to drought
Genomic studies using high-throughput genotyping assays like KASP have made it possible to genotype large populations at various loci within a very short time [44]. Several recent studies used KASP markers to identify the allelic variation of functional genes in wheat cultivars from China [44], United States [71], and Canada [72]. In our study the clear association of TEF-7A and TaCwi-4A with GY and Dreb1 with root surface area (RoSuAr) under SA conditions indicated the usefulness of deploying these markers in wheat breeding for drought tolerance.
Dehydration responsive element binding proteins, Dreb1, have been shown to be induced by water stress, low temperature and salinity [73]. In this study, TaDreb1 was associated with increased root surface area, root volume and delayed end of senescence which indicated multitrait effects of this transcription factor which were not previously reported. The TaCwi-4A marker was correlated with GY, root shot ratio (Ro:Sh ratio) and PH under SA conditions. It was previously reported that storage carbohydrate accumulation in drought susceptible and tolerant cultivars depends on the expression of gene for cell wall invertase (TaCwi) in anthers [74]. The effect of drought on pollen fertility is irreversible and may cause grain loss or yield reduction under drought conditions. Since these genes tightly control sink strength and carbohydrate supply, deployment of favorable alleles of these genes could maintain pollen fertility and grain number in wheat. The drought tolerance and correlation of yield-related traits in CIMMYT winter wheat germplasm was strongly associated with TaCwi-4A which ultimately increased the grain sink size during drought stress. The selection of favorable allele for TaCwi-4A gene can enhance drought adaptability in marker-assisted breeding. Similarly, the association of several important traits like SenEnd, GY, HD, and GA with the flowering time related gene PRR73-A1 is interesting and indicated the expanded role of this plant development gene in senescence timing. Previously, Zhang et. al. [36] identified that PRR73-A1 was associated with plant height and explained up to 7.5% of the total phenotyping variability in Chinese wheats. This gene also showed association with plant height in this experiment in 2018. It was previously observed that flowering time related genes are very important for wheat adaptability in target environments, and these genes are correlated with several yield component traits [38]. Our results provided a set of target genes which could be manipulated to further finetune the expression of important drought-tolerance traits.

Conclusion
In the Mediterranean environment wheat is grown mostly under semiarid conditions and frequent drought affecting wheat yield severely. The strategy of developing drought-tolerant wheat varieties depends on understanding and identifying below-ground and above-ground traits for drought tolerance together with use of marker-assisted selection. In this experiment we used high-throughput root phenotyping techniques like shovelomics for characterizing root system architectural and RGB imaging based vegetation index for canopy senescence dynamic traits. We conclude that higher number of crown roots per shoot was a key trait for yield increase under drought conditions. Use of the RGB-based vegetation index to characterise the canopy green area dynamics could save time and increase precision in selection. The strong correlation of green area index with GY at flowering was encouraging and indicated this index can be used as tool for early stage selection for higher GY. In this study we have evaluated five established functional genes for drought-tolerance traits in field conditions that were previously reported elsewhere and for the first time we have validated them in Turkish wheat breeding lines. The genetic marker TaCwi.4A, responsible for drought tolerance was associated with higher GY in drought conditions in these experiments and could be used for future breeding. Our results also provided new insight on effects of root system architecture traits, for example, the importance of root angle under irrigated conditions and roots per shoot under drought for increasing grain yield which could be important for developing drought-tolerant cultivars.
Supporting information S1 Table. Code