Climate change in a typical Apulian region for table grape production: spatialisation of bioclimatic indices, classification and Future Scenarios

The present research focused on the characterisation of climate evolution in a typical Apulian region for table grape production under the protected geographical indication, “Uva di Puglia I.G.P.”Two thirty-year time window period (TW) were analysed: 1961-1990 and 1991-2020. Georeferenced maps for both TWs were produced to delimit homogeneous zones and to evaluate the climate variability within the investigated area by means of the two bioclimatic indices, Heliothermal Index (HI) and Winkler Index (WI). Spatial analysis of HI and WI was performed using the regression-kriging (RK) interpolation method and the Digital Elevation Model/DEM (10 x 10 m) as a prediction attribute.An increase in both the minimum and maximum temperatures was observed, and locations above 300 m a.s.l. shifted from HI+1 “temperate warm” to HI+2 “warm” according to the Geoviticulture Multicriteria Climatic Classification System. WI values similarly increased between the periods 1961–1990 and 1991–2020, shifting all the sites grouped in the Elevation Classes defined as being below 300 m a.s.l. from Region IV to Region V of the Winkler Classification.According to HI and WI, presumed maturity was calculated as being reached 9 to 15 (HI) and 12 to 28 days (WI) earlier in 1991–2020 than in 1961–1990, taking into account the heat requirements of cv. Italia table grape (representative of Apulian table grape production), were set at 2200 for both indices on the basis of literature data.Moreover, three table grape vineyards, located in the three main producing provinces of Apulia (Bari, Taranto and Barletta-Andria-Trani (BAT)), were considered for future scenarios analysis on the basis of two different Representative Concentration Pathways (RCPs), 4.5 and 8.5, and classified according to the Geoviticulture Multicriteria Climatic Classification System (MCC). Future scenarios scored WI values that exceeded the threshold of 2700 in the BAT and TA provinces in the 2061–2090 time window period for RCP 8.5. In contrast, RCP 4.5 led to a mitigating effect, which was not noticeable until 2040, with a consequent reclassification of the investigated areas on the basis of HI and Cool Night Index (CI).These findings suggest that in order to prevent or overcome heat stress, it will be necessary to implement strategies, such as vineyard relocation to unexplored elevations or latitudes and/or the exploitation of new table grape varieties able to fulfill the optimal maturity parameters, even when the duration of the phenological phases is shorter.


INTRODUCTION
The existence of climate change is now incontrovertible, having clear consequences for agriculture and viticulture (Koufos et al., 2018). Bioclimatic indices are an effective tool for the analysis and precise evaluation of climatic conditions in a given grape growing area. Among them, the Heliotermal Index (HI) (Huglin, 1978) and the Winkler Index (WI) (Amerine and Winkler, 1944) (Table 1) are frequently applied in the viticultural sector, as they correlate well with quantitative and qualitative parameters, like yield, chemical composition, sugar and acidity content and accumulation of secondary metabolites (Tonietto and Carbonneau, 2004). Other bioclimatic indices which rely on temperature have been employed in geoviticultural climatic characterisation; for example, Diurnal Thermal Range (DTR) (Koufos et al., 2013), Cool Night Index (CI) (Tonietto, 1999), Biologically Effective degree-day index (BEDD) (Gladstones, 1992) and Growing Season Temperature (GST) (Jones, 2006).
The Dryness Index (DI) (Riou et al., 1994) differs from these indices as it is based on rainfall and evapotranspiration with the aim of classifying a geographical area in relation to its hydric components.
In climate research, a sufficiently large and representative dataset needs to be consideredcovering at least 30 years -and preferably referred to the specific historic period of , as defined by the World Meteorological Organization (WMO). Recently, the European Centre for Medium-large Weather Forecasts (ECMWF) have started to consider the ten year time window period of 1991-2020 as the main reference period, as previously unavailable satellite data can now be used (C3S Climate Bulletin, 2021).
In macroclimate, mesoclimate and local climate research, data are normally spatially represented and visualised by means of GIS software (Neethling et al., 2019). However, climatic stations that are far apart -which is particularly the case in macroclimate analysis -can result in uncovered  Amerine and Winkler, 1944;b from Tonietto and Carbonneau, 2004. spaces which have to be filled by spatial data products; these are obtained by the interpolation of existing long-term, quality-controlled data sources using numerous techniques . Universal kriging (UK), kriging with external drift (KED) and regression-kriging (RK) are the most suitable and effective interpolation methods for reducing standard error in the spatialisation of bioclimatic indices in the following cases: when elevation is considered as a single independent variable (Moral et al., 2014;Honorio et al., 2018), in comparative studies (Hengl et al., 2003), and compared to other interpolation methods like localised regression models employed to produce climatic maps (Attorre et al., 2007;Blasi et al., 2007). These interpolations, defined as "hybrid" (McBratney et al., 2000), involve non-static geostatistical methods which use auxiliar predictors with spatially correlated data in order to perform univariate or multivariate regression models, and thus effectively representing a regression analysis method (Wackernagel, 2003 . RCPs are trajectories of future greenhouse gas concentrations based on potential radiative forcing. In RCP 4.5, gas emissions, partially curbed by the application of green technologies, will continue to increase until 2040, then they will decline and stabilise by 2100 at 4.5 W/m 2 in comparison to the year 1850 (Thomson et al., 2011). In contrast, in the RCP 8.5 pathway, no gas emission reduction policies for climate mitigation are applied, causing a radiative forcing at 8.5 W/m 2 in comparison to the year 1850 by 2100 (Riahi et al., 2011).
In recent years, much research has dealt with evaluating climate change and its impact on viticulture, mainly with regard to wine grapes Li and Zhou, 2014;Sgubin et al., 2019;Koufos et al., 2020;Piña-Rey et al., 2020;Santos et al., 2020).
Little attention has been given to table grape, which represents one of the most important agricultural sectors in the world. However, the climatic characterisation of Apulia has been undertaken from an economic (Sardaro et al., 2021) and environmental (Polemio and Lonigro, 2014) perspective, and the effect of future climate on the reduction of both grape yield and wine production in the first half of the 21st century in Apulia has been investigated (Lionello et al., 2014). Worldwide Table grape production stood at 27 million tons in 2018, with Italy in seventh place with a production of 1.1 million tons, around 4 % of worldwide production, (OIV, 2019). Moreover, in 2020, ISTAT reported that the Apulia region produced 60 % of total Italian table grape, corresponding to 0.63 million tons mostly concentrated in Bari, Taranto and the BAT provinces. Apulian table grape production is therefore protected by the geographical indication, "Uva di Puglia I.G.P." (UE Reg. n.680/2012, 2012, which delimits the production areas at elevations below 330 m a.s.l. and is restricted to Italia, Regina, Victoria, Michele Palieri and Red Globe varieties. The aim of this study was to characterise and describe the climate evolution of the Apulia region, the most important table grape producing area in Italy, by means of bioclimatic indices based on a historic meteorological series and two time window periods (TW): 1960-1990 (TW 1 ) and 1991-2020 (TW 2 ). A spatial analysis of bioclimatic indices was performed to better represent climate changes over the last sixty years. Moreover, three table grape producing vineyards, located in the three main producing provinces of Apulia (Bari, Taranto and BAT), underwent a future scenarios analysis based on RCPs: 4.5 and 8.5, and they were classified according to the Geoviticulture Multicriteria Climatic Classification System (MCC).

Climate data
A historic series of the monthly Tmax and Tmin on 30 Apulian weather stations was downloaded from the Civil Protection site at https://protezionecivile.puglia.it for two thirty-year time window periods (TW): 1961-1990 (TW 1 ) and 1991-2020 (TW 2 ).
In order to evaluate changes in the Heliothermal Index (HI) and Winkler Index (WI) across both TWs (Tables 4 and 5), the 30 sites were grouped into five Elevation Classes (EC) depending on their 100 m range within an elevation of 0 to 500 m a.s.l. Differences in monthly and total HI and WI were evaluated using the non-parametric Kruskal-Wallis test. Moreover, linear regression equations between cumulative HI or WI with Julian Days (JD) were calculated, with the aim of evaluating the number of JD ± SE needed by a representative "Uva di Puglia I.G.P." table grape cultivar, like Italia, to reach 2200 HI or WI for each of the five groups within each of the two TWs. The value 2200 was chosen based on the requirements reported by Fregoni (2005) for late or very late table grape cultivars.

Bioclimatic Indices
The Bioclimatic Indices and relative classes of viticultural climate are listed in Table 1.
The Winkler Index was calculated as follows: ETP (Penman, 1948) was calculated by means of CROPWAT 8.0 (FAO, 1992), which provides an Estimated ETP on the base of Tmax and Tmin when other climatic data are missing.
In the occurrence of outliers, medians were taken into account in order to obtain monthly PR (ISPRA, 2013).

Geostatistical interpolation techniques
By applying the geostatistical technique, regression-kriging (RK), it was possible to combine the linear or multiple regression analysis of auxiliary information and subsequently the kriging algorithm for residuals. The independent or predictor auxiliary variable is the digital elevation model (DEM), with a resolution of 10 m in the UTM WGS 84 zone 33N projection system (Tarquini and Nannipieri, 2017). The elevation predictor produces maps with detailed results and high accuracy of the spatial prediction, since its values, which are extracted from a raster format image, are regularly positioned (grid), are available in every point of the regionalised variable and possess a high resolution 10 x 10 m (Hengl et al., 2003;Hengl, 2007).
The validation of the predicted bioclimatic indices was processed by Leave-One-Out Cross Validation (LOOCV) in both TW 1 and TW 2 for each weather station, providing an independent comparison of the calculated and interpolated HI and WI indices. The Mean Error (ME) and Root Mean Square Error (RMSE) of the comparison were determined to evaluate the interpolation error, and therefore the reliability/accuracy of the kriging model.
The open-source SAGA GIS software (System for Automated Geoscientific Analyses) version 2.3.2 and QGIS version 3.10 were used to perform the interpolation procedure by means of kriging and the cartographic rendering in GRID format with a resolution of 10 x 10 m respectively.

Future scenarios
The DEAR-Clima database (http://meteo3.geo. auth.gr:3838/) was employed to retrieve Tmax, Tmin and PR of three representative table grape vineyards: Farm BAT area, Farm BA area and Farm TA area. Geographical and elevation data are given in Figure 1.
The time series were based on high horizontal resolution simulations produced by the Regional Climate Model (RCM) from the Coordinated Regional Downscaling Experiment (CORDEX) research program. RCM data processed in the DEAR-Clima database have a high spatial resolution ( HI and WI were calculated for the historic and RCP scenarios for each of the three vineyards, together with linear regression equations between cumulative HI or WI with Julian Days (JD) and Multicriteria Climatic Classification (MCC).
Trends in Tmin and Tmax in the historic and RCP scenarios for the three vineyards were graphically plotted using weighted moving averages within a 5-year period.
The predicted JD ± SE required to reach 2200 HI or WI on the three vineyards and in the three different scenarios were also calculated.

RESULTS
Thirty Apulian sites ( Figure 1) were studied in two thirty-year time window periods (TWs): 1961-1990 (TW 1 ) and 1991-2020 (TW 2 ). The Heliotermal Index (HI) and the Winkler Index (WI) were calculated for each of the two TWs in the Apulian area defined by the Protected Geographical Indication "Uva di Puglia I.G.P.". It comprises all the provinces of Apulia except Lecce province, and limits table grape production to elevations below 330 m a.s.l., which in practice excludes almost the whole province of Foggia. Therefore, the research focused on the central part of Apulia dedicated to table grape production in an area that includes the BAT, Bari, Taranto and Brindisi provinces.  (Figure 1), was highly correlated with the bioclimatic indices. As can be seen in Table 2, HI and WI hare highly correlated; FIGURE 2. Spatialisation of Heliotermal index (HI) and of Winkler index (WI) in a typical Apulian territory for table grape production in two time window periods: 1961-1990 (TW 1 ) and 1991-2020 (TW 2 ).
this is because they are based on temperature (Moral et al., 2014;Honorio et al., 2018). Table 2 also shows that elevation has a significant negative linear correlation with HI and WI in both Tws, in particular WI TW2 (r = -0.79); meanwhile, no significant correlation was found between longitude/latitude and the two bioclimatic indices. The spread of the interpolation error within the study area, and therefore the reliability of the model, was confirmed by the Leave-One-Out Cross Validation (LOOCV) for both TW 1 and TW 2 , for each weather station ( Table 3). The RK yielded the lowest root-mean-square error (RMSE) with respect to other interpolation methods (data not shown), especially for HI in both TWs, indicating an effective spatial variability representation. The two indices calculated from interpolation (predicted data) based on the 30 Apulian weather stations, were very close to the observed one in both TWs, which was confirmed by the low Mean Error (ME) value and the high Pearson correlation coefficient ( Figure 3). Table 4 and Table 5 report the comparison of HI and WI between TW 1 and TW 2 relative to the five ECs. The number of sites included in each EC ranged from a minimum of 4 in EC V to a maximum of 9 sites in EC I. A total of 21 out of the 30 sites were located between 0 and 300 m. a.s.l., an altitude suitable for table grape cultivation. Given the small number of sites falling into each of the EC, the non-parametric Kruskall-Wallis test was preferred to a parametric ANOVA.
Differences in HI between the TWs were found in EC I (0-100), EC III (201-300) and EC IV (300-400), while no significant differences were found in the remaining two ECs, and HI even showed an increase of 196     In order to analyse the historic  and future (2021-2050 and 2061-2090) climate scenarios (Representative Concentration Pathway -RCP 4.5 and RCP 8.5), monthly Tmin, Tmax and Precipitation (PR) of three table grape vineyards located in three different Apulian provinces (Farm BAT, Farm BA and Farm TA; Figure 1) were retrieved from http://meteo3.geo.auth.gr:3838/. WI, HI, Cool Night Index (CI) and Dryness Index (DI) were calculated for all the scenarios ( Table 6).
The differences in terms of Tmax and Tmin between the three vineyards for all the scenarios are shown in Figure 4. In the historic scenario, the recorded temperatures on Farm BA were lower than the other two vineyards by around 1 °C for both Tmax and Tmin and this affected the bioclimatic indices. MCC in both RCP 4.5 and RCP 8.5 shifted Farm BA to CI-2 (warm nights); this is the same class as the other two vineyards, for which no classification change occurred, excepted in terms of HI, which shifted to HI+3 (very warm). Table 6 shows linear equation regressions (R 2 > 0.97) between HI or WI and JD for all the scenarios and time series for all three vineyards. Based on these equations, once HI or WI were equal to 2200, a box-plot of the JD ± SE needed to reach the presumed harvest date for the Italia variety was calculated ( Figure 5). An early presumed maturation was found for all the vineyards, gradually becoming more intense over the years according to the different scenarios. In the historic scenario, 273 JD were necessary to reach 2200 HI and thus the presumed maturity in Farm BA, which corresponded to 30 September; this is in line with the actual period of maturation for Italia cv., a late table grape that is usually harvested in the third week of September or, at the latest, the first week of October.  Instead, Farm TA and BAT, being located in warmer Apulian areas and at lower elevations, needed 256 and 259 JD respectively for maturation, corresponding to 13 and 16 September, two weeks earlier than Farm BA. Similarly, Farm BA achieved 2200 WI in 317 JD (13 November), which is later than expected, while 295 JD (22 October) were needed by the other two vineyards.
A shift to earlier presumed maturity was evident in the 2021-2050 period for all three vineyards in terms of HI, ranging from 10 JD for Farm BAT and 11 JD for Farm TA to 12 JD for Farm BA in RCP 4.5, with differences when compared to RCP 8.5 of 2 to 4 days. WI values showed the same trends as HI, with maturity being 20 JD, 24 JD and 19 JD earlier in RCP 4.5 for Farms BAT, BA and TA respectively, with differences when compared to RCP 8.5 of 3-5 days. As previously stated, differences between RCP 4.5 and 8.5 were found to be more substantial in the following TW, 2061-2090. According to HI, by comparing the historic scenario with RCP 4.5 and RCP 8.5, both Farm BAT and Farm TA showed an early presumed maturation of 19 JD and 28 JD respectively.  Farm BA was calculated to be even earlier, with a 22 and 33 JD difference when the historic scenario was compared to RCP 4.5 and 8.5 respectively.
A shift to 2200 WI was likewise reached in 2061-2090, with values ranging from 28 JD in Farms BAT and TA to 34 JD in RCP 4.5. Large differences were found in the historic scenario compared to RCP 8.5, with presumed maturity occurring 39 JD early in Farm BAT and TA, and, even more remarkably, 48 JD early in Farm BA.
However, the box-plot in Figure 5 shows that overall, when HI ± SE was calculated from the linear regressions, the ranges of historic and RCPs relating to the time window period 2021-2050 and RPC 4.5 of the 2061-2090 substantially overlap. On the other hand, a lower JD value was evident in RCP 8.5 in the period 2061-2090, in which the temporal range of the estimated JD value ± SE did not overlap with the previous scenarios, suggesting a real worsening of climate conditions which will affect grape phenology in the future. The box-plot relative to the WI values is in some ways even more emblematic in 2061-2090 for both RCP 8.5 and RCP 4.5, indicating that the countermeasures envisaged by RCP 4.5 may not even be sufficient to reduce the effects of climate change on the phenology of grapevine now and in the coming decades.

DISCUSSION
The present research aimed to characterise the changes to climate over the last decades in the Apulia region by means of bioclimatic indices spatialisation. The comparison of two thirty-year historic series allowed us to understand how climate changed in Apulia between 1961 and 2020 and to predict what will happen in the future from the recorded temperatures, which mainly result in a shift to earlier phenological events in grapevine (Koufos et al., 2013).
Our results are consistent with Cotecchia (2014), who reported an average annual temperature in Apulia of 15-16 °C and higher average values in the Ionian-Salento area and Tavoliere, where rainfall values are usually recorded as being lower than the 600 mm/year regional average. In this regard, we found the same zones to accumulate the highest values of HI and WI in both TW 1 (1961( -1990( ) and TW 2 (1991. These results are in line with Santos et al. (2012), who analysed macroclimate and viticultural zoning in Europe in the period 1950−2009 -which overlaps both our TWs -and found significant positive trends in WI and HI, particularly for southern Europe.
On the other hand, the areas with the lowest HI values were those located along the entire Murgian arch; this plateau crosses Apulia from north to south up to the hinterland of the provinces of Taranto and Brindisi, with altitudes often exceeding 300 m. a.s.l., and therefore poorly suitable in the past for the production of table grapes.
When comparing HI and WI, clear differences were found between the two bioclimatic indices in terms of the JD needed to reach the optimum of 2200. In particular, more JD were needed to reach at least 2200 WI than to reach 2200 HI, which is in line with Navrátilová et al. (2021), but with more fluctuations. These differences are mainly related to the periods that the formulas were based on: HI was based on the period of 1st April to 30th September, while WI was based on the period 1st April to 31th October, with the period of calculation ending after harvest. In addition, other differences can be explained by the fact that WI considers the vine to stop its development below 10 °C, while the base temperature threshold for accumulation is likely cultivar specific and is possibly as low as 5 °C (García de Cortázar-Atauri et al., 2009). Furthermore, HI takes into account the average of maximum and the mean daily temperatures, and also introduces a coefficient K which takes into account the latitude and, therefore, the different durations of daylight (Piña-Rey et al., 2020).
Overall, the generalised increase in average temperatures throughout the study area calls into question the maintenance of the elevation limit of 330 m a.s.l. for table grape production in Apulia. The concept of vineyard relocation, in fact, is now widely under discussion (Santos et al., 2020), which focuses mainly on heat and drought escape by combining two different but complementary strategies: the exploitation of wider latitudinal bands not restricted to conventional 30 to 50 degrees latitudes in both hemispheres (Karvonen, 2014), and a shift to cooler sites characterised by higher elevations and lower solar radiation (Fraga et al., 2017). It is known, in fact, that the required temperature for vines during the growing season is between 12-13 °C and 22-24 °C; Apulian temperatures normally fall within this range during the growing season (April-October), and occasionally late frost occurs. High temperatures and heatwaves, however, can increase the average seasonal temperatures to above 22-24 °C; this can induce heat stress on vine, particularly in dry climates (Schultz and Jones, 2010), and can often result in earlier phenological stages, thereby affecting fruit metabolism and in turn grape yield and quality (Venios et al., 2020).
The prediction of JD needed to reach optimum HI and WI, which was based on linear regression equations, yielded results of a systematic earlier shifting of grapevine phenology. The requirement of 2200 HI or WI corresponds to the Italia variety (chosen as an example of table grape cultivar grown in the Apulia region) reaching presumed maturity; in the present study it was obtained 9 to 28 days earlier in TW 2 respect to TW 1 , depending on the bioclimatic index and the elevation of the site. Our results confirmed those stated by other authors (Jones, 2007;Alikadic et al., 2019), who reported values of presumed maturity to be 6 to 25 days early for different grapevine varieties in Europe, with an average of 3 to 6 days response per 1 °C of warming during the last 30-50 years.
However, it is important to note that the most common vine training system for table grape, especially in Apulia, is the overhead tendone trellis (Puglia type). Tendone involves the distribution of the canopy on a continuous horizontal plane approximately 2 m above ground level. The canopy is supported by a trellis system and a grid of steel wires, on which shoots are trained to create an upper layer of leaves that overlaps and shadows the underside productive band of bunches (Pascuzzi, 2013).
In addition to this, table grape vineyards are normally coverd with plastic sheets, in order to protect grapes and to allow a harvest delay. Depending on the technical properties of the material from which it is made, the covering modifies the microclimate under the tendone by lowering solar radiation and air temperature, increasing RH with effects on vine water status, and preserving TSS concentration and TSS/TA ratio, polyphenol, flavonoids and volatile compounds, such as linalool, typical of 'Italia' muscat aroma (de Palma et al., 2019). Moreover, late table grape cultivars like Italia are normally covered late (at veraison) with plastic sheets to protect the clusters from rain and hailstorm, thus delaying harvest by 30 to 90 days.
That said, it is well-known that the technological maturity of table grapes tends to diverge and, like Greek producers (Koufos et al., 2020), Apulian growers do not regularly record observations of phenological phases and harvest dates, except for a few varieties and only sporadically. Therefore, as it was difficult to find vineyards in the region with recorded data, in the present study three vineyards representative of the production of table grapes in Apulia were identified in order to perform future projections for climate change.
As expected, RCP 4.5 and RCP 8.5 did not show any large differences until around 2040. In the long term (i.e., from 2061 onwards), the effects of the countermeasures foreseen by the RCP scenario 4.5 became effective in containing the growing trend of both Tmin and Tmax. This is in line with Lionello et al. (2014), who reported an acceleration in temperature increase and a decrease in precipitation in Apulia in the second half of the century. In the latter case, a huge decrease in DI would be expected in our future scenario, since it is based on rainfall. In any case, rainfall amount was not predicted to change much, with just a few millimetres difference from one scenario to another.
Indeed, a weakness of DI is the use of total monthly rainfall in the formula and the fact that it does not take into account the occurrence of extreme events; these may lead to abnormal daily precipitation rates that contribute to the monthly total, but is in fact water which is not available for vineyards.
In terms of Tmin, while all the vineyards were projected to become CI-2 (Warm nights) in 2061-2090, the effect of elevation in preventing a Tmax increase was evident in Farm BA (262 m a.s.l.), which had a milder RCP 8.5 future scenario than Farm BAT and Farm TA, which are located at sea level, or, in any case below 100 m a.s.l. (Figure 1) and showed an HI of > 3000 and a WI of > 2700. In contrast, neither the historic nor future scenario showed differences in terms of DI for any of the three vineyards; they even produced negative values which increased from one thirty-year series to the next.

CONCLUSIONS
This study focused on the spatialisation and variation over time of bioclimatic indices as an indicator of past and future climate change in a typical Apulian region for table grape production. In this region, temperatures have risen in the last sixty years and will continue to do so in the future at a different intensity according to the scenarios. An increase in heat accumulation will have the obvious consequence of earlier harvest dates, which in turn will impact both the production and management aspects of the vineyard. Targeted strategies in viticulture should thus be developed -and in some cases this has already started -to contain or overcome heat stress in the future by i) the adoption of technical expedients like vineyard sheet-covering to mitigate the microclimate in the vineyard, ii) the use of cultivars whose grapes can reach the optimal technological parameters, even when the duration of the phenological phases is shortened, and iii) vineyard relocation to elevations or latitudes unexplored until today.