Identification of quantitative trait loci associated with canopy temperature in soybean

A consistent risk for soybean (Glycine max L.) production is the impact of drought on growth and yield. Canopy temperature (CT) is an indirect measure of transpiration rate and stomatal conductance and may be valuable in distinguishing differences among genotypes in response to drought. The objective of this study was to map quantitative trait loci (QTLs) associated with CT using thermal infrared imaging in a population of recombinant inbred lines developed from a cross between KS4895 and Jackson. Heritability of CT was 35% when estimated across environments. QTL analysis identified 11 loci for CT distributed on eight chromosomes that individually explained between 4.6 and 12.3% of the phenotypic variation. The locus on Gm11 was identified in two individual environments and across environments and explained the highest proportion of phenotypic variation (9.3% to 11.5%) in CT. Several of these CT loci coincided with the genomic regions from previous studies associated with canopy wilting, canopy temperature, water use efficiency, and other morpho-physiological traits related with drought tolerance. Candidate genes with biological function related to transpiration, root development, and signal transduction underlie these putative CT loci. These genomic regions may be important resources in soybean breeding programs to improve tolerance to drought.

. Planting date and weather data including maximum temperature (MaxT), minimum temperature (MinT), No. of days without rainfall, estimated vapor pressure deficit (VPD), and estimated soil moisture deficit at the time of canopy temperature measurements for the field experiments conducted at Pine Tree (PT) and Rohwer (RH) in 2017, 2018, and 2019. a AE denote averaged across environments. b h 2 narrow sense heritability. www.nature.com/scientificreports/ Phenotypic evaluation. CT measurements were made multiple times (two to five times) in all environments between 1200 and 1400 h when the sun was unobstructed by clouds. All measurement dates occurred once the canopy was completely closed to avoid the CT measurements being overwhelmed by the heat signature from the soil that occurred prior to canopy closure. Dates for CT measurement were targeted when the estimated soil moisture deficit was > 45 mm 67 and soybean development was between full bloom and beginning seedfill 9,14 . Data ultimately used for QTL analysis for each environment were chosen for those dates in which the estimated soil moisture deficit was the greatest as previous research had determined that genotypic differences in CT were not apparent under water replete conditions 9 . Canopy temperature was determined by using aerial infrared image analysis. A drone (DJI Phantom 4 Pro, www.dji.com/phant om-4-pro) equipped with an infrared (IR) camera (Model: FLIR Tau 2, www.flir.com) with a sensor size of 640 × 512 pixels, 25 mm focal length, sensitivity of 0.05 °C, and an accuracy of ± 5%, was flown at a height of 120 m above the ground at full canopy coverage. A digital video recorder recorded the video stream from the camera. The IR camera is not calibrated. That is, values of CT from the IR range from 0 (cool) to 255 (hot) and cover a range of 12.5 °C, but the specific temperature of the canopy is not determined. Herein, we report values directly from the IR camera as a measure of CT.
The IR images were extracted from the video file using VLC video player (www.video lan.org), and 6 to 17 images representing each plot multiple times were selected manually. Selected IR images were processed using FieldAnalyzer software (www.turfa nalyz er.com/field -analy zer) to extract CT values from the average values of 400 to 2000 pixels from the center portion of the IR image of each plot and was used as a measure of CT. There were multiple CT values of each plot extracted from multiple selected images. The final CT values used for analysis was the average CT values determined from analyzing multiple images and after removing values that were more than ± 2 standard errors from the mean. Statistical analysis. The phenotypic analysis of CT was performed using SAS v.9.2 software (SAS Institute, 2013). The normality of CT in an individual environment was checked using a Q-Q plot of residuals and the Shapiro-Wilk test 70 . The presence of statistical differences between parents for CT was estimated using a t test. Pearson's correlation analysis between environments and analysis of variance (ANOVA) averaged over environments were performed using PROC CORR (α = 0.05) and PROC MIXED procedures, respectively. Heritability was estimated from the variance components calculated with restricted maximum likelihood (REML) method of VARCOMP procedure. Narrow sense heritability (h 2 ) was calculated using the following formula 71 : where σ 2 G , σ 2 GE , σ 2 e are genotypic variance, genotype-by-environment variance, and error variance, respectively, and E and R are the number of environments and replications, respectively. Because F 5 -derived RILs were used in this research, σ 2 G was composed entirely of additive variance and additive × additive epistasis variance, with negligible variance associated with other components of dominance variance. As the result, this heritability should be considered as narrow sense heritability. BLUPs (best linear unbiased predictions) were calculated using PROC MIXED procedure for individual environments and averaged across environments, considering all factors in the model as random. Environment was considered a fixed effect in the combined data analysis. QTL analysis of CT was performed using BLUP values to reduce the environmental variations. QTL analysis. The genetic map for the KS4895 × Jackson mapping population was previously described by Hwang et al. 54 and used for QTL analysis in the present study. Briefly, the linkage map consists of 38 simplesequence repeat (SSR) and 510 single-nucleotide polymorphism (SNP) markers with a map size of 2089 cM. The list of markers and their physical position on specific chromosomes is provided in Supplementary file 1. The BLUP estimates calculated for individual environments and averaged across environments were used for QTL analysis. The QTL analysis was performed with WinQTL Cartographer v.2.5 72 using composite interval mapping (CIM) and multiple interval mapping (MIM) methods. CIM was performed using Model 6 (standard model) of the Zmapqtl program module 73 . Markers as co-factors to control background noise were selected with forward and backward stepwise regression methods with a walk speed and window size of 1 cM. A significant LOD (log of odds) threshold score was determined by a permutation test with 1000 runs and with a genome wide type I error of 5% 74 . The most likely position of QTLs and an estimate of the magnitude of their additive effects were determined using the CIM method 75 . The confidence interval for putative QTL positions was determined by one-LOD drop on either side of the LOD peak.
Multiple interval mapping (MIM) is a stepwise model procedure 76 in which presence of significant QTLs and QTL × QTL interactions are detected using QTLs identified in the CIM method as an initial MIM selection model. This pre-selected MIM model was optimized for identified QTLs, search for new QTLs, and QTL × QTL interactions by using the 'optimize QTLs positions' , 'search for new QTLs' , and 'QTL interactions' options, respectively. The MIM model was determined with the minimum Bayesian information criterion (BIC), c(n) = In(n), and with genome walk speed and window size of 1 cM. The criterion used to declare coincident QTLs between Across environments: www.nature.com/scientificreports/ environments was based on at least a 10 cM overlap in QTL intervals on the linkage map. In this study, a QTL that explained more than 10% of phenotypic variation was considered a major QTL.

Candidate genes identification.
Candidate genes for genomic regions underlying putative QTLs for CT identified in each environment and across environments were searched using the genome browser of Soybase (www.soyba se.com). The candidate genes falling within ± 1 MB from the nearest marker of putative QTLs were selected according to the Glyma2.0 assembly in Soybase (www.soyba se.org) with consideration for those genes having biological function related to transpiration, canopy temperature, rooting, and plant water relations.

Results
Phenotypic evaluation. On the dates for CT measurement, the maximum temperature ranged from 29 (PT17) to 37 °C (PT19) and minimum temperature ranged from 15 (PT17) to 25 °C (RH17) ( Table 1). The estimated VPD was ≥ 2.3 kPa on day of CT measurements in all environments except for RH18 in which VPD was 1.6 kPa. There had been no rainfall from 4 to 13 days prior to CT measurements. The estimated soil moisture deficit ranged from 49 mm to more than 75 mm on the day of CT measurements. Irrigation is recommended if estimated soil moisture deficit exceeds 37 mm for silt loam soils present at Pine Tree and 50 mm for silty clay soils present at Rohwer 67 , indicating drought-stress conditions on the day of CT measurements.
The CT values had a wide range in all environments with RIL means that ranged from 50 to 64 ( Fig. 1, Table 2). Jackson had lower CT values than KS4895 in all environments except RH17 (Table 2), and parents were significantly (P < 0.05) different in PT17, RH19, and across environments (AE) ( Table 2). Averaged over environments, IR values were 12 units less for Jackson than KS4895, which indicate a CT that was approximately 0.59 °C cooler for Jackson. The distributions of CT values were approximately normal (P < 0.001) except for PT18 and RH18, which were slightly skewed left as indicated by the Shapiro and Wilk test (data not shown, Fig. 1). There were weakly significant correlations for CT between environments for PT17 and PT18 (r = 0.18*), PT17 and RH19 (r = 0.26**), RH17 and RH19 (r = 0.15*), PT19 and RH18 (r = − 0.16*), and PT19 and RH19 (r = 0.49***). Across environments, ANOVA indicated that there was a significant effect (P < 0.05) of RILs, environment, and

QTL analysis.
Analysis of QTLs associated with CT in individual environments identified seven QTLs present on five chromosomes using the CIM method. Of these seven QTLs, one QTL was identified in RH17, PT18, and RH18 and two QTLs were identified in PT19 and RH19. No QTLs were detected in PT17. These QTLs had additive effects that ranged from − 0.08 to − 1.17 and explained 5.7% to 12.3% of the phenotypic variation ( Table 3). The QTLs identified in PT18 (1), PT19 (2), RH18 (1), and RH19 (1) were also identified by the MIM method. Two additional QTLs present on Gm13 (at 31355907 bp) in RH18 and on Gm16 (at 26745906 bp) in PT19 were identified by the MIM method but were not identified by the CIM method (Table 3). Across environment (AE) QTL analysis of CT identified five QTLs present on Gm02 (1), Gm11 (1), and Gm18 (3) with additive effects ranging from − 0.39 to − 0.51 and explaining phenotypic variation from 5.3 to 9.3%. Three out of five AE QTLs were common between CIM and MIM methods, and one new QTL on Gm15 (at 51379618 bp) was identified only by the MIM method ( Table 3). The QTL present on Gm11 (at 10350509 bp) was common among PT19, RH19, and AE. The QTL on Gm18 (at 50727159 bp) was identified in RH17 and AE, and the QTL on Gm18 (at 56161047 bp) was common for RH19 and AE. All other QTLs were environmentally specific ( Table 3). The favorable alleles for deceasing CT for all the QTLs were from KS4895 except for the QTLs on Gm10 and Gm16. QTL × QTL interactions were identified between the QTLs present on Gm11 and Gm16 Table 3. The QTLs associated with canopy temperature identified by composite interval mapping (CIM) and multiple interval mapping (MIM) in a RIL population of KS4895 and Jackson which were evaluated at Pine Tree and Rohwer in 2017, 2018, and 2019. a Glycine max chromosome on which QTL was present. b Environment in which specific QTL was identified. Prefixes PT and RH denotes Pine Tree and Rohwer, respectively followed by 17 (2017), 18 (2018), and 19 (2019) for years. AE denotes averaged across environments. c QTL position in base pairs on respective chromosomes according to Glycine max genome assembly 2.0 (Glyma.Wm82.a2; www.soyba se.com). d Nearest marker to the QTL. e Additive effect of the QTL. f Proportion (%) of phenotypic variation explained by specific QTL. g Allele that decreases CT values was considered as the favorable allele; Positive sign indicates that favorable alleles (decreasing CT) were from Jackson and negative sign indicates the KS4895 allele. h Flanking markers (FM) present near or at 95% confidence interval (CI) of the maximum likely QTL positions. The LOD values with ± 1 declination was used to estimate the 95% confidence interval. Candidate gene identification. There were more than 1092 candidate genes present within ± 1 MB of the nearest markers for putative QTLs with a range from 53 to 206 genes for individual QTLs. A list of these genes (Glyma2.0 ID) with their biological annotation is provided in Supplementary file 2. Out of these 1092 genes, those having biological function related to stomatal complex morphogenesis, regulation of stomatal movement, response to water deprivation, response to abscisic acid (ABA), ABA mediated signaling pathway, ABA transport, root hair elongation, root hair cell differentiation, primary and adventitious root development, root morphogenesis, water transport, response to osmotic and oxidative stress, signal transduction, and response to different hormones stimulus were considered to play a potential role in controlling CT in response to different soil moisture conditions (bold text in Supplementary file 2).

Discussion
Remote sensing with thermal infrared imagining provides a powerful approach to measure and compare temperature differences of plant canopies of a large number of genotypes for field scale experiments and has been widely used in various crops 3,[13][14][15]43,77 . The present study investigated the genetic control of CT in a population derived from a cross between KS4895 and Jackson, which was evaluated across six environments. Under replete soil moisture and aerial environmental conditions, plants continue to transpire through open stomata. In contrast, as soil moisture becomes limiting, plants close stomata as a preventive mechanism, transpiration decreases, and there is an increase in CT 31 . Those genotypes that have access to soil moisture continue transpiration during drought stress, resulting in a cooler canopy. High soil moisture deficit and VPD at the time of CT measurements resulted in drought stress in all environments ( Table 1). The PT19 and RH19 environments had greater soil moisture deficit and VPD, resulting in greater differences among RILs and higher heritability of CT as compared to other environments. For all environments there was transgressive segregation among RILs with lower and greater CT than the parents, indicating the distribution of favorable alleles for cooler canopy temperature in both parents. CT is highly influenced by environmental conditions (soil moisture availability, vapor pressure deficit, air temperature) and plant morphology (canopy and root architecture conditions) 34,78 , resulting in significant genotype × environment interactions. While the h 2 of CT was 35% when estimated across six environments, and the range of h 2 among environments was from 7 to 62% ( Table 1). The two environments in which h 2 was > 50% were noted for a severe soil moisture deficit (> 75 mm), maximum temperatures ≥ 36 °C, and VPD ≥ 3.3 (Table 1). From this it appears that high h 2 did not appear to be solely due to soil moisture deficit as PT18 and RH18 had severe soil moisture deficit but had h 2 ≤ 11%. Environments with high h 2 (PT19 and RH19) both had severe soil moisture deficit (≥ 75 mm) combined with high temperature (≥ 36 °C) and high VPD (≥ 3.3 kPa). Along with making measurements when the canopy is completely closed and the sun is completely unobstructed, other conditions for increasing h 2 for CT include making measurements when there is a high soil moisture deficit, VPD exceeds 3.3 kPa, and when maximum temperatures are high. Additional research is needed in order to define better the conditions that would optimize CT measurements and increase h 2 . Kaler et al. 14 reported a broad sense heritability (H) of 19% for CT in a diverse panel of soybean accessions evaluated in different environments using GWAS. Likewise, low to moderate broad sense heritability of CT/CTD was reported in wheat (H = 16% to 38%) [79][80][81] , in sugarcane (H = 9% to 44%) 82 , and in rice (H = 21% to 30%) 83 . The low heritability of CT/CTD in soybean and other crops indicate that this trait is highly influenced by environmental conditions 84 .
The complexity of CT in soybean is further indicated by detection of multiple QTLs by CIM and MIM methods. The QTL analysis for CT using CIM and MIM methods identified a total of nine QTLs in individual environments and six QTLs when data were averaged across environments ( Table 3). Most of these QTLs were environmentally specific except the QTLs present on Gm11 (at 10350509 bp) and Gm18 (at 50727159 bp and 56161047 bp), which were common in at least one individual environment and across environments ( Table 3). The QTLs present at Locus 5 on Gm11 explained phenotypic variation that ranged from 9.3% to 11.5% and is considered a major QTL ( Table 3). The markers linked with these QTLs have potential utility using marker assisted selection or genomic selection in a breeding program. The inconsistency of QTLs among environments might be due to different environmental factors such as field moisture status, soil temperature and depth, solar radiation, and VPD. Although QTLs were generally environmentally specific, most of these QTLs were detected by both CIM and MIM methods, which increases the confidence of these results. Of particular interest are environmental conditions that could improve consistency and increase heritability of CT. More research is needed to increase the utility of markers linked with QTLs identified in specific environments in selecting genotypes with a cooler canopy and to determine environmental conditions that optimize heritability.
Considering the overlapping confidence interval of QTLs identified in individual environments and across environments, we found 11 loci on eight chromosomes (Table 3, Fig. 2). The flanking markers position of each loci is provided in the Supplementary file 3. Even though Jackson tended to have lower CT among environments (Table 2), there were nine favorable loci from KS4895 and two favorable loci from Jackson. Hwang et al. 54 found that KS4895 contributed many of the favorable alleles for slow canopy wilting that were identified in multiple biparental populations. Bai and Purcell 9 found a positive correlation (0.37 < r < 0.62) between cooler CT and slow canopy wilting. There is a possibility that KS4895 has both slow canopy wilting and cooler canopy temperature alleles in the KS4895 × Jackson population.
For comparative analysis of CT loci with QTLs associated with plant water relations and drought tolerance related traits, the putative CT loci were aligned on the soybean linkage map in Soybase (www.soyba se.com). The QTLs previously mapped within the 95% confidence interval of putative QTLs in this study were considered to be Scientific Reports | (2020) 10:17604 | https://doi.org/10.1038/s41598-020-74614-8 www.nature.com/scientificreports/ present in the same genomic regions. The CT Loci 1 and 2 (on Gm02) coincide with a canopy wilting QTL found in a population derived from KS4895 and Jackson (09705KJ population) that is different from the population used in the present research and in a population derived from Benning and PI 416937 54 . The favorable alleles for slow canopy wilting at this locus are from KS4895 and PI 416937, consistent with the present research that the KS4895 allele is associated with cool CT. These loci localized in a genomic cluster found in meta-analysis of canopy wilting QTLs using multiple biparental populations 55 . The loci were also mapped in an association panel in the genomic regions associated with canopy wilting 56 , canopy coverage 61 , and δ 13 C 60 . The co-localization of these drought-related QTLs with CT indicates the strong relationship among transpiration, WUE, canopy wilting, and CT. Locus 3 (on Gm03) and Locus 4 (on Gm10) coincide with QTLs for CT 14 and canopy wilting 56 , respectively, that were identified in GWAS. Locus 5 (on Gm11) maps to the same genomic region associated with CT 14 , canopy wilting 56 , and canopy coverage 61 identified in GWAS analysis conducted using a diverse panel of soybean accessions. Locus 5 also coincides with a QTL for canopy wilting 54 and for yield (unpublished results, Bazzer and Purcell) identified in population with KS4895 as a parent; the favorable alleles for slow canopy wilting in both these populations were from KS4895 as they were for CT in the present study.
Locus 6 (on Gm13) overlaps with a QTL for δ 13 C found by GWAS 60 . Nguyen et al. 85 mapped a QTL for root area and root length at the same genomic location in a population derived from Tachinagaha × Iyodaizu. Rooting depth and area are drought avoidance mechanisms that may increase water availability 86 . The coincidence of CT and root morphology QTLs may point to a root system that extracts water from deeper soil horizons and results in cooler canopy during drought. In wheat, CT QTLs have been linked with rooting traits that allow extraction of more water from soil under drought 87 .
Locus 7 (on Gm15) and Locus 8 (on Gm16) coincide with canopy wilting 56 and canopy coverage 61 that were identified by GWAS. In addition, Locus 8 also overlaps with δ 13 C 60 . Earlier canopy coverage helps to decrease the water loss by soil evaporation relative to transpiration and improve WUE 88 . Locus 9 (on Gm18 at 17323638 bp) coincides with a WUE QTL mapped in a Young × PI416937 population with the 'Young' allele increasing WUE 89 . Locus 9 and Locus 10 (on Gm18 at 50727159 bp) also fall in the genomic region harboring QTLs for CT 14 , canopy wilting 56 , and δ 13 C 59,60 identified by GWAS. In addition, Locus 10 overlaps with the oxygen isotope ratio 60 , which is a proxy for measurement of transpiration and is associated with stomatal conductance 90 . www.nature.com/scientificreports/ Locus 11 (on Gm18 at 56161047 bp) coincides with δ 13 C identified in the same population as the present research 58 . The favorable allele for δ 13 C at this locus was from Jackson, while KS4895 provided the favorable allele for CT at this locus. The coincidence of δ 13 C and CT QTLs illustrates a shared genetic relationship between these two physiological traits. The co-localization of putative CT loci with QTLs associated with other morphophysiological traits such as WUE, canopy wilting, canopy coverage, and transpiration may be a pleiotropic effect of the same genes controlling these traits or that the genes are spaced closely together on specific chromosomes.
The candidate gene search underlying putative CT loci identified genes with functions related to plant water relations, root morphology, and transpiration. These include genes involved in stomatal complex morphogenesis, regulation of stomatal movement, response to water deprivation, response to abscisic acid (ABA), ABA mediated signaling pathway, ABA transport, root hair elongation, root hair cell differentiation, response to oxidative stress, and signal transduction.
The upregulation of root morphology (lateral root formation, root hair elongation, root development response to ABA) related genes during drought may result in extracting residual soil moisture that maintains primary growth and developmental processes. In wheat, deeper root development in response to drought stress resulted in a cooler canopy and an increase in yield 25,91 . Aquaporin-related genes were also found underlying CT loci and these membrane proteins allow movement of water throughout plant in response to stress 92 . The co-localization of CT loci with QTLs associated with drought tolerance related traits and with underlying candidate genes with biological function related to transpiration, stomatal conductance, and plant water relations increases the probability that putative CT loci are associated with variation in CT in the present research. Additional research is needed to confirm the canopy temperature QTLs in different populations and in different environments to increase the efficiency of these genomic regions in selecting genotypes with a cooler canopy and with drought tolerance.

Conclusions
In the present study, we identified genomic regions associated with CT in a recombinant inbred population derived from KS4895 and Jackson that was phenotyped in six different environments. These results represent the first QTLs for CT identified in soybean using a biparental population. The population segregated for CT in all environments and was used for QTL analysis. The heritability of CT was relatively low as compared to other morpho-physiological traits due to greater influence of different environmental factors. Eleven genomic loci present on eight chromosomes for CT were identified across several environments. The CT locus present on Gm11 explained phenotypic variation more than 10% and was considered as a major QTL. The favorable allele for cooler canopy for most of the loci were from KS4895, which were also coincident with canopy wilting QTLs identified in multiple biparental populations by Hwang et al. 54 . The identified CT QTLs coincided with QTLs associated with drought tolerance related traits mapped in previous studies, and genomic regions underlying these putative CT have the genes with biological function related to transpiration, stomatal conductance, and plant water relations. More research is needed to confirm these QTLs in different genetic backgrounds and in multiple/different environments to evaluate the efficiency of these QTLs to use in soybean breeding for selecting genotypes with a cooler canopy and drought tolerance.