European Grapevine Moth and Vitis vinifera L. Phenology in the Douro Region: (A)synchrony and Climate Scenarios

: The European grapevine moth ( Lobesia botrana ; Denis and Schiffermüller, 1775) is considered a key pest for grapevine ( Vitis vinifera L.) in the Douro Region, Portugal. The phenology of both the grapevine and the pest has changed in the last decades due to the increase in temperature. Here, we assess the potential impact of climate change on the (a)synchrony of both species. The results show that the phenological stages (budburst, ﬂowering and veraison) undergo an advancement throughout the region (at an ~1 km resolution) under a climate change scenario (Representative Concentration Pathways, RCP8.5) for the period 2051–2080, with respect to the historic period (1989–2015). For cv. Touriga Nacional and Touriga Franca , the budburst advances up to 14 days, whereas for ﬂowering and veraison the advancements are up to 10 days (mainly at low elevations along the Douro River). For the phenology of Lobesia botrana , earliness was also veriﬁed in the three ﬂights (consequently there may be more generations per year), covering the entire region. Furthermore, the third ﬂight advances further compared to the others. For both varieties, the interaction between the third ﬂight (beginning and peak) and the veraison date is the most relevant modiﬁcation under the future climate change scenario (RCP8.5, 2051–2080). The aforementioned outcomes from the phenology models help to better understand the possible shifts of both trophic levels in the region under future climate, giving insights into their future interactions.


Introduction
Plants and insects are dependent on the accumulation of heat units for their development. Plant phenology is directly influenced by weather and climate, with the following parameters being considered: temperature, photoperiod, relative humidity, precipitation and CO 2 [1]. Furthermore, pests are also influenced directly by climatic factors, although other features such as habitat structure, overwintering, food quality, and the length of the growing season also influence their development [1,2].
In recent years, changes in climatic conditions (particularly changes in temperature and the seasonal pattern of precipitation distribution) have already influenced the interactions between the pest and the host. In this way, it is expected that future climate change could have a strong impact on phenological stages, population dynamics, development Generalized sigmoid models (GSM) were developed to predict three phenological stages of the grapevine, which are budburst (Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie scale, BBCH 07), flowering (BBCH 65), and veraison (BBCH 81) [32]. These models were calibrated and validated with a large number of varieties (51, including cv. Touriga Nacional and Touriga Franca) from four Portuguese wine regions: Douro, Dão, Vinhos Verdes and Lisboa, as described by Reis et al. [5]. Phenological stages were recorded (overall period: 1990-2018) through site observations based on the Baggiolini scale [33] and assessed when at least 50% of a pre-defined and time-invariant (homogeneous) subsample of plants reached the corresponding stage. Lastly, cv. Touriga Nacional and Touriga Franca were selected, as these varieties are considered to have high expression and economic importance in the region.

LB Flight
Degree-day models were developed based on data collected over 20 years in the DDR. The flight model of LB, as described by Carlos et al. [17], is based on the degreedays required for the occurrence of its main flight events (beginning and occurrence of 50% of captured males, i.e., peak flight in pheromone traps). The lower and upper development thresholds were defined at 7.3 • C and 33 • C, respectively, according to Carlos et al. [17]. Lastly, in this study, the accumulation of degree-days since 1 January was used for each flight.

Robustness and Predictability
The models applied to predict the LB flights and the phenological stages of the grapevine (budburst, flowering and veraison) show high robustness and effectiveness based on data observed in the field. According to the literature, phenology models may operate with different inputs (temperature, precipitation, relative humidity, solar radiation and others) [5,11,17,[34][35][36][37] based on the scale of interest and specific objectives. The robustness of the models is particularly important, especially in the reliability of projections of climate change impacts [38]. However, a greater number of atmospheric elements may improve the model performance and its predictability. Elevation at an ~1 km grid resolution is also displayed.

Climate Dataset for the Historic Period and Future Scenarios
E-OBS originally (available at: https://www.ecad.eu/download/ensembles/download. php) (accessed on 29 September 2021) consisted of gridded daily minimum (T N ), maximum (T X ) and mean (T G ) temperature values and daily precipitation totals at a resolution of approximately 25 km [39]. Later, the mean sea level pressure was added as a grid variable [40]. The data were developed as part of the ENSEMBLES project to be used in climate change studies and for the validation of regional climate models [41]. In this context, a new high-resolution grid (~1 km) was used for the daily T G , T N and T X covering the DDR, based on the original dataset (~25 km, 1950-2015) [42]. In this way, from the set of observational data, a period of 27 years (1989-2015) was extracted.
The CORDEX project (Coordinated Regional Climate Downscaling Experiment; http: //wcrp-cordex.ipsl.jussieu.fr/) (accessed on 1 October 2021) [43] aims to provide an internationally coordinated framework for improving climate scenarios [44]. EURO-CORDEX is the European branch of this project, with 29 working groups [45]. A large set of simulations is available, depending on the combinations of GCM (Global Climate Model) and RCM (Regional Climate Model) [46]. For our purpose, the combinations used (GCM-RCM chains) were CNRMALADIN, ICHECDMI, IPSLINERIS and MPICLM. More details can be found in Reis et al. [22]. Daily T G , T N and T X were extracted from these climate model experiments.
For the future scenarios, two periods were considered: 2021-2050 (short-term) and 2051-2080 (medium-long term), under RCP4.5 and RCP8.5. In addition, the raw simulated data were available at an~11 km grid resolution, covering the entire region. However, the abovementioned simulated data were converted to a~1 km grid resolution to be compared with the historic period and used in the phenological models. In this way, three steps were undertaken to obtain high-resolution temperature datasets. First, the climate model simulated data, previously bias-corrected using E-OBS as a baseline [47], were subtracted from the E-OBS dataset on a daily timescale and at an~11 km resolution (daily anomalies). Second, a bilinear interpolation (in the latitude × longitude grid) was applied to interpolate the anomalies from~11 km to~1 km grid resolution. Third, the~1 km grid resolution anomalies were summed to the corresponding variable for the historic baseline period , obtained from a high-resolution dataset in Portugal [42].

Impact of Future Climate Scenarios
To predict the impact of climate change on LB flights, degree-day models were used for each flight (three flights). The starting point of the accumulation of degree-days was from 1 January. Lastly, the lower and upper development thresholds used were 7.3 • C and 33 • C, respectively [17].
The generalized sigmoid models (GSM) [5] for each phenological stage (budburst, flowering and veraison) were used for cv. Touriga Nacional and Touriga Franca, to evaluate the impact of climate scenarios on phenology. These models incorporate the F-forcing parameter for each variety and phenological stage. Moreover, this parameter is the critical thermal forcing (sum of units) required to complete a given phenological stage. According to Reis et al. [5], the following F-forcing parameters are considered: a) cv. Touriga Nacional: F = 46 (January 1st-budburst); F = 36 (budburst-flowering); F = 65 (flowering-veraison); b) cv. Touriga Franca: F = 50 (January 1st-budburst); F = 39 (budburst-flowering); F = 58 (flowering-veraison). GSM was applied with daily T G for the two future periods (2021-2050 and 2051-2080) and under both RCP4.5 and RCP8.5. In addition, the same method was applied for flight phenology models, using daily T N and T X .
The output comprised the simulated days (DOY) for the beginning and peak of each flight, in the case of LB, and the dates of the three grapevine phenological stages for the two varieties under study. Spatio-temporal maps were drawn up with the outputs obtained from LB and grapevine (cv. Touriga Nacional and Touriga Franca) at an~1 km resolution and also considering the temporal average for each period (future and historic periods).
When comparing the periods (future and historic) and assessing the statistical significance of the corresponding differences, we assumed two independent samples: X with m observations (2021-2050; 2051-2080) from the future period and Y with n observations (1989-2015) from the historic period. Averages (X and Y) and variances (S 2 x and S 2 y ) are shown below: We considered unequal variances, so the statistical test for unequal variance t-test is: For use in significance tests, the statistic test distribution has an approximate Student's t distribution, with the number of degrees of freedom given by Satterthwaite's approximation [48][49][50].
Lastly, the following null hypothesis was considered: The value obtained, if 0, indicates that the t-test does not reject the null hypothesis at the 5% significance level (p ≤ 0.05), even if equal variances are not assumed.

Interaction of Trophic Levels
We subtracted the Julian days (DOY) from both trophic levels for the historic  and future periods (2051-2080; RCP8.5) at an~1 km grid resolution over the entire region (for the development of the high-resolution grid see Section 2.3).
The difference obtained took into account only the third flight (beginning and peak) of the LB with the veraison date (cv. Touriga Franca and Touriga Nacional), according to the available phenology models (Section 2.2), since this is indeed the most relevant interaction for the Douro wine region in both the historic and future periods (2051-2080; RCP8.5).
We assumed that the synchrony between LB and the grapevine is when the pest can be hosted on the host. In more detail, it is when LB's third flight synchronizes with the veraison's date, i.e., the difference in their timings is small ((∆ DOY) ∼ = 0). For asynchrony, it is when there are differences, in Julian days, between the third flight and the veraison's date ((∆ DOY) = 0) [4].
Lastly, we defined for the entire region elevations below or equal to 400 m as low elevation areas. Between 401 and 499 m is an average elevation. Above or equal to 500 m is high elevation.

Historic Period
In the historic period , the days obtained by the phenology models (budburst, flowering and veraison) were validated with visual observations in the field. In this way, slight differences in phenology (period average DOY) were observed between cv. Touriga Nacional and Touriga Franca in the three phenological stages evaluated (Figure 2A-F).
For budburst, it was verified that cv. Touriga Nacional (DOY begin 70) occurs earlier than cv. Touriga Franca (DOY begin 75) at low to average elevations, mainly along the Douro river. At the high elevations, in both varieties, budburst is observed later, since the daily thermal accumulation is smaller (DOY ≥ 110).

Future Period: Climate Scenarios
We examined the difference between the averages of the future (RCP4.5 and RCP8.5) and historic periods. The future period average consists of two sub-periods, namely 2021-2050 and 2051-2080. The average of each sub-period is based on the ensemble of climate models previously mentioned. In the analysis of the difference between periods, the alternative hypothesis was followed, which is represented: H1: Future ( ) -Historic ( ) ≠ 0. In this way, we assessed the statistical significance of the differences in the future sub-periods compared to the historic period.
For budburst in both varieties, there is advancement at low and average elevations in DS and BC with respect to the historic period (up to 14 days). In addition, at high For flowering, the same earliness behaviour is observed for cv. Touriga Nacional and Touriga Franca at low to average elevations (DOY begin 131 and DOY begin 137, respectively). At the high elevations, it is observed for both varieties that flowering occurs later, with DOY ≥ 160, particularly for cv. Touriga Franca.
For veraison, an inversion of the observed dates of the two varieties is observed. Therefore, cv. Touriga Franca tends to be earlier than cv. Touriga Nacional at low to average elevations, along the Douro river (DOY begin 197 and DOY begin 199, respectively). At high elevations, the veraison happens later in both varieties (DOY ≥ 227). However, in some sites, it may not happen, likely due to insufficient daily thermal accumulation, resulting in difficulties in grape ripening.

Future Period: Climate Scenarios
We examined the difference between the averages of the future (RCP4.5 and RCP8.5) and historic periods. The future period average consists of two sub-periods, namely 2021-2050 and 2051-2080. The average of each sub-period is based on the ensemble of climate models previously mentioned. In the analysis of the difference between periods, the alternative hypothesis was followed, which is represented: H 1 : Future X -Historic Y = 0. In this way, we assessed the statistical significance of the differences in the future sub-periods compared to the historic period.
For budburst in both varieties, there is advancement at low and average elevations in DS and BC with respect to the historic period (up to 14 days). In addition, at high elevations in the CC, the advancement is smaller (the daily thermal accumulation is smaller), and, therefore, it is roughly between 3 to 6 days ( Figure 3A,B).
For flowering, there is a weaker advancement. In DS, especially at low and average elevations, this will be up to 10 days. Furthermore, in BC and CC, at low to average elevations (including along the Douro river), the projected advancement will be roughly 1 to 7 days. Lastly, at high elevations, in CC and BC, there is no advancement. This aspect is mainly observed for cv. Touriga Nacional, since it has a larger area covered ( Figure 3C,D).
For veraison, at low and average elevations in BC and DS, it was found that the advancement is projected to be up to 10 days for both varieties. In CC, at high elevations, there is slight advancement of 2 to 7 days ( Figure 3E,F). Overall, it can be concluded that there will be a greater advancement in budburst with respect to flowering and veraison for both varieties. In Figure S1 (Supplementary Material), the spatial-temporal maps with the ensemble-mean DOY (averaged over the ensemble of climate models) of each sub-period and climate scenario are presented.

Historic Period
We examined the average of the DOY period of the beginning and peak of the LB flight. Furthermore, this analysis was considered for the three flights under study ( Figure 4A-F).
For the first flight, at low and average elevations (warmer areas), it was observed that the beginning occurred earlier than at high elevations (cooler areas). Therefore, in those elevations, the general DOY begin range is 77-92, mainly along the Douro river. At high elevations, the daily thermal accumulation is smaller and, consequently, there is a delay equal to or greater than 28 days (DOY begin ≥ 105). However, there are also locations with high elevations, in the range of DOY 93 to 105. At the peak, more zones are observed near the earliest date, once that it covers more area in five intervals of the map (between DOY 100 and 121). In this way, at low and average elevations, the general DOY peak range is 100-121. Lastly, we observed a difference of 40 days or more at some high elevations (DOY peak ≥ 140).
In the second flight, the beginning is similar to the first flight in terms of earliness throughout the region. Hence, mainly along the Douro river, at low and average elevations, the general DOY begin range is 151-166. In areas farther away from the river, the phenology delays at some high elevations, with the difference being equal to or greater than 29 days (DOY begin ≥ 180). The peak is earlier at low and average elevations because its general DOY peak range is 164-180. At higher elevations, there is a delay equal to or greater than 31 days (DOY peak ≥ 195).
In the third flight the earliest date reached is observed (DOY begin = 197), mainly along the river. Additionally, it is also noteworthy that, at some high elevations towards the DDR limits, dates close to the earliest were observed. However, at some high elevations, there is a difference in the earliest date, at least 38 days (DOY begin ≥ 235). Consequently, the flight may not happen, because the thermal accumulation is insufficient, as the diapause is induced by the photoperiod [51]. The peak, at low and average elevations, reaches the earliest date (DOY peak = 211). At high elevations, there is a difference to the earliest date of elevations (including along the Douro river), the projected advancement will be roughly 1 to 7 days. Lastly, at high elevations, in CC and BC, there is no advancement. This aspect is mainly observed for cv. Touriga Nacional, since it has a larger area covered ( Figure 3C,D).
For veraison, at low and average elevations in BC and DS, it was found that the advancement is projected to be up to 10 days for both varieties. In CC, at high elevations, there is slight advancement of 2 to 7 days ( Figure 3E,F). Overall, it can be concluded that there will be a greater advancement in budburst with respect to flowering and veraison for both varieties. In Figure S1 (Supplementary Material), the spatial-temporal maps with the ensemble-mean DOY (averaged over the ensemble of climate models) of each subperiod and climate scenario are presented. Overall, on the one hand, we examined that there are different dynamics in the phenology of each flight, as it depends at least on elevation and thermal accumulation to reach a certain stage. On the other hand, this dynamic means that there are different timings of the phenology in the region (e.g., at low elevations, especially near the river, the peak of the first flight occurs, while LB is still at the beginning of the first flight at high elevations).

Historic Period
We examined the average of the DOY period of the beginning and peak of the LB flight. Furthermore, this analysis was considered for the three flights under study ( Figure  4A-F). .

Future Period: Climate Scenarios
We analyzed the difference of the averages of the future and historic periods for the flight phenology. A similar statistical methodology was applied to the grapevine phenology. It was verified that there are statistically significant changes with respect to the historic period in the three sub-regions, mostly in the sub-period 2051-2080 in the two climate scenarios (RCP4.5 and RCP8.5). However, we analyze herein with more detail the most severe climate scenario (RCP8.5), leaving RCP4.5 results for Figure S2 (Supplementary Material).
In the first flight, it was observed that the advancement about the historic period is up to 14 days. This advance is observed at low and average elevations, mainly in BC and DS. However, in CC and BC, the advancement is smaller, at high elevations, but the beginning advances further regarding the peak of the flight.
In the second flight, in general, the peak advances more near the beginning. At the peak, the advancement of phenology regarding the historic period is up to 12 days. This advancement covers practically the entire region (BC, CC and DS) regardless of the elevation. At the beginning of the flight, it advances less, especially at average and high elevations, in CC and BC, with this advancement being roughly between 1 to 5 days.
On the third flight, there was greater advancement compared to the other flights. Another aspect considered, in general, is the peak that advances further regarding the beginning of the flight. At peak, the advancement can reach 22 days, in BC, CC and DS (the latter in a smaller area). However, at the beginning of the flight, the advancement observed is lower in the three sub-regions, which can reach up to 20 days.
In general, it can be concluded that there is generalized advancement, covering the entire DDR for the three LB flights, adding that the third advances further compared to the others ( Figure 5A-F). Finally, in annexes S3 and S4, the spatial-temporal maps with the DOY average (set of climate models) of the evaluated sub-period and climate scenarios (RCP4.5 and RCP8.5) are available.

Historic Period
The analysis of the interaction between the trophic levels (LB and grapevine) is important to examine whether there is synchrony between them, which underlies likely interactions between the LB and grapevine. The most interesting interaction was found between the third flight (beginning) and the veraison date in the historic period . Synchrony between the beginning of the third flight and the veraison date of both varieties (cv. Touriga Nacional and Touriga Franca) is observed. In greater detail, the beginning of the third flight is synchronized with the two varieties when the differences in their corresponding timings are small ((∆ DOY) ∼ = 0). That phenological interaction is mainly observed in two sub-regions for cv. Touriga Franca: CC and DS ( Figure 6B; white color on the map), whereas it is observed in the three sub-regions for cv. Touriga Nacional ( Figure 6A; white color on the map).
On the one hand, concerning asynchrony, this is observed when the beginning of the third flight occurred later than the date of the veraison in cv. Touriga Nacional in the three sub-regions BC (≤+10 (∆ DOY) ≤ +1), CC and DS (≤+4 (∆ DOY) ≤ +1), regardless of the elevation ( Figure 6A; light to dark orange color on the map). The outcomes for cv. Touriga Franca in the three sub-regions are as follows: BC (≤+14 (∆ DOY) ≤ +1), CC (≤+4 (∆ DOY) ≤ +1) and DS (≤+8 (∆ DOY) ≤ +1) ( Figure 6B; light to dark orange color on the map). On the other hand, asynchrony was also observed when the beginning of the third flight occurred earlier than the date of the veraison on cv. Touriga Nacional in the three sub-regions BC (≤−4 (∆ DOY) ≤ −1), CC and DS (≤−10 (∆ DOY) ≤ −1), irrespective to the elevation ( Figure 6A; light to dark blue color on the map). For cv. Touriga Franca, the results in two sub-regions are the following: CC and DS (≤−8 (∆ DOY) ≤ −1) ( Figure 6B; light to dark blue color on the map).

Future Period
Anticipating the impact of climate change on the two trophic levels (LB and grapevine) is of high importance. Therefore, predicting possible changes in both trophic levels from simulations obtained by phenology models can help to understand this complex mechanism. Therefore, the future sub-period (2051-2080) in a more severe climate scenario (RCP8.5) was analyzed in greater detail. The results obtained suggest that there may be synchrony ((∆ DOY) ∼ = 0) between the beginning of the third flight and the veraison date (cv. Touriga Nacional and Touriga Franca), at both average and high elevations, mainly in BC and DS ( Figure 7A,B; white color on the maps). For asynchrony, the beginning of the third flight is projected to occur earlier than the veraison date in both varieties over most of the region (Figure 7A,B; light to dark blue color on the maps). The difference will be ≤−14 (∆ DOY) ≤ −1, except in the BC sub-region, which will occur later, mainly at high elevations (≤+4 (∆ DOY) ≤ +1) ( Figure 7A,B; light orange color on the maps). Consequently, at the peak of the same flight, there will be synchrony, mainly at low elevations (along the Douro river) in the CC and DS sub-regions and both varieties ( Figure 7C,D; white color on the maps). Lastly, asynchrony is observed over the three sub-regions BC (≤+18 (∆ DOY) ≤ +1), CC (≤+10 (∆ DOY) ≤ +1) and DS (≤+14 (∆ DOY) ≤ +1). In this case, the peak of the flight will occur later about the date of the veraison of the two varieties ( Figure 7C,D; light to dark orange color on the maps).
Agronomy 2022, 12, x FOR PEER REVIEW 12 of 19 (the latter in a smaller area). However, at the beginning of the flight, the advancement observed is lower in the three sub-regions, which can reach up to 20 days. In general, it can be concluded that there is generalized advancement, covering the entire DDR for the three LB flights, adding that the third advances further compared to the others (Figure 5A-F). Finally, in annexes S3 and S4, the spatial-temporal maps with the DOY average (set of climate models) of the evaluated sub-period and climate scenarios (RCP4.5 and RCP8.5) are available. In addition, the Douro river and main tributaries are also outlined (thick blue line). hand, asynchrony was also observed when the beginning of the third flight occurred earlier than the date of the veraison on cv. Touriga Nacional in the three sub-regions BC (≤−4 Δ DOY ≤ −1), CC and DS (≤−10 Δ DOY ≤ −1), irrespective to the elevation ( Figure 6A; light to dark blue color on the map). For cv. Touriga Franca, the results in two sub-regions are the following: CC and DS (≤−8 Δ DOY ≤ −1) ( Figure 6B; light to dark blue color on the map).

Future Period
Anticipating the impact of climate change on the two trophic levels (LB and grapevine) is of high importance. Therefore, predicting possible changes in both trophic levels from simulations obtained by phenology models can help to understand this complex mechanism. Therefore, the future sub-period (2051-2080) in a more severe climate scenario (RCP8.5) was analyzed in greater detail. The results obtained suggest that there may be synchrony (Δ DOY ≅ 0) between the beginning of the third flight and the veraison date (cv. Touriga Nacional and Touriga Franca), at both average and high elevations, mainly in BC and DS ( Figure 7A,B; white color on the maps). For asynchrony, the beginning of the third flight is projected to occur earlier than the veraison date in both varieties over most of the region (Figure 7A,B; light to dark blue color on the maps). The difference will be ≤−14 Δ DOY ≤ −1, except in the BC sub-region, which will occur later, mainly at high elevations (≤+4 Δ DOY ≤ +1) ( Figure 7A,B; light orange color on the maps). Consequently, at the peak of the same flight, there will be synchrony, mainly at low elevations (along the Douro river) in the CC and DS sub-regions and both varieties ( Figure  7C,D; white color on the maps). Lastly, asynchrony is observed over the three sub-regions BC (≤+18 Δ DOY ≤ +1), CC (≤+10 Δ DOY ≤ +1) and DS (≤+14 Δ DOY ≤ +1). In this case, the peak of the flight will occur later about the date of the veraison of the two varieties ( Figure  7C,D; light to dark orange color on the maps).
In conclusion, the simulations indicate that the rhythm of LB may be different in comparison to the grapevine phenology under a future climate scenario.

Discussion
This study aimed to analyze the impact of climate change on the cycle of grapevine and LB in DDR. Our results indicate that the phenological stages will occur earlier across the region (at a ~1 km resolution) for 2051-2080 under RCP8.5. In this way, the budburst advances up to 14 days, flowering and veraison up to 10 days (mainly at low elevations, along the Douro river). The generalized advancement is valid for both target varieties, cv. Touriga Nacional and Touriga Franca. This analysis is in line with several recent studies, In conclusion, the simulations indicate that the rhythm of LB may be different in comparison to the grapevine phenology under a future climate scenario.

Discussion
This study aimed to analyze the impact of climate change on the cycle of grapevine and LB in DDR. Our results indicate that the phenological stages will occur earlier across the region (at a~1 km resolution) for 2051-2080 under RCP8.5. In this way, the budburst advances up to 14 days, flowering and veraison up to 10 days (mainly at low elevations, along the Douro river). The generalized advancement is valid for both target varieties, cv. Touriga Nacional and Touriga Franca. This analysis is in line with several recent studies, which highlight the impact of climate change on phenology, namely earlier phenological stages and shorter intervals [15,27,30,[52][53][54]. Moreover, earlier phenological stages will result in heterogeneity. Earlier budburst and flowering can result in a substantial increase in frost damage risks [30]. Lastly, the climate projections of the grapevine phenological stages provide valid information for medium-long term planning/management [55].
In the three sub-regions, at high elevations, thermal accumulation and phenology advancement is weaker. Thus, for both varieties (cv. Touriga Nacional and Touriga Franca), there was less advancement in the budburst, mainly at high elevations of the CC sub-region. For flowering, no advancement was found for BC and CC (cv. Touriga Nacional with the largest area covered). Finally, for the veraison, there was lower advancement in phenology, mainly at high elevations of the CC sub-region in both varieties. According to Caffarra and Eccel [56], the increase in temperature can make more areas suitable for viticulture at higher elevations, where the climate is currently too cold. However, the suitability of the viticulture areas must be carefully evaluated. In temperate zones, the temperature variations will affect the development of grapevine, which depends at least on local climatic conditions, the cultivated variety and cultural practices.
For LB phenology, a similar statistical methodology was applied to the grapevine. For the first flight, the advancement is up to 14 days for 2051-2080 under RCP8.5 (mainly at low and average elevations in BC and DS). For the second flight, the flight peak advances more than its beginning. The advancement of the peak of the flight is up to 12 days, almost covering the whole region (BC, CC and DS), regardless of the elevation. The third flight experiences the strongest advancement when compared to the other flights. At peak, the advancement can reach 22 days in BC, CC and DS. However, the advancement is smaller at the beginning of the flight, for which it can reach 20 days at low elevations in the three sub-regions. In general, it can be concluded that there is generalized advancement throughout the DDR for the three LB flights, though the third flight undergoes a more significant advancement. The earliness of LB flights in this climate scenario is in line with other previous studies [2,23,57]. According to Reis et al. [22], it was observed that LB phenology advances in the locations evaluated in the DDR. The earliness observed in the three flights in the 2021-2080 period was 7 to 12 days, under RCP4.5, and 15 to 24 days, under RCP8.5, when compared to the historic period (2000-2019). Therefore, that study suggests that a fourth complete flight is likely to occur in the future. However, the number of days with excessively high temperatures (i.e., above the upper threshold for development (>33 • C)), is projected to increase. This may decrease LB populations (in the total number of male catches in the traps/percentage of bunches attacked on the second and third flights/generations) [22] in some warmer zones in DDR.
The simulations obtained by the phenology models (LB and grapevine) suggest that the warm climate will prevent the two trophic levels from developing at the same rhythm in a climatic scenario. The most interesting results regarding the interactions between both trophic levels were between the third flight and the veraison. The simulations suggest that there is synchrony between the beginning of the third flight and veraison in both varieties (Cv. Touriga Nacional and Touriga Franca) for 2051-2080 under RCP8.5. The observed synchrony is mainly in the BC and DS sub-regions at average and high elevations. Another predicted difference compared to the historic period is the peak of the same flight with the veraison date. In this case, the synchronization can occur mainly at low elevations, in the CC and DS sub-regions in both varieties. The projections reveal that a change in the phenology between the two trophic levels can happen under climate change. Our simulations are consistent with other preceding studies [4,[58][59][60].

Conclusions
This study provides important insights regarding the shifts in phenology in DDR. This evolution is related to (a)synchrony of LB and grapevine in DDR. In addition, this study provides new knowledge at a regional scale and with a medium-long term projection (2051-2080). The projection mainly takes into account the RCP8.5 climate scenario, but the results for RCP4.5 are also provided.
The simulations highlight generalized advancement of the phenological stages (budburst, flowering and veraison), covering the entire region (at an~1 km resolution) and for both selected varieties. The earliness is greatest in budburst in the sub-period 2051-2080 and under RCP8.5. Regarding the LB phenology, for the same climate scenario and sub-period, there is also generalized advancement of the three flights (beginning and peak), covering the entire region. Furthermore, it is important to note that the third flight is projected to undergo the strongest advancement compared to the others. Therefore, the most interesting interactions in future climate scenarios (RCP8.5) are the third flight (beginning and peak) with the veraison date for both varieties.
Lastly, simulations with phenology models can help to understand the impact of climate change on the phenology of both trophic levels (LB and grapevine). Therefore, future findings should study interactions in more detail at different space and time scales, with the aim of better adaptation of the species.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/agronomy12010098/s1, Figure S1: Phenology maps for the future period (2051-2080, RCP8.5) covering the entire Douro Demarcated Region. A color scale represents the Julian days (DOY) obtained by the phenology models for each phenological stage, namely budburst (A) and (B), flowering (C) and (D) and veraison (E) and (F). The varieties chosen were cv. Touriga Nacional (left) and Touriga Franca (right). In addition, it is integrated into the Douro river (thick blue line) along with the three sub-regions. Figure S2: Difference between future and historic (RCP4.5) phenology maps covering the entire Douro Demarcated Region. The color scale shows the difference between future and historic ((∆ DOY)) for each beginning (left) and peak (right) of the flight, namely first (A) and (B), second (C) and (D) and third (E) and (F). In addition, it is integrated into the Douro river (thick blue line) along with the three sub-regions. Figure S3: Phenology maps for the future period (2051-2080, RCP4.5) covering the entire Douro Demarcated Region. The color scale represents the Julian days (DOY) obtained by the phenology models for each beginning (left) and peak (right) of the flight, namely the first (A) and (B), second (C) and (D) and third (E) and (F). In addition, it is integrated into the Douro river (thick blue line) along with the three sub-regions. Figure S4: Phenology maps for the future period (2051-2080, RCP8.5) covering the entire Douro Demarcated Region. The color scale represents the Julian days (DOY) obtained by the phenology models for each beginning (left) and peak (right) of the flight, namely the first (A) and (B), second (C) and (D) and third (E) and (F). In addition, it is integrated into the Douro river (thick blue line) along with the three sub-regions.
Portuguese Foundation for Science and Technology, under the project UIDB/04033/2020. Fátima Gonçalves is also grateful to the European Social Fund (ESF) through the Regional Operational Program North 2020, within the scope of the Program "Contratação de Recursos Humanos Altamente Qualificados", Norte-06-3559-FSE-000188 and to the FCT for financial support by national funds FCT/MCTES to CIMO (UIDB/00690/2020).