Association mapping and genetic dissection of drought-induced canopy temperature differences in rice.

Drought-stressed plants display reduced stomatal conductance, which results in increased leaf temperature by limiting transpiration. In this study, thermal imaging was used to quantify the differences in canopy temperature under drought in a rice diversity panel consisting of 293 indica accessions. The population was grown under paddy field conditions and drought stress was imposed for 2 weeks at flowering. The canopy temperature of the accessions during stress negatively correlated with grain yield (r= -0.48) and positively with plant height (r=0.56). Temperature values were used to perform a genome-wide association (GWA) analysis using a 45K single nucleotide polynmorphism (SNP) map. A quantitative trait locus (QTL) for canopy temperature under drought was detected on chromosome 3 and fine-mapped using a high-density imputed SNP map. The candidate genes underlying the QTL point towards differences in the regulation of guard cell solute intake for stomatal opening as the possible source of temperature variation. Genetic variation for the significant markers of the QTL was present only within the tall, low-yielding landraces adapted to drought-prone environments. The absence of variation in the shorter genotypes, which showed lower leaf temperature and higher grain yield, suggests that breeding for high grain yield in rice under paddy conditions has reduced genetic variation for stomatal response under drought.

Introduction

The increasing variation in temperature, precipitation, and their interaction, resulting from global climate change is predicted to increase the variability in global crop yield by >30% (Ray ). Among the cereals, rice is especially sensitive to water limitation and heat stress, particularly during the reproductive stage (Jagadish ; Sandhu ). Climate change and the increasing probability of both prolonged and intermittent periods of drought are therefore likely to seriously affect rice production, particularly in rain-fed lowland farmlands which account for >30% of the world’s total rice cultivation area (Bailey-Serres ). Thus, plant breeders aim to develop varieties with improved yield performance under both favourable and water-limited conditions (Kumar ).

Crop germplasm collections stored in gene banks worldwide represent a large and potentially valuable reservoir of favourable alleles that can be used to develop new crop varieties that provide yield stability under both favourable and stressful environments (Tester and Langridge, 2010; Huang and Han, 2014). Over the last 10–20 years, rapid improvements in the throughput and cost-effectiveness of sequencing and genotyping have made it possible to generate extensive information about plant genetic variation at the genome level. This genomic information can be combined with phenotypic data for genetic analyses. The development of phenotyping tools has not progressed as rapidly, resulting in a ‘phenotyping bottleneck’ (Furbank and Tester, 2011; Cobb ) which limits the genetic dissection of complex traits such as tolerance to drought stress. However, new, non-destructive, non-invasive, image-based approaches to phenotyping in both the field and controlled environments are increasingly available, greatly enhancing the potential to phenotype large populations (Furbank and Tester, 2011; White ; Cobb ; Reynolds ). The use of indirect ‘proxy’ indicators for stress can be particularly useful and a powerful resource for field-based phenotyping (Jones, 2014). Among them, canopy temperature, measured by thermal imaging, has already proven to be a good indicator of drought stress in the field, as it indirectly measures stomatal conductance (Leinonen ; Munns ), one of the main physiological traits involved in the regulation of water loss (Schroeder ).

Several recent field studies successfully utilized infrared thermography to measure genotypic variation in stomatal conductance in a large number of genotypes (Jones ; Rebetzke ; Zia ; Prashar ; Rutkoski ). Critical to the success of such studies was the use of appropriate normalization techniques to overcome the environmental fluctuations (air temperature, humidity, wind speed, and incident radiation) that induce variation in canopy temperature during the process of imaging. The same studies also suggest that thermal image analysis of crop canopies is maximally effective in water-limited environments, as the genotypic differences in stomatal conductance are maximized under these conditions. As a consequence, thermal imaging provides a potentially useful phenotyping strategy for the selection of drought-tolerant genotypes (Jones ; Prashar ).

In the present study, we assessed the effectiveness of thermal imaging to quantify genetic variation in canopy temperature/stomatal conductance in tropical rice, using a population of 293 indica accessions grown in the field under control and drought conditions at the International Rice Research Institute (IRRI) in the Philippines. Statistical analyses revealed a relationship between canopy temperature during flowering, plant height, and grain yield. We also report what is, to our knowledge, the first genome-wide association (GWA) analysis of leaf temperature in rice, demonstrating that there is genetic variation for this trait, and pinpoint genomic loci and a priori candidate genes that underlie this variation.

Materials and methods

Description of the field experiment

A population consisting of 293 accessions of Oryza sativa subsp. indica was used in a field trial experiment at the IRRI, Los Baños, Philippines (14°11'N, 121°15'E; elevation 21 m above sea level) during the 2014 dry season. The field trial was located at the IRRI upland farm where the soil profile does not represent a constraint for deep penetration of the roots. The accessions are largely the same as those in the PRAY indica panel (http://ricephenonetwork.irri.org) which includes traditional and improved indica rice lines originating from rice-growing countries in tropical and subtropical regions around the world. The same panel was recently used in studies where a number of diverse traits were phenotyped as the basis for genome-wide association studies (GWAS) (Qiu ; Al-Tamimi ; Rebolledo ; Kadam , 2018; Kikuchi ). The field trial was carried out in two separate fields, one that served as control and the other for the drought stress treatment. Each field comprised three replicates of the population (A, B, C for control and D, E, F for drought) arranged in a serpentine design (see Supplementary Fig. S1 at JXB online). Each replicated accession consisted of 48 plants covering a 2.5×0.8m area and arranged as four rows of 12 plants each. To manage the differences in flowering phenology, the accessions were sown at 7 d intervals and transplanted to create subgroups that allowed us to synchronize flowering. Eight subgroups were created according to the number of days required to reach 50% flowering. Each group was progressively sown and transplanted into the field with an interval of 7 d between each group. Inside each subgroup, accessions were transplanted from the shortest to the tallest to minimize the positioning of short and tall genotypes next to each other (plant height data collected during the dry season 2013). The imposed drought stress treatment consisted of 14 consecutive days of water withholding applied only to the stress field at the reproductive stage (targeting 50% flowering). Weather data were collected during the entire experiment by a weather station located in the middle of the two fields. Twenty-six tensiometers were randomly distributed over the stress field to record soil water potential. At the end of the stress period, the stress field was re-watered until all accessions reached the maturity stage for harvest. The control field was constantly kept in a flooded condition (paddy field). At harvest (on average 30 d after re-watering), the following traits were scored for all replicates: plant height (cm), grain yield (g m–2), shoot biomass (g m–2), and harvest index (the ratio between grain yield and total biomass). The dates of initial flowering, 50% flowering, and 100% flowering were also recorded for replicated trials of each accession under both treatments. A more detailed description of the experiment, including the description of the same field trial conducted during the 2013 dry season, can be found in Kadam .

Thermal imaging

A FLIR B660 (FLIR systems, USA) infrared camera was used for taking both infrared and visual images. The thermal camera is assembled with a focal plane array (FPA) uncooled microbolometer that operates in the spectral range of 7.5–13 μm with a resolution of 640×480 pixels. The thermal camera is also equipped with a digital camera with a resolution of 3.2 Megapixels. All pictures were taken from 3.5 m height (Jones ) with each image covering ~50 m2 (Fig. 1; Supplementary Fig. S1). The distance between the camera and the centre point of the field in the image was kept constant, resulting in a camera angle of ~20° from the ground. Thermal pictures were taken, during the morning, 8, 9, and 10 d after the stress was applied (from 2 to 4 April 2014). We collected images to fully cover ‘Rep B’ (control), ‘Rep E’, and ‘Rep F’ (drought). For each replicate, we collected images on two consecutive days at two different times during the morning period. To image an entire replicate, it took on average 45 min. ‘Rep B&E’ were covered by 18 pictures and ‘Rep F’ by 27–28 (Supplementary Fig. S1; Supplementary Table S1A, B). Camera settings were kept constant during the entire process of imaging with atmospheric temperature set to 30 °C and emissivity set at 0.95 according to Jones and Prashar and Jones (2014, 2016).

Fig. 1.

Plot identification and picture analysis. Example of a digital picture with the corresponding thermal picture area as indicated by a rectangle (A). Thermal picture used for plot identification (B). Selection of specific polygonal areas for the quantification of the genotypes’ temperature (C). Temperature of the polygonal areas selected in (C) (D).

Plot identification and picture analysis

Plot identification was achieved following the experimental design and by the use of three T-boards placed at known positions in each image (Fig. 1). In addition, every plot in the field trial was marked by a stick placed between two consecutive plots. The stick was characterized by a relatively higher temperature than the plant canopies and it was visible in the thermal images (Fig. 1B, C). Temperature quantification was performed by loading the images into the ThermaCAM Researcher Professional 2.10 software (FLIR systems), selecting a rectangular area for each plot canopy, and using the mean temperature of the pixels in the enclosed rectangular area as representative for the specific plot (Fig. 1C, D). The temperature of each T-board reference surface (black and white) was determined in the same way.

Plot image normalization

We considered three methods to normalize plot temperatures. In the first method (‘image mean’), plot temperatures in each image were multiplied by the ratio between the mean temperature of all the plots in the replicate and the mean temperature of all the plots in the image. The normalized temperature of plots occurring in two consecutive images was calculated as the mean of the two resulting values (Prashar ). The second and third methods are based on the same procedure but using the mean temperature of the reference surfaces (‘white reference’ and ‘black reference’) to calculate the ratios. Correlation analysis between normalized data was performed to evaluate which normalization method produced the highest degree of reproducibility for the same replicate imaged over two consecutive days.

Statistical analysis

Statistical analysis of the data was conducted by using R statistical software (version 3.4.3; The R Foundation for Statistical Computing). Correlation analysis and graphical matrices were produced using a modified function of the ‘corrplot’ R package. Box-Cox transformation of traits that were not normally distributed (Shapiro-Wilk’s P-value <0.05) (replicates B03, E02, and F04) was calculated using the ‘forecast’ R package. Single nucleotide polymorphism (SNP)-based principal component analysis (PCA) was performed using the prcomp function in the ‘stats’ R package.

GWA analysis with the 45K SNP map

GWAS of normalized trait values (Supplementary Fig. S2) were performed using a linear-mixed model in EMMAX (Kang ), which corrects for population structure by including a kinship matrix (IBS matrix) as covariate. EMMAX also provides an estimate of the phenotypic variance (heritability estimated by markers, h2) explained by the IBS matrix. Of the 293 accessions in the field experiment, 271 matched the original panel (indica PRAY panel of 339 accessions) that was used to generate a 47K SNP map using genotyping by sequencing (GBS). The 47K SNP GBS map displayed 8.75% missing data that were imputed by Fast Phase Hidden Markov Model (Scheet and Stephens, 2006) as reported by Kadam , 2018), thus reducing the missing rate to 0%. The reduced number of accessions (271) altered the minor allele frequency (MAF) threshold of the 339 accessions panel, originally set at 0.05. To exclude rare alleles from the present study (n=271 accessions), the 47K SNP map was re-filtered for MAF >0.05, resulting in 45 505 SNPs available for GWAS. PCA based on the 45 505 SNPs was conducted to quantify subpopulation structure (Fig. S3). The main component (PC1) explained only 8.72% of the genetic variation, and a combination of the first three components failed to clearly separate groups of accessions. Therefore, no PC covariates were added to the linear-mixed model to correct for subpopulation structure, following the approach used in McCouch . GWAS results are presented as Manhattan and Quantile–Quantile plots using the ‘qqman’ R package.

To avoid Type 1 errors (false positives), the significance threshold used to identify marker–trait associations was set at P<0.00001 (i.e. –log10P>5.0). This significance threshold was higher than –log10P>4.0, commonly adopted in other GWAS on rice using SNP maps of similar density (Zhao ; Dimkpa ; Kadam et al., 2017, 2018), but lower than the Bonferroni threshold (–log10 P>6.0, for α=0.05), which was considered too stringent because of the increasing occurrence of Type 2 errors (false negatives). Using a threshold of –log10P>5.0 for individual analyses, discovery of the same marker–trait associations across experiments provided additional support for the associations.

GWA analysis using a high-density imputed SNP map

Imputation of the 47K SNP map was conducted using the Rice Imputation Server (RIS) following Wang . The map of the 47K SNPs (.hmp format) was first converted to Plink format (.ped/.map) and then to Oxford format (.gen/.sample) before being uploaded as a compressed folder (.tar.gz) in the RIS (http://rice-impute.biotech.cornell.edu). The RIS-imputed map was downloaded as a single Plink file (.bed/.bim./.fam format) and divided into 12 individual chromosomes. Focusing only on chromosome 3, missing SNPs were imputed with Beagle version 4.1 (Browning and Browning, 2007). Finally, the Beagle-imputed map of chromosome 3 was filtered at MAF >0.05, resulting in a set of 186 012 SNPs available for mapping on this chromosome (4039 SNPs with MAF >0.05 were present on chromosome 3 in the 45K SNP map). GWA analysis using the chromosome 3 imputed map was conducted as described above for the 45K SNP map. The IBS matrix of kinship used as a model covariate was calculated based on the 45K GBS SNP map (Wang ).

Linkage disequilibrium (LD) analysis and a priori candidate gene selection

The local pairwise LD pattern near the significant SNPs was calculated and graphically represented by the ‘snp.plotter’ R package (Luna and Nicodemus, 2007). The annotations of genes located within LD blocks were obtained from the MSUv7 rice genome database (http://rice.plantbiology.msu.edu/). The exact localization and functional annotation of significant SNPs was conducted using SNPEff version 2.05 (Cingolani ) with MSUv7 as the reference genome.

Results

Thermal imaging and data normalization

Fluctuations in environmental conditions are the main obstacle to the use of thermal imaging to reliably analyse plant canopy temperature. Supplementary Table S1B shows the changes in air temperature, humidity, wind speed, and solar irradiance measured on the days and in the time windows during which the thermal imaging was performed. To reduce the overall effect of these factors on the analysis of canopy temperature, we imaged the field replicates only during the mornings on three consecutive days. Mornings were selected for imaging due to the sharp increase in wind speed experienced every afternoon in the field location, and previous reports describing wind as a major factor strongly impacting stomatal conductance values (Jones, 1999; Maes and Steppe, 2012). Nevertheless, our data documented variation caused by environmental fluctuations between images (Supplementary Fig. S4) and therefore normalization was needed (Prashar and Jones, 2014). We applied three different procedures to reduce the variation caused by environmental fluctuations (‘image mean’, ‘white reference’, and ‘black reference’) (see the Materials and methods). By comparing the same field replicate imaged over two consecutive days, we found that ‘image mean’ normalization produced higher correlation values (B03–B04, from 0.075 to 0.37; E02–E03, from 0.23 to 0.69; F03–F04, from 0.25 to 0.65) than normalization based on white and black references (Supplementary Fig. S5). Hence, ‘image mean’ normalization was used in all subsequent analyses.

Canopy temperature

Drought stress induced a strong increase in canopy temperature. The mean value of the normalized stressed replicates E and F (EF mean) was 2.27 °C higher than that of the control B replicate (B mean) (Fig. 2). Together with the treatment, the time of day strongly impacted leaf temperature. Control replicate B03, which was imaged in the early morning, showed a lower canopy temperature range (mean difference –0.90 °C) than B04 which was imaged in the late morning (~2 h later) the following day. As with the control replicates, drought replicate E02 showed a lower temperature (mean difference –1.85 °C) than E03. The temperature difference between the two E replicates was due to the combined effect of time of day (E02 was measured earlier in the morning than E03) and an additional day of drought stress. Tensiometer readings showed that the soil water potential of the drought field decreased sharply during the days of imaging, moving from an average of –34 kPa on 2 April to –53 kPa on 4 April (Supplementary Fig. S6). Of the two F replicates, F03 was imaged half an hour earlier than F04 but the latter was exposed to one more day of stress, resulting in a higher temperature for F04 (mean difference 0.33 °C) than for F03 (Fig. 2).

Fig. 2.

Canopy temperature of the different field replicates. Boxplots representing the mean temperatures of control (B mean) and drought (EF mean) replicates (on the left), and of the six single replicates (on the right), together with the information on the days and time windows of picture taking (in the middle). White and grey boxplots represent control and stress values, respectively.

Relationships between canopy temperature and agronomic traits

Canopy temperature was measured in the second of two years (2013 and 2014 dry seasons) in which a field experiment was conducted to collect information on phenotypic trait performance of all 293 rice accessions evaluated under well-watered and drought stress conditions (Kadam ). Plant height (PHT), grain yield (GY), shoot biomass (Shoot bio.), and harvest index (HI) were among the agronomic traits recorded at the time of harvest in both years. A ‘flowering’ variable (FLW) was calculated during the stress period of the 2014 field trial only, by subtracting the date of 50% flowering for every genotype in each replicate from the date of thermal imaging. Drought stress similarly affected all traits in both years (Supplementary Fig. S7) and particularly reduced GY and HI, but minimally affected Shoot bio. and FLW (2014 only). To investigate the relationship between canopy temperature at the time of stress exposure and plant traits at harvest time, we conducted a Spearman correlation analysis between normalized temperature and agronomic trait values (Fig. 3). Canopy temperature showed significant (P<0.001) and negative correlations with GY and HI, and positive correlations with PHT, Shoot bio., and FLW, under both conditions, but with higher correlation coefficients under stress than control; FLW was not significantly (P>0.001) correlated with canopy temperature under drought (Fig. 3A, B). Remarkably, almost identical correlations were found between canopy temperature measured in the 2014 dry season and the agronomic traits scored during the field trial conducted in the 2013 dry season (Fig. 3C, D). In both years, the highest correlation coefficients were found between the mean values of canopy temperature under drought stress and PHT (r=0.56 in 2014; r=0.58 in 2013), GY (r= –0.48 in 2014; r= –0.38 in 2013), and HI (r= -0.53 in 2014; r= –0.54 in 2013). For these three traits, the percentage of variance explained by the linear models (R2) associating temperature and trait was equal to 34% for PHT in both years, and 25% and 16% for GY, and 30% and 28% for HI in 2014 and 2013, respectively (Supplementary Fig. S8). Overall, these results show that, under drought stress, thermal imaging of rice canopies at flowering time can detect canopy temperature differences that correlate with plant performance at the time of harvest. Furthermore, the almost identical correlations between the agronomic traits scored during the 2013 field trial and canopy temperature measured in 2014 indirectly validate the robustness of the temperature results across two seasons of field trials.

Fig. 3.

Correlation analysis between mean values of canopy temperature and agronomic traits. Mean canopy temperature values selected under control (CON) and drought stress (DRO) conditions were determined in the 2013 dry season, while agronomic traits were determined in both field trials conducted during the 2013 and 2014 dry seasons. Correlation values were determined by Spearman’s correlation analysis. All the correlations coloured in blue and red are significant at P<0.001. Correlations coloured in white have P>0.001. PHT, plant height; GY, grain yield; Shoot bio., shoot biomass; HI, harvest index; FLW, flowering; ‘B mean’, mean temperature of control field replicates; ‘EF mean’, mean temperature of drought field replicates.

GWAS and LD analysis using the 45K SNP map

The results described above demonstrate the effectiveness of thermal imaging in detecting quantitative differences in canopy temperature. Hence, we decided to try to use canopy temperature as a trait for association mapping analysis. Considering the strong influence of time of day and the day itself on canopy temperature (Fig. 2), the values of the separate field replicates were analysed, in addition to the mean temperature values of control and drought replicates. GWA mapping was conducted using a 45K SNP map and a stringent threshold of genome-wide significance (–log10P>5.0) to detect only highly significant marker–trait associations. Quantile–Quantile plots relative to the GWA analyses are reported in Supplementary Fig. S9.

Mean temperature variation of control replicates (B mean) as well as the temperature of separate control replicates (B03 and B04) was characterized by a low fraction of heritability estimated by markers (h2=0.17–0.24) (Supplementary Table S2) and, indeed, GWA analysis did not find marker–trait associations above the threshold of significance (Supplementary Fig. S10). Low heritability levels and absence of significant marker–trait associations suggest that no major genetic determinants are responsible for canopy temperature variation under control conditions in this panel.

Mean canopy temperature values under drought (EF mean) showed higher marker heritability (h2=0.5) than under control conditions (Supplementary Table S2). Marker heritability was also higher for canopy temperature evaluated on separate drought replicates compared with control replicates (Supplementary Table S2). GWA analysis identified two significant markers associated with ‘EF mean’ (Fig. 4A). The markers (SNP_12262251 and SNP_12529189) were located 267 kbp apart on chromosome 3, showed a similar level of significance (–log10P>5.0), and minor alleles at both loci were associated with higher canopy temperature than the major alleles (Supplementary Table S3). GWA analysis for F03, the field replicate with the highest heritability (h2=0.45), identified five markers significantly associated with canopy temperature (Fig. 4B). Both SNP_12262251 and SNP_12529189 were significantly associated with canopy temperature for both ‘EF mean’ and F03 (Fig. 4A, B), and minor alleles at all the five markers associated with F03 were associated with higher canopy temperature (Supplementary Table S3).

Fig. 4.

Manhattan plots of the significant GWA mapping results using the 45K SNP map. Manhattan plots of the GWA mapping results using the 45K SNP map for mean canopy temperature under drought (EF mean) and for the single field replicate F03. The red dashed line indicates the genome-wide threshold for significant associated markers (–log10P>5.0). SNPs above the red line are highlighted in red.

To exclude the possibility that the significant marker–trait associations detected under drought might have been influenced by genome-wide subpopulation structure among the 271 indica accessions, we re-conducted GWA mapping, adding the main marker-based PC (PC1, which explains 8.72% of the total genetic variation) as a covariate in the linear-mixed model (Supplementary Fig. S3). The addition of the covariate did not affect the outcome of the GWAS (Supplementary Fig. S11). These results suggest that genome-wide subpopulation structure among the accessions has a negligible influence on the chromosome 3 marker–trait association detected for ‘EF mean’ and F03. We therefore report results for the model run without covariates.

The higher heritability of canopy temperature under drought suggests a stronger genetic control for the trait under this stress condition. Mapping results indicate that the marker–trait associations detected on chromosome 3 are robust, with an estimated mean allele effect on ‘EF mean’ canopy temperature of +0.61 °C for SNP_12529189 and +0.57 °C for SNP_12262251 (Supplementary Table S3). The effect of these two alleles is consistent across all the individual drought field replicates (Supplementary Fig. S4).

Local LD was examined across an ~1.2 Mbp region encompassing ~500 kbp upstream and downstream of the two significant markers (SNP_12262251 and SNP_12529189). The pairwise LD estimates (r2) of the 102 SNPs in this region revealed that the two significant markers map to different LD blocks (r2>0.6–0.7) (Supplementary Fig. S12A). It is noteworthy that, using the 45K SNP map, SNP_12262251 maps to its own LD block, with the closest upstream marker (SNP_11994173) located ~270 kbp away; these two markers show a very low value of pairwise LD (Supplementary Fig. S12B). A low density of GBS markers around the two significant SNPs does not allow a precise determination of the LD configuration, and leads to a likely underestimation of the size of the LD block containing SNP_12262251.

Fine mapping of chromosome 3 using a high-density imputed SNP map

To improve the mapping resolution, GWA mapping of ‘EF mean’ canopy temperature values was performed again only for chromosome 3 using a high-density SNP map generated by haplotype-based genetic imputation using the RIS (Wang ). Following imputation, the number of SNP markers available for mapping of chromosome 3 increased 46×, from 4039 (on the 45K SNP map) to 186 012 on the high-density imputed map, and the imputed map drastically improved the mapping resolution. The number of significant SNPs (–log10P>5.0) associated with the ‘EF mean’ canopy temperature increased from 2 (45K map) to 65 (Fig. 5A; Supplementary Table S5). All newly imputed markers localized in the same region as the original significant markers (SNP_12262251 and SNP_12529189) identified using the GBS map (Supplementary Table S5).

Fig. 5.

Manhattan plots of the GWA mapping results for chromosome 3 using the imputed map. Manhattan plot of GWA mapping results for chromosome 3 using 186K SNPs for ‘EF mean’ canopy temperature. SNPs above the red line are highlighted in red (A). Zoom in on the region of the significant SNPs. The red arrows indicate the two QTLs (qCT1 and qCT2). SNPs of the GBS map are highlighted in light blue (B). In both figures, the red dashed line indicates the genome-wide threshold for significant associated markers (–log10P>5.0).

Zooming into the region of significance (12.1–12.7 Mbp) on the high-density map showed two distinct quantitative trait loci (QTLs), qCT1 and qCT2 (Fig. 5B). Previously, qCT2 could not be clearly detected due to the lack of markers around SNP_12262251 (Supplementary Fig. S12B). The significant SNPs of both QTLs were characterized by minor alleles associated with higher canopy temperature values (Supplementary Table S5). In the imputed high-density map, the top SNP associated with qCT1, SNP_12523765 (–log10P=6.26), was 5.4 kbp away from the top SNP identified using the GBS map (SNP_12529189, –log10P=5.19). In qCT2, a series of 15 SNPs showed the same significance (–log10P=5.57), and SNP_12267148 was closest (4.9 kbp away) to the top SNP on the GBS map, SNP_12262251 (–log10P=5.10) (Fig. 5B; Supplementary Table S5). Thus, the imputed map supports the location of the QTL identified on chromosome 3 using the 45K SNP map, but provides improved resolution, making it possible to differentiate two, closely linked genomic associations.

To determine if the two QTLs were independently associated with canopy temperature, we again performed GWA mapping, but this time we fitted the most significant SNP of qCT1, and subsequently of qCT2, as covariates in the linear-mixed model. Fixing either SNP_12523765 (qCT1) or SNP_12267148 (qCT2) as a model covariate yielded a similar loss of signal from all markers in the region (–log10P<3), causing both QTLs to disappear (Supplementary Fig. S13). This suggests that a lack of recombination between the two QTLs in our panel of indica varieties may limit our ability to genetically dissect this complex QTL region.

Haplotype analysis of qCT1–qCT2 and their relationships to agronomic traits

To test this hypothesis, we analysed haplotype variation across the qCT1–qCT2 region on chromosome 3 to determine whether we could identify recombinant genotypes, and, if so, whether the recombinants were significantly associated with canopy temperature. We were also interested to determine whether recombinants were significantly associated with two agronomic traits, PHT and GY, that were, themselves, significantly correlated with canopy temperature (Fig. 3).

To conduct the haplotype analysis, we selected SNPs from the high-density map that showed the most significant association with canopy temperature for ‘EF mean’ (3 SNPs for qCT1 and 15 SNPs for qCT2) plus the two most significant GBS SNPs (SNP_12529189 and SNP_12262251, in qCT1 and qCT2, respectively). Among the 271 accessions in the panel, we identified three major recombinant haplotype groups (present in >5% of accessions) (Fig. 6A; Supplementary Table S6). The most common haplotype, Haplotype I, was found in 206 accessions (~76% of the panel) and carried major alleles at all 20 SNPs across the two QTLs. These accessions displayed a mean canopy temperature under drought of 34.24 °C (Fig. 6A). Haplotypes II and III were found in 22 (~8%) and 18 (~7%) accessions, respectively. Haplotype III was the mirror image of Haplotype I; it carried minor alleles at all 20 SNPs, and accessions carrying this haplotype displayed a mean canopy temperature under drought of 35.05 °C (Fig. 6A). Haplotype II was a recombinant between Haplotype I and Haplotype III; it carried minor alleles for the SNPs in qCT1 (like Haplotype III) and major alleles for the SNPs in qCT2 (like Haplotype I). Accessions carrying Haplotype II displayed an intermediate mean canopy temperature under drought of 34.85 °C (Fig. 6A).

Fig. 6.

Haplotype variation across the qCT1–qCT2 region and relationship with agronomic traits. Haplotypes for qCT1 and qCT2 identified among the accessions of the panel. (A) The 20 most significant SNPs for the two QTLs are listed together with the information relative to the minor and major allele (min./maj.), physical position (bp), and significance (–log10P). Minor alleles are coloured in green and major alleles in brown. For each haplotype (I, II, and III) are reported the number of lines carrying the haplotype and their mean value of canopy temperature (‘B mean’ and ‘EF mean’), plant height (PHT), and grain yield (GY) under control and drought stress conditions. Boxplots representing the range of variation of mean canopy temperature (A), plant height (B), and grain yield (C) for the accessions of the three haplotypes (I, II, and III) under control and drought stress conditions. White and grey boxplots represent control and stress values, respectively. *, **, *** indicate the level of significant difference (P<0.05, P<0.01, P<0.001) between the groups.

We next compared the phenotypic performance of accessions carrying Haplotype II with the performance of those carrying Haplotype I and Haplotype III for canopy temperature, PHT, and GY (Fig. 6B, C, D). These comparisons were performed using Welch’s t-tests to accommodate the unbalanced sample sizes of the haplotype groups. A highly significant (P<0.001) canopy temperature difference under drought (EF mean) was detected between accessions carrying Haplotype I and accessions carrying Haplotype II (mean difference 0.61 °C) and Haplotype III (mean difference 0.81 °C), while no significant difference was detected between accessions carrying Haplotype II and Haplotype III (mean difference 0.20 °C) (Fig. 6B). Given that Haplotype I accessions differ from Haplotype II and Haplotype III accessions at qCT1, but not at qCT2, this suggests that the difference in canopy temperature is driven by gene(s) underlying qCT1. Similar results, but less significant (P<0.01), were found for canopy temperature under control conditions (B mean) (Fig. 6B). Under control conditions, the three haplotype groups displayed lower temperature differences (mean of Haplotype I=32.02 °C, Haplotype II=32.39 °C, and Haplotype III=32.29 °C) than under drought stress (Fig. 6A). We detected an almost identical pattern of significant differences between the three haplotype groups for PHT under both treatments and GY under drought (Fig. 6C, D). To further explore the relationship between haplotype groups, PHT, and GY under drought, we examined the mean PHT of accessions carrying different haplotypes and observed that those carrying Haplotype I were ~20 cm shorter than accessions carrying Haplotype II and Haplotype III (Fig. 6A).

Among the 48 accessions shorter than 100 cm (based on PHT in control conditions), 47 carried Haplotype I and only a single accession carried Haplotype II (cv Binulawan from Philippines); while none carried Haplotype III (Supplementary Table S6). The short, semi-dwarf accessions were essentially fixed for Haplotype I and thus carried only major alleles at qCT1 (associated with lower canopy temperature), while tall accessions carried all three haplotypes. Finally, we considered the differences in GY for accessions carrying different haplotypes (Fig. 6D). Under control conditions, GY for accessions carrying Haplotype I (mean 475 g m–2) was marginally greater than for those carrying Haplotype II (mean 382 g m–2; P<0.05) but was not significantly different from accessions carrying Haplotype III (mean 417 g m–2) (Fig. 6D). Under drought stress, GY differences between Haplotype I (mean 296 g m–2) and the other two groups (Haplotype II mean 169 g m–2 and Haplotype III mean 187 g m–2) became highly significant (P<0.001) (Fig. 6D). This suggests a negative effect of minor alleles at qCT1 for GY under stress in our experimental conditions.

LD analysis of qCT1 and a priori candidate gene identification

To identify candidate genes underlying qCT1, we first explored the extent of LD surrounding the most significant SNPs associated with ‘EF mean’ canopy temperature (despite the fact that we determined that the difference in canopy temperature is driven by qCT1 we provide the same LD analysis for qCT2 in Supplementary Fig. S14).

Using the high-density map, we first identified SNPs showing values of association of –log10P>2 located 125 kbp upstream and downstream of the qCT1 most significant marker, SNP_12523765 (Fig. 7A). Pairwise LD estimates (r2) of the 418 SNPs in this region showed that all the SNPs in qCT1 form a single, large LD block characterized by r2≥0.5 (Fig. 7A). While there are many markers in the region that do not meet the significance threshold (–log10P>5) for association with ‘EF mean’, all map within a sub-block characterized by high pairwise LD (r2>0.8). Notably, the 20 most significant SNPs with ‘EF mean’ (shown in red in Fig. 7A) are localized between two highly significant markers, SNP_12523765 and SNP_12460502. Thus, we considered the region delimited by these two markers (~42 kbp; indicated by the dashed black line in Fig. 7A) as the most interesting target for the identification of a priori candidate genes.

Fig. 7.

Localized linkage disequilibrium analysis of qCT1. Manhattan plots displaying the level of significance (y-axis) over genomic positions (x-axis) in a window of 125 kbp upstream and downstream of the marker (12523765) most significantly associated with canopy temperature of ‘EF mean’ and located on chromosome 3 in qCT1 (A). Localized region (zoom in, black dashed triangle in A) showing the genes (black arrows) underlying the most significantly associated (–log10P-value >5.0) markers’ loci (B). Different colours are used to represent the pairwise LD estimates (r2) for each genomic location. Genomic locations of the most significant markers are projected on the LD matrix and on gene positions (only in B) by red lines. SNPs of the GBS map are highlighted in dark blue.

Seven genes are included in this region (Fig. 7B; Supplementary Table S7). Among these seven, two were considered interesting candidate genes for their possible role in physiological processes regulating stomatal function. One gene is known to be responsive to abiotic stress (GOSlim terms) and encodes a putative mitochondrial fumarate hydratase (LOC_Os03g21950). This gene is located between SNP_12545362 and SNP_12560502, two of the three most significant markers in qCT1 (Fig. 7B). Fumarate hydratase (fumarase) is responsible for the conversion of fumarate to malate, a solute involved in the mechanism of stomatal opening/closure. The other interesting gene in this region (LOC_Os03g21890) encodes a plasma membrane high-affinity potassium (HAK) transporter (Bañuelos ) and is located in close proximity to the most significant marker (SNP_12523765) within qCT1 (Fig. 7B). HAK transporters are involved in guard cell K+ flux that controls stomatal opening/closure (Jezek and Blatt, 2017). Finally, we considered the predicted function of each of the 20 significant SNPs (–log10P>5) in this region to determine whether any of them could result in an amino acid change or a putative regulatory change affecting a specific gene candidate (Supplementary Table S8). None of the SNPs was associated with predicted non-synonymous mutations that could point towards a particular candidate, but many are located upstream of the gene coding region and thus potentially associated with changes in regulation of gene expression, including the most significant one, SNP_12523765, which is located 711 bp upstream of the HAK transporter.

Discussion

Normalization and physiological implication of canopy temperature results

The main determinants of canopy temperature in plants include genetic components affecting stomatal aperture and canopy structure (which may also affect aerodynamic resistance and radiation interception) and a range of environmental factors. The ability to reduce environmental fluctuations in humidity, irradiance, and wind speed is key to screening for true genetic variation in stomatal conductance (Prashar ). Jones demonstrated that the variation in canopy temperature between different rice genotypes can be detected in field experiments by thermal imaging with the use of appropriate normalization techniques. In this study we took this approach to the next level by screening a diversity panel consisting of 293 indica accessions. To date, only a few studies have followed a similar extensive approach using thermal imaging of plant canopies under field conditions (Prashar ; Zia ; Rutkoski ). Normalization of the raw canopy temperature data by ‘image mean’ reduced the influence of environmental factors. This is clearly shown by the increased correlation between the normalized canopy temperature data of the same field measured over two consecutive days (Supplementary Fig. S5A, B). Our results show that drought stress strongly increased canopy temperature. Other factors that may significantly affect canopy temperature are the time window and day of imaging and, for the stressed replicates only, the drought exposure time (Fig. 2). These results confirm that canopy temperature determined by thermal imaging is a reliable field proxy for stomatal conductance (Jones, 2014; Prashar and Jones, 2014), and that this trait is characterized by a very dynamic response to changing environmental conditions (Vico ; Drake ). The results of this study also suggest that water limitation reduced this dynamic response, as evidenced by the higher correlation coefficients under drought between the temperatures for the same field replicate imaged at different moments (Supplementary Fig. S5B).

Canopy temperature and agronomic traits at harvest time

Genotypic variation in stomatal conductance in rice may be responsible for differences in photosynthesis, even under optimal growing conditions (Ohsumi ; Ouyang ). Reduction of stomatal conductance is a well described physiological response to drought stress in rice (Centritto ; Ji ). The resulting limitation in leaf CO2 diffusion has been shown to cause grain yield reduction in rice genotypes grown under water-limited conditions with the stress imposed at the flowering stage in the field (Centritto ; Lauteri ). There is evidence—obtained with a limited number of genotypes—that canopy stomatal conductance monitoring by thermal remote sensing at the flowering stage could be an effective criterion for the selection of high-yielding rice genotypes (Horie ). In the present study, canopy temperature measurements were used to screen 293 rice accessions for stomatal conductance differences at anthesis and to explore correlations between the genotypic variation in canopy temperature and several plant agronomic traits.

Canopy temperature was strongly and positively correlated with plant height, and an equally strong negative correlation was found with grain yield and harvest index, particularly under stress conditions (Fig. 3; Supplementary Fig. S8). It is interesting to note that very similar correlations were documented between canopy temperature and agronomic traits scored in both years of field trials (Fig. 3), though canopy temperature was measured in only one year. The consistency of these correlations indirectly reinforces the effectiveness of our normalization procedure in reducing the influence of environmental factors on canopy temperature results. It can therefore be assumed that the effective detection of true genotypic differences in canopy temperature in one year could be similarly correlated with the performance of the accessions for agronomic traits, such as plant height and grain yield, over the 2 years of field trials.

The fact that the correlations between canopy temperature and agronomic traits under drought were stronger than under control conditions supports the idea that stomatal conductance is more important for plant performance under stress than under optimal conditions. Similar results have been reported in previous large field studies in other C3 and C4 cereals. For example, a similar negative correlation between canopy temperature in the reproductive stage and grain yield was described for segregating bi-parental wheat populations (Saint Pierre ). Zia also found a negative correlation between canopy temperatures at anthesis, and grain yield in 150 maize single cross-hybrids under water-limited conditions.

On the other hand, the positive correlation between plant height and thermal data found in this study (Fig. 3) is contrary to what has been observed for other crops, including wheat (Giunta ; Rebetzke ) and potato (Prashar ). In those studies, a negative correlation was observed between canopy temperature and plant height, and it was interpreted in terms of an atmospheric temperature profile where an increased aerodynamic resistance in the shorter genotypes was responsible for their higher leaf temperature (Rebetzke ). It is unclear why this does not apply in our rice trial, but it might be that the taller genotypes really do have more closed stomata than the shorter genotypes, with this effect over-riding any aerodynamic effect of height. This hypothesis would need to be tested using direct stomatal conductance measurements (i.e. porometer measurements), but it is indirectly confirmed by the outcomes of the genetic analysis (discussed in detail below). The tendency for shorter rice genotypes to have more open stomata may be linked to the fact that the new, high-yielding semi-dwarf varieties of rice were selected under irrigated conditions (Pingali, 2012; Kumar ) without considering water as a limiting factor. The presence of standing water in the rice paddy system differentiates it from most wheat production environments.

Another aspect that may confound the interpretation of the thermal imaging data is the effect of flowering on the canopy temperature measurements. Our data set allowed the quantification of the contribution of differences in flowering stage to canopy temperature results. The accessions of the panel were sown and transplanted to the field in a staggered way to synchronize phenology with the aim of imposing stress at 50% flowering in all the varieties. Despite the good synchronization of flowering observed among genotypes (Kadam ), a perfect synchronization is difficult to achieve with such a large and diverse panel, grown under varying conditions. Indeed, drought affects flowering time and in many cases accelerates it, a phenomenon referred to as drought escape (Zhang ). The significant correlation between canopy temperature and flowering stage in control plants (Fig. 3A) is in agreement with a 2 °C higher temperature in unstressed wheat canopies with spikes, compared with unstressed canopies without spikes (Hatfield ), and suggests that the quantification of plant canopy temperature can be significantly affected by differences in flowering. However, flowering time differences in this study did not significantly affect the canopy temperature under drought (Fig. 3B). This result may be explained by the fact that under water-limited conditions the rise in leaf temperature may be larger than the increase resulting from the presence of panicles.

Association mapping and QTL identification

In the present study we used the observed phenotypic variation for rice canopy temperature under control and drought conditions in a GWA mapping experiment to identify the genetic factors contributing to this variation. The genotypic differences observed in canopy temperature under the two treatments show that there is substantial genetic variation, which is especially visible under drought. The results indicate a low marker heritability under control conditions both for the single replicates and for their averaged canopy temperature (Supplementary Table S2). This may indicate either that the environmental noise masks the genetic factors under control conditions, or that the absence of stress did not trigger their expression. The absence of significant marker–trait associations for mean canopy temperature values under control conditions (Supplementary Fig. S10) indirectly confirms this.

In accordance with other reports (Jones ), thermal data under drought stress showed good heritability, substantially higher than under control conditions (Supplementary Table S2). This suggests that stress maximizes the genotypic differences in canopy temperature, which therefore can be more effectively detected by thermal imaging. A main marker–trait association was identified on chromosome 3 using the mean temperature values of the single stress replicates and the 45K SNP (GBS) map (Fig. 4). Subsequent mapping of chromosome 3 using the imputed high-density SNP set increased the signal strength of the marker–trait associations, and even more the mapping resolution (Fig. 5), supporting the idea that imputation is a quick and cost-effective tool for adding value to existing genotyped panels (Wang ). The increased resolution provided by the high-density map allowed the identification of two distinct, neighbouring QTLs (Fig. 5B). This distinction could not be resolved using the GBS map because of the low density of markers across the region (Supplementary Fig. S12). Haplotype analysis of recombinant accessions carrying either major or minor SNP alleles across the two QTL regions revealed that qCT1 was primarily responsible for canopy temperature variation (Fig. 6A, B). This suggests that the detection of qCT2 was due to the extensive LD across the region, where a majority of genotypes carried haplotypes with blocks of major or minor alleles across the two QTLs (e.g. Haplotype I and Haplotype III in Fig. 6A). Minor alleles for the SNPs defining qCT1 were associated with higher canopy temperature and occurred almost exclusively in taller, low-yielding genotypes (Fig. 6A; Supplementary Tables S5, S6) whereas major alleles across qCT1 were widely distributed across the panel, found in ~75% of accessions, but notably were almost completely fixed in the shorter (plant height <100 cm), high-yielding genotypes (Supplementary Table S6). This finding may support the hypothesis that rice genetic variation for stomatal conductance (here indirectly determined by canopy temperature measurements) was reduced as a result of selection for short-statured, high-tillering, and productive genotypes in flooded environments. The taller accessions of the population are mostly low-tillering and low-yielding landraces, many of which were selected for drought-prone, rain-fed environments (Kumar ) where alleles responsible for reduced transpiration are preferred, even if this negatively affects grain yield under non-stressed conditions (Passioura, 2012). Fixation of the major alleles at qCT1 in short-statured genotypes suggests that, in rice breeding, selection for high grain yield under flooded conditions reduced the genetic variation available for traits related to more conservative water use.

The marker–trait association detected on chromosome 3 using the mean temperature values of stress replicates was also identified using the temperature values of replicate F03 alone (Fig. 4B), but it fell below significance using F04 (Supplementary Fig. S10, the other drought replicate that shared the field location with F03 (Supplementary Fig. S1). This difference may be explained by the different severity levels of the imposed stress between the two days of imaging, in combination with other environmental variation. F04 was imaged at almost the same time of day, but during a time window characterized by higher mean solar radiation and wind speed (Supplementary Table S1B). Furthermore, F04 was imaged 1 d later than F03, such that the water limitation further increased (Supplementary Fig. S6), resulting in higher canopy temperature values than observed in F03 (Fig. 2). It is likely that, due to the increased severity of the stress in F04 (Fig. 3), the stomata closed in a larger number of genotypes, thus reducing the variability and sensitivity to detect genetic differences in stomatal closure that were still detectable the day before. This confirms the very dynamic response of stomata to changing environmental conditions (Vico ; Drake ) and suggests the need for reducing the time window during which the thermal imaging is performed to increase the number of replicates imaged per day. This target can be achieved by assembling thermal cameras on unmanned aerial vehicles (Shi ).

Candidate genes

A region of ~42 kbp was identified inside the qCT1 LD block by considering the most significant SNPs associated with mean canopy temperature under drought (Fig. 7B). This region contains seven genes (Supplementary Table S7), of which two were targeted as interesting a priori candidate genes, a mitochondrial fumarase (LOC_Os03g21950) and a plasma membrane HAK transporter (LOC_Os03g21890). Plant guard cells accumulate solutes such as K+, Cl–, and malate during stomatal opening and release/metabolize them during stomatal closure. During these processes, solute flux through the plasma membrane of guard cells is highly active, with K+ intake driving stomatal opening with the involvement of different types of transporters, including HAK-type transporters (Jezek and Blatt, 2017). Gago and Santelia and Lawson (2016) recently reviewed the role of guard cell and adjacent mesophyll cell metabolism in stomatal movement, highlighting the importance of malate (negatively charged) intake as a direct counter ion for K+ in the vacuole of guard cells or as a signalling molecule involved in the activation of guard cell vacuolar Cl– channels during stomatal opening. Fumarase is a mitochondrial enzyme involved in the production of malate, through the hydration of fumarate, in a critical step of the tricarboxylic acid (TCA) cycle (Sweetlove ). In transgenic tomato plants, inhibition of fumarase resulted in a reduction in TCA cycle activity. This reduction had little effect on leaf metabolism but markedly reduced plant biomass because of a deficiency in stomatal function that resulted in reduced stomatal conductance (Nunes-Nesi ). The co-location of fumarase and of a HAK transporter, both important for the mechanism of stomatal opening, in the QTL region of highest interest associated with canopy temperature/stomatal conductance variation reinforces our mapping results.

Finally, our analysis of the predicted function of the SNPs within the ~42 kbp region of highest significance within qCT1 did not highlight specific variants responsible for amino acid changes (Supplementary Table S8) but identified many upstream gene variants, suggesting that changes in gene regulation may explain the phenotypic variation associated with this region. This hypothesis should be further investigated by sequencing the candidate genes and their promoter regions in a subset of contrasting lines carrying the major and minor alleles at the significant markers within qCT1, and testing gene expression differences in response to stress. This would help to pinpoint the functional nucleotide polymorphisms and to assess how they impact stomatal conductance and response to water stress.

The significant SNPs across qCT1 may also be of interest to breeding programmes aiming to develop new varieties adapted to drought-prone agroecosystems. These SNPs represent a useful target for either marker-assisted selection or genome editing using CRISPR/Cas9. Additionally, this should help to determine whether introducing minor alleles at qCT1 into semi-dwarf high-yielding varieties equips them with new genetic potential for a more conservative water use strategy under stress. This trait is currently not available in this germplasm. It will be of great interest to see whether improved water use efficiency in these high-yielding varieties can be accomplished without negatively impacting their productivity (in the presence of a non-pleiotropic effect of qCT1 on grain yield) when water is abundant.

Conclusions

Physiological profiling of plant traits combined with genetic analysis has the potential to greatly accelerate crop improvement (Reynolds and Langridge, 2016). The present study shows that changes in stomatal conductance, an important physiological response to water limitation, can be indirectly measured by thermal imaging and that the latter technique can be used to quantitatively screen a large panel of rice accessions. Canopy temperature during stress is a good predictor of grain yield performance and, therefore, thermal imaging represents an effective tool that can be used to accelerate physiological selection in plant breeding. In addition, association mapping of thermal data revealed the presence of genetic variation controlling canopy temperature under stress. The a priori candidate genes that were identified as underlying this genetic variation suggest that differences in the regulation of genes involved in guard cell solute intake affect stomatal behaviour, which we detected as canopy temperature differences. Finally, our analysis shows that the major donors of genetic variation for canopy temperature/stomatal conductance are the tall accessions of the panel. These old varieties and landraces present in crop germplasm collections represent a strategic reserve of genetic variation that can be tapped for developing new understanding of stress response and new varieties that are physiologically adapted to highly variable, water-limited environments.

Supplementary data

Supplementary data are available at JXB online.

Fig. S1. Field trial at the IRRI during the dry season 2013–2014.

Fig. S2. Density of distribution of the standardized residuals for all the imaged field replicates.

Fig. S3. Principal component analysis plots for the 271 indica rice accessions.

Fig. S4. Mean temperature of thermal pictures relative to the E03 drought field replicate.

Fig. S5. Correlation matrices of non-normalized and normalized field replicates.

Fig. S6. Soil water potential of the drought field.

Fig. S7. Boxplots representing the range of variation for the recorded agronomic traits.

Fig. S8. Scatterplots between drought stress mean temperature values and plant height, grain yield, and harvest index scored during the 2013 and 2014 dry seasons.

Fig. S9. Quantile–Quantile plots of expected versus observed P-values for the GWA mapping results.

Fig. S10. Manhattan plots of the GWA mapping results of control replicate B mean, B03, and B04 and of the separate drought replicates E02, E03, and F04.

Fig. S11. Comparison between Manhattan plots of the GWA mapping results using the 45K SNP map without the marker-based PC or including PC1 as a covariate in the linear-mixed model for ‘EF mean’ and F03.

Fig. S12. Manhattan plots displaying the level of significance over genomic positions of the chromosome 3 region of the two markers’ loci significantly associated with canopy temperature of ‘EF mean’.

Fig. S13. Manhattan plots of the GWA mapping results for the QTL region of chromosome 3 using the imputed map for ‘EF mean’ canopy temperature.

Fig. S14. Localized linkage disequilibrium analysis of qCT2.

Table S1A. Detailed information on imaged field replicates.

Table S1B. Weather station data of the days of imaging.

Table S2. Heritability estimated by markers (h2) of temperature results for all the field replicates.

Table S3. Significant SNPs identified by GWA mapping of canopy temperature under drought stress using the GBS SNP map.

Table S4. Allelic effect on canopy temperature for markers of the most significant SNPs, located on chromosome 3, in the different drought field replicates and for their mean values.

Table S5. Significant SNPs identified by GWA mapping of ‘EF mean’ canopy temperature using the imputed map.

Table S6. List of the 246 accessions carrying one of the three haplotype groups, identified considering the most significant SNPs of qCT1 and qCT2, and their phenotypic performance.

Table S7. Genes included in the localized region delimited by the most significantly associated SNPs with canopy temperature of ‘EF mean’ and located inside qCT1.

Table S8. Effect of the 20 SNPs present in the qCT1 region and significantly associated with ‘EF mean’ canopy temperature.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.