Predictions of future grazing season length for European dairy, beef and sheep farms based on regression with bioclimatic variables

Summary

Grazing season length (GSL) on grassland farms with ruminant production systems can influence farm economics, livestock disease transmission, environmental impact, milk and meat quality, and consumer choice. Bioclimatic variables are biologically meaningful climate variables that may enable predictions of the impact of future climate change on GSL on European farms. The present study investigated the spatial relationship between current GSL (months) measured by EUROSTAT on dairy, beef and sheep farms in 706, 774 and 878 regions, respectively, and bioclimatic variables. A stepwise multiple regression model revealed a highly significant association between observed GSL and bioclimatic variables across Europe. Mean GSL was positively associated with the mean temperature of the coldest quarter and isothermality, and negatively associated with precipitation in the wettest month. Extrapolating these relationships to future climate change scenarios, most European countries were predicted to have a net increase in GSL with the increase being largest (up to 2·5 months) in the north-east of Europe. However, there were also predictions of increased variability between regions and decreases in GSL of up to 1·5 months in some areas such as the west of France, the south-west of Norway and the west coast of Britain. The study quantified and mapped the potential impact of climate change on GSL for dairy, beef and sheep farms across Europe.

INTRODUCTION

There is currently a scientific consensus that anthropogenically induced climate change is occurring (Oreskes 2004; Cook et al. 2013). This climate change presents many potential future challenges to agriculture worldwide, not least changes to grazing season length (GSL) and to disease on ruminant livestock farms (Gale et al. 2009; Morgan & Wall 2009). Grazing season length is an important parameter when defining ruminant production systems. A long grazing season can increase the annual proportion of grazed grass in ruminant diets, which can reduce feed monetary costs and can thereby increase farm profitability in many ruminant production systems (Peeters 2009; Finneran et al. 2012). Dillon et al. (2005) showed that the cost of milk production across eight different countries (Australia, Denmark, Germany, Ireland, Netherlands, New Zealand, USA and the UK) was negatively related to the proportion of grazed grass in the diet.

Grazing season length can also change the environmental impact pathways of ruminant production systems. Webb et al. (2005) found that increasing GSL by 1 month resulted in 7–9% lower ammonia (NH₃) emissions from cattle farms in England and Wales, but that most of this conserved nitrogen (N) was lost in increased nitrate (NO₃ ⁻) leaching. Grazing season length is also important for product marketing as grazing of cattle and sheep is generally positively perceived by consumers when compared with indoor feeding of cattle and sheep (Font i Furnols et al. 2011).

Grazing season length is an important variable in the transmission of helminth parasites such as nematodes (roundworms) and trematodes (fluke) in grazing livestock, as these parasites are generally only transmitted when the livestock are grazing (Bowman et al. 2003). Therefore, potential future changes to GSL in response to climate change need to be considered when evaluating the potential impacts of climate change on helminth parasite transmission in grazing livestock (Morgan & Wall 2009). An international European research project (http://www.gloworm.eu) has been investigating the potential impact of climate change on helminth parasites of grazing livestock (Vercruysse et al. 2014), but future predictions of GSLs are required for these investigations.

The scientific literature to date contains no quantified predictions of future GSLs across Europe. It might be assumed that GSL will respond directly to changes in grass growth and therefore be predicted through grass growth modelling. However, even though farmers control the GSL, they generally do not measure grass growth and they can be influenced by other factors such as land trafficability after heavy rainfall or animal welfare in cold and hot weather (Ekesbo 2009; Creighton et al. 2011). Furthermore, predicting GSL based on grass growth modelling is problematic because grass growth models can rely on farm management data inputs (Hurtado-Uria et al. 2013). An alternative approach is to treat GSL as a management response to climate, investigate current relationships between climate variables and GSL and potentially extrapolate such relationships to climate change scenarios.

Bioclimatic variables are biologically meaningful variables that are derived from monthly, seasonal and annual combinations of temperature and precipitation data (Beaumont et al. 2005; Hijmans et al. 2005). These bioclimatic variables can be limiting or influential to species biology and are therefore typically used for ecological niche modelling to investigate the effect of climate on species distribution or phenology (Pearson & Dawson 2003). However, due to their biological implications, they may also be useful for predicting current and future farm management such as GSL on grazing livestock farms.

The objectives of the present study were to investigate the spatial relationship between current regional GSL and bioclimatic variables for European dairy, beef and sheep farms and, if a significant relationship exists, use it to predict potential future changes to GSL under the most recent IPCC AR5 (Intergovernmental Panel on Climate Change, Fifth Assessment Report) climate change projections.

MATERIALS AND METHODS

Geographical regions

The geographical area of this study included 32 European countries that provided GSL data to EUROSTAT (the statistical office of the European Union): Austria, Belgium, Bulgaria, Croatia, Cyprus, Czech Republic, Denmark, Estonia, France, Finland, Germany, Greece, Hungary, Iceland, Ireland, Italy, Latvia, Lithuania, Luxemburg, Malta, Montenegro, Netherlands, Norway, Poland, Portugal, Romania, Slovakia, Slovenia, Spain, Sweden, Switzerland and the UK. For all results, the observational units were geographical regions identified in the ‘nomenclature of units for territorial statistics’ classification system (NUTS; European Commission 2011). This is a geocode standard for referencing the sub-divisions of European countries for statistical purposes. Within the NUTS system, the highest level of resolution is NUTS level 3 and the lowest is NUTS level 1. The present study used NUTS level 3 regions for all countries except for Germany (NUTS level 2 regions due to availability).

Farm type

Within each region, three specialist grazing livestock farm types defined by the European Commission (2009a ) were identified: (i) specialist dairying (dairy farms), (ii) specialist cattle rearing and/or fattening (beef farms) and (iii) sheep, goats and other grazing livestock (sheep farms). These farm types were based on the proportional contribution from each livestock type (dairy cows, other cattle or sheep/goats) to the economic standard output of each farm (total monetary output at farm-gate price).

Observed current grazing season length

Current GSL data for each region were sourced from the results of the EUROSTAT Survey on Agricultural Production Methods (SAPM). The SAPM was initially conducted by national statistical bodies in each participating country and then collated by EUROSTAT. The reference year for the questions asked in the SAPM was 2010 in all countries except Spain and Portugal (reference year 2009). It was collected as part of an agricultural census (all known farms) in Austria, Bulgaria, Czech Republic, Estonia, France, Italy, Lithuania, Luxemburg, Malta, Montenegro, Netherlands, Portugal, Romania and Slovakia, whereas it was collected from a sample survey in the following countries with sample size (proportion of all farms) in parentheses: Belgium (0·03), Croatia (0·08), Cyprus (0·18), Denmark (0·34), Finland (0·22), Germany (0·26), Greece (0·06), Hungary (0·06), Ireland (0·19), Latvia (0·23), Norway (0·20), Poland (0·12), Slovenia (0·11), Spain (0·07), Sweden (0·10), Switzerland (0·26) and the UK (0·10). Full details of the SAPM methodology are described in European Commission regulation 1200/2009 (European Commission 2009b ).

Within the SAPM questionnaires, GSL was defined as ‘the number of months for which animals have been grazing during the reference year’ and was asked in two separate questions for land owned, rented or otherwise allocated to the farm (farmland) and land used by the farm but not allocated to it (commonage). The results were collated by EUROSTAT and provided in the following categories for farms that grazed their livestock: (i) 1–2 months, (ii) 3–4 months, (iii) 5–6 months, (iv) 7–9 months and (v) 10 months or more. Each category contained the number of farms, the grazed land areas, the numbers of the different livestock types (including age classes) and the total grazing livestock units (LU; one LU = feed requirement equivalent of one standard dairy cow; European Commission 2009b ). From the above SAPM categorical results provided by EUROSTAT, a weighted mean regional GSL was calculated from the number of LU in each GSL category using the central GSL within each category (1·5, 2·5, 3·5, 5·5, 8·0 and 11·0 months for categories (i)–(v) above, respectively). Before access to the data were granted, EUROSTAT used the following data screening in order to preserve confidentiality and remove dominance effects: any farm that accounted for at least 0·85 of a value for one region was excluded, any regional value that was calculated from less than five farms was excluded and all final regional values were rounded to the closest multiplier of 10. Weighted means calculated from fewer than 30 farms were considered unreliable and also excluded. In total, there were 706, 774 and 878 regions for dairy, beef and sheep farms, respectively. The weighted mean GSL for each geographical region and farm type was then calculated using the following equation:

$$GSL = \sum {\left[ \matrix{\left( {\displaystyle{{{F_{(n)}}} \over {\sum {F_{(n1,\; n2,...,n5)}}}} \times N} \right) \cr + \; \left( {\displaystyle{{{C_{(n)}}} \over {\sum {C_{(n1,\; n2,...,n5)}}}} \times N} \right)\; \hfill \cr + \; \left( {\displaystyle{{{B_{(n)}}} \over {\sum {B_{(n1,\; n2,...,n5)}}}} \times N} \right) \hfill} \right]_{n1,n2,....n5}} $$

where F was the number of LU grazing only farmland, C was the number of LU grazing only commonage, B was the number of LU grazing both farmland and commonage, N was the mean number of months in each GSL category and n was the nth mean GSL category.

Bioclimatic variables

Spatial data on bioclimatic variables (biologically meaningful combinations of temperature and precipitation) were provided by the WORLDCLIM dataset (Hijmans et al. 2005; http://www.worldclim.org/). In total, there were 19 bioclimatic variables available: annual mean temperature (°C), mean diurnal temperature range (°C), isothermality (°C; mean diurnal temperature range/annual temperature range), temperature seasonality (°C), maximum temperature of warmest month (°C), minimum temperature of coldest month (°C), annual temperature range (°C), mean temperature of wettest quarter (°C), mean temperature of driest quarter (°C), mean temperature of warmest quarter (°C), mean temperature of coldest quarter (°C), annual precipitation (mm), precipitation of wettest month (mm), precipitation of driest month (mm), precipitation seasonality (mm), precipitation of wettest quarter (mm), precipitation of driest quarter (mm), precipitation of warmest quarter (mm) and precipitation of coldest quarter (mm).

These data had a spatial resolution of 30 arc-seconds (c. 1 km²) and were provided for the current temporal baseline (interpolations of observed data between 1950 and 2000) and for future projections for the years 2050 (mean for 2041–2060) and 2070 (mean for 2061–2080) under the global climate model HadGEM2-ES for the four representative concentration pathways (2·6, 4·5, 6·0 and 8·5 W/m²) as predicted by phase five of the Coupled Model Intercomparison Project (CMIP5) that was used in the most recent (fifth) assessment report by the IPCC (Jones et al. 2011; Taylor et al. 2012).

In order to match the spatial resolution of the GSL data, mean regional bioclimatic variables were obtained for NUTS 3 (NUTS 2 for Germany) regions using the zonal statistics tool in ArcMap (10·1, ESRI, Redlands, CA, USA) overlain with polygons of the 2010 NUTS shapefiles available from EUROSTAT (2014).

Statistical analyses and mapping

The current relationship between observed GSL and bioclimatic variables was analysed in a multiple linear regression model with farm type included as a class variable. In total, there were 706, 774 and 878 regions for dairy, beef and sheep farms, respectively. Stepwise multiple linear regression was used through the GLMSELECT procedure in SAS (9·3, Cary, NC, USA) using the methodology advised by Cohen (2006). This methodology randomly selected 0·5 of the observations for training the model and the remaining 0·5 for validation testing. The ‘stepwise = validate’ statement specified that the stepwise process finish when adding or removing any variable resulted in an increase in the mean squared error in the validation data. All bioclimatic variables, their square functions and interactions were included in the analysis. The bioclimatic variables selected by the model were checked for multicollinearity with variance inflation factors and Pearson correlation coefficients using the PROC CORR and PROC REG procedures (respectively) as described by Ngo (2012).

Future predictions of GSL were calculated using the final regression equation from the above model and the predicted bioclimatic variables for the years 2050 and 2070 under the global climate model HadGEM2-ES for the four representative concentration pathways (2·6, 4·5, 6·0 and 8·5) predicted by CMIP5. Grazing season lengths were only predicted for regions that had a current observed GSL and were capped at the lower and upper limits of the observed GSL dataset (1·5 and 11·0 months, respectively). Within each country and farm type, the observed current GSL, predicted current GSL and all predicted future GSLs were compared in an analysis of variance (ANOVA) using the PROC GLM statement in SAS (9·3, Cary, NC, USA) with the Bonferoni adjustment for pairwise comparisons. Changes in predicted future changes to GSL were calculated relative to the predicted current GSL and results were mapped using ArcMap (10·1, ESRI, Redlands, CA, USA).

RESULTS

Relationship between observed current grazing season length and bioclimatic variables

Observed current GSL was highly significantly associated with bioclimatic variables (final R ² = 0·66, final mean square error = 1·248 months, P < 0·001, Table 1). The mean temperature of the coldest quarter was the variable most closely associated with GSL, accounting for approximately half of the variation (R ² = 0·52), being the first step selected in the model (P < 0·001; Table 1) and being positively associated with GSL. Following this, the class variable of farm type was the next variable selected in the model with GSL being a mean 1·24 and 0·60 months less on dairy and beef farms than on sheep farms, respectively (P < 0·001; Table 1). Isothermality and the precipitation of the wettest month were the next variables kept in the model being positively and negatively associated with GSL, respectively (P < 0·001; Table 1). In the next step of the model, the precipitation of the wettest month had an interaction with isothermality (P < 0·001; Table 1). This interaction resulted in the negative association between GSL and precipitation in the wettest month being greatest at low isothermality (Fig. 1(a)) and the positive association between GSL and isothermality being greatest at high precipitation of the wettest month (Fig. 1(b)). The final step in the model included the square of isothermality being negatively associated with GSL (P < 0·001; Table 1). This final step resulted in a slight curvilinear relationship between isothermality and GSL at high values (Fig. 1(b), ○).

Fig. 1. Modelled effect (a) of precipitation of the wettest month (mm) on predicted grazing season length (GSL) for the 5th (○) and 95th (●) percentiles of isothermality (°C) and modelled effect (b) of isothermality on predicted GSL for the 5th (○) and 95th (●) percentiles of precipitation in the wettest month.

Table 1. Step results of a stepwise linear multiple regression model testing the prediction of regional GSL from regional bioclimatic variables across Europe

Predicted grazing season length

When the observed current GSL values were plotted against the predicted current GSL values (Figs 2(a), 2(c) and 2(e)), the intercepts were not significantly different to 0·0 (0·11; P = 0·339) and the slopes not significantly different to 1·0 (1·01; P = 0·345). However, plots of the residuals revealed a trend of increasing negative residual error with increasing GSL (Figs 2(b), 2(d) and 2(f)). As Fig. 3 shows, this was due to large variation in observed GSL between some neighbouring regions in the Mediterranean and Eastern Europe. However, the model did predict the broad trend of longer GSLs in the west and south of Europe in comparison to the east and north, as well as the shorter GSLs in the alpine mountains (Fig. 3).

Fig. 2. Regression between current (2010) observed and predicted grazing season length for dairy (a), beef (c) and sheep (e) farms (y = 1·01x − 0·11, R ² = 0·66, P = 0·001 for all). Corresponding residuals are shown for dairy (b), beef (d) and sheep (f) farms.

Fig. 3. Maps of grazing season length (GSL; months) in 2010 as observed for dairy farms (a), predicted for dairy farms (b), observed for beef farms (c), predicted for beef farms (d), observed for sheep farms (e) and predicted for sheep farms (f). Observed and predicted were based on EUROSTAT results and regression with bioclimatic variables, respectively. (Colour online)

For each country, differences between observed current GSL, predicted current GSL and all predicted future GSLs for dairy, beef and sheep farms are shown in Tables 2–4, respectively. For dairy farms, the predicted current GSL was longer than the observed current GSL for Belgium, Estonia, Germany, Hungary and the Netherlands but shorter for Bulgaria, France, Latvia and Lithuania (P < 0·05; Table 2). For beef farms, the predicted current GSL was longer than the observed for Belgium, Ireland, Netherlands, Poland and Romania, but shorter for Bulgaria, Czech Republic, Denmark, Estonia, Finland, France, Hungary, Latvia and the UK (P < 0·05; Table 3). For sheep farms, predicted current GSL was longer than observed for Ireland, Lithuania, Netherlands, Norway, Poland, Romania, Spain and Switzerland but shorter for Croatia, Czech Republic, Denmark, Estonia, France, Germany and the UK (P < 0·05; Table 4).

Table 2. Analysis of variance results for mean dairy farm GSL (months) between that observed in 2010 and predicted from bioclimatic variables for 2010 and the four representative concentration pathways (2·5, 4·5, 6·0 and 8·5 W/m²) of the HadGEM2-ES climate change scenarios in 2050 and 2070

The number of regions in each country (n) and standard error (s.e.) are also presented.

Table 3. Analysis of variance results for mean beef farm grazing season length (months) between that observed in 2010 and predicted from bioclimatic variables for 2010 and the four representative concentration pathways (2·5, 4·5, 6·0 and 8·5 W/m²) of the HadGEM2-ES climate change scenarios in 2050 and 2070

The number of regions in each country (n) and standard error (s.e.) are also presented.

Table 4. Analysis of variance results for mean sheep farm grazing season length (months) between that observed in 2010 and predicted from bioclimatic variables for 2010 and the four representative concentration pathways (2·5, 4·5, 6·0 and 8·5 W/m²) of the HadGEM2-ES climate change scenarios in 2050 and 2070

The number of regions in each country (n) and standard error (s.e.) are also presented.

Relative to predicted current GSL, predicted future GSL for 2050 and 2070 increased under climate change scenarios for most countries for all three farm types (P < 0·05, Tables 2–4). However, there were no differences in GSL (P > 0·05) between predicted current and future scenarios for dairy farms for Croatia, Greece, Norway, Spain and the UK; beef farms for Croatia, Greece, Norway, Portugal and Spain or sheep farms for Belgium, Cyprus, Greece, Iceland, Malta and Portugal (Tables 2–4). Relative to predicted current GSL, predicted future GSL did not decrease in any country under any climate change scenario for 2050 and 2070 (Tables 2–4). However, comparisons could not be made in Cyprus, Iceland, Luxemburg, Macedonia, Malta or Montenegro for dairy and beef farms or in Cyprus, Luxemburg and Montenegro for sheep farms due to there being only one region in each of these country–farm-type combinations (Tables 2–4).

For each region, the mapped predicted future changes in GSL for representative concentration pathways 2·5 and 8·5 W/m² and the years 2050 and 2070 are shown in Fig. 4 (dairy farms), Fig. 5 (beef farms) and Fig. 6 (sheep farms). Both future increases and decreases in GSL were predicted within some countries such as France, Norway, Germany, Italy, Spain and the UK. For most regions, increases in GSL were predicted, with the increase being greatest in the Alps and the northeast of Europe, and greatest in 2070 under the highest representative concentration pathway (Figs 4–6(d)). However this latter scenario also predicted the greatest variation in predicted changes to GSL with decreases in GSL between 1·0 and 1·5 months in some areas such as the west of France, the south-west of Norway and the west coast of Britain (Figs 4–6(d)).

Fig. 4. Maps of predicted future changes to grazing season length (GSL; months) on dairy farms under the CMIP5 HadGEM2-ESclimate change scenarios for 2050 (a, b) and 2070 (c, d) under representative concentration pathways 2·5 (a, c) and 8·5 (b, d). (Colour online)

Fig. 5. Maps of predicted future changes to grazing season length (GSL; months) on beef farms under the CMIP5 HadGEM2-ES climate change scenarios for 2050 (a, b) and 2070 (c, d) under representative concentration pathways 2·5 (a, c) and 8·5 (b, d). (Colour online)

Fig. 6. Maps of predicted future changes to grazing season length (GSL; months) on sheep farms under the CMIP5 HadGEM2-ES climate change scenarios for 2050 (a, b) and 2070 (c, d) under representative concentration pathways 2·5 (a, c) and 8·5 (b, d). (Colour online)

DISCUSSION

Association between grazing season length and bioclimatic variables

The positive correlation between GSL and the mean temperature of the coldest quarter was a logical result. Temperature is one of the main determinants of grass growth, which is generally assumed to start only when the 10-day moving mean temperature is >3–5 °C, assuming sufficient soil moisture is present (Schapendonk et al. 1998). Therefore, for many regions across Europe, winter temperature is most likely to be the main limiting factor to grass growth. However, winter temperature may also result in livestock being housed to protect them from cold conditions. Most cattle and sheep have the ability to acclimatize to cold conditions, with lower critical temperatures of –13 °C being reported for cattle with winter (acclimatized) coats in calm, dry conditions (Ekesbo 2009). However, this lower critical temperature increases rapidly in wet or windy conditions, or when livestock have not acclimatized, and can be as high as 15 °C for cattle with wet summer coats in windy conditions (Ekesbo 2009). For these reasons, combined with the effect of the mean temperature of the coldest period described above, the negative association between GSL and precipitation in the wettest month, the positive association between GSL and isothermality, and the interaction between the latter two can generally be interpreted as livestock being housed in winter in cold, wet and/or changeable conditions.

It appears surprising that bioclimatic variables associated with high temperature or dry conditions were not associated with GSL. Previous studies have suggested that heat stress generally presents a greater risk for farmed ruminants than cold stress (Hemsworth et al. 1995) and the negative effects of low soil moisture on grass growth are well documented (Volaire et al. 2009). However, in the SAPM results, many regions in southern Europe which experienced hot and dry conditions had quite long observed GSLs. In such conditions, grazing managers may respond to the climatic stresses with lower stocking densities (particularly on sheep farms) or feed/water supplementation at pasture (particularly on dairy farms), rather than housing the livestock. Furthermore, it should be noted that the present study only used farms that had a GSL to measure (i.e. practiced grazing). Many farms in these regions (particularly large dairy farms) may house livestock all year round and have no GSL to measure. The negative effect of precipitation in the wettest month on GSL, rather than drought conditions, may also be due to land trafficability in many regions. For example, a survey of Irish dairy farmers found that c. 60% of them identified soil condition as the main limiting factor for GSL (Creighton et al. 2011) and modelling studies have previously identified wet soil as a limiting factor for GSL in the UK (Rounsevell et al. 1996; Parsons et al. 2001). This effect is most likely a preventative measure by farmers, as treading damage on wet soil can substantially reduce subsequent grass growth (Phelan et al. 2013).The linear regression methodology approach of the present study was useful for determining bioclimatic predictor variables for GSL. It avoided many of the assumptions of management (both inputs and responses) and complexity that would have been required if a grass growth model was used to predict GSL. However, it was applied over a broad geographical scale (all of Europe in the present study). At a more regional level, individual studies exploring that region's threshold levels of bioclimatic variables for GSL might also be useful in predicting that region's response to climate change. For example, the proportion of farmers that house livestock when precipitation or temperature reaches certain thresholds might be extrapolated to future climate change predictions for a specific region.

Bioclimatic variables are often used in ecological studies of wild species distributions (Beaumont et al. 2005). However, the highly significant association between bioclimatic variables and a farmer-controlled variable (GSL) in the present study shows that bioclimatic variables can be useful for other applications such as agricultural studies.

Limitations

One of the main limitations to the present study is the fact that observed GSL data were only available for 2010 (2009 for Spain and Portugal), whereas the baseline for bioclimatic variables used in the study was the period between 1950 and 2000. Global climate has already changed slightly within this period (Stocker et al. 2013) and this may have increased errors in the model. Furthermore, the fact that the GSL data were only available for 1 year may have increased the likelihood of unusual weather in that year (in comparison to the 50-year mean), influencing the observed GSL. However, the present study area had large spatial variation in bioclimatic variables in comparison with any recent temporal or between-year variation. For example, the mean linear decadal rate of increase in winter temperature in Europe between 1901 and 2000 has been reported as 0·1 ± 0·07 °C (Luterbacher et al. 2004). In contrast, winter temperature (mean temperature of the coldest quarter) in the spatial dataset used in the present study varied by 26 °C (from –13·2 °C in Lappi, Finland, to 12·8 °C in Melila, Spain).

While the association between GSL and bioclimatic variables was highly significant, there were variations between some neighbouring regions that were not predicted. This may have been due to various other influences such as localized production systems/producers, land type, land availability, etc. There did appear to be greater variation between neighbouring regions for GSL in Mediterranean and Eastern European regions, potentially suggesting greater variation where summer droughts are experienced, although this was not apparent in the present study's results.

Another limitation of the available data is that GSL definition was quite broad. It was defined as ‘the number of months for which animals have been grazing during the reference year’ by EUROSTAT (European Commission 2009b ), but may have varied slightly when collected by the national statistical bodies in each country (due to translation or terminology differences). This could cause definition issues where livestock are fed indoors, but have access to a small pasture for health and welfare. Some farmers may define this as grazing and others as not grazing. In such situations, the proportion of grazed grass in the livestock's diet may vary considerably between farms of similar recorded GSL. Furthermore, the fact that GSL data were calculated from weighted means resulted in GSL ranging from 1·5 to 11 months rather than from >0 to 12 months, which reduced the range of values in GSL available for regression. However, the above limitations would most likely have increased error, rather than bias, in the multiple regression model and resulted in lower R ². Despite this, a significant correlation, high R ² and low mean square error was output by the model, most likely due to the overwriting effect of the size of the dataset and the large spatial variation in bioclimatic variables and GSL. Finally, the partitioning of the data into separate training and validation subsets (0·5 for each) reduced the likelihood of model overfitting and increased the suitability of applying the regression results to other (future) scenarios (Picard & Cook 1984).

Predicted future grazing season length

Aside from any assumptions made in the HadGEM2-ES future climate change projection models (Jones et al. 2011), the present study assumes that the current relationship between bioclimatic variables, farm type and GSL will apply in the future. However, future technological adaptations, management practices or economic changes may increase or decrease the importance of bioclimatic variables in determining GSL. For example, the increasing use of robotic milking systems may reduce the amount of grazing by dairy cows due to perceived difficulties of including such technology in grazing systems, regardless of climate change (Kristensen et al. 2010). Farm expansion and lack of availability of grazing land for that expansion may also result in greater use of feeds other than grazed grass, regardless of changes to the climate.

The predictions of the present study also assume that farm management will respond to bioclimatic variables over time in the same manner that they currently respond to them over space (between regions). This assumption has been used elsewhere, most notably in the climate change analogue tool developed by the Consultative Group on International Agricultural Research (CGIAR) and the Research Program on Climate Change, Agriculture and Food Security (CCAFS) (Ramirez-Villegas et al. 2011). However, farm management may be slower to respond to climate change than anticipated. For example, in European arable farming systems, drilling, tilling and harvesting dates have already been found to be correlated with temperature over time (mean r = –0·53; Menzel et al. 2006), but the response of these activities was found to be much smaller (mean advance of 0·4 days per decade) when compared with the observed phenological changes in plants (mean advance of flowering and leaf unfolding by 2·5 days per decade; Menzel et al. 2006). This might be due to cultural tradition or a lack of perception of climate change until it is large enough to be recognized by the farmer. Therefore, there may be a lag time in the response of farm management to future changes to the climate.

Despite the above assumptions, the significant association between current GSL and bioclimatic variables found in the present study and the scientific consensus on likely future changes to these bioclimatic variables do indicate that GSL will be altered by climate change in Europe, with most regions predicted to have increased future GSL. However, it is notable that GSL was not predicted to increase in some regions and may decrease by up to 1·5 months in a minority of regions due to greater winter precipitation and lower isothermality. A previous study in the UK that modelled land suitability classes in response to climate change also found that the positive effect of increasing temperature on the distribution of grassland suitability may be completely offset by reduced land trafficability due to increasing precipitation (Rounsevell et al. 1996).

Increased GSL on grazing livestock farms can reduce feed costs (Dillon et al. 2005; Läpple et al. 2012). For example, Läpple et al. (2012) analysed national statistical farm data in Ireland and found that the cost of milk production was negatively related to GSL, with a reduction of €0·016/l for each 10-day increase in GSL. Increased GSL may also result in changes to the environmental impacts of grazing with a trend towards lower greenhouse gas (GHG) emissions, but increased nitrate leaching (Webb et al. 2005; Murphy et al. 2013). However, increased GSL may also result in increased exposure to helminth parasites, which may be associated with reduced milk yields (Charlier et al. 2005; Verschave et al. 2014), and decreased production efficiency, counterbalancing the benefits from reductions in GHG emissions (Kenyon et al. 2013). Furthermore, changes to timing of the grazing season may result in changes to the timing of helminth parasite burdens for which current farm management practices are unprepared. For example, it may result in increased exposure to liver fluke (Fasciola hepatica L.) in late autumn or earlier exposure to gastrointestinal nematodes in spring (Van Dijk et al. 2010). However, these effects may also depend on the parasites’ response to climate change (Rose et al. 2015). A holistic approach to predicting changes to parasite epidemiology under climate change is therefore needed, which integrates effects on parasite biology and those on relevant management variables such as GSL.

CONCLUSION

Current observed GSL on European dairy, beef and sheep farms was highly associated with bioclimatic variables. Grazing season length was positively correlated with the mean temperature of the coldest quarter and isothermality and negatively correlated with precipitation in the wettest month. Extrapolating these correlations to future climate change scenarios, most European countries were predicted to have a net increase in GSL, with the increase being largest in the north-east of Europe. However, there were also predictions of increased variability between regions and decreases in GSL between 1·0 and 1·5 months in some areas such as the west of France, the south-west of Norway and the west coast of Britain. The present study quantified and mapped predictions of the potential impact of climate change on future GSL across Europe for the first time and may be used to explore climate change implications for farm feed costs, livestock disease and the environmental impacts of farming.