The annual EUROSTAT crop yield statistics at the NUTS 2 level (i.e., basic regions for the application of regional policies) for the period 1975–2009 were obtained from the European Commission (http://appsso.eurostat.ec.europa.eu/nui/show.do?dataset=agrracsh&lang=en, downloaded September 2014). This dataset documents all primary crop types, including wheat and grain maize. The wheat yield data documented in EUROSTAT are mixed with both winter and spring wheat yield observations. Planting areas of winter wheat exceed 83% of the total planting area of wheat crops (both spring and winter wheat) over TE (Sacks et al., 2010). In this study we focus on wheat, which has relatively continuous spatiotemporal coverage (Wu et al., 2014). In the EUROSTAT database, crop yield is reported as the amount of dry matter suitable for consumption (Moors et al., 2010). We first examined the length of available wheat yield records for each NUTS 2 region, and discarded crop yield time series <10 years in length. We then rasterized the crop yield data from the regions at NUTS 2 level into gridded data with a spatial resolution of 0.5° matching the spatial resolution of the gridded climate and land cover data. The rasterizing process followed that of Wu et al. (2014). We calculated the gridded winter wheat yield, Y¯, using Equation (1):

where A_k is the area fraction of the k^th NUTS 2 region meeting in one grid, Y_k is the wheat yield of the k^th region recorded in EUROSTAT and n is the number of regions meeting within one grid. In cases where extreme values were present in Y_k (i.e., above or below three standard deviations of the crop yield series for kth region during 1975–2009), we assigned the corresponding A_k as zero to minimize possible bias introduced by unrealistic yield records prior to the rasterizing process (Wu et al., 2014).

The gridded dataset (E-OBS V11.0) of daily maximum and minimum temperature and precipitation, with a spatial resolution of 0.5°, during 1950–2014 was obtained from the European Climate Assessment project (http://www.ecad.eu/download/ensembles/download.php) (Hofstra et al., 2009). The WFDEI climate forcing for driving the ORCHIDEE_crop model (for detail see section Calculation of Growing Degree Days (GDD) and VDD) was obtained from the ISIMIP project (https://www.pik-potsdam.de/research/climate-impacts-and-vulnerabilities/research/rd2-cross-cutting-activities/isi-mip/about/isi-mip2).

GDD for each grid was calculated using the daily minimum (T_min) and maximum (T_max) temperatures, as shown in Equation (2):

where, t and n are days from planting and the total number of days from planting to harvest for winter wheat, respectively. Degree days (DD) was calculated by Equation (3):

where T_t is the mean daily temperature, calculated as (T_min + T_max)/2. In this study, we set a base temperature of 0 °C for winter wheat (McMaster and Wilhelm, 1997). We determined the growing season (from planting to harvest) of winter wheat for each grid over TE, based on an observational gridded crop calendar dataset (the SAGE dataset, http://nelson.wisc.edu/sage/data-and-models/crop-calendar-dataset/index.php) (Sacks et al., 2010). This globally available crop calendar dataset has been widely used in global/regional crop yield simulations and statistical analyses of crop yield–climate interactions (Lizumi et al., 2013; Ray et al., 2015). However, this static crop calendar dataset did not allow us to consider the possibly changing planting and harvesting dates for winter wheat over TE.

Two different kinds of vernalization effectiveness functions of daily temperature (VF) were used in this study to calculate the VDD during the vernalization period (from planting to start of anthesis) and to validate each other. One VF followed a quadratic function of daily temperature used in the ORCHIDEE_crop model (Supplementary Figure S1)—the same function used in the STICS model, in detail see Brisson et al. (2008) and Wu et al. (2016)—as shown in Equation (4):

where T_t, T_opt, and T_amp are the mean daily temperature, optimal vernalization temperature and thermal semi-amplitude of vernalizing effect, respectively. T_opt and T_amp for winter wheat are parameterized in this study as 6.5 and 10.0 °C, respectively, according to Brisson et al. (2008). Another VF was derived from the SIMTAG model (Stapper, 1984), which is a piece-wise linear function of mean daily temperature and similar to that used in the CERES-Wheat model and that of Reinink et al. (1986). The response of VF to T_t in the SIMTAG model is illustrated in Equation (5) (Supplementary Figure S1):

We then calculated the VDD for winter wheat in each grid by multiplying T_t and the corresponding VF factors during the vernalization period. Unfortunately, spatially detailed observations of the start of anthesis for winter wheat over TE are lacking. We alternatively used a recently developed process-based agro-land surface model ORCHIDEE_crop (Wu et al., 2016) to simulate the spatial patterns in the anthesis date for winter wheat over TE during 1979–2009, driven by WFDEI climate forcing and the SAGE gridded planting dates for winter wheat.

In the ORCHIDEE_crop model, the vernalization requirement for winter wheat growth is determined by the genotype-dependent number of vernalizing days (JVC) and the vernalizing value of a given day (JVI), which are counted since germination (IGER). A minimum number of vernalizing days (JVC_min) is required for winter wheat—in our study it was parameterized as 7 according to Brisson et al. (2008). Then, the progress in winter wheat vernalization, RFVI, increases gradually after IGER until it reaches the value 1, following Equation (6):

where JVI_i is the vernalizing value of a given day i since IGER, I is the day when RFVI reaches the value 1—for detail see Brisson et al. (2008). The ORCHIDEE_crop model has been validated for both phenology (e.g., harvest date) and yields of winter wheat across seven sites covering a variety of climate conditions over Europe (Wu et al., 2016). In this study, we compared the simulated anthesis dates with the limited field observations in Europe (Supplementary Table 1). Our comparison showed that the ORCHIDEE_crop model can capture the spatial pattern in anthesis dates of winter wheat across different climate conditions in Europe (y = 1.11x – 15.4, R² = 0.69, p < 0.05) (Supplementary Figure S2). Thus, we used the simulated spatial distribution in anthesis dates from the ORCHIDEE_crop model during 1975–2009 to calculate the VDD over TE.

Multivariate linear regression (MLR) was applied to each grid with available winter wheat yields and corresponding climate data to calculate the interannual sensitivity of winter wheat yields to variations of different climate factors over TE for the period 1975–2009, as shown in Equation (7):

where Y is the detrended time series of annual crop yields, γ_i is the regression coefficient for the i^th predictable detrended climate variable X_i; n and ϵ is the number of independent variables and residual error, respectively. GDD, VDD and total growing-season precipitation (PRE) were included in the MLR analyses. A quadratic polynomial function of time was fitted to the wheat-yield series to account for possible non-steady changes in technical improvement, agricultural investment and human management (e.g., fertilization and pesticide use) (Wu et al., 2014). Then, Y was calculated by subtracting the fitted values from the raw yield series. Note that the non-linear effects of technical improvement and human management on winter wheat yields could not be completely eliminated even though we filtered the annual yield series (Lobell and Field, 2007; Wu et al., 2014). The GDD, VDD, and PRE were also detrended using the same method as for crop yields. The detrended time series allowed us to address the effects of interannual variations in climate factors on the interannual variations of winter wheat yields (Wu et al., 2014). We did not constrain a positive contribution of GDD to wheat yields as was done by Butler and Huybers (2013). Their study introduced another variable considering the effects of extremely high temperature on crop development and grain yields, whereas we calculated GDD as the total thermal accumulations during the whole growing season.

Solar radiation, especially during grain filling, is reported to be a very important climate factor regulating winter wheat yields (Murchie et al., 2009; Xiong et al., 2012). However, an inclusion of detrended solar radiation (derived from the ISIMIP project, see section Calculation of Growing Degree Days (GDD) and VDD) summed from the start of anthesis to harvest did not markedly improve the regression model performance, with mean R² = 0.41 (range 0.28–0.88), indicating that the yields of winter wheat over TE were primarily regulated by variations in thermal and water conditions. To simplify our model structure we did not include solar radiation in our final analyses. To further verify the results of MLR we also performed independent ridge regression analyses to investigate the interannual sensitivity of winter wheat yields to climate variations. Ridge regression has been demonstrated to be effective for addressing collinearity among independent variables in multivariate regression (Hoerl and Kennard, 1970).

In this study, the effects of some other factors, such as photoperiod (Brooking and Jamieson, 2002), carbon dioxide fertilization (Lobell and Field, 2008) and fertilizer application (Godden et al., 1998) on winter wheat yields were not considered. The sensitivity of winter wheat yields to climate factors in six major Köppen–Geiger climate zones was also investigated. The six major climate zones, including arid (grouping BWk, BWh, BSk, and BSh), warm temperate humid (grouping Cfa, Cfb, and Cfc), warm temperate arid (grouping Csa, Csb, Csc, Cwa, Cwb, and Cwc), cold humid (grouping Dfa, Dfb, Dfc, and Dfd), cold summer dry (grouping Dsa, Dsb, Dsc, and Dsd) and tundra, were reclassified based on Köppen–Geiger climate classification (Kottek et al., 2006) over TE (for detail see Supplementary Figure S3).

Exponential quantile regressions at 25th, 50th, and 75th percentiles were applied to investigate the spatial relationships between γ_VDD and the multi-year mean VDD using the R package of quantreg (V5.21, https://cran.r-project.org/web/packages/quantreg/index.html) (Ruppert et al., 2012). Analyses were performed for data derived from different categories of regions, including all cropland of TE and cropland with significant (p < 0.05) fitted MLR models. A global land cover type data product MCD12Q1 (http://glcf.umd.edu/data/lc/) at a spatial resolution of 0.5° was used to identify cropland over TE. This dataset identifies 17 land cover classes based on the International Geosphere Biosphere Programme vegetation classification scheme (Loveland et al., 2000), which includes 11 natural vegetation classes, three developed and mosaicked land classes, and three non-vegetated land classes. We only retained grids with cropland fraction ≥5% for our final analyses.

The MLRs between winter wheat yields and GDD during growing season (from planting to harvest), VDD (based on VF function in the ORCHIDEE crop model) from planting to start of anthesis, and total growing-season PRE for the period 1975–2009 were significant (p < 0.05) in ~61% of grids with sufficiently long records for crop yields over TE. The goodness-of-fit of MLRs had mean R² = 0.39 (range 0.27–0.88) (Supplementary Figure S4A). The goodness-of-fit of MLR was much smaller (mean R² = 0.28, range 0.15–0.67) when VDD was excluded from the MLR (Supplementary Figure S4B).

Generally negative interannual sensitivity of winter wheat yields to GDD (γ_GDD) was observed in most of TE, with a mean of −2.30 ± 1.31 kg ha⁻¹ GDD⁻¹ (Figure 1A). This result indicates that higher GDD accumulation during the growing season can lead to lower winter wheat yields over TE. There was more negative γ_GDD in parts of temperate-humid regions and cold-humid regions (e.g., eastern France and eastern Europe) of TE (Figure 1A, Supplementary Figure S3). The spatial pattern in the interannual sensitivity of winter wheat yields to total growing-season precipitation (γ_PRE) showed great variation across TE (Figure 1C). In general, there was negative γ_PRE in western TE (mainly France) and positive γ_PRE in parts of eastern TE (Figure 1C).

In addition to these findings, we focused on γ_VDD during 1975–2009. Our results showed a generally positive γ_VDD (based on VF function in the ORCHIDEE_crop model) in most of TE, with a mean of 2.8 ± 1.5 kg ha⁻¹ VDD⁻¹ (Figure 1B, Table 1). Higher γ_VDD was observed in continental TE (Figure 1B). To verify the spatial patterns in γ_GDD, γ_VDD and γ_PRE derived from MLR, we performed an independent ridge regression analysis to get a more “accurate” estimation for regression coefficients of predictor variables (Hoerl and Kennard, 1970). The independent ridge regression analysis showed consistent spatial patterns in the interannual sensitivity of winter wheat yields to variations of climate factors (Supplementary Figure S5). Notably, ridge regression analyses produced consistent spatial patterns in the interannual sensitivity of winter wheat yields with two different VDD calculations based on different response functions between vernalization effectiveness and daily temperature (for detail see section Materials and Methods) (Supplementary Figures S5, S6).

The fitted 25th, 50th, and 75th non-linear quantile regressions consistently revealed significant (p < 0.05) negative exponential relationships between γ_VDD of winter wheat yields and the multi-year mean VDD for two different calculations of VDD during 1975–2009 in cropland grids with a significant fitted MLR model (y = 0.19e^−0.005x and y = 0.07e^−0.00083x for 50th quantile exponential fits for VDD calculated based on VF from the SIMTAG and ORCHIDEE crop models, respectively) (Figure 2). There was a higher γ_VDD of winter wheat yields in cropland grids with lower multi-year mean VDD (Figure 2, Supplementary Figure S7).

Consistent patterns were also found for data derived from all cropland grids (including grids with and without a significant MLR model) over TE for two different VDD calculations (Figure 3). This verified that robustness of the increasing γ_VDD along with a decreasing VDD was neither susceptible to analysis regions (i.e., grids with significant MLR model vs. all grids within croplands of TE) nor to the different calculations of VDD (Figures 2, 3).

There was great spatial difference in the linear trends in VDD calculated by VF from the ORCHIDEE_crop model during 1975–2009 in TE (Supplementary Figure S8). Our analyses revealed decreasing VDD in parts of southern and continental TE during 1975–2009. Investigation of the relationships between the linear trends of VDD and multi-year mean VDD (based on VF function in the ORCHIDEE_crop model) during 1975–2009 revealed a more rapid decrease in VDD in regions with higher VDD over TE (Figure 4).

The observed negative γ_GDD over TE is consistent with previous findings that winter wheat yields in this region are susceptible to climate warming (Porter and Semenov, 2005; Wu et al., 2014) and warmer temperatures during the growing season lead to a stagnation, or worse, decreases in crop yields (Brisson et al., 2010; Semenov and Shewry, 2011; Wu et al., 2014). We identified a much more negative γ_GDD in parts of temperate-humid and cold-humid regions (e.g., eastern France and eastern Europe) of TE, where crops are suffering more severe high temperature stress (Ciais et al., 2005; Wu et al., 2014). Such discernible negative effects of climate warming on the widely observed stagnation or collapse in yields of wheat have also been observed in many other regions worldwide (Calderini and Slafer, 1998; Lobell and Field, 2007; Brisson et al., 2010; Ray et al., 2012), despite the confounding effects of socio-economic and agricultural policy (Lobell, 2012). It was recently reported that the risk of high temperatures to winter wheat yields will increase further in the future, particularly around sensitive development stages, such as flowering (Semenov and Shewry, 2011).

The γ_PRE of winter wheat yields showed a clear spatial pattern—predominately positive and negative in regions characterized with temperate continental and oceanic climate, respectively (Figure 1C). Winter wheat grown in temperate continental regions with drier spring and summer suffers water stress during the growing season (Ciais et al., 2005; Maracchi et al., 2005; Wu et al., 2014). Thus, winter wheat yields exhibit a positive response to variations in water supply, and this conclusion was confirmed by the positive relationship between winter wheat yields and total growing-season precipitation in temperate continental regions (Figure 1C). In contrast, the positive γ_PRE disappeared in humid regions of TE, probably because of higher precipitation and the great spatial variations in soil water conditions among climate zones. Given the low irrigation intensity over TE, particularly for winter wheat, changes in precipitation amounts and seasonality accompanying climate warming will influence the magnitude of drought stress on winter wheat yields (Olesen and Bindi, 2002). Current trends and future scenarios depict an increase in drought stress for TE (Spinoni et al., 2015), which will potentially increase competition for water between agriculture and other sectors and uses (Lavalle et al., 2009).

Despite the great spatial variations in γ_PRE, it exhibited an increasing interannual sensitivity of winter wheat yields to water availability over TE during recent decades (Lavalle et al., 2009; Wu et al., 2014), illustrating an increasing effect of drought stress on winter wheat yields. Increase in water use efficiency caused by increasing carbon dioxide concentration could partly compensate for some of the negative effects of increasing water limitation. However, lower harvestable yields, higher yield variability and reduction in suitable areas for winter wheat are expected in these areas (Maracchi et al., 2005; Lavalle et al., 2009).

Positive γ_VDD of winter wheat yields in most of TE illustrates a crucial effect of vernalization fulfillment on winter wheat yields. We identified higher positive γ_VDD in regions with lower multi-year mean VDD in croplands of TE, as revealed by robust analyses of non-linear quantile regressions for two kinds of calculations of VDD. These statistical findings qualitatively point out that winter wheat yields are sensitive to variations in VDD, particularly in regions with low VDD. It can be reasonably expected that winter wheat yields will dramatically decline in response to even a subtle decrease in VDD (Figure 2). Our results highlight the possible vulnerability of winter wheat yields to warming-mediated interannual variations in vernalization fulfillment in regions with low VDD, in particular accompanying much higher warming rates in winter and spring (Cohen et al., 2012; IPCC, 2013). Warmer temperatures during the vernalization period can reshape the crop vernalization fulfillment process across diverse bioclimatic regions. Vernalization variations affect the floral initiation time, tiller and spikelet formation, onset of reproductive growth, grain numbers and consequently crop yields (Brooking, 1996; Fischer, 2011). Sufficient vernalization fulfillment before anthesis could partially compensate for the adverse effect of global warming on winter wheat heat stress (Rezaei et al., 2015; Stratonovitch and Semenov, 2015).

The predicted future warming over TE with increases in climate variability (Schär et al., 2004; Seneviratne et al., 2006; Wu et al., 2014) and in the frequency of high temperature and heat waves (Ciais et al., 2005; Perkins et al., 2012; Reichstein et al., 2013) could affect vernalization variation, especially in regions with low VDD (e.g., southern and continental TE). This in turn would reshape the sensitivity of winter wheat yields to variations of VDD in such regions. However, whether and when projected further climate warming will trigger a failure of vernalization (i.e., being unable to meet the minimum vernalization requirement) for current winter wheat cropping systems, particularly in the southern limit of winter wheat planting areas of TE, remains poorly understood. From gene perspective, vernalization process is dominantly regulated by the interactions between genes Vrn 1 and Vrn 2 (Yan et al., 2004). Previous study reported that loss of function of Vrn 2, whether by natural mutations or domestication, could result in reduction in vernalization requirement for winter wheat (Yan et al., 2004). Nevertheless, there still remains a big challenge for building an adaptive support system for sustainable agriculture, given that the breeding and adopting of new varieties is not likely to keep pace with rapid climate warming (Challinor et al., 2016). A northward shift of planting areas for vernalization-sensitive crops (e.g., winter wheat) is reasonably expected in response to the predicted future climate warming in TE (Olesen et al., 2007; IPCC, 2013). However, such shifts may not be sufficiently rapid to compensate for the negative effects of climate change on wheat yields, since virgin land for agriculture expansion is limited. Instead, the increasing γ_VDD in such regions could inevitably affect current food production. Note that the effects of VDD changes on winter wheat yields may have already been countered by local adaptation measures of individual farmers, such as changing sowing dates or introducing new cultivars. Any such gradual changes would not be captured by the statistical models using detrended data (for detail see section Materials and Methods). Thus, the interannual sensitivity of winter wheat yields to VDD variations reported here can be viewed as a “baseline” of the impacts of warming-mediated VDD variations on winter wheat yields.

Undoubtedly, the responses and adaptation of winter wheat systems to historical and future warming-induced changes in VDD also depend on human management practices, particularly the availability and levels of applied water and adaptive planting schemes (Wang et al., 2015). Although irrigation is seldom applied to winter wheat in Europe, it is locally used to avoid drought stress throughout the whole growing season (Rodríguez Díaz et al., 2007; Zhao et al., 2015), particularly in water-limited regions (see Figure 3A in Zhao et al., 2015). In addition, the possible difference in cultivars planted in rain-fed and irrigated croplands of TE could also bias the sensitivity of winter wheat yields to changes in VDD. However, the lack of sufficient observations of planted cultivars during recent decades across TE does not allow us to estimate the sensitivity of yields for different winter wheat cultivars to warming-induced VDD variations. Thus, to what extent reactions in such human management can modify the sensitivity of winter wheat yields to VDD changes remains unclear. Caution should be also exercised on potential biases in our conclusions that may originate from error propagation from different datasets and data processing (e.g., rasterizing of NUTS 2 crop yields). More accurate and detailed datasets for both crop yields and human management are urgently needed to constrain such uncertainty.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.