Thermography as a Tool to Assess Inter-Cultivar Variability in Garlic Performance along Variations of Soil Water Availability

: Climate change entails increasingly frequent, longer, and more severe droughts, especially in some regions, such as the Mediterranean region. Under these water scarcity conditions, agricultural yields of important crops, such as garlic, are threatened. Finding better adapted cultivars to low water availability environments could help mitigate the negative agricultural and economic impacts of climate change. For this purpose, plant phenotyping protocols based on remote-sensing technologies, such as thermal imaging, can be particularly valuable since they facilitate screening and selection of germplasm in a cost-e ﬀ ective manner, covering a wide range of temporal and spatial scales. In this study, the use of a thermal index known as the crop water stress index (CWSI) was tested as a predictor of bulb biomass and for the assessment of inter-cultivar variability of ﬁve garlic cultivars in response to a gradient of soil volumetric water contents (VWCs). Three experimental assays, one in the 2018 season and two in 2019, covering a wide range of water availability levels were carried out. Di ﬀ erent linear models were developed, with CWSI and VWCs as continuous predictors of bulb biomass, and the factor cultivar as a categorical predictor. The results support the existence of inter-cultivar variation in terms of sensitivity to water availability. The most productive cultivars under favorable conditions were also the most sensitive to water availability. In contrast, the cultivars with lower bulb production potential displayed lower sensitivity to water availability and higher stability across experimental assays. The results also support that CWSI, which was sensitive to inter-cultivar variability, is a good predictor of garlic bulb biomass. Therefore, CWSI can be a valuable tool for garlic phenotyping and cultivar screening. and


Introduction
Erratic climatic conditions brought about by climate change entail a great instability in crop yields throughout crop seasons [1][2][3], which causes economic losses [4]. These economic and yield damages could be buffered with more resilient agroecosystems based on crop biodiversity [5,6]. In this sense, finding sustainable, resilient, and more yield-stable cultivars under different environmental conditions becomes crucial, especially under climate change scenarios [7][8][9][10][11].
Crop phenotyping is currently in the search for improved genotypes/cultivars against abiotic and biotic stresses [12,13]. Thus, intraspecific characterization of functional and yield responses to the environment should be a key approach in crop breeding and selection programs [14]. For this purpose, crop phenotyping and cultivar selection must be performed in a cost-effective manner with practical implementation [15][16][17]. In this sense, the use of new technologies, along with enhanced

Study Site and Weather Conditions
Three field assays were carried out during two consecutive seasons, 2017-2018 and 2018-2019 (hereafter referred as "2018" and "2019", respectively). During 2018, we conducted one of the assays (Alb18) at the "Centro de Investigación Agroforestal de Albaladejito" (40 • 04 31.436N, 2 • 12 18.061W, altitude 902 m.a.s.l) in Cuenca, inland Spain. The climate of this region is classified as a Mediterranean-continental climate, characterized by cold winters (mean minimum temperature) and hot and dry summers, with strong thermal oscillations with mean (minimum-maximum) temperatures of −1.5 • C in January to 30.1 • C in July [52]. Historical average annual precipitation and temperature are ca. 500 mm and 11.8 • C, respectively [52]. The soil of Albaladejito is a sandy silt loam. For additional details of its physicochemical characteristics, see [53].
In 2019, two additional assays were placed in two locations. The first one was placed at Albaladejito (Alb19) in an experimental plot close to that of Alb18, and the other one (Isl19) was placed at IMIDRA's experimental station, "La Isla" (40 • 18.75 N; 3 • 29.89 W; 528 m.a.s.l), in Arganda del Rey, Community of Madrid, 140 km away from Albaladejito. "La Isla" is located in the central region of the Iberian Peninsula, within the domain of the Mediterranean-temperate climate, with a Mediterranean-continental tendency. Historical average annual rainfall is between 440 and 490 mm and the average annual temperature is around 13.8 • C [52]. The soil in "La Isla" is a sandy silt loam with the following chemical characteristics: pH of 8.5; 1.3% organic matter; and 0.3% nitrogen, 42 mg phosphorus (P)/kg, and 434 mg potassium (K)/kg. Soil analyses were performed by "Laboratorio de Suelos IMIDRA" (Instituto Madrileño de Investigación y Desarrollo Rural, Agrario y Alimentario, Madrid, Spain). Soil analyses were carried out before the establishment of the assay, but determinations from previous years show very constant inter-annual values.
Meteorological data at Albaladejito were recorded with a micro-meteorological station, Vantage pro 2 model 6152EU, Davis Instruments, Diablo, US, in 2018 and 2019. At "La Isla" these data were collected from a weather station placed in the experimental farm named "Arganda", which belongs to the national network of the agroclimatic information system for irrigation, managed by the MAPA (http://eportal.mapa.gob.es/websiar/Inicio.aspx). Mean temperatures (Tmean), sun radiation, wind speed, and vapor pressure deficit (VPD) varied across seasons and locations. Relatively mild weather conditions were recorded in Alb18 while more severe climatic conditions were recorded in the assays of 2019, especially in Isl19 ( Figure 1).

Plant Material, Experimental Design, and Watering Treatment
The five studied cultivars were selected according to their different genetic background and different origins with contrasted climatic conditions ( Table 1). The studied cultivars were "Gardacho", GAR (introduced "American type" commercial cultivar); "Purple of "Pedroñeras", PED (popular and appreciated cultivar from Castile-La Mancha recognized by the European Union as a Protected Geographical Indication, PGI); "Fino de Chinchón", CHI (traditional cultivar from the south of Madrid); Cbt00089, DRYII (traditional cultivar from Tenerife Island); and Port07990, RAIN (traditional cultivar from the north-west of Portugal). The bulbs for the experiment were obtained from the Bank of Plant Germplasm of IMIDRA (FAO code ESP198), the garlic cooperatives ("San Isidro el Santo" and "Coopaman S.C.L", "Las Pedroñeras", Cuenca, Spain), the Agricultural Biodiversity Conservation Centre of Tenerife (CCBAT), and the Portuguese Bank of Plant Germplasm.  The five garlic cultivars were planted on field plots across a gradient of soil water availability. In Alb18, 4 cloves of each of the five cultivars were planted within 11 plots. In Alb19 and Isl19, three cloves x cultivar were planted within nine plots at each location. Before planting, the cloves were treated with 1.5% sodium hypochlorite and the fungicide Prelude (prochloraz 20%, BASF SE, Ludwigshafen am Rhein, Germany). The biggest and healthiest cloves of each cultivar were selected and manually planted within each plot. The distance between plants within a line was 0.15 m and the distance between lines was 0.40 m. The cloves were randomly placed within each plot and watered after planting to help plant establishment.
The plots were watered when needed to maintain the soil volumetric water content (VWCs) within the pre-stablished target range (Table 2). In 2019, to reach lower VWCs, the plots were covered with rainfall cover structures during rain events to avoid uncontrolled water supply. VWCs was measured at hourly intervals with FDR probes (ECH2O EC-5, Meter group, Inc., Pullman, WA, USA) connected to dataloggers (Em50, Meter group, Inc., Pullman, WA, USA) on each plot during the whole crop cycle.

Crop Water Stress Index (CWSI), Bulb Biomass, and Bulb Diameter Measurement
Ground-based thermal images (Fluke Ti300 camera, Fluke Corporation, USA, 7.5-14 µm, 240 × 180 pixels, ε = 0.96) were taken to remotely assess leaves' temperature. The procedure described by Jones [54] was followed to include reference surfaces. Leaf references were attached to a cork board, where the 'wet' reference (Twet) was obtained by spraying fresh water on two detached garlic leaves while the 'dry' reference (Tdry) was obtained by coating detached garlic leaves with petroleum jelly. Zenithal thermal images, at a 1-m distance from the canopy, were taken so that they covered the whole plant and surface references. Thermal images were taken for every plant within each plot at least two (Alb18) or three times (Alb19 and Isl19). All the measurements were carried out throughout the bulbing phase [55]. The measurement interval was based on the growth development of each assay. The first infrared images were taken when most of the plants reached ca. 7-8 leaves while the second and third infrared images were taken when the plants achieved the peak of growth (ca. 10-14 leaves). Finally, the software, SmartView (Fluke Smartview 4.3, Fluke Corporation, Everett, WA, USA), was used for the analysis of the thermal images. To reduce the sensitivity of the analysis to variation in the leaf angles, average leaf temperatures were calculated by selecting representative areas of the youngest fully expanded leaves (Tleaf) and references using the function polygon markers of the software SmartView ( Figure 2). The mean ± S.D of these Tdry and Twet references for the whole experiment were 34.8 ± 4.4 • C and 20.1 ± 3.5 • C respectively. Finally, processed data were used to calculate the crop water stress index (CWSI) as: Bulbs were harvested at the end of the growth cycle, when the three youngest leaves were dry. Afterwards, bulbs were stored in a ventilated warehouse and weighed until constant weight to the nearest mg. At the same time, bulb diameter was measured with a caliper.

Statistical Analysis
At the end of each growing season, the VWCs data recorded on the dataloggers were collected and the integrated VWCs over the whole period for each plot were calculated using a macro ('area below curves') implemented in SigmaPlot 11.0 (Systat Software Inc., San Jose, CA, USA). Then, the average VWCs per day was calculated as the ratio between the area below the curve of VWCs and the number of days measured on each plot. The representative CWSI value used in the models was calculated as the average of CWSI measurements for each cultivar within the same plot. The representative bulb biomass used in the models was calculated, similarly, as the average bulb biomass of cultivar replicates on each plot. Regression linear models were used to test for the fixed effects of the continuous predictors VWCs per day and CWSI, and the categorical factor cultivar on bulb biomass. In addition to the main effects of these variables, the interaction term between cultivar and VWCs and/or CWSI were also included. Shapiro-Wilk's, Breusch-Pagan, and variance inflation factor tests were used to test for normality, homogeneity of variances, and multicollinearity, respectively. Bulb biomass was log-transformed to meet the assumptions of normality and homoscedasticity. Additional mixed linear regression models were fitted to test for the effects of the random factor (assay). ANOVA analyses were also performed to test for the significance of the regression coefficients and to test for differences between assays. Based on the predictions of the best model, estimated marginal means of bulb biomass were plotted across selected levels of the continuous predictors and Tukey test comparisons of the slopes among cultivars were computed to identify homogeneous groups. All the models and subsequent analyses were performed using R (version 1.2.5019; RStudio Team, 2019).

Models' Verification and Variable Reliability as Bulb Biomass Predictors
Alternative models differing in complexity and the number of factors were evaluated (see Table 3). The simplest models were those with only one of the continuous predictors (Cp) included, VWCs or CWSI. Both variables have been previously reported as predictors of yield and biomass gain in other crops [56][57][58][59]. The results of this study confirm that bulb biomass of garlic was positively correlated with VWCs while negatively correlated with CWSI (Supplementary Materials Table S1), with both Cp being significant ( Table 4). The lower the soil water availability, the lower the bulb biomass, as displayed in previous studies in garlic [37,53,60]. Table 3. Description and variables included in each of the tested models.
Model 4, which included both continuous predictors and the interaction of cultivar with either of them, improved the AIC and R2adj, resulting in the best models tested. Additional linear mixed models including the random factor assay (model 5 in Table 3) did not improve these models (see model 5 in Table 4). Despite CWSI usually being highly correlated with water availability variables, such as volumetric water content and evapotranspiration coefficients [30,65,66], in this study, none of the models that included both predictors showed multicollinearity (variance inflation factor test). Besides, the increased goodness of fit and explained variance of the models that included both predictors in comparison with those models that included only one of the predictors suggest that CWSI accounts for a significant portion of variability that VWCs cannot explain. Moreover, the Model4Cult*CWSI, which also included the interaction of cultivar with CWSI, resulted in the best fit. CWSI is tightly related with leaf/stem water potential [67][68][69], which has been found to vary among cultivars grown under the same soil moisture conditions [70][71][72]. Besides, CWSI includes the effect of other weather conditions besides plant water status [68]. In fact, other climatic and biotic factors have been proven to induce changes on CWSI values, such as VPD, radiation, wind speed, etc. [67,[73][74][75][76].
In other horticultural crops, CWSI has already been used for irrigation scheduling [40,41] and as a yield predictor [58,76,77]. However, the determination of irrigation requirements is a complex issue that needs further research [66]. Overall, the great deal of interest shown in the last decade on thermal imaging is mainly related with its rapidity of determination, non-invasive nature, and robustness. These characteristics makes thermal imaging a promising technique for modern crop phenotyping [16]. Under this context and based on the findings of this study, thermal indicators, such as CWSI, can be a valuable complementary tool for garlic crop management (e.g., soil moisture monitoring, irrigation scheduling, cultivar selection, and yield prediction).

Inter-Cultivar Variability Analysis on the Sensitivity of Bulb Production to VWCs and CWSI Gradients
Bulb biomass production was cultivar dependent. In particular, the bulb biomass of "Gardacho" and "Pedroñeras" was consistently higher than that of the local cultivars Cbt0089 and Port07990 along the gradients of both VWCs and CWSI ( Figure 3A,B), in agreement with the findings of Sánchez-Virosta and Sánchez-Gómez [78].
Interestingly, this inter-cultivar variability was more pronounced in plots with higher VWCs or lower CWSI, while inter-cultivar differences were less pronounced under low water availability levels. The pattern of reduced inter-cultivar variability under more limiting conditions has also been observed in other crops, such as cotton [79] or field pea [80]. In fact, it has been stated that the interaction between yield and the genetic variance in grain crops tends to decrease as the stress intensity increases (see [81] and references therein). This agrees with the results observed for garlic in this field study and in the pot-based experiments of Sánchez-Virosta and Sánchez-Gómez [78].
Accordingly, the significant interaction terms of the models that resulted in the best fits (model 4) confirmed the existence of this inter-cultivar variability in the response to VWCs and CWSI (Table 4, Figure 4; coefficient estimates and further detailed interpretation are provided in Table S1 and Box S1). Previous studies also showed this inter-cultivar sensitivity variation with CWSI in other herbaceous crops, such as sunflower [82], lentil [28], or maize [83]. In this study, the more water-limited the environments, the less inter-cultivar differences on bulb biomass-estimated marginal means (B-EMMs) among cultivars. For example, significant differences on B-EMMs were not found among cultivars under VWCs below 15%. In contrast, significant differences were found in VWCs above 15% or CWSI below ca. 0.60 ( Figure 4). The observed inter-cultivar differences on bulb biomass under higher VWCs levels were likely associated with underlying physiological mechanisms related to water consumption and growth potential (see Sánchez-Virosta and Sánchez-Gómez [78]). These findings agree with other studies in wheat [84,85], which found that a constitutive higher yield potential can involve a decrease in absolute yield stability, and inversely, lower growth potential can be associated with higher yield stability. Despite not all cultivar-pairs comparisons were significant (see table below Figure 4), as the models reflected variability in the sensitivity to water availability. Local cultivars Port07990 and especially Cbt00089 showed the lowest slopes, in agreement with the findings of Sánchez-Virosta and Sánchez-Gómez [78]. This suggests that these local cultivars might have physiological mechanisms that provide homeostasis along a wide range of environmental conditions, yet lower growth and yield potential. Such functional performance could be beneficial under more severe drought conditions [86]. In contrast, "Gardacho" and "Pedroñeras" showed steeper slopes and higher yield potential (Figure 4). This is in agreement with the results of Sánchez-Virosta and Sánchez-Gómez [78], where greater yield potential was coupled with increased capacity for water and resource uptake. However, this increased capacity for resource acquisition is sometimes associated with an increased frequency of stress experiences [86] and lower stress resilience [10], which are important aspects that should be taken into account under climate change scenarios.

Climatic Conditions and Bulb Biomass Production Cross Experimental Assays
Other climatic conditions besides water availability varied across the different assays and were especially contrasted between years ( Figure 1). Here, 2019 was one of the years of the last decades with higher temperatures in the Mediterranean region, which negatively affected summer crop yields [87]. Accordingly, in this study, the most limiting conditions occurred in 2019. In contrast, wetter than usual conditions in the study area in 2018 were also reflected in higher than average yields [88]. All this environmental variability could also influence the final bulb biomass and diameter.
Bulb biomass showed great variability and dispersion within the same levels of VWCs, especially at Alb18, where the most favorable conditions occurred. However, this could be explained, as previously mentioned, by the highest inter-cultivar differences of garlic performance under more favorable conditions. Nonetheless, other environment factors different from soil water availability and inter-cultivar genetic variability can also affect garlic bulb biomass gain and growth [49,89,90]. For example, above-optimal temperatures (>25 • C), as those occurred in 2019, have been proved to impair photosynthesis and bulb biomass in garlic [91]. In fact, when different assays were compared at a similar soil moisture, bulb biomass differences between assays were still noticeable ( Figure 5). Significant differences were found between Alb18 and Alb19 within the 20% > VWCs > 17% target (F (1/14) = 17.3, p = 0.000) and between Alb19 with Isl19 within the 17% > VWCs > 15% target (F (1/14) = 9.6, p = 0.008). These results indicate that other uncontrolled environmental factors different from VWCs and CWSI that differed among assays also had an impact on bulb biomass gain. Despite, the potential effect of uncontrolled factors, the models explained a large proportion of the variation of the data (68%), identifying VWCs, CWSI, and cultivar as the main determinants of bulb biomass in garlic. This study serves as a basis for future improved models that account for inter-cultivar variation in the response patterns of biomass gain to environmental factors. It also contributes to a better understanding of yield stability of the crop under different environmental conditions, and to evaluate its adaptive potential to more water-limiting scenarios [2,[92][93][94].

Conclusions
Garlic bulb biomass was negatively affected by water availability. The variables VWCs, CWSI, and cultivar were the main determinants of bulb biomass. The inclusion of CWSI significantly improved the models and the variance explained. Thus, we consider that monitoring of the canopy temperature of the crop can be a valuable complement to soil water content monitoring for irrigation and crop management in garlic. In addition, the existence of inter-cultivar variability of garlic sensitivity to water availability was confirmed in this study. In general, the more water limited the environments, the smaller the inter-cultivar differences on bulb biomass production. In contrast, under favorable environments, greater differences among cultivars were found. However, the most productive cultivars under non-limiting environments were also the ones that displayed the steepest model slopes and, henceforth, the highest sensitivity to water availability, which should be considered in future selection programs for climate change mitigation. Overall, this study suggests that the use of thermal imaging methodologies, and in particular, specific thermal indices, such as CWSI, can become valuable tools for practical applications on garlic farming and cultivar selection.