Three varieties of strawberry, i.e., Xiangye, Hongyan, and ZhangJi, introduced from Japan and widely cultivated in China were used in this study, which differ in drought-tolerance according to the preliminary screening by field survey. The experiments were carried out in a semi-controlled greenhouse (L×W×H: 10 m×5 m×4 m) in January and April of 2022 in Huai’an (119°01′E, 33°35′N), Jiangsu Province, China. One seedling at the four-leaf stage was transplanted per pot (upper diameter 16 cm, lower diameter 13 cm, height 18 cm, 1.5 L in total) with profile porous ceramic substrate (PPC, diameter: 0.2 mm, pH: 5.5 ± 1, porosity: 74, and CEC: 33.6 mEq/100 g). Seven replicates (three for destructive sampling) for each genotype were set in a completely randomized design on the high-throughput physiological phenotyping platform, also known as ‘Plantarray’ ( Figure 1A , Plant-Ditech, Israel, Halperin et al., 2017). Water and nutrient supplements (Yamazaki nutrient solution) were controlled by the automatic irrigation system of “Plantarray” (Xu et al., 2015).

As shown in Figure 1B , the irrigation regimes were scheduled as sufficient irrigation across the whole period (CK), water recovery following progressive drought at the critical threshold (θ_cri) when Tr starts to reduce significantly (WR_θ_cri), and water recovery following progressive drought at severe drought (WR_SD), where the wilting leaves were observed. Experiments of WR_θ_cri and WR_SD included three phases, i.e., sufficient irrigation, progressive drought and water recovery. The water supply of the WR_SD treatment was stopped 5 days earlier than WR_θ_cri (according to the preliminary experiment) to guarantee the same day of water recovery; thus, the functional physiological traits of different treatments could be compared under the same environmental conditions. At the sufficient irrigation and water recovery phases, nutrient solution was provided by irrigation for 240 seconds (oversaturated) at 0:00, 1:00, 2:00, and 3:00 hours, and no nutrient solution was supplied during the progressive drought phase.

Air temperature (T_air, °C) and relative humidity (RH, %), photosynthetically active radiation (PAR, μmolm⁻² s⁻¹) above the canopy, and soil volumetric water content (VWC, m³ m⁻³) in each pot were measured by the VP-4 sensor (Decagon Devices, Pullman, Wash, USA), 5TM sensor (Decagon Devices, Pullman, Wash, USA) and PYR solar radiation sensor (Decagon Devices, Pullman, Wash, USA), respectively. Diurnal variations in T_air, RH and PAR during the spring growth season are exhibited in Figure 2 . The vapor pressure deficit (VPD, kPa) was calculated according to Halperin et al. (2017).

The physiological parameters for plant water relations were obtained from the “Plantarray” system. This system provides a simultaneous measurement and detailed physiological response profile for plants in each pot in the array over time periods ranging from a few minutes to the entire growing season under normal, stress and recovery conditions and at any phenological stage (Li et al., 2020). The whole-plant transpiration rates (Tr, g min^-1), daily whole-plant transpiration (E, g d^-1), and daily fresh plant growth (PG, g d^-1, only under overirrigated conditions) were calculated using the SPAC software implemented in the “Plantarray” system according to the continuously monitored system weight. Briefly, Tr was computed by multiplying the first derivative of the measured system weight by -1. E was calculated by the difference in the measured system weight at predawn (averaged between 4:00 and 4:30) and that at evening (averaged between 19:00 and 19:30). During the well irrigation phase, the PG was calculated as the difference in the predawn system weights on two adjacent days. Difference between measurement data were compared through analysis of variance and that between groups with Q test.

Daily whole-plant midday transpiration (Tr_m, averaged over the period between 12:00 and 14:00) was fitted to a two-piecewise linear function of the corresponding soil volumetric water content (VWC) during the dynamic period of water deficit, which was described in our previous study (Wu et al., 2021). To offset the influence of daily environmental variation, Tr_m is normalized to VPD (Tr_m,VPD). The dynamic water response function is then determined using the following equation:

where θ_cri is the critical soil water content quantified by the VWC in this study, a is the maximum Tr_m,VPD, k is the decreasing slope of Tr_m,VPD when VWC is below θ_cri.

New fully expanded leaves of strawberry variety Xiangye of CK, WR_θ_cri and WR_SD treatments were sampled at approximately 12:00 am the day before water recovery around the end of progressive drought (12 h before water recovery, i.e., CK 12 h-, WR_θ_cri 12 h-, WR_SD 12 h-), the day of water recovery (12 h after water recovery, i.e., CK 12 h+, WR_θ_cri 12 h+, WR_SD 12 h+), and the day after water recovery (36 h after water recovery, i.e., CK 36 h+, WR_θ_cri 36 h+, WR_SD 36 h+), respectively. Three biological replicates (pots) were set for each treatment. Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA). Five micrograms of total RNA from each sample was used to construct an RNA-Seq library by using the TruSeq RNA Sample Preparation Kit according to the manufacturer’s instructions (Illumina, San Diego, CA). A total of 27 libraries were constructed and then sequenced on the Novaseq 6000 platform. The accession number for the sequence data is PRJNA891140 (https://www.ncbi.nlm.nih.gov/search/all/?term=PRJNA891140).

The raw reads were filtered and trimmed using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. The clean reads were combined and aligned to the strawberry reference genomes (https://www.rosaceae.org/Analysis/9642085). The count of the mapped reads from each sample was derived and normalized to fragments per kilobase of transcript length per million mapped reads (FPKM) for each predicted transcript using Cufflinks (v2.2.1).

Pairwise comparisons were made between samples collected at the same timepoint and different phases. To reduce noise and false results, only genes having an FPKM ≥ 1 in at least 3 of the 6 samples of a comparison (3 replicates each sample) and a coefficient of variation (CV)< 0.2 in the 3 replicates were considered. The genes exhibiting a difference of at least twofold change with the corrected P value after adjusting with false discovery rate (FDR) ≤ 0.05 were considered significantly differentially expressed.

Gene functional categories (GO enrichments) were analyzed using AgriGO V2.0 (Tian et al., 2017) under a q-value threshold of 0.01 for statistical significance. Unigenes were annotated to corresponding KEGG gene/KEGG ortholog descriptions and their associated pathways using the KEGG Automatic Annotation Server (KAAS, http://www.genome.jp/tools/kaas/) (Moriya et al., 2007). The list of gene identifiers and log2FC values from each treatment were imported into MapMan software version 3.5.1R2 and assigned to functional categories using the strawberry mapping file obtained by Mercator 4.

The midday Tr normalized to VPD (Tr_m,VPD) and soil VWC averaged between 12:00 and 14:00 were well fitted with a piecewise function ( Figure 3 , R² ranged from 0.82 to 0.91 for the dataset obtained from spring season). The transpiration variation caused by environmental fluctuations, as shown in Figure 2 , was largely smoothed by dividing by VPD, where Tr_m,VPD remained relatively stable when there was no severe water deficit ( Figure 3 ). The inflection point of the fitting curve was determined as the critical soil water content, i.e., θ_cri in Eq. 1, below which Tr_m,VPD decreased linearly. θ_cri, which represents the sensitivity of stomatal closure (Halperin et al., 2017), varied within different varieties of strawberry. The θ_cri of Hongyan was highest (0.1944 ± 0.0087, Figure 3A ), followed by Xiangye (0.1410 ± 0.0031, Figure 3B ), and Zhangji (0.1327 ± 0.0087, Figure 3C ), indicating that within the three strawberry varieties, Hongyan was most sensitive to water deficit and reduced Tr the earliest to survive under drought. In addition, the descending slope, i.e., k in Eq. 1 also differed within varieties, where Xiangye (4.6011 ± 0.4396, Figure 3B ) was significantly different from the other two (Hongyan: 2.3915 ± 0.3178, Figure 3C ; Zhangji: 2.7851 ± 0.7534, Figure 3A ), suggesting a faster process of stomatal closure. There was no significant difference between θ_cri of a particular variety obtained from the winter season and spring season, but there was little difference for k. That is, θ_criis a relatively conservative, genotype-independent parameter that characterizes the dynamic pattern of Tr in response to the soil water content.

The midday Tr (Tr_m) and daily transpiration (E) of treatments with three irrigation schedules, i.e., well irrigation (CK), water recovery at θ_cri (WR_θ_cri), and water recovery under severe drought (WR_SD), were compared, as shown in Figure 4 for Hongyan, Xiangye, and Zhangji. There was no significant difference in Tr_m or E during the period before progressive drought treatment within all three varieties ( Figure 4 ). With the experiment, WR_θ_cri was close to or slightly slower than CK (no significant difference in two-way ANOVA), where Tr_m and E were both reduced by ~6% averaged for three varieties across the progressive drought phase (blue line in shaded area of Figure 4 ). However, the Tr_m and E of WR_SD both decreased significantly compared to CK (P<0.01), by ~37% and ~35%, respectively (red line in shaded area of Figure 4 ). Within the three varieties, the Tr_m of Hongyan, which is highly sensitive to water deficit according to θ_cri, decreased the most during the progressive drought period in the WR_SD treatment (~42%, Figure 4A ), followed by Xiangye (~37%, Figure 4B ) and Zhangji (~31%, Figure 4C ), but was different for E (Zhangji: reduced by ~41%, Figure 4D ; Hongyan: reduced by ~34%, Figure 4E ; Xiangye: reduced by ~31%, Figure 4F ). After water recovery, the Tr_m and E of WR_θ_cri recovered rapidly to levels similar to those of CK, where Tr_m and E only decreased by ~2% averaged for the three varieties across the water recovery phase ( Figure 4 ). In accordance with Tr_m and E during the progressive drought phase, the recovery of Tr_m was poorest for Hongyan (reduced by ~29%, Figure 4A ), followed by Xiangye (~18%, Figure 4B ) and Zhangji (~16%, Figure 4C ), while there was the largest gap between the E values of WR_SD and CK for Zhangji (reduced by ~32%, Figure 4F ), followed by Hongyan (~17%, Figure 4D ) and Xiangye (~16%, Figure 4E ).

Daily plant growth (fresh weight, g) was also detected nondestructively by the “Plantarray” system during the water recovery phase, and water use efficiency (WUE) during the 5-day window after water recovery was also calculated by accumulated plant growth and E. The daily plant growth of WR_SD during the water recovery phase was significantly lower than WR_θ_cri and CK in all three varieties, mainly due to the daily transpiration, as there was no significant difference in WUE either between varieties or treatment ( Table 1 ).

Physiological traits before and after water recovery indicated that there is an essential difference for strawberry growth between WR_θ_cri and WR_SD. To gain insights into the gene regulatory profiles under the respective conditions, gene expressions 12 h before (end of progressive drought, 12 h-), 12 h after (12 h+), and 36 h after (36 h+) water recovery in the WR_θ_cri and WR_SD treatments were compared with CK based on RNA-Seq data. An overview of pairwise comparisons is shown in Figure 5 . The number of both upregulated and downregulated DEGs of WR_SD at 12 h before water recovery compared to CK (WR_SD vs. CK 12 h-) were greatest, while the opposite was true for WR_θ vs. CK 12 h-.

Previous studies have reported that water recovery under moderate drought will not reduce and even increase crop growth and yield, but the threshold of deficit irrigation is unclear or determined empirically according to different ratios of ETc (Williams et al., 2009; Kapur et al., 2018). Our study determined the inflection point (θ_cri) of the Tr_m,VPD-VWC curve fitted by a piecewise function as the critical soil water content under drought ( Figure 3 ) and explored the potential of θ_cri as the threshold of deficit irrigation. Midday Tr, daily E before and after water recovery, and daily plant growth and WUE after water recovery obtained from the “Plantarray” phenotyping system indicated that there was no significant difference in physiological traits between water recovery at θ_cri (WR_θ_cri) and CK, but they were significantly reduced when water recovery occurred under severe drought except WUE ( Figure 4 , Table 1 ). In modeling studies, Tr/VPD is scaled to biomass accumulation (Steduto et al., 2012; Vadez et al., 2014); therefore, θ_cri also represents the threshold when water stress significantly affects biomass assimilation. GO and KEGG analyses also supported that photosynthesis-related processes and carbohydrate pathways were downregulated in the WR_SD vs. CK set but not significant in the WR_θ_cri vs. CK set ( Figures 6 , 7 ). Our study indicated that θ_cri was a genotype-independent parameter characterizing dynamic plant water use behavior. With the rapid development of phenomics, high-throughput screening of θ_cri has been realized with low cost, thus making it possible to make irrigation schedules individually according to the θ_cri of various varieties (Li et al., 2020). As no irrigation was applied during the progressive drought phase of WR_θ_cri, this method was obviously more water-saving than the ET_c methods depending on daily reference evapotranspiration.

GO terms and KEGG pathways both indicated that mineral absorption is a vital physiological process influenced by drought, not only in the WR_θ_cri vs. CK set but also in WR_SD vs. CK ( Figures 6 , 7 ). A series of DEGs visualized in the metabolism overview also included YSL, NRAMP and VIT, which are widely known as Fe transporters. Iron plays an irreplaceable role in alleviating the stress imposed by drought (Tripathi et al., 2018; Rai et al., 2021). Many studies have demonstrated that the application of Fe nutrition to plants under drought conditions enhances tolerance, yield and yield components (Di Caterina et al., 2007; Pourgholam et al., 2013; Mirjahanmardi and Ehsanzadeh, 2016). Recently, Xu et al. (2021) found that iron homeostasis is disrupted by drought stress in sorghum roots using time-series root RNA-Seq data and that loss of a plant phytosiderophore iron transporter impacts microbial community composition, leading to significant increases in the drought-enriched lineage Actinobacteria. Mineral absorption-related genes were upregulated in WR_θ_cri vs. CK 12 h- and WR_θ_cri vs. CK 12 h+, but entire WR_SD vs. CK set, suggesting that iron homeostasis was repaired in WR_θ_cri treatment 36 h after, while not recovered at severe drought treatment even 36 h after water recovery ( Figure 7A, C ). Moreover, necroptosis and ferroptosis, which are related to programmed cell death, were significantly enriched in upregulated DEGs of the entire WR_SD vs. CK set but were not retrieved after water recovery in WR_θ_cri treatment, indicating the reversible injury of WR_θ_cri treatment ( Figure 7A, C ). This is also why Tr and plant growth were not recovered in the WR_SD treatment ( Figure 4 , Table 1 ).

Drought stress triggering the biosynthesis of secondary metabolites has been widely reported (Yadava et al., 2019). Flavonoids are low-molecular-weight polyphenolic secondary metabolic compounds and play important roles in drought tolerance (Yadava et al., 2019). Ma et al. (2014) found that flavonoid pathway gene expression and the accumulation of flavonoid compounds were closely related to drought tolerance in wheat. Yang et al. (2020) found that B. chinense adapts to drought stress through physiological changes and the regulation of flavonoids in different plant tissues. Meng et al. (2021) demonstrated a role for the CcCIPK14-CcCBL1 complex in drought stress tolerance through the regulation of flavonoid biosynthesis in pigeon pea. Li et al. (2021) researched an overly insensitive drought mutant and indicated that flavonoids improve the drought tolerance of maize seedlings by regulating the homeostasis of reactive oxygen species. Flavonoids prevent stomatal closure by reducing the H₂O₂ level in guard cells and alleviate drought-induced ROS damage. Interestingly, our results indicated that DEGs related to flavonoid biosynthesis were upregulated in the entire WR_θ vs. CK set and WR_SD vs. CK 36 h+ but downregulated in WR_SD vs. CK 12 h- and WR_SD vs. CK 12 h+. Although it is more common for genes enriched in flavonoid biosynthesis to be upregulated under drought, key enzyme genes in flavonoid biosynthesis, e.g., F3H, 4CL, IFS, and DFR, have also been observed to be downregulated under severe drought (Gu et al., 2020; Yang et al., 2020). We hypothesize that moderate drought triggers flavonoid biosynthesis to defend against the negative effects of water stress but is inhibited under severe drought until 36 h after water recovery. Shojaie et al. (2016) observed that the total flavonoid concentrations in roots and shoots both reached their highest levels under mild drought stress (Ψw = –0.5 MPa) at the end of the drought experiment, while that of severe drought (Ψw = –0.9 MPa) significantly dropped after it peaked in the early stage, which is in accordance with our results.