The experiments were conducted over three consecutive seasons, from 2017 to 2019, in two vineyards located in the Negev Highlands in Israel: the Mitzpe Ramon (MR) vineyard (30°38’48.6”N 34°47’24.5”E, 850 m asl) and the Ramat Negev (RN) vineyard at the Desert Agro-Research Center (30°58’43.4”N 34°42’31.6”E, 300 m asl). Both vineyards shared the same experimental setup, comprising 20 red wine cultivars ( Figure 1 ), grafted onto 140 RU rootstock; both vineyards were planted in 2012 in a randomized block design, as described earlier (Gashu et al., 2020). The RN vineyard experienced a slightly higher daily mean temperature (with a difference of 1.5°C) and vapor pressure deficit (VPD) compared to MR over the course of the three seasons; Supplementary Figure 1 shows the differences in VPD between locations over the three seasons. The climatic conditions between the two vineyards have been extensively described in Gashu et al. (2020).

Daily degree days (DDD) were calculated following the method described by Jones et al. (2010); the sum of DDD from véraison to harvest was computed for each cultivar to assess the DDD effect on fruit quality parameters, as follows:

Where, Tmax is the daily maximum temperature and Tmin is the daily minimum temperature.

To evaluate the effect of maximum and minimum temperatures on fruit metabolism, we introduced accumulated degree hours indices expressing heat stress (Hs) and relaxation (Relx). Hs and Relx count, for each cultivar, the duration from véraison to harvest during which the hourly air temperatures were higher than 30°C and lower than 20°C, respectively.

Where, Temp is hourly temperature.

During each season, at véraison and harvest, berries were sampled for metabolite extraction and berry indices. At véraison, each cultivar was sampled when berries reached approximately 50% color change. At harvest, berries were sampled from each cultivar approaching the 23 ± 1°Brix level. Samples were collected in the morning from four biological replicates at each location from each cultivar. In each sampling, at least 30 berries per replicate were pooled from five different vines in each block on the east side of the vine (six berries per vine were sampled from the top, middle, and bottom of the bunch), and immediately snap-frozen in liquid nitrogen. Berries were peeled while still frozen, by carefully separating the skin from the pulp, and the seeds were removed. The skin was kept at -80°C until further analysis.

Grape skin samples were lyophilized and ground under liquid nitrogen using a retsch‐mill (Retsch, Haan, Germany) with pre-chilled holders and grinding beads. For metabolite extraction, 40 mg of frozen skin powder was weighed and extracted in a 1-ml pre-cooled methanol:chloroform:water extraction solution (2.5:1:1 v/v) with ampicillin (1 mg ml⁻¹ in water) and corticosterone (1 mg ml⁻¹ in methanol) as internal standards as described by Hochberg et al. (2013) and Degu et al. (2014). Skin extracts were filtered (0.22 μm Millipore, MA, USA) and transferred to glass vials for analysis using ultra-performance liquid chromatography coupled to a quadrupole time‐of‐flight mass spectrometer (UPLC QTOF‐MS; Waters, MA, USA).

Chromatographic separation and the MS conditions were exactly as described previously (Hochberg et al., 2013).

MassLynxTM version 4.1 (Waters) was used for system control and data acquisition. Metabolites were annotated based on fragmentation patterns searched against those in the Chemspider metabolite database (http://www.chemspider.com/), and the consistency of their retention times with those of identified metabolites was compared with the data in the scientific literature.

Frozen ground skin tissue from the same sample used for metabolite extraction were used for gene expression analysis as previously described by Degu et al. (2014). RNA extraction was performed from three biological replicates of 50 mg of berry skin tissues sampled at véraison and harvest in the above mentioned five cultivars, grown at RN and MR in 2019 season. Given that the climate conditions in 2019 resemble the typical Negev desert climate, samples collected during the 2019 season were used for RNA extraction. Total RNA from grape skin tissues was extracted using Quick-RNA™ Miniprep Kit (https://www.zymoresearch.com), following the manufacturer’s instructions, with minor modifications. Briefly, RNA was extracted using buffer containing CTAB (2%), 150 mM trisodium citrate, 20 mM Na₂EDTA, 6M guanidine hydrochloride, PVP40 (2%) and 1% β-ME The extracted RNA was subsequently precipitated using isopropanol and then subjected to column purification. The cDNA was synthesized using a qScript™ cDNA Synthesis Kit (https://www.quantabio.com) using 250 ng of purified RNA following the instruction manual. The qPCR was performed with an ABI 7500 Real-Time PCR system (https://www.thermofisher.com) using Power SYBR^® Green (https://www.thermofisher.com) and data acquisition was done by 7500 software v2.3 (Applied Biosystems, USA). All reactions were performed in triplicates from three sets of biological replicates. Relative quantification was performed using the Actin gene as an internal reference gene, and RNV as a calibrator sample by the 2^-ΔΔCT method (Livak and Schmittgen, 2001).

All the statistical analyses were performed using ‘R’ version 3.6.0 (R Core Team, 2017). A three-way factorial analysis was used to assess the effects of cultivar (C), location (L), and growing season (Y), and the interactions between them, using the built-in aov function. The differences between locations for each cultivar were tested using the Wilcox.test function. Clustered heatmaps were created using Complexheatmap (Gu et al., 2016). Clustering of samples was calculated by Euclidean distances and the Ward.D2 clustering method in the functions get_dist and hclust, and the built‐in “dendextend” and “factoextra” packages (Galili, 2015). Correlation-based network analyses were conducted using the MetScape application and the NetworkAnalyzer tool, available in Cytoscape version 3.7.2. All correlation analyses and preparation for network visualizations were generated in ‘R’ using the built‐in cor function with the “Pearson” algorithm. Correlations were incorporated into the network if they were statistically significant (p value< 0.05) and their correlation coefficient was higher than 0.3 or lower than −0.3. Principal component analysis (PCA) was plotted using JMP^® version 13 (SAS Institute Inc., Cary, NC, 1989–2007). The norm of reactions that represents phenotypic changes (y-axis) due to location and seasonal climate changes (x-axis), was calculated using linear regression models for each trait computing the slope. The slope represents phenotypic plasticity.

Metabolite profiles obtained from berries, sampled at véraison during the 2017, 2018 and 2019 seasons, showed distinct patterns between cultivars and seasons ( Supplementary Figure 2 ). At véraison, HCL analysis highlighted the separation between seasons ( Supplementary Figure 2 ). Inter-seasonal differences were likely attributed to a lower anthocyanin content in 2017 samples. In contrast, 2018 samples grouped among 2017 and 2019 clusters showing an in-between pattern of change in metabolite content. Within each seasonal cluster, cultivars were grouped by location. Among the cultivars, only Sangiovese, Pinot noir and Grenache noir, which undergo véraison earlier than others, did not separate neither by season nor by location ( Supplementary Figure 2 ).

At harvest, considerable varietal differences were evident compared to véraison regardless of the vineyard location ( Figure 2 ). Through HCL analysis, we identified metabolites that exhibit a predominant cultivar effect. The lowest anthocyanin content was shown in Pinot noir, Sangiovese, and Grenache Noir ( Figure 2 ; Supplementary Table 1 ) irrespective of the season. Exceptions were cyanidin-3-glucoside (157.66 and 142.45 µg g^-1 DW at MR and RN, respectively) highest in Sangiovese, and peonidin-3-glucoside (3.64 and 4.84 µg g DW^-1 at MR and RN, respectively) in Pinot noir ( Figure 2 ; Supplementary Table 1 ). Among the stilbenes, higher content was measured in Pinot noir and Touriga nacional cultivars; where the relative content of trans piceid and cis piceid were 10-fold higher than in Cabernet Franc and Cabernet Sauvignon. Regardless of the location, Pinot noir, Malbec, and Zinfandel exhibited relatively high levels of naringenin chalcone-4-glucoside, taxifolin and astilbin. In contrast, the lowest content of these metabolites was observed in Pinotage, Tempranillo and Tinta cao cultivars ( Figure 2 ).

When values were expressed as fold-change to visualize the differences in metabolite content between locations ( Figure 3 ), three groups of metabolites were of particular interest at véraison. Cluster (1) comprised of metabolites showing minor differences between locations regardless of the seasons, e.g. phenolic acids, flavan-3-ols, stilbenes, and quercetin and quercetin-3-glucuronide, among flavonols ( Figure 3A ; Supplementary Figure 3 ); cluster (2) included mainly flavonols, flavanonols, and aromatic amino acids such as tryptophan and phenylalanine, measured at higher content in MR samples, particularly in the 2017 season while differences between the sites were cultivar dependent in the 2018 and 2019 seasons; cluster (3) was composed of anthocyanins, myricetin and its aglycons, syringetin-3-glucoside and isorhamnetin-3-glucoside, that markedly accumulated in MR during the study with the exception of season 2018 ( Figure 3A ), when high spring temperature likely triggered the biosynthesis of anthocyanins in véraison berries at the warmer RN.

At harvest, cluster (1) included metabolites displaying inter-seasonal variability between the vineyard sites, e.g. stilbenes (trans piceid, cis piceid and D-viniferin) markedly accumulated at MR in the 2017 and 2018 seasons but not in 2019; cluster (2) included metabolites displaying minor differences between locations, e.g. phenolic acids, amino acids, flavan-3-ols, a few coumaroylated anthocyanins, and quercetin and myricetin aglycones; and cluster (3) was characterized by metabolites showing predominant location effect ( Figure 3B ), mainly flavonols and anthocyanins, generally lower at warmer RN with very few exceptions among cultivars, such as Ruby Cabernet, Merlot, Malbec, and Petit Verdot.

For the purpose of visualizing the pace of metabolite pattern of change from véraison to harvest in each location, phenylpropanoid profiles obtained from berries at harvest were divided to the values at véraison (harvest/véraison), and metabolite response values were log transformed prior to HCL analysis. Three clusters were of particular interest at both locations ( Supplementary Figures 4A, B ). Cluster (1) included metabolites which decreased towards harvest, and included flavan-3-ols, phenolic acid, quercetin and its aglycons, rutin, kaempferol-3-glucuronide, taxifolin, naringenin ch-4-glu and astilbin. The extent of decrease in metabolite content from véraison to harvest was greater at warmer RN, where seasonal fluctuation ( Supplementary Figure 4 , cluster 1) were also observed. In contrast, at MR the extent of decrease remained consistent across the three seasons. Cluster (2) and cluster (3), comprising mainly of flavonols, amino acids, stilbenes, and anthocyanins, which generally increased during berry ripening at both sites.

Differences between the vineyards in respect to the relations between metabolites and climate factors (such as DDD, Hs, and Relx) were explored by correlation-based network analysis. First, we constructed networks for RN and MR samples separately ( Figures 4A, C ). The correlation threshold level for all networks was set at p-value< 0.05 and ‘r > 0.4 or r< −0.4’. The results showed that the RN network was characterized by a slightly higher number of edges (222) and nodes (56) than MR (216 and 53 edges and nodes, respectively) ( Supplementary Table 2 ).

When networks were constructed containing edges and nodes unique to each site (i.e. differential graphs), results showed a similar number of edges at both sites ( Figures 4B, D ). However, the relation between anthocyanin and flavonols was different between the locations. In the MR network, glycosylated, acylated, and p-coumaroylated anthocyanins were negatively correlated with flavonols (quercetin, quercetin-3-glucuronide, and isorhamnetin-3-glucoside) and phenolic acids such as ferulic and caftaric acids ( Figure 4B ). In contrast, in the RN network, anthocyanins were positively correlated with flavonols, with the exception of the negative correlation of quercetin-3-O-glucoside and quercetin-3-glucuronide with Malvindin‐3‐O‐(6″‐acetyl‐glucoside) and malvidin‐3‐O‐glucoside, respectively ( Figure 4D ). In both networks, only a few metabolites were linked to environmental factors. In the MR network, only myricetin-3-O-glucoside was positively correlated with DDD, while none of the metabolites were correlated with Hs and Relx. In the RN network, leucine was positively correlated with DDD, whereas astilbin was negatively correlated with Hs and accumulated radiation. It should be noted that the network analysis was performed combining the data of all cultivars, which may conceal the response of individual cultivars to environmental elements.

To compare skin polyphenol metabolites in the berries, a Pearson’s correlation matrix was created using UPLC-QTOF-MS generated metabolite profiles at harvest ( Figure 5A ), revealing a strong variation between cultivars regardless of the location and season, with the exception of cv. Tempranillo, which was unequivocally affected by the location factor. The correlation values were used as distance coefficients to build an HCL dendrogram, which described the berry metabolome in greater depth ( Figure 5B ). Still, the HCL dendrogram showed significant varietal influence, whereas the vineyard location and season had no significant impact ( Figure 5B ).

Principal component analysis was subsequently used to unravel the location and season effects specifically on Tempranillo berries ( Figure 5C ), and the result revealed a strong sample separation according to location. In addition, Tempranillo samples at each location grouped according to the seasons ( Figure 5C ). When correlation-based network analysis was performed on Tempranillo harvest data ( Figures 5D, E ), the results showed that major anthocyanins at MR were negatively correlated with cumulative minimum temperature (Relx) ( Figure 5D ). In RN network, only a few flavan-3-ols (procyanidin B1, B3), flavonols (kaempferol-3-glucoside, kaempferol-3-glucuronide), and astilbin were negatively correlated with Relx ( Figure 5E ). Interestingly, most of the flavonols and flavanonols at the warmer RN were positively correlated with DDD ( Figure 5E ). Of all cultivars, Tempranillo seemed to be the most location- and season-sensitive cultivar; and cumulative minimum temperature the most influencing environmental factor.

Next, we ranked cultivars according to their plasticity to measured location and seasonal temperature changes, by calculating the norm of reaction for all metabolite traits (Pigliucci, 2005). In the norm of reaction plot, the slope represents phenotypic plasticity (Toubiana et al., 2020); thus, the higher the slope, the greater the phenotypic plasticity of a particular genotype for a particular trait in response to the location and season ( Figure 6A ). The more moderate slope for Tinta Cao and Malbec indicates that the differences between MR and RN sites were mainly due to the interaction between season and location ( Figure 6B ). In contrast, Tempranillo displayed the steepest slope ( Figure 6C ). These results are based on the changes across all metabolites, however, a given cultivar may have high plasticity for one trait and low plasticity for another. For a few metabolites, the observed differences between locations were specific to the variety, while for others it was the interaction between cultivar, season and location the main factor explaining their variance ( Supplementary Figure 5A ). For example, among metabolites, rutin exhibited the steepest slope ( Supplementary Figure 5C ), indicating that the differences between MR and RN sites were largely attributed to a consistent climate variation between location across the seasons, while other metabolites such as malvidin 3’’6’’-caffeoyl-glucoside) and malvindin‐3‐O‐(6″‐acetyl‐glucoside) displayed more moderate slopes, ( Supplementary Figure 5B ), suggesting that differences between MR and RN sites were due to the contributions of cultivar and season.

To examine if the observed differences between sites in phenylpropanoid metabolites are related to the expression of key genes in the phenylpropanoid metabolism, we performed targeted gene expression analysis on the berry skin of five cultivars: Cabernet franc, Merlot, Malbec, Sangiovese, and Tempranillo. These cultivars were selected for their phenological and metabolic responses to environmental conditions. Tempranillo was chosen for its sensitivity to climate, Sangiovese for its consistent metabolic response to locations and seasons, and Cabernet Franc, Merlot and Malbec for their significant inter-seasonal phenological variation at the warmer RN vineyard (Gashu et al., 2020). The results showed that the effect of environmental conditions at the gene expression levels was developmental and varietal dependent. For example, comparing varieties, the transcript levels of both the upstream and downstream phenylpropanoid pathways were more strongly downregulated in Tempranillo than in other varieties at both sites, particularly during harvest ( Figure 7 ). The reduction in transcript levels in Tempranillo berries from véraison to harvest was higher at the warmer RN site than at the cooler MR site, suggesting that this could be a reason for the lower flavonol and anthocyanin accumulation at RN compared to MR. However, there was not always a direct link between the gene expression and the metabolic data for other varieties, probably due to the involvement of non-linear reaction kinetics mechanisms, enzyme activity, substrate affinity, and post-translational modifications (Zelezniak et al., 2014; Chu et al., 2019).

The genes involved in flavonol biosynthesis, such as VviF3H, VviFLS and VviLAR were expressed at similar levels between locations in véraison berries, with the exception of Tempranillo in which VviFLS gene was significantly reduced at MR. In contrast, in fully matured berries the expression of VviFLS was significantly reduced in Cabernet Franc, Merlot and Sangiovese at MR compared to RN, indicating that the differences between sites in gene expression were developmental and varietal dependent. The transcript levels of UDP-glucose:flavonoid 3-O-glucosyltransferase (VviUFGT), anthocyanin-O-methyltransferase (VviOMT) and the transcription factor MYB (VviMYB14) were generally down regulated at harvest without clear differences between sites. An exception to this was VviMYBPA1 which was strongly downregulated at cooler MR site. Our results are supported by prior works that show the UFGT gene expression increase observed in the initial 7 days after véraison under low temperature, which remains consistent in both high and low temperature treatments 21 days post véraison (Mori et al., 2007; Movahed et al., 2016).

The mechanism behind the suppression of anthocyanin accumulation in berry skins during high-temperature ripening conditions remains unclear. In certain red-pigmented plants, such as apples, the less intense red color observed in fruits ripening at higher temperatures has been linked to the downregulation of genes associated with anthocyanin biosynthesis (Lin-Wang et al., 2011). However, such an effect was not observed in most of the wine grape varieties including Cabernet Sauvignon berries ripening under control and higher-temperature conditions, suggesting that anthocyanin levels may be reduced in ripening berry skins by chemical and/or enzymatic degradation (Mori et al., 2007). Hence, the slower rate of anthocyanin accumulation in berries exposed to high temperatures after véraison is more likely attributed to the direct inhibition of enzymes within the phenylpropanoid biosynthesis pathway rather than the synthesis of the corresponding mRNAs.