Monitoring of drought stress and transpiration rate using proximal thermal and hyperspectral imaging in an indoor automated plant phenotyping platform

Abstract

Background

Thermography is a popular tool to assess plant water-use behavior, as plant temperature is influenced by transpiration rate, and is commonly used in field experiments to detect plant water deficit. Its application in indoor automated phenotyping platforms is still limited and mainly focuses on differences in plant temperature between genotypes or treatments, instead of estimating stomatal conductance or transpiration rate. In this study, the transferability of commonly used thermography analysis protocols from the field to greenhouse phenotyping platforms was evaluated. In addition, the added value of combining thermal infrared (TIR) with hyperspectral imaging to monitor drought effects on plant transpiration rate (E) was evaluated.

Results

The sensitivity of commonly used TIR indices to detect drought-induced and genotypic differences in water status was investigated in eight maize inbred lines in the automated phenotyping platform PHENOVISION. Indices that normalized plant temperature for vapor pressure deficit and/or air temperature at the time of imaging were most sensitive to drought and could detect genotypic differences in the plants’ water-use behavior. However, these indices were not strongly correlated to stomatal conductance and E. The canopy temperature depression index, the crop water stress index and the simplified stomatal conductance index were more suitable to monitor these traits, and were consequently used to develop empirical E prediction models by combining them with hyperspectral indices and/or environmental variables. Different modeling strategies were evaluated, including single index-based, machine learning and mechanistic models. Model comparison showed that combining multiple TIR indices in a random forest model can improve E prediction accuracy, and that the contribution of the hyperspectral data is limited when multiple indices are used. However, the empirical models trained on one genotype were not transferable to all eight inbred lines.

Conclusion

Overall, this study demonstrates that existing TIR indices can be used to monitor drought stress and develop E prediction models in an indoor setup, as long as the indices normalize plant temperature for ambient air temperature or relative humidity.

Introduction

Effects of water deficit (WD) on plant physiology have been widely studied, as they can have a major impact on plant performance and yield [1]. Thermography is a popular tool to monitor changes in plant water-use behavior, especially those in response to WD. The amount of emitted thermal infrared (TIR) radiation of a plant relates to its temperature, which in itself depends on incoming radiation and transpiration rate (E) because of the latent heat of vaporization at the leaf surface [2]. Transpiration is important in regulating plant temperature for enzymatic processes and in providing nutrients from the soil via the xylem. It is driven by the evaporative demand in the atmosphere and controlled by stomatal conductance (g_s). Under water-limiting conditions, plants will reduce their E and water loss by closing their stomata, resulting in an increase in plant temperature (T_p), which can be detected within leaves or at the whole leaf, plant and canopy level. As a consequence, the difference between plant and air temperature (T_a) will be less negative in plants exposed to WD as compared to well-watered (WW) plants, allowing the detection of drought stress.

Thermal infrared sensors have been used in a wide range of studies on abiotic and biotic stresses, such as plant WD, salinity, and pathogen infection [3,4,5,6,7,8,9]. The performance in detecting stress strongly depends on stress-induced changes in g_s and E. In the case of plant WD, thermography has been implemented for automating irrigation [10, 11]. Thermography has also become a promising tool for plant phenotyping, because T_p correlates with both water-use behavior and physiological traits, such as E, g_s, relative water content, water potential (ψ) and non-photochemical quenching [12,13,14,15,16]. Thermal infrared sensors have been used to screen for stress-tolerant genotypes [5, 17] and to monitor the dynamic stomatal responses to abiotic stress [7, 18].

Thermography can be applied on different spatial scales from organ to field level [4, 5, 18, 19]. Most phenotyping studies have focused on field conditions, while the usability of thermography in indoor phenotyping platforms is much less investigated. The latter allows for the comparison of treatments or genotypes under controlled environmental conditions, which can improve the accuracy of thermal imaging methods, as these are strongly affected by light intensity, temperature, vapor pressure deficit (VPD) and wind speed. The setup requires nevertheless temperature normalization and the use of reference surfaces, as indoor TIR detection is still influenced by all infrared radiation emitting surfaces, including greenhouse or growth chamber infrastructure. Most studies that have implemented TIR imaging in indoor phenotyping platforms have done direct comparisons of T_p or have normalized plant temperature by calculating TIR indices [7, 19,20,21,22].

The simplest T_p normalization is the Canopy Temperature Depression (CTD) index, which subtracts T_a from T_p to correct for ambient temperature conditions. This index has been applied to plant disease, heat tolerance and drought stress studies of crops or transgenic plants [8, 23, 24], as well as for automated irrigation scheduling of wheat in the form of Stress Degree Days (SDD, [25]). An alternative approach to standardize plant temperature is by normalizing it against the temperature of a non-transpiring and fully transpiring leaf or plant. This approach was implemented in the popular Crop Water Stress Index (CWSI), in which CTD was normalized against its maximum (non-transpiring) and minimum (fully transpiring) values. Three methods have been described to calculate the CWSI: the analytical (or energy balance), empirical (or baseline) and direct approach [26]. The latter two are more popular, as they require less meteorological information. The CTDs of non- and fully transpiring leaves are estimated with an empirical baseline, a temperature histogram, or with wet and dry reference leaf temperatures (T_wet and T_dry, respectively, [27]). The CWSI has been correlated with leaf water potential, g_s and evapotranspiration [9, 12, 13, 19, 23, 28, 29], but is not linearly related to g_s in the direct approach [26]. Consequently, the CWSI has mainly been applied to detect stress and not to estimate g_s or E. The latter can be done using the energy balance model that describes the energy exchange between a leaf and its environment [2]. This mechanistic approach simplifies and improves the biological interpretation of the relationship between T_p and g_s or E, but is based on many meteorological and physical variables that can be difficult to accurately estimate at the appropriate spatial scale and increase the complexity of data processing [18, 30]. Simplifications of the energy balance model have been developed to reduce the amount of required data by incorporating T_wet and/or T_dry [30]. The reference materials used to estimate T_wet and T_dry need to have the same optical and thermal properties as the leaves and have to be exposed to the same environmental conditions [31]. Incorporating this in an imaging setup can be technically challenging, consequently the energy balance approach has not been widely applied in phenotyping studies.

In phenomics, TIR imaging is often combined with fluorescence and hyperspectral imaging, as these systems can capture complementary information on photosynthetic efficiency, leaf anatomy, pigment and water content [32,33,34]. The main reason for combining imaging systems is to simultaneously monitor multiple physiological traits, even though the integration of different data types may produce more robust trait prediction models. The combination of multi/hyperspectral, TIR and/or fluorescence data has been investigated in disease, drought and yield detection in the field [15, 20, 35,36,37,38,39,40,41]. An increased prediction or classification accuracy was obtained in the combination of TIR and hyperspectral indices in stepwise multiple linear regression, Partial Least Squares Regression (PLSR) or Support Vector Machine Classification models [15, 37, 38].

In this study, commonly used TIR drought detection and E prediction methods developed under field conditions were investigated in a maize drought stress experiment using PHENOVISION, a high-throughput phenotyping platform in semi-controlled greenhouse conditions [42, 43] with the following aims:

• The drought sensitivity of existing and new TIR indices was compared using eight maize inbred lines.

• The advantages of combining TIR, environmental and hyperspectral data in E prediction models was evaluated, as the performance of this approach has only been investigated in field applications [15, 37, 38]. This question was tackled by developing and comparing empirical models created using different data combinations and modeling approaches, such as linear, Least Absolute Shrinkage and Selection Operator (LASSO), PLSR and Random Forest (RF) regressions.

• Empirical models are often setup specifically, limiting the transferability to other genotypes, species and imaging setups/platforms. To evaluate setup specificity, the empirical models trained on one genotype were applied on seven additional maize inbred lines.

• The best performing empirical model was also compared to a simplified version of the energy balance approach, which is a setup independent mechanistic model, to determine the most effective approach to estimate E.

Methods

Thermal imaging setup and environmental monitoring in the PHENOVISION plant phenotyping system

PHENOVISION is a plant phenotyping platform that automatically irrigates and images 392 plants (built by SMO, Eeklo, Belgium,[44]). Three imaging systems are available: a multi-view red–green–blue (RGB) imaging system, a top view TIR camera and a hyperspectral imaging system [42, 43]. The TIR camera is positioned in an enclosed imaging cabin to eliminate outside radiation. The thermal images are captured using a FLIR SC645 with a 24.6-mm lens (25° × 18.8° FOV) (FLIR Systems Inc., Belgium), which has an accuracy of 2% for objects with a temperature range of -20 to 150°C. The wavelength band ranges from 7.5 to 13 µm. The camera has a fixed top-view (close to nadir) position. Plant height is measured by means of a light curtain at the entrance of the cabin to adjust the position of plants (lift height) to the focus distance of the camera lens. At a distance of 3.5 m, the pixels in the 640 X 480 array have a spatial resolution of 5.76 mm². Inside the imaging cabin, light racks containing light emitting diode (LED) lamps are attached to the ceiling and to one wall to provide consistent illumination. The LED lamps provide a photosynthetic active radiation (PAR) light intensity of 74 µmol photons m⁻² s⁻¹ and emitt no radiant heat. A Lambertian surface made of aluminum foil and a black metal reference plate are positioned within the field of view of the camera. The aluminum foil is used to measure the reflected temperature, which is the energy emitted by the surroundings of the plant and reflected by the plant. It is required to correct the measured thermal energy of the plant and consequently the T_p. The black reference plate is connected to a thermocouple to monitor deviations in the camera accuracy over time (‘drift’), which is further explained in the image processing section. Besides this, the camera is allowed to perform automatic non-uniformity correction (NUC), which corrects for detector drift due to camera frame heating during operation. In this study, the imaged temperature of the black plate was also used as a surrogate for the temperature of a non-transpiring plant (T_dry, Additional file 1, [4, 27]). Thermal data were combined with environmental measurements collected at three different locations in the greenhouse, namely in the growth zone (gz), and inside (in) and outside (out) the cabin (Fig. 1).

Fig. 1

Location of environment monitoring positions in PHENOVISION. A indicates the growth zone. B shows the monitoring position outside the imaging cabin. This is the location where the plants are waiting before imaging. C indicates the position inside the TIR imaging cabin

Experimental setup

Plants were grown in a semi-controlled environment with adjustable T_a and relative humidity (RH). During seedling establishment, T_a was set at 23 °C during daytime and 22 °C at night and vapor pressure deficit was maintained at 1.5 kPa by means of RH adjustments. A diurnal gradient in environmental conditions was established from the V4 (four fully developed leaves) stage onwards by gradually increasing T_a from 22 °C at night to 28 °C in the afternoon, which resulted in VPD values ranging from 1.5 kPa during the night to 3 kPa in the afternoon. A 16/8-h day/night cycle was maintained in the greenhouse using high-pressure sodium vapor lamps. The light intensity (photosynthetically active radiation, PAR) was 280 µmol m⁻² s⁻¹ on average. Additionally, RH, T_a, VPD, PAR, and black sphere temperature (T_BS) were measured in the growth zone by four weather stations containing a SKH 2053 humidity and temperature sensor, a PAR SKL 2625 sensor (Skye Instruments, UK) and a metal black sphere (Testo BE, Belgium) with a MWTC/MWTC-D thermocouple inside (OMEGA Engineering Ltd., UK). The black sphere temperature represents radiant heat, which may be higher than T_a under influence of solar radiation or heat produced by lamps. Slight variability in PAR and T_a was observed within the growth zone which had a negligible effect on the phenotype (Additional file 2: Table S1). Consequently, the gz measurements were averaged. Outside and inside the imaging cabin, RH and T_a were monitored using El-USB-2 loggers (Lascar electronics, UK). Maize plants were sown in 7-l plastic pots filled with peat-based soil containing Osmocote® fertilizer (N. V. Van Israel, Belgium). The plants were fertilized weekly with 40 ml of a 200-ppm N solution of Peters® Excel CalMag Grower (ICL Specialty Fertilizers, The Netherlands).

Two drought stress experiments were performed to determine the effectiveness of TIR imaging for drought stress detection and E predictions. The first small drought experiment (DR) investigated the effect of WD on the transpiration rate and stomatal conductance of maize genotype B104, while the second experiment tested the transferability (TF) of thermal image-based phenotyping to other maize genotypes.

During the DR experiment, B104 maize plants were positioned randomly on the platform. WD treatments were applied based on the gravimetric measurement of soil water content and automated weighing and watering in the phenotyping platform. The DR experiment included four WD treatments with a soil water content of 0.8 (soil water potential < − 1500 kPa), 1.0 (< − 1500 kPa), 1.4 (± − 100 kPa) and 1.6 g H₂O g⁻¹ dry soil (± -40 kPa), and one WW treatment with a soil water content of 2.4 g H₂O g⁻¹ dry soil (± -10 kPa). Twenty plants were randomly assigned to each treatment and were surrounded by border plants, which received a mild WD treatment (1.8 g H₂O g⁻¹ dry soil, ± − 25 kPa). At sowing, a fixed amount of water was added to each pot to ensure germination. After this, plants were weighed on a daily basis and water was added according to the target soil water content until the plants reached V10-11 (10 to 11 fully developed leaves, the end of the experiment).

The TF experiment contained a total of 306 plants belonging to eight inbred lines, namely B104, H99, MS71, NC358, OH43, TX303, Tzi8 and W153R, and 80 border plants (B104). The inbred lines contained stiff stalk (B104), non-stiff stalk (H99, MS71, OH43 and W153R), tropical-subtropical (NC358 and Tzi8) and tropical-subtropical mixed (TX303) genotypes [45] that varied in drought sensitivity, morphology and developmental timing. For each genotype, 19 plants were assigned to a WD treatment and 19 to a WW treatment. All plants were irrigated to a WW soil water content of 2.4 g H₂O g⁻¹ dry soil (soil water potential of -10 kPa) until they reached the V5 (five fully developed leaves) stage, after which water was withheld from WD plants until a soil water content of 1.4 g H₂O g⁻¹ dry soil (soil water potential of ± -100 kPa) was attained. Once the WD soil water content was reached, plants were irrigated daily to maintain this soil water content.

Thermal and hyperspectral image processing

Imaging of individual plants was performed daily between 7.30 and 14.00 during the DR and TF experiments (14,744 images of eight genotypes). These images were supplemented with 288 afternoon images that were captured on TF sampling days for physiological traits. These took place at two timepoints: the V5 stage + 10 days and the V13 stage ± 2 days (Additional file 2: Table S2). The images collected during the DR experiment were used (together with the TF dataset) to create E prediction models, while the TF images were also used to evaluate the drought detectability of thermal indices and to investigate the transferability of indices and models to multiple genotypes.

Thermographic data were processed using the ‘raw2temp’ function of the ‘Thermimage’ R package [46], which converts raw radiation values to temperature using standard equations applied in thermography. This function includes a correction for background long-wave emission (reflected temperature), which was estimated using crumpled and flattened aluminum foil (Lambertian surface) with an emissivity of 1. The ‘raw2temp’ function requires the distance between the plant and camera (m), T_a (°C), RH (%), an estimate of leaf emissivity, which was set to 0.96 [47], reflected temperature (°C), Infrared Window Temperature = 20 (°C), Infrared Window transmission = 1, and the constants PlanckR1 = 16,963.094, PlanckB = 1435.3, PlanckF = 1, PlanckO = − 4328 and PlanckR2 = 0.014514672 (Fig. 2). Deviations in the camera accuracy (camera ‘drift’) were determined using a black plate with thermocouple. If a consistent difference higher than 1°C between the temperature estimated by the camera and the thermocouple was observed, the imaged temperature was adjusted using the median difference in temperature measured by the two devices. After image pre-processing, the plant was segmented out of the background by aligning the thermal image with the corresponding segmented RGB image. RGB images were collected using an Allied Vision Technologies Prosilica GE4000C 11-megapixel camera (Allied Vision Technologies GmbH, Germany) equipped with a Canon EF 24 mm f/1.4L II USM lens (Canon Inc., Japan). RGB plant segmentation was performed with a convolutional neural network (CNN) model developed in Pytorch (Fig. 2). The CNN model was composed of a pre-trained DenseNet M161 backbone with added bilinear upsampling layers after each ‘dense’ layer, where each new layer is upsampled to the same resolution as the closest higher-resolution DenseNet output and then combined using a weighted sum. This is done until a feature map at the original image resolution is created. The training and validation of the CNN model was performed on 1MP foreground-containing crops of 255 manually annotated image-segmentation pairs. A subset of 10% of the images were randomly chosen for validation. Optimization was done for 7300 steps (using a batch size of 2) with binary cross-entropy loss and weight updates calculated using stochastic gradient descent with momentum (learning rate = 0.01 and momentum = 0.95). Additionally, several standard image augmentation, class weighting and image weighting techniques were used during training. More details on the CNN model can be found in De Meyer et al. [48]. The segmented and raw RGB images were downscaled from 2672 × 4008 to 1336 × 2004 pixel resolution for the thermal image analysis. Both thermal and RGB cameras were first calibrated using a metal plate with holes or a chessboard (Matlab, Simon Donné, IMEC—IPI—Ghent University), which resulted in intrinsic and extrinsic camera matrices. Due to a slight shift in the RGB camera position, the extrinsic matrix could not be utilized. Instead, the translation parameters of the extrinsic matrix were expressed in function of the distance between the RGB camera and the plant. This was done because the alignment of the TIR and RGB images was affected by the position (height) of the plant. The transformation matrices were used to align the thermal and RGB images (Fig. 2). This required the distance between the plant and the cameras, which was not available for each pixel. Instead, the distance between 2/3 of the plant height and the cameras was used. Using 2/3 of the plant height as an average depth value for all plant pixels did not perform optimally for large plants because the leaves were positioned further apart. To resolve this issue, an erosion was applied on the segmentation filters (binary image) to remove remaining background pixels (‘erode’ function of the ‘mmand’ R package [49]). The median temperature of each plant was calculated and used for drought detection and E predictions (Fig. 2).

Fig. 2

Schematic representation of the image processing and transpiration rate prediction model development pipeline. The top light gray square summarizes the image processing for the Red–Green–Blue (RGB: green), Thermal Infrared (TIR: orange) and hyperspectral visible and near-infrared (VNIR) and shortwave-infrared (SWIR) (blue) imaging systems. The RGB image processing consisted of a plant segmentation using a convolutional neural network (CNN). The TIR image processing included a radiation to temperature (T) conversion, TIR plant segmentation by aligning RGB and TIR images and the calculation of median plant temperature (T_p) and TIR indices. The hyperspectral processing consisted of a radiometric calibration and plant segmentation using the Red-edge Normalized Difference Vegetation Index (Re-NDVI) and a random forest (RF) model for VNIR and SWIR, respectively. After the hyperspectral segmentation, a correction of the distance between the white reference plate and the plant was performed followed by a brightness classification of the VNIR data and a VNIR-SWIR alignment to extract illumination classes from the data. Indices and reflectance were calculated for the intermediate light class [42]. The bottom gray square summarizes the different modeling approaches applied to develop transpiration rate (E) prediction models. Mechanistic energy balance and empirical models were developed by combining different data types including environmental data (env: pink), TIR data and hyperspectral data (hyp). The empirical modeling algorithms were Least Absolute Shrinkage and Selection Operator (LASSO), Random Forest (RF), linear model with stepwise selection, Partial Least Square Regression (PLSR) and linear models combining one TIR index with environmental data and hyperspectral indices

Hyperspectral data were collected by two pushbroom line-scanners in a second imaging cabin in PHENOVISION. The visible and near-infrared (400–1000 nm) and shortwave-infrared (970–2500 nm) scanners had a spectral resolution of 3 and 6.3 nm, respectively. The specifications of this imaging system have been described in Mertens et al. [42]. The pre-processing included image calibration, plant segmentation and illumination classification (Fig. 2) as detailed in Mertens et al. [42]. The illumination classification was not able to remove all the light effects, because the distance between the white reference plate, which is used to calibrate the hyperspectral images, and the top of the plant was too big for large plants. This issue was solved by correcting the reflectance distribution of a plant (pixel level) to the distribution it would have had at the white reference height tile in the PHENOVISION imaging cabin (at 1.2 m from the imaging system). Three plants were therefore imaged at different lift heights and the relationship between average reflectance and the distance between plant and camera was determined (Fig. 3). This relationship was then used to adjust the reflectance distribution on a pixel level with the following formula.

Fig. 3

Relationship between reflectance and the distance between the plant and the camera. Average plant reflectance was scaled between zero and one to remove biological variation between the plants. Three plants were used to determine this relationship, indicated with gray, yellow and blue dots. The average relationship and the 95% confidence interval are indicated by a black line and gray shading, respectively

$${\rho }_{adj i}={\rho }_{i}*\frac{{\rho }_{scaled white i}}{{\rho }_{scaled current i}}$$

(1)

In this equation, ${\rho }_{adj i}$ is the adjusted reflectance of wavelength i, ${\rho }_{i}$ corresponds to the unadjusted reflectance, ${\rho }_{scaled white i}$ represents the scaled average reflectance of the plant positioned at the white reference and ${\rho }_{scaled current i}$ is the scaled average reflectance at the actual distance from the camera. The reflectance distribution correction was performed before illumination classification and reflectance averaging (Fig. 2). To simplify the analysis, only one light class (intermediate light class) was selected based on the percentage of pixels per plant it contained during all experiments and the ability to detect drought [42].

Physiological trait measurements

Gas exchange measurements were collected once per week between 7.30 and 14.00 in the DR experiment. Five plants per treatment (total of 25) were measured 5–10 min before imaging (Additional file 2: Tables S2 and S3). The gas exchange measurements were collected using a portable LICOR 6400-XT infrared gas analyzer (LI-COR Biosciences, USA). A steady-state CO₂ level of 400 µmol mol⁻¹ was maintained in the leaf chamber, while temperature and PAR were adjusted to the greenhouse temperature (25–31°C) and PAR (50–700 µmol photons m⁻² s⁻¹) conditions at the time of measurement.

The same measurement approach and schedule was applied in the TF experiment, in which fluorescence and gas exchange measurements were collected on three plants per genotype-treatment combination (total of 48 plants). Here, additional physiological measurements were collected in the afternoon (13.00–18.00) at V5 + 10 days (acute drought phase) and V13 ± two days (13 fully developed leaves) for nine plants per genotype-treatment combination (total of 144 plants, Additional file 2: Tables S2 and S3). Non-destructive fluorescence and gas exchange measurements consisted of effective quantum yield of photosystem II (ϕPS₂), energy harvesting efficiency by oxidized PSII (F_v’/F_m’), E and g_s of H₂O. In addition, leaf ψ of a top leaf, which was clearly visible in the image, was destructively measured 5–10 min after imaging using a PMS model 1000 pressure chamber (PMS Instrument Company, USA). Six plants per genotype-treatment combination (96 plants in total) were selected to monitor ψ during the TF experiment. Each individual plant was measured every two weeks resulting in 48 measurements per week. Additional ψ measurements were also collected at V5 + 10 days and V13 ± two days (144 plants, Additional file 2: Table S2 and S3).

Indices

Thermal infrared indices were calculated to detect plant WD and to estimate E for different maize genotypes. The indices combined T_p with T_dry (black plate visible in the image) or environmental measurements collected in the growth zone, and inside and outside of the imaging cabin (Fig. 1). Six existing and two new indices, T_BS-T_p and temperature ratio index (TRI), were evaluated (Table 1). The new T_BS-T_p is an adaptation of T_dry-T_p that was tested, because T_dry as well as T_BS measured the combined effect of incoming radiation and T_a on an object’s temperature. The temperature ratio index (TRI) was evaluated because it could detect WD treatments and because it correlated with g_s and E. The ICWSI and CWSI indices were calculated using the baseline approach. This required the estimation of the maximum temperature or minimum CTD of a non-transpiring plant and the minimum temperature or maximum CTD of a fully transpiring plant. The maximum temperature was assumed to be T_a + 3°C based on the maximum CTD observed in this study and literature [50]. The T_p or CTD of a fully transpiring plant was determined with a baseline that related T_p or CTD to VPD for each genotype (12 WW plants, 38 ± 4 images per genotype). The baseline models were created using the ‘lm’ function of the ‘stats’ R package [51]. For both the ICWSI and CWSI, a development-corrected baseline was developed, while for the CWSI, a second baseline model was developed that corrected for non-constant environmental conditions (Additional file 2: Fig. S1). Baselines were therefore created for young and mature plants separately (ICSWI, CWSI_dev), or by adding T_a as a predictor to the baseline models (CWSI_Ta, [26]).

Table 1 Thermal infrared and hyperspectral reflectance indices evaluated in this study

Hyperspectral indices were evaluated by determining their relationship with E and g_s. Ten indices were selected for this analysis, of which six were correlated with E and four with water content or ψ (Table 1). Correlations of indices and environmental data with E and g_s were calculated using the ‘cor.mtest’ function of the ‘corrplot’ R package [52]. Indices and environmental data that showed strong correlations were subsequently used to create linear prediction models for E.

Transpiration rate modeling

Different combinations of datasets, including data on T_p and TIR indices, environmental and hyperspectral data, and modeling approaches were evaluated (Fig. 2). Linear models were created to relate T_p, indices and environmental data with E. Index-based linear models had E as dependent variable and an index as independent variable, while correlation-based linear models combined one strongly correlating index with environmental data (‘lm’ function of ‘stats’ R package). The significant contribution of different predictors to the linear model was evaluated using the Chi-square test of the ‘anova’ function (‘stats’ R package). To investigate if combining multiple indices with T_p and environmental data could improve prediction accuracy, RF, LASSO and stepwise selection models were created using the ‘randomForest’, ‘glmnet’ and ‘olsrr’ R packages, respectively [63,64,65]. The RF model hyperparameters were fine-tuned by leave-one-out cross validation (LOOCV) and the optimal number of predictors was determined with the ‘VSURF’ R package [66]. To combine thermal, environmental and hyperspectral data, correlation-based linear models, that combine TIR and hyperspectral indices, RF, LASSO, stepwise selection and PLSR models were evaluated [67]. The PLSR models were created with the PLSR function of the ‘pls’ R package as described in Mertens et al. [42]. The E prediction models were developed using ordered quantile normalized data (‘orderNorm’ function, ‘bestNormalize’ R package, [68]). The prediction accuracy of the models was evaluated by calculating the ‘out-of-bag’ test Mean Absolute Percent Error (MAPE, formula 2), the Root Mean Square Error (RMSE) and R-squared (R²) (‘postResample’ function, ‘caret’ R package, [69]) of 100 bootstrap samples. A final model was created using the whole B104 dataset. This model was subsequently used to test the transferability of the models to different genotypes.

$$MAPE=\frac{1}{n}\sum_{t=1}^{n}\left|\frac{{A}_{t}-{F}_{t}}{{A}_{t}}\right|*100$$

(2)

The description and results of the mechanistic energy balance model that was evaluated in this study can be found in the Additional file 1.

Drought detection and statistics

The effectiveness of TIR indices in detecting the occurrence of plant WD was evaluated by comparing when significant differences between the WW and WD treatments were first detected and for how many days. To take variation in the start of the WD treatment in the TF experiment into account, time was expressed as days after V5 (start drought). Treatment differences were tested for each day and for each genotype (start experiment: n = 19, end: n = 5) with a factorial ANOVA using the ‘lsmeans’ R package [70]. Factorial ANOVA was also used to test treatment differences in E, g_s, ψ and ϕPS₂ for every group of four days and genotype (n = 3, can be higher on sampling days). The days were grouped together because measurements at each timepoint were spread over a few consecutive days. The physiological data were split into morning (< 13.00) and afternoon (> 13.00) measurements to take the time-of-day effect on physiology into account. The TIR index data were not split into morning and afternoon, as the majority of the images were collected in the morning. All the p-values were corrected using the sidak method (‘MHTmult’ R package, [71]).

Genotypic differences in the sensitivity of TIR indices and physiological measurements to drought stress effects were tested for each day after V5 using the ‘contrast’ function in the ‘lsmeans’ R package [70]. The genotype analysis of the physiological traits could only be performed at 10 days after V5 (n = 8–10 for each treatment). This is because the sample size was too small before that timepoint because not all plants had reached their target WD. The p-values of this analysis were corrected using the sidak method (‘MHTmult’ R package, [71]).

Results

Drought affects the water-use behavior of B104 maize

In the TF experiment, a WD treatment was applied on B104 maize from the V5 developmental stage onward. It consisted of two stages, namely an acute drought stage, in which water was withheld until the lower soil water content was reached (V5 + maximum 10 days), and a steady-state drought during further development, when the plants were grown continuously at this lower soil water content. Plant physiology was most strongly affected during and shortly after the acute drought stage, in which a decrease in the g_s, E, ϕPS₂ and ψ was observed in the WD treatment (Fig. 4A). The treatment differences were only significant for leaf ψ and g_s measured in the afternoon at 10 days after the onset of drought. The absence of significant differences in the other physiological traits and timepoints may be related to the limited number of available datapoints. The difference between WD and WW plants became less pronounced during the steady-state drought as the WD plants adapted to the lower soil water content (Fig. 4A). This adaptation included a decrease in biomass, as the WD B104 plants had on average 39% lower fresh weight at silking compared to WW plants (WW: 615.3 ± 70.9 g, WD: 375.3 ± 34.2 g), resulting in lower water requirements during the steady-state drought.

Fig. 4

Responses of physiological traits and thermal infrared indices to drought. (A) Drought effects on plant physiology during the TF experiment. The B104 genotype and four physiological traits, i.e. stomatal conductance (g_s, mol H₂O m⁻² s⁻¹), effective quantum yield of photosystem II (ϕPS₂), transpiration rate (E, mmol H₂O m⁻² s⁻¹) and leaf water potential (ψ, MPa), were selected for this analysis. (B) Responses of thermal infrared (TIR) indices to drought. B104 maize plants (n_WD: 19, n_WW: 19) were imaged daily and TIR indices were calculated using the formulas described in Table 1. Well-watered (WW) and water-deficit (WD) maize plants were monitored from V4 until the silking stage. The average trends of the WW and WD treatments are indicated by blue solid and red dashed lines, respectively. The 95% confidence interval of the average is represented by gray shading (A) and the standard deviation by blue and red shading for the WW and WD treatments, respectively (B). Individual measurements of the plants are visualized by blue dots (WW) and red circles (WD). The black vertical line indicates the start of the WD treatment. The days on which significant treatment differences were observed are marked with a light gray vertical shading behind the average trend and dots (P < 0.05)

Thermal infrared indices can detect drought

All indices that normalized T_p using environmental variables were calculated for each of the three environment monitoring positions to determine which one was the most suitable for drought detection. They showed slight differences in VPD and T_a. During the TF experiment, VPD and T_a were the highest outside the imaging cabin (VPD_out: 2.2 kPa, T_{a_out}: 27 °C), followed by the positions inside the cabin and in the growth zone (VPD_in: 1.9 kPa, T_{a_in}: 25.53 °C, VPD_gz: 1.8 kPa, T_{a_gz}: 25 °C) (Additional file 2: Fig. S2). Consequently, indices produced slightly different values depending on the monitoring position with those using VPD and/or T_a measured inside the imaging cabin performing better (Additional file 2: Fig. S2; Table 2). The indices that included VPD and/or T_a collected outside the cabin or in the growth zone were excluded to simplify the drought detection analysis.

Table 2 Effect of the environmental measurement location on the drought sensitivity of the thermal infrared indices

Thermal indices also showed stronger differences between treatments during the acute drought period and a recovery during the steady-state drought. The CSI demonstrated the strongest drought stress effects with significant treatment differences from 4 until 23 days after the onset of drought (Fig. 4B). It was followed by the indices TRI, CTD and T_dry-T_p in which the drought stress effects lasted from 4 until 20 days after the onset of drought, and finally ICWSI in which the first drought detection occurred at 8 days and lasted until 23 days. The indices CWSI_dev, CWSI_Ta and T_BS-T_p were less or not sensitive to drought as significant effects were absent (CWSI_Ta) or only visible on a few days (CWSI_dev and T_BS-T_p, Fig. 4B). The thermal indices were not only influenced by soil WD, but also showed a trend during the morning and early afternoon with lower CSI and TRI values for WW plants around noon compared to the morning (Fig. 5). These diurnal changes corresponded with an increase in VPD and T_a in the greenhouse and a decrease in CTD. The decrease was less pronounced for WD plants, which often resulted in larger treatment differences around noon than in the morning. Trends were also visible in the physiological measurements, where a decrease in ψ was observed during the morning, which was more pronounced for WD compared to WW plants (Fig. 5G). Overall, TIR indices were able to detect drought stress effects in an automated phenotyping platform as long as environmental data, such as T_a and VPD, were used to normalize T_p.

Fig. 5

Daytime trend in environmental data, thermal infrared indices and leaf water potential of B104 plants. For the environmental data and indices, one day during the acute drought period was selected that corresponded with about 5 days after the onset of drought. For leaf water potential (ψ), all measurements collected during the experiment were used because of the limited number of B104 measurements per day. Average and individual environmental measurements (A–C) are shown with black lines and black dots, respectively. The average index and ψ values of the well-watered (WW) and water-deficit (WD) treatment are represented by a blue line and a red dashed line, respectively (D–G). The individual plants are visualized with blue dots (WW) and red circles (WD). Daytime patterns of meteorological data: (A) PAR in the growth zone (PAR_gz), (B) VPD inside the cabin (VPD_in), and (C) air temperature inside the cabin (T_{a_in}). Morning patterns of TIR: (D) CSI_in index, (E) TRI_in index, and (F) CTD_in index. (G) Daytime trend in leaf ψ

Environmental data and thermal indices can predict transpiration rate

Besides drought detection, TIR indices can be used to monitor g_s and E. In the PHENOVISION setting, relatively strong correlations with g_s and E were observed for T_dry-T_p (r_E = 0.61, r_gs = 0.50, P < 0.05), CWSI_Ta (r_E = − 0.61, r_gs = − 0.58, P < 0.05) and CTD (r_E = − 0.63, r_gs = − 0.64, P < 0.05) (Fig. 6B). In the case of CWSI_Ta and CTD, the highest correlations were found when environmental data measured outside the imaging cabin were used. This is also the location where the gas exchange measurements were performed. Consequently, also E and g_s showed the highest correlations with environmental data collected at this position (Fig. 6A). PAR and T_BS, which were only measured in the growth zone, were also significantly correlated with E and g_s (r_E-PAR = 0.54, r_E-TBS = 0.58, r_gs-PAR = 0.54,r_gs-TBS = 0.54, P < 0.05). The correlated environmental data and indices were combined to create empirical models that predict E.

Fig. 6

Correlation of stomatal conductance and transpiration rate with measured independent variables. The independent variables included (A) environmental variables: temperature of a non-transpiring plant (T_dry), photosynthetically active radiation (PAR), black sphere temperature (T_BS), air temperature (T_a), relative humidity (RH) and vapor pressure deficit (VPD) measured in the growth zone (gz), and inside (in) and outside (out) the cabin, (B) thermal infrared indices (Table 1), and (C) relative reflectance at 532, 992, 1407 and 1881 nm and hyperspectral indices (Table 1). Significant (P < 0.05) correlations are indicated in blue (positive) and red (negative). g_s stomatal conductance, E transpiration rate

The indices with the strongest correlations, CWSI_{Ta_out}, CTD_out and T_dry-T_p (Fig. 6B and Additional file 2: Fig. S3), had the highest prediction accuracy of all the index-based models with an RMSE of 0.68, 0.66 and 0.61 and an R² of 0.39, 0.44 and 0.41, respectively (Table 3, Fig. 7). The test prediction accuracy of these models was improved by combining the indices with environmental data, such as VPD_out, PAR and T_BS (RMSE_{CWSI_TaOut-env} = 0.57, RMSE_{CTD_Out-env} = 0.60, RMSE_Tdry-Tp-env = 0.56, R²_{CWSI_TaOut-env} = 0.52, R²_{CTD_Out-env} = 0.57 and R²_Tdry-Tp-env = 0.55, Table 4, Fig. 7). VPD_out was the only environmental variable that improved the prediction accuracy of all three index-based models. PAR and T_BS were added to the index-based models because they relate to the incoming radiation, which may additionally increase T_p independently of T_a.

Table 3 The prediction accuracy of normalized thermal infrared indices

Fig. 7

Relationship between measured and predicted transpiration rate (E). E was measured with a portable LICOR infrared gas analyzer. (A) The first row visualizes the test prediction accuracy of thermal infrared (TIR) index-based models (CTD_out, CWSI_{Ta_out}, and T_dry-T_p), while the second row shows the performance of these indices when they are combined with environmental data not included in the index (Table 4). The third row contains two figures that show the accuracy of models that combine one TIR index (CTD_out or CWSI_{Ta_out}) with environmental data and hyperspectral indices (RMP_1483/1430, WBI). This type of model was not created for the T_dry-T_p index, as hyperspectral indices did not significantly contribute to this model (Chi-square test, P < 0.05). (B) Prediction accuracy of the LASSO, RF and PLSR models that combine multiple TIR indices with hyperspectral wavelengths and/or environmental data. The semi-transparent black dots represent the predicted and measured E of the 100 bootstrap samples, while the black line shows the one-to-one relationship. The test RMSE and R² are added to each scatterplot

Table 4 Prediction accuracy of transpiration rate prediction models that use one thermal infrared index

To evaluate the impact of combining different indices and environmental data on the prediction accuracy of E, stepwise selection, LASSO and RF models were developed and compared with the single index models. Including multiple indices and environmental variables could indeed improve the prediction accuracy even further. The LASSO model had the highest accuracy of all TIR models with a median test RMSE = 0.47 and test R² = 0.63 (Table 5, Fig. 7). Interpreting the variables of a LASSO model is difficult as it removes redundant variables with potential biological relevance. The variables of the RF model (second-best) were therefore evaluated (Table 5). RF selected four environmental variables (T_{a_out}, VPD_gz, PAR_gz and T_{BS_gz}) and 11 indices (CTD_out, CSI_in, CSI_out, ICWSI_in, CWSI_{dev_in}, CWSI_{dev_out}, CWSI_{dev_gz}, CWSI_{Ta_gz}, T_BS-T_p, TRI_in and TRI_gz). Most of these variables were significantly correlated with E, except for the CSI indices, ICWSI_in and the indices using gz data. The RF model incorporated environmental data of all three monitoring positions directly or indirectly (through indices). This may indicate that all three monitoring positions capture relevant information to predict E in PHENOVISION. The value of the different monitoring positions may be explained by the fact that the plants were not stabilized and acclimatized to the environment of the cabin during TIR imaging. In PHENOVISION, maize plants (and E) can adjust to the growth zone environment and the waiting area outside the imaging cabin, as they remain at these positions for some time, while in the imaging cabin they are immediately imaged to avoid acclimatization to an irrelevant environment. Consequently, the transportation of the plants around the platform increases the need of monitoring the environment at multiple locations.

Table 5 Accuracy of transpiration rate prediction models based on multiple thermal infrared indices

Hyperspectral data improve transpiration rate prediction accuracy

Combining TIR and hyperspectral data may improve the accuracy and robustness of E prediction models, as significant correlations between hyperspectral data and g_s or E were observed (42; Fig. 6C). The wavelength regions that correlated with E were located around 523 nm, 992 nm, 1407 nm, 1881 nm and 2346 nm. All wavelengths, except for 523 nm, are related to water content and leaf anatomy [32, 72]. Six out of 10 selected indices (Table 1) correlated significantly with E (Fig. 6C, P< 0.05). The strongest negative relationship was present between E and the WBI index (r = − 0.52), while the most prominent positive correlation was observed for RMP_1483/1430 (r = 0.53). Adding hyperspectral indices (WBI and RMP_1483/1430) to TIR index-based models significantly improved the CTD_out and CWSI_{Ta_out} models with about 11% and  9%, respectively (Table 4, Fig. 7, Chi-square test, P < 0.05). The absence of improvement in the T_dry-T_p model may be related to the limited amount of hyperspectral information (wavelengths) contained in these indices. This constraint was reduced by training machine learning models (PLSR, RF, LASSO, step-wise selection) on the complete hyperspectral, environmental and TIR datasets. The LASSO and RF models with hyperspectral data were compared to the previously developed LASSO and RF E prediction models to evaluate the contribution of the hyperspectral imaging system. The prediction accuracy was improved to a small extent: the test RMSE of the RF model was slightly reduced by 10%, while the test R² of both models was improved by 8 ± 2% (Table 5, Fig. 7). The best performing model was the LASSO model, which had a RMSE of 0.49 and a R² of 0.69. All algorithms selected environmental and thermal data, as well as hyperspectral information. The T_BS, CTD_out and reflectance around 500 nm were selected by all model algorithms, while PAR, VPD_out, T_dry–T_p, ICWSI_in and reflectance around 2500 nm were predictors in three out of four models. All the variables selected by the models are listed in Table 5. The wavelengths selected by the E models were mainly related to water content and leaf anatomy (around 2500 nm) or photosynthesis/pigment content (around 500 nm). The hyperspectral models also included environmental data from all three monitoring positions directly or indirectly. Including environmental data into the E prediction models was pivotal as LASSO models created using only image-based data (hyperspectral and/or T_p) had an approximately 30% lower prediction accuracy than the LASSO model that incorporated all available data (Table 5).

Validating thermal indices and models on other maize genotypes

All empirical E prediction models developed so far were trained on the B104 genotype. The transferability of these models to other genotypes (H99, MS71, NC358, OH43, TX303, TZi8 and W153R) was uncertain, as genotypes may differ in drought sensitivity, water-use behavior, leaf anatomy and reflectance. In this study, all genotypes showed a similar drought response in measured physiological traits with larger treatment differences in the period of acute drought. The magnitude of these drought stress effects differed among genotypes (Fig. 8). Leaf ψ and ϕPS₂ showed the largest genotypic differences with more pronounced treatment effects in NC358, OH43, TZi8 and W153R compared to B104, H99, MS71 and TX303. However, only TX303 differed significantly from NC358, OH43 and W153R at 10 days after the onset of drought. Drought stress effects on E and g_s were less variable between genotypes. Treatment differences were only slightly larger in B104, MS71, NC358, TX303 and TZi8 compared to H99, OH43 and W153R (Fig. 8). This difference was mainly caused by higher E and g_s values in WW plants, while the treatment differences in ψ and ϕPS₂ resulted from decreases in the WD plants. The most drought-sensitive indices in B104 (CSI_in and TRI_in) were able to detect drought stress effects in all eight genotypes and showed differences between genotypes. TX303 and TZi8 had the strongest drought stress effects, as significant treatment differences were observed during 28.5 ± 1.5 and 25.5 ± 1.5 days of the TF experiment, respectively. They were followed by B104, MS71 and NC358 with 20.5 ± 1.5, 20.5 ± 0.5 and 20.5 ± 0.5 days, respectively. The weakest drought stress effects were observed in H99, OH43 and W153R, for which only 15.5 ± 0.5, 8.5 ± 0.5 and 10 ± 1 days of the experiment showed significant treatment differences (Fig. 9).

Fig. 8

Genotypic differences in the responses of physiological traits to drought stress. The genotypes H99, MS71, NC358, OH43, TX303, TZi8 and W153R and four physiological traits, i.e. stomatal conductance (g_s, mol H₂O m⁻² s⁻¹), effective quantum yield of photosystem II (φ_PS2), transpiration rate (E, mmol H₂O m⁻² s⁻¹) and leaf water potential (ψ, MPa), were selected for this analysis. Physiological traits of well-watered (WW) and water-deficit (WD) maize plants were monitored from V4 (four fully developed leaves) until the silking stage. The average trends of the WW and WD treatments are indicated by blue solid and red dashed lines, respectively. The gray shading around the lines represents the 95% confidence interval of the average. Individual measurements of the plants are visualized by blue dots (WW) and red circles (WD). The black vertical line indicates the start of the WD treatment. The days on which significant treatment differences were observed are marked with a light gray vertical shading behind the average trends and dots (P < 0.05)

Fig. 9

Responses of thermal infrared indices to drought. The H99, MS71, NC358, OH43, TX303, TZi8 and W153R maize plants were imaged daily and thermal infrared indices, including canopy temperature depression (CTD_in), canopy stress index (CSI_in), simplified stomatal conductance index (T_dry-T_p), temperature ratio index (TRI_in), Idso's crop water stress index (ICWSI_in), crop water stress index (CWSI_{Ta_in}, CWSI_{dev_in}) and T_BS-T_p, were calculated using the formulas described in Table 1. The averages of the well-watered (WW) and water-deficit (WD) treatments are indicated by blue and red dashed lines, respectively. The blue and red shading represents the standard deviation around the mean. The black vertical line indicates the start of the WD treatment. The grey vertical shading indicates the days on which significant treatment differences were observed (P < 0.05, n=19)

The relationship between E and TIR indices across genotypes influenced the prediction accuracy of the E prediction models. Most models were able to make reasonable predictions of E for multiple genotypes, except for the LASSO model trained on TIR and environmental data, and the stepwise selection model that combined the environmental information with both hyperspectral and TIR data (Figs. 10 and 11). These models had small R² (< 0.1), high RMSE (> 1) and MAPE (> 90) values for certain genotypes. The high MAPE value indicated that the prediction error was higher than 90% of the measured value. The LASSO and stepwise selection model had an acceptable prediction accuracy for B104, suggesting that they are not transferable to other genotypes. When comparing the remaining models, the RF, LASSO and PLSR algorithms performed better than the other modeling approaches for most genotypes (Figs. 10 and 11). As the RF algorithm was trained with and without hyperspectral data, the contribution of this data type to the robustness and transferability of E prediction models could be evaluated. Adding hyperspectral data slightly improved the prediction accuracy for five out of eight genotypes, suggesting that the additional information contained in these data increases the robustness of E predictions (Figs. 10 and 11). The RF model that combined both imaging data (TIR and hyperspectral) was further used to compare the prediction accuracy of different genotypes. The R² of all seven genotypes, other than B104, was lower than the R² of B104 (Fig. 10); however, the RMSE and MAPE values were more similar (Fig. 11). The lower R² values may result from the larger E range in B104 compared to most other genotypes (Fig. 12). The genotypes could be subdivided into three groups based on the performance of the RF E model. The first group contained two genotypes, NC358 and W153R, for which the E model performed similar to B104 with an error of approximately 34% of the measured value (MAPE: ± 34% and RMSE > 0.59, Figs. 11 and 12). The second group was more variable containing both low and high prediction errors. H99 had a low MAPE of only 22%, but a high RMSE of 0.69, which may be caused by the lower prediction accuracy of the higher E values (Fig. 11). In contrast, TZi8 and OH43 had a relatively high MAPE value (TZi8: 48%, OH43: 57%) combined with a relatively low RMSE value (TZi8: 0.62, OH43: 0.61). Figure 12 shows that the E model slightly overestimated the measured E values for TZi8 and OH43. The last group contained MS71 and TX303 and showed the highest variability in the scatterplot comparing predicted to measured E values (Fig. 12). The prediction accuracy of these genotypes was low with an RMSE of 0.78 and a MAPE of 29% for MS71 and an RMSE of 0.65 and MAPE of 57% for TX303.

Fig. 10

R-squared accuracy of transpiration rate prediction models for different genotypes. Eight genotypes were compared: B104, H99, MS71, NC358, OH43, TX303, TZi8 and W153R. The test prediction accuracy of the 100 bootstrap samples was used for B104, while for the other genotypes, the accuracy was determined with a test set that was not included in the training of the models. The shading in this figure represents the R² value with darker blue colors corresponding to higher accuracies

Fig. 11

Accuracy of transpiration rate prediction models for different genotypes. Eight genotypes were compared: B104, H99, MS71, NC358, OH43, TX303, TZi8 and W153R. The test prediction accuracy of the 100 bootstrap samples was used for B104, while for the other genotypes the accuracy was determined with a test set that was not included in the training of the models. The prediction accuracy shown in this figure are (A) the root-mean square error (RMSE) and (B) the mean absolute percent error (MAPE). The shading in this figure represents the value of the accuracy measure with darker red colors corresponding to higher values and lower accuracies

Fig. 12

Relationship between measured and predicted transpiration rate for different genotypes. The Random Forest models that combined thermal, hyperspectral and meteorological data were used to predict transpiration rate (E). Genotypes were grouped based on the prediction accuracy of the model. The accurate group contained three genotypes (B104, NC358 and W153R) and is indicated by a light blue header. The moderately accurate group has an orange header and contains H99, OH43 and TZi8, while the less accurate group is marked with a red header and contains MS71 and TX303. The dots represent the predicted and measured E of individual plants, while the black line shows the one-to-one relationship. For B104, leave-one-out cross-validation predictions are shown, while the predictions of the other genotypes were created using a test set

Discussion

Thermal indices detect drought in an automated phenotyping platform

The use of TIR indices to detect drought is well established under field and orchard conditions, as this approach has been applied in irrigation/drought studies of potatoes, sesame, cotton, coffee, wheat, maize, forage grasses, pistachio, olive, peach and apple [7, 12, 16, 19, 73,74,75,76,77,78,79]. The few studies that used TIR indices in greenhouse experiments [7, 19, 77, 80] showed that both the direct and empirical CWSI approaches, the stomatal conductance index (I_g) and CTD indices were able to detect drought in sesame, wheat, maize and peach trees. So far, the use of these indices is less common in indoor whole-plant phenotyping platforms. Here, the CWSI index was not able to detect mild drought differences, which may be related to the lower accuracy of the non-water stressed baseline, as more variation was observed around the regression of VPD with CTD compared to the one with T_p (adj. R²_{CTD_in-baseline}: 0.51, adj. R²_{Tp_in-baseline}: 0.66, Additional file 2: Fig. S1). The baseline relating VPD to T_p was used to calculate the ICWSI index, which in contrast to the CWSI, detected treatment differences eight days after the onset of drought. Increasing the sample size used to develop the CWSI baseline may improve its performance. In addition to the lower accuracy of the baseline, inaccurate estimations of the CTD of a non-transpiring plant (3°C) may reduce the drought sensitivity of the CWSI. Alternative methods to estimate this CTD value include measuring T_a and T_p of a leaf covered with petroleum jelly or calculating CTD based on the difference between the saturated vapor pressure at T_a and at the temperature of a non-transpiring plant, which can be estimated by adding T_a to the intercept of the relationship between VPD and CTD of a fully transpiring plant [7, 26, 81].

The drought sensitivity of the CTD_in index was confirmed in this study, as this index detected significant drought stress effects four days after the onset of WD. Nevertheless, it was outperformed by the CSI_in and TRI_in indices, which showed significant differences on the second day. The performance of the CSI_in index may be related to the diurnal changes in VPD that were created in the greenhouse, and that induced a larger range in E. As a consequence, T_p values could not be normalized by T_a alone. The CSI_in index incorporated the required correction for RH by normalizing CTD with VPD [50]. It was originally developed to detect drought stress in wheat fields and was shown to capture physiological responses both at the single-plant leaf and population-field level [54].

Transpiration rate prediction based on thermal indices or an energy balance approach

Plant temperature and TIR indices have been related to different water-use behavior traits, such as g_s, E, relative water content, and stem and leaf water potential [12, 13, 15, 19, 23, 28, 29, 77, 82]. One of the drought-sensitive indices in this study, CSI, has been correlated with g_s and E in the field [54], but this was not confirmed in this study under greenhouse conditions. CSI_in was not correlated with the gas exchange measurements, while CSI_out was non-linearly related to these traits. A negative relationship between CSI_out and g_s and E was observed when the index values exceeded -1.5, which corresponds with the results observed by Rodriguez et al. [54]. Points that deviated from this negative relationship (CSI_out < − 1.5) had E values < 2 mmol m⁻² s⁻¹ and were collected in the morning when stomata were not yet completely open (Additional file 2: Fig. S3). The opening of stomata is triggered by light [83] and it requires some time for all stomata to open and for E to reach its maximal value. During this time period, the relationship between VPD and T_p is less strong or even absent [50], which may explain the lower performance of the CSI_out index during these timepoints. The TIR indices CTD, CWSI_Ta and T_dry-T_p (r_{CTD_out}: − 0.63–− 0.64, r_{CWSI_Ta_out}: − 0.58–− 0.61, r_Tdry-Tp: 0.5–0.61) did have a relatively strong correlation with the gas exchange measurements. Negative relationships between CTD or CWSI and g_s have been observed in olive trees, sesame, strawberry, grapevine, both in field and greenhouse setups [13, 19, 23, 28, 84]. Maes and Steppe [26] simulated the relationship between g_s, CTD and the direct CWSI approach for different environmental conditions and observed a non-linear relationship between these indices and g_s. Nevertheless, several studies have described a linear or almost linear relationship of CTD and CWSI with g_s [19, 23, 57, 84]. Some of these more linear relationships may have resulted from a low g_s range, because g_s is almost linearly related with T_p and CWSI, as long as its value is smaller than < 0.5 mol m⁻² s⁻¹ [18]. The other strongly correlated index T_dry-T_p had a positive relationship with g_s and E. It is a simplified version of the stomatal conductance index (I_g,[57]) and describes an almost linear relationship with g_s, which was also visible in this study (Additional file 2: Fig. S3).

The relationship between physiological traits and TIR indices can be used to develop empirical prediction models, but this was done in only a few studies for g_s, and stem and leaf ψ [12, 27, 85, 86]. Environmental factors such as T_a, PAR, wind speed and VPD, affect the relationships between water-use behavior traits and TIR indices [26]. Their variation within and between days under field conditions reduces the reliability of these models. In indoor phenotyping platforms, environmental conditions are semi-controlled, making the use of these types of models more feasible. Here, TIR index-based models had higher prediction accuracies when indices were combined with environmental data such as VPD and PAR. The best performing model was the RF model that incorporated TIR indices and environmental variables collected directly or indirectly at the three monitoring positions. Environmental monitoring outside the growth zone of the platform is not required in phenotyping platforms with mobile imaging systems in which the cameras move to the plants. In these platforms, imaging is performed in a less controlled light environment and plants are often grown in a canopy structure to simulate field conditions, which will increase the influence of shading on T_p and E predictions. Within-canopy illumination variation will induce leaf temperature variation that is unrelated to changes in g_s and E, and leads to scatter in the canopy temperature-physiology relationships. Sunlit pixels, which have a more pronounced slope between CTD and g_s [26, 87], can be extracted from the images [88, 89] to reduce the non-biological variation in temperature. The potentially higher variability in T_p and g_s may nevertheless mask subtle differences in the plant water status [89, 90]. More research is needed to elucidate the effects of illumination variation in indoor thermography applications and to develop optimal approaches that can remove these effects.

Empirical models are platform-specific and cannot be transferred to the field, where meteorological conditions are more variable [91].The energy balance approach does not have this limitation. This mechanistic model has been evaluated in field studies [30, 31, 92], but its use in indoor phenotyping platforms is less common because of the difficulty of distinguishing longwave radiation from the surrounding environment [93]. The need of measuring net isothermal radiation can be eliminated using a simplified version of the energy balance approach that incorporates the temperature of a non-transpiring leaf (dry reference surface). Here, this approach had a lower accuracy compared to the empirical models due to incorrect model assumptions (Figs. 10 and 11). Additional file 1 provides a detailed discussion of the issues observed in the implementation of the energy balance approach and suggestions for optimization. One important aspect of the energy balance approach is the steady-state assumption, which states that the environmental conditions are stable during imaging and that leaf temperature is at equilibrium with its environment. This assumption is difficult to maintain in the field because the environmental conditions are constantly changing. Indoor phenotyping systems with semi-controlled environments seem more appropriate in this sense, but here as well, factors such as plant transportation, changing light levels, and diurnal temperature and VPD trends, influence the acclimatization of the plants to the imaging environment and to reach the steady-state condition. Vialet-Chabrand and Lawson [93] propose a method that predicts g_s and E under dynamic environmental conditions by applying an energy balance model with a built-in dynamic model of g_s, which is fitted on observed temperature measurements. The method derives g_s under a fluctuating light environment and detects temporal response variations within and between wheat leaves [18, 93].

Hyperspectral data to improve transpiration rate prediction models

Hyperspectral reflectance in the near-infrared (NIR, 700–100 nm) and short-wave infrared (SWIR, 1000–2500 nm) region has been related to leaf cuticle thickness, water content and anatomy [72, 94]. Cuticular transfer, intercellular space, g_s and mesophyll cell wall conductance determine total leaf conductance and therefore directly or indirectly affect E [2]. Gas exchange methods are often used to measure E and g_s, also in this study, and actually measure total leaf conductance, which is most strongly influenced by g_s, but also by the aforementioned factors. Hyperspectral data may be able to capture information on these other factors and consequently improve E and g_s prediction models. The relationship between g_s or E and hyperspectral reflectance was confirmed in this study and in Mertens et al. [42]. Both studies showed significant relationships with blue-green (523–532 nm) and the NIR water-absorption trough (976–992 nm) reflectance. The number of significant correlations was higher in this study, as additional wavelengths located in the SWIR water absorption troughs (1407, 1881 and 2346 nm) also correlated with E and g_s. The diurnal relationship between red reflectance and E was the only one described by Mertens et al. [42] that was not detected here. The reason may be the limited number of E measurements collected within one day and/or the interaction between drought and diurnal effects, which increased the changes in red reflectance and decreased E variation. Overall, several studies have used the relationships between reflectance and E to develop hyperspectral indices that can be used to monitor this trait [42, 58, 95]. The WBI index was originally created to monitor water content in plants, but Marino et al. [95] observed correlations with g_s and whole-plant E, which were also found in this study (r_WBI-E = -0.52). This correlation differs from the results observed by Mertens et al. [42], where WBI was not significantly related to E. These differing results may be caused by dissimilarities in the time between imaging and measurements, which ranged up to 1 h in the study of Mertens et al. [42] and was limited to 10 min in this study. The other two indices that had relatively high correlations with E were SR_1440/1460 (r = − 0.49) and RMP_1483/1430 (r = 0.53). These indices included (SWIR) water absorption regions, suggesting that these parts of the spectrum may contain information about E.

The combination of TIR and hyperspectral indices in empirical models has mainly been used to predict yield, chlorophyll content and relative water content [15, 37]. These studies observed a slight improvement in R² when TIR and hyperspectral indices were combined in a PLSR or stepwise multiple linear regression model. Adding hyperspectral information to predict E was only beneficial here for the single TIR indices-based models that used CWSI_Ta and CTD. In these models, the water content WBI and RMP_1483/1430 indices replaced the VPD variable, suggesting that they may contain additional information about the relationship between leaf water content, E and VPD. When hyperspectral reflectance was combined with multiple TIR indices and environmental data, a slight improvement was visible; however, this combination performed better for some of the other inbred lines, especially W153R. This improvement suggests that hyperspectral data can capture genotypic differences in water-use behavior, but more research is needed to determine the contribution of this data type to E prediction models. The most robust approach to predict E was the energy balance model (Additional file 1), which had similar, albeit low (R₂ = 0.32, RMSE = 0.68 and MAPE = 45%) prediction accuracies for all eight inbred lines. The accuracy of this approach can be improved by refining the T_dry estimate, which may make this a viable method to monitor E in automated phenotyping platforms.

Thermal indices detect genotypic differences in drought sensitivity

Genotypic differences in T_p or TIR indices have been observed in several heat tolerance or drought sensitivity studies in rice, wheat and maize [17, 24, 96,97,98]. They noted larger CWSI differences between WW and WD treatments in drought stress-sensitive genotypes compared to tolerant ones. In this study, genotypic differences in drought sensitivity were observed in both TIR indices and physiological traits. The strongest physiological differences were present in the leaf ψ measurements, in which the drought-induced decrease was significantly less pronounced in TX303 than OH43, NC358 and W153R. This partially corresponded with the available drought sensitivity data, in which plant height and leaf length of OH43 and W153R were more affected by drought compared to TX303 (data not shown). These genotypic differences diverged from the CSI_in and TRI_in results, in which TX303 was more affected by drought compared to OH43, W153R and H99. The index results seemed to correspond with the g_s and E measurements, which showed slightly larger (but not significant) treatment differences in TX303 compared to the other three genotypes. These results suggest that TX303 might close its stomata more quickly at a less negative leaf ψ, which corresponds with a water-conservative type of drought response. The limited number of significant differences observed in the genotype comparison analysis may have resulted from confounding factors. Measurements and images were grouped based on how long the plants had received the WD treatment, which started at the V5 stage. This procedure was chosen to facilitate the comparison of genotypes that differed in developmental timing and thus the date at which the drought started. This resulted, however, in combining and comparing data collected on different measurement days. The potential difference in environmental conditions may have created additional variation in the E and g_s measurements. To make this analysis more robust, a larger number of physiological E and g_s measurements are required, or the start of the WD treatment should be synchronized for all the genotypes. Consequently, more research is needed to elucidate the genotypic differences in physiology and TIR indices of these eight inbred lines.

Conclusions

Thermal infrared imaging has been a popular tool to detect drought and estimate g_s or E in the field for many years. The use of thermography in indoor automated phenotyping platforms has been less investigated and is mainly limited to comparing T_p or CWSI between drought-stressed and well-watered plants. In addition, the advantages of combining different imaging systems, such as TIR and hyperspectral imaging, to predict E has not been evaluated before in indoor phenotyping setups. In this study, the accuracy of TIR indices developed under field conditions were investigated in an indoor automated phenotyping platform. The added value of combining TIR and hyperspectral data in empirical E prediction models was determined and compared to mechanistic models. The results demonstrated that TIR indices that corrected T_p for VPD and/or T_a (CSI, TRI), were most sensitive to drought and were able to detect genotypic differences but were not strongly correlated with E and g_s. Instead, E and g_s were correlated with the commonly used indices CWSI and CTD, which could be used to develop empirical E prediction models. Model performance was the highest when TIR indices were combined with environmental data and hyperspectral indices or wavelengths. Empirical models had the highest prediction accuracies for the genotype on which they were trained (B104), but their performances were inconsistent for other genotypes, indicating that the genotypic differences should be considered during model development. In addition, it is important to monitor environmental data at multiple positions on a phenotyping platform with fixed camera positions, as transportation will influence the acclimatization of the plant to its environment and E prediction accuracy.