Modeling heat stress under organic dairy farming conditions in warm temperate climates within the Mediterranean basin

We studied the effect of heat stress on milk quality in Spanish organic dairy farms using published milk productivity equations. We collected data from 23 weather stations and 14,424 milk test-days for milk yield and milk fat and protein content for the period July 2011 to June 2013. As an indicator of heat stress, we used the maximum daily temperature–humidity index (THI) from 2 days before the milk test date. We fitted the data using hierarchical regression models stratified by farm, cow parity and monthly test-day milk records. The effect of THI was deemed low on biological costs through milk yield. However, the known negative relationship between milk yield and milk quality (protein and fat content) became even steeper when the THI increased, suggesting a significant negative correlation between heat stress and milk quality. Therefore, although the milk yield of cows in the organic farming systems analyzed appeared resilient to heat stress conditions, milk quality, a major selling point for organic dairy products, was negatively affected. The model presented here could be used to predict the potential impacts of different climate change scenarios on dairy farming, and to delineate adaptation strategies within organic systems.


Introduction
During the past few decades, there has been a dramatic intensification of animal production in response to the rising global demand for meat and milk. In dairy cattle systems, intensification https://doi.org/10.1007/s10584-020-02818-y Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10584-020-02818-y) contains supplementary material, which is available to authorized users. of production can decrease the thermoregulatory capacity of the animals, and as a result heat stress has become a major challenge facing the dairy sector today (Polsky and von Keyserlingk 2017), with impacts on cattle health, welfare, and productivity (Özkan et al. 2016;Vitali et al. 2015).
The issue of heat stress is of particular importance in the context of climate change, which is expected to cause large increases in temperature in regions such as the Mediterranean basin (e.g., Segnalini et al. 2013). Studies in farm animals have shown that across most of this area, animals are at increasing risk of suffering heat stress (Bouraoui et al. 2002;Segnalini et al. 2011) as the frequency and intensity of heat waves continue to rise (Pasqui and Giuseppe 2019). These changes are reflected in an increase across Europe in the number of days when the temperature-humidity index (THI) exceeds the thermal comfort threshold for livestock, and necessitate the development of adaptation strategies to minimize negative effects on farm animals, particularly in the Mediterranean basin (Segnalini et al. 2013).
Ongoing climate change can be expected to have different effects on different parts of the dairy sector, since there is considerable diversity in how cows are housed and managed in Europe (EFSA 2015). Mikovits et al. (2019) showed that, in Central Europe, the indoor climate of livestock confinement systems (often referred to as industrial or intensive farming systems) are more sensitive to heat stress than systems where animals are kept outdoors. This sensitivity is caused primarily by high indoor temperatures resulting from the heat load produced by the animals, with such temperatures often higher than those outdoors even in the summer. In addition, animals in intensive systems are producing at a higher throughput level than those in more extensive systems, resulting in an elevated internal heat load and thermal stress (Polsky and von Keyserlingk 2017).
Numerous studies have been undertaken to establish thresholds for tolerance to heat stress in livestock on the basis of THI values (Segnalini et al. 2013). The heat stress threshold has been set at daily THI values of 68-72 Polsky and von Keyserlingk 2017), although lower values have been found for temperate areas (Lambertz et al. 2014). Previous work conducted in regions close to and in the same years as the present study (Carabaño et al. 2014;Cerqueira et al. 2016) showed that under intensive livestock production, THI had a detrimental effect on milk production parameters. For Holstein Friesian cows in southern Spain, Carabaño et al. (2014) used modeling to describe heat stress response for production traits. At 25 and 34°C the estimated rates of decay for milk, fat, and protein yield were 36 and 170, 3.8 and 3.0, and 3.9 and 8.2 g/day per degree Celsius, respectively. In non-organic commercial farms in northern Portugal, Cerqueira et al. (2016) focused on physiological indicators and milk production losses and reported that the step of the THI of 72 over 78 meant an increase of 1.3°C and 37.3 breaths per minute in dairy cows. THI above 78 corresponded to a loss of milk production of 1.8 kg/cow/day (with THI surpassing the value of 72 on 14% of the days). In addition to effects on milk yield, there is strong evidence for an influence of heat stress on milk composition, with studies finding that increasing THI was associated with a decline in total protein and total fat content (Bouraoui et al. 2002;Bernabucci et al. 2015;Hammami et al. 2015;Hill and Wall 2015). For example, Lambertz et al. (2014) found a significant decrease in daily fat content when THI levels increased from 60-65 to above 65. Compared to intensive systems, different heat stress scenarios could be expected in alternative systems such as organic dairy production systems which are characterized by high-roughage diets based on home-grown feedstuffs, access to outdoor pasture, and restrictive use of antibiotics (EC 834/2007). Specifically, cows raised organically may be less susceptible to heat stress due to their less intensive management. In the context of increased awareness regarding climate change and environmental issues, organic farming is often suggested as a solution to the negative environmental effects of current food production (Smith et al. 2019). However, while thermal stress in dairy cows has been extensively studied in high-yielding cows (Gauly et al. 2013), there has been very little discussion in the scientific literature regarding heat stress under extensively managed cattle systems, and the extent to which heat stress affects the ability of cows to adapt to an organic production environment has not been well explored. One challenge for studies of organic farming systems is that such systems include a range of production environments, varying from very extensive small-scale farming to large intensive systems more similar to conventional production (Wallenbeck et al. 2018;Blanco-Penedo et al. 2019).
There is, as yet, little discussion or research on the challenges that organic livestock production systems face at the production unit (farm) scale in the Mediterranean basin. In the light of this gap in current knowledge, the objective of this study was to determine the effect of heat stress on milk production and quality by studying the effects of THI on milk yield and milk fat and protein content in organic dairy farms in Spain.

Scenario and farm information
The study collected data from 26 organic dairy farms located in seven regions in northern (Asturias, Basque Country, Cantabria, Catalonia and Galicia) and central (Madrid) Spain (Fig. 1). These regions correspond to areas of high dairy cattle concentration. The climatology of these areas is quite diverse and has been classified according to the Köppen-Geiger climate classification (http://koeppen-geiger.vu-wien.ac.at/) as warm temperate, but with some diversity within this classification, with precipitation ranging from no dry season (north) to dry summer (central), and temperatures ranging from warm to hot summer.
The selection of the farms was based on the criteria of the European Commission Seventh Framework Programme (FP7) project IMPRO (http://www.impro-dairy.eu/) to ensure that the sample was representative of organic dairy production in Spain. The inclusion criteria were as follows: (1) time under organic conversion (a minimum of 1 year), (2) availability of official milk recording scheme records, (3) intention to continue in organic production for at least 5 years, and (4) a herd size typical for the country. The surveyed farms comprised approximately 33% of the total official census of organic dairy farms in Spain. The 26 selected farms contained in total 1486 dairy lactating cows from one to 12 parities (average parity 3.4).
A data set was retrieved from the Spanish Holstein Friesian National Confederation (CONAFE) containing 14,424 test-day records, including the parameters of milk yield (kg/ day), fat and protein (g/100 ml milk), milk somatic cell count (SCC), calving dates, number of calves, and days in milk (DIM), from July 2011 to June 2013. During this period, monthly measurements of lactation were taken for each of the same cows on herd records.
Weather data were obtained from 23 weather stations located in the proximity of the farms within a distance of between 1.9 km and 61 km (median distance = 14.3 km, mean distance = 16 (±11.2) km). Each farm was associated with the nearest weather station. Meteorological data from weather stations were acquired from the State Meteorological Agency (AEMET) in Spain. Meteorological data were used to calculate the THI, assumed to be the best available indicator of heat stress in dairy cows (Berman 2005). Following previous studies ; Polsky and von Keyserlingk 2017) THI thresholds for heat stress of 65 and 72 were focused on in this study. The lower of these two values is a little below the lower level (68) estimated by these authors, in order to account for studies suggesting that a THI of 65 and above can affect milk quality parameters (Lambertz et al. 2014).
General statistics from the data used in the analysis are presented in Table 1. Basic farm structural information was obtained from farmers through face-to-face questionnaires. The average herd size in the sample was 50.1 cows (range: 11.5-312.2) with Holstein Friesian as the predominant breed (pure and crosses) in 85% of the farms, and 15% of other predominant breeds (Fleckvieh or Brown Swiss). Half of the herds included crossbreeds such as Fleckvieh (35.7%), Holstein Friesian (26.6%), Brown Swiss (27%), and Jersey (10.7%). Fig. 1 Map of the sampled farms and nearest weather stations. The red dots indicate the 23 weather stations from which daily temperature and relative humidity data were acquired, which were utilized to calculate the maximum and minimum temperature-humidity index Unfortunately, the dataset did not facilitate the linking of data on breeds to the milk recording data for individual cows, meaning that it was not possible to differentiate the milk records of Holstein Friesian from those of crossbreed individuals-although in all cases crossbreeds had Holstein Friesian genes. As a result, breed could not be introduced as a parameter in the model. However, the modeling considered individual performance variability as it followed each individual animal's parity.
The average annual milk yield per cow was 5658 kg (range: 3400-8371) on 305-day lactation. The average daily milk yield of the herds was 26.9 kg (range: 4.2-65.1). In the study herds, all cows had access to pasture during the whole year, comprising 330 days (range 200-365) in total. Cows had an average of 13.4 (range: 5-22) grazing hours per day. Although characterized at a general level for farms, grazing activity data did not capture daily variations. As grazing varies greatly from day to day depending on weather conditions, it is likely to correlate with THI. Therefore, nesting the model parameters by farm was used as a proxy to represent between-farm differences in grazing and housing-related factors in subsequent modeling.
Milk/concentrate (kg/kg: milk yield in kg per kg of concentrate) was 5.27 (range: 2.2-28.3), and average concentrate per hectare (kg/ha) per year was 0.32 (range: 0-1.28). Housing across the sampled farms was diverse, with 56% of farms having outdoor climate buildings. Buildings on all farms were naturally ventilated. Ventilators were placed in the milking parlor if present. According to organic standards (checked by audits under the organic certification scheme), animals should have sufficient shade under hedges and trees and adequate shelters, in addition to guaranteed free access to water.

Data management
The THI was used to characterize the climate conditions on the farms. THI was calculated according to the National Research Council guidelines (1971), assumed to be the best available environmental indicator of heat stress in livestock (Berman 2005). THI is then described as: where T is the daily maximum air temperature in degrees Centigrade, and RH is the relative humidity at the hottest moment of the day (reached at 1300 hours). All measurements were taken 2 days before milking, following Bertocchi et al. (2014) and Saizi et al. (2019). All temperatures were corrected according to the difference in altitude between each farm and its corresponding weather station, decreasing 6.49°C per 1000 m increase (ICAO 1993). The maximum daily THI was chosen (see Fig. 2) because it represents the highest risk conditions for heat stress experienced by the cows. It also fits well with results by Brügemann et al. (2012), who indicated higher sensitivity of test-day milk yield to extreme values of maximum THI compared with daily average THI values.

Biological and statistical models
The effect of THI on milk yield and milk quality was assessed with two main sets of analysis within the Bayesian framework. (In Bayesian reasoning, uncertainty is attributed to the parameters, while the sampled data are regarded as a fixed quantity once collected. All parameters are modeled by distributions Dohoo et al. 2010). A Bayesian approach was chosen because of its underlying assumption that all parameters are random quantities, the flexibility to model nonlinear relationships, and the fact that the model's uncertainty permeates to all modeled layers. In Bayesian analysis, a parameter is characterized by an entire distribution of values instead of one fixed value as in classical frequentist analysis. The significance of a variable's effect with respect to the model is then determined by evaluating the a posteriori distribution of the parameter of interest, in magnitude and spread (typically, a parameter should not overlap 0 to have an effect in the model). Apart from the parameter magnitude, the overall model performance was measured using the deviance information criterion (DIC) (Spiegelhalter et al. 2002), a hierarchical modeling generalization of the Akaike information criterion (AIC). This is particularly useful in Bayesian model selection problems to compare models of increasing complexity.

Effects of THI on milk yield
We developed a Bayesian hierarchical model that builds upon existing knowledge of milk productivity, models given days in milk (DIM), and accounts for variability and uncertainty at three different hierarchical levels: farm (f), individual cow parity (p), and monthly test-day milk records (m). We use the milk productivity curve proposed by Wood (1967), because it has been shown to be the most parsimonious predictor for milk yield (Koçak and Ekiz 2008). We stratified variability by parity to fit a separate lactation model to each cow on each parity. Milk yield (Y m, p, f ) is then assumed to follow a normal probability distribution with mean μ m, p, f and variance σ so that where γ f sets the initial value independently by farm, β p,f controls the initial rate of increase in yield after the first days and the maximum production independently by parity, and ρ sets the rate of decay in yield as a continuous function of DIM m,p,f . We tested and included for the effect of parity number in the parameter β p,f as.
We then tested for the effect of THI on the general productivity level, so that The THI recorded 2 days before the milking date was used for each milk record (Bertocchi et al. 2014;Saizi et al. 2019). We tested THI as a class variable with different critical values, but none improved the model.
Effects of THI on milk fat and protein Preliminary analysis of the relationship between THI and the amounts of fat-and protein-corrected milk (FPCM) revealed no relationship (see Fig.  2c). We exploited the natural negative correlation between milk yield and fat and protein content (Neitzel 1967) and tested for the effects of THI on the respective regression parameters. In this way, we did not segregate the data by parity but rather included parity number as a covariate of each measurement. The milk parameter Z m,f (fat or protein content) then followed a normal distribution with mean μ m,f and variance σ so that where Y m,f is milk yield, and A and B are farm-level intercept parameters that follow a normal distribution with mean μ A and μ B , and variance σ A and σ B , respectively.
In all models, all parameters not otherwise specified were defined by uninformative a priori probability distributions centered at zero when applicable. A Bayesian approach via JAGS was used to obtain solutions for all the unknowns and statistics for Bayesian model assessment (Plummer 2012) until convergence, and sampling parameters over 10,000 extra sampling iterations.

Results
Across all farms, milk yield, fat, and protein yield averaged (± SD) 21.14 (±8.3 kg), 3.7 (± 0.87 kg), and 3.2 (± 0.40 kg), respectively, over the whole study period. During the research period, THI exceeded the threshold of 65 on 31.3% of the days, and THI values exceeded 72 on 10.6% of the total days (Fig. 2). Dairy cows were exposed to heat stress defined by the THI thresholds during summer months on 43% of test-days for the lower THI threshold and 40% for the higher threshold. However, these thresholds were also exceeded in spring (15% and 6% of test-days), autumn (50% and 44% of test-days), and even some days in winter (1% of all measurement above THI 65 occurred in February). Overall, the THI threshold of 65 was exceeded as early as February and as late as November, with the 72 threshold exceeded as early as May and as late as October.

Effects of heat stress on milk yield
No direct evidence was found for a negative effect of heat stress on milk productivity (adjusted milk yield means) in the tested organic dairy farms (i.e., the parameter estimate for THI overlapped 0). However, the addition of THI to the model improved the fit of the model to the data at model comparison [reduced deviance information criterion (DIC); Table 2]. The overfitting generated by a correlation between the initial yield value (γ f ) and the THI parameter effect (δ), which was estimated as positive, contributed to the model improvement. This is just an approximate estimation of the parameters, as the sampling algorithms of the model including THI did not converge (the algorithms struggled to find a parsimonious solution).
Regarding inter-and intra-farm variability, our results show that each farm has a unique milk productivity level-initial yield value and variability-and this is highly dependent on the cow's parity (individual animal response) (Fig. 3a) becoming the most important factor influencing milk yield (Fig. 3b).

Effects of heat stress on milk content
The most parsimonious model included both milk yield and the effect of THI on milk yield. In addition, the quadratic effect of parity also significantly improved the model fit (Table 2). Increasing THI steepened the slope of the known positive relationship between milk yield and milk quality (milk quality declined with milk yield) for all the farms. The regression analysis of  Wood: Lactation curve according to Wood (1967). β p,f : stratified parameter for initial rate of increment in yield. Parity: effect of parity number on β p,f . THI: effect of THI on the yield level the test-day records confirmed that both fat and protein content were positively correlated with milk yield, although this effect was smaller for fat than for protein (Figs. 3a and 4a). Thus, apart from the direct effect of yield on fat and protein, we also found an indirect effect of THI on this relationship, that is, the higher the THI, the greater the negative effect of yield on fat and protein.

Discussion
This study aimed to investigate the extent to which milk production and milk quality from cows on organic farms are affected by heat stress. Our results indicated no clear relationship between milk production in organic systems and THI; THI levels reached thresholds expected to cause heat stress without an apparent impact on milk yield. This discrepancy relative to previous studies undertaken in nearby areas (Carabaño et al. 2019) might suggest differences in heat stress tolerance between production systems.

Milk yield and THI on organic farms
Previous studies suggest a number of potential mechanisms that could explain why milk yield on organic farms is not affected by THI (stratifying the data did not show any effect in our study to establish a relationship between THI and milk yield) in the Fig. 4 Correlation between days in milk (DIM), parity and milk yield, and milk yield and milk parameters. The curve in the first plot is the fitted curve same way as it appears to be in conventional systems. These include genetic differences in herds between the systems, and environmental differences arising from differences in management and infrastructure. With regard to animal genetics, Dunn et al. (2014) linked differences in the threshold values of THI that lead to heat stress to genotype and production level. High-yielding cows have been shown to be more susceptible to heat stress than lowyielding cows (Bohmanova et al. 2007;Collier et al. 2006). Kadzere et al. (2002) proposed that this could be linked to changes in the physiology of thermoregulation of cows with intensive genetic selection towards high milk production. In addition, the increase in yield variability as temperature increased in the present study can be explained by the fact that in a restrictive environment, issues such as high heat loads, may only impair the expression of the full genetic potential of high-producing animals (Carabaño et al. 2016). Thus, the lesser impact of THI on milk yield in organic farms may be due to organic herds being less strongly selected for high milk production (Wallenbeck et al. 2018;Blanco-Penedo et al. 2019). In fact, the heterogeneity of variance might be considered in line to be consistent with the maintenance of genetic variance and variability of production in the herd. In organic farming, local and adapted breeds are preferred over high yielding breeds for their greater resistance to diseases and their heat tolerance, given that under organic management, more emphasis is placed on the adaptability of the animal to conditions than on the use of inputs to shield highly specialized animals from environmental variability (Carabaño et al. 2019). The resulting genetic diversity in organic herds may therefore offer organic farming systems a degree of resilience to climate change impacts.
With regard to management, technological solutions relating to barn design and environmental controls applied to tackle heat stress in intensive indoor production systems (Ji et al. 2017) may not be feasible for more extensive organic systems. However, environmental conditions within organic farm systems offer opportunities to both protect cattle and enable them to adapt to conditions. These management actions and biological processes extend from the scale of hours (i.e. providing shade, possibility of faster recovery outdoors at night, and moving grazing schedule to low-heat-stress hours), days (heat accumulation, hormonal responses, and change in feeding strategies), and seasons (cows spending more time standing in shade or eating more during the night, biological acclimation of animals, and lactation peak), to years (composition of herds by age, sex, body weight, and resistance to heat stress) and beyond (genetic selection through breeding) (Silanikove 2000;Kendall et al. 2006;Renaudeau et al. 2012). Advanced planning with an understanding of animal responses to thermal challenges and the ability to provide management options is needed to prevent or mitigate the adverse consequences of heat stress (Nienaber and Hahn 2007). However, the extent to which conditions in organic systems affect heat stress effects is still relatively unknown. Implementation of research on the effects of changes, such as the introduction of innovative housing solutions to better protect livestock against climate change-related weather extremes, is essential. More broadly, the effectiveness and efficiency of actions to tackle heat stress might be improved by an exchange of ideas between organic and conventional farming communities in order to identify options that have been developed for or within one community but which might benefit or be adapted for use in the other.
The importance of parity to milk yield in the present findings has not been described previously. It is well known that organic herds differ from conventional herds in age distribution (Frank et al. 2019;Sorge et al. 2016;Villar et al. 2016) (older animals at first calving and a higher proportion of animals with a higher number of lactations). This systemrelated difference might explain the greatest source of variation in milk yield compared to studies in non-organic herds.

Milk quality and THI on organic farms
Our results corroborate those of previous studies showing that cows producing less milk produce milk with lower fat and protein content. Additionally, our study reveals that the correlation between THI and lower fat and protein content of milk is more evident if cows produce less milk. This effect is particularly important in organic farms producing milk for systems that reconstitute milk solids or convert milk into other dairy products. Milk quality has been an important selling point for organic farming, with organic milk fat being found to have significantly higher levels of phospholipids compared with conventional milk fat (e.g. Ferreiro et al. 2015). Organic milk has also been found to contain high levels of conjugated linoleic acid (CLA), which has been shown to be altered by heat stress (Liu et al. 2017;Liu et al. 2020).
Again, the impact of THI on milk quality in organic systems may be related to both genetic and management factors in these systems; for example, high phospholipid content in organic relative to conventional milk has been attributed to differences in diet, season, stage of lactation, and breed between organic and conventional systems (Schwendel et al. 2015). A major issue in untangling the effects of these factors has been that the role of heat stress is often confounded by seasonal changes in feeding patterns (Liu et al. 2020). Differences in diet, including access to grazing in extensive organic systems, represents an important difference between organic and intensive conventional systems in Spain (Orjales et al. 2019). The extensive management of the farms in our study (grazing and nights outdoors in summer, where cows can recover better from the heat load in the day) may therefore have been a major contributor to our findings, together with the system-related farm practices mentioned above. Other parameters that may affect the relationship between THI and milk quality include exposure to wind and the availability of pasture with shade trees or shade shelters (Lee and Hillman 2007). Providing shade and shelter to the animals in pasture is part of the compulsory control scheme in organic farms (EC 834/2007) and may allow animals to reduce at least 30-50% of their total heat load (West 2003). However, under heat load, breed and genetics might be largely responsble for variation in milk quality, since the response is highly dependent on cow-factors (different heat stress susceptibility). Decreased feed intake has been shown to account for 35% of reduced milk production in hot conditions, while the remaining 65% is attributable to direct physiological effects of heat stress (Rhoads et al. 2009); however, lower intake was probably less important in our study, and other aspects might have been more important here. Apart from lipid synthesis in the mammary gland, diet is also a key determinant of the phospholipid content of raw milk (Ferreiro et al. 2015). The decreased milk protein content due to heat stress is the result of the specific downregulation of mammary protein synthetic activity (and not an artifact of the overall milk yield reduction). The mechanisms involved remain largely unknown. However, it is suggested to be less variable than milk fat composition, where the major contributor for milk protein synthesis is (non-diet-dependent) the limited amino acid supply to the mammary gland (Gao et al. 2017) (Table 3).
With regard to genetics, our results show (Fig. 3) that responses of milk quality to THI vary between individual animals; although much of this variation could be explained by parity, genetically based differences in susceptibility to heat stress have been previously demonstrated (Schwendel et al. 2015). Further studies are needed to investigate the relationship between milk quality and THI by breed, including breeds chosen to maximize milk solids production in organic systems in places such as Ireland and New Zealand (Rodríguez-Bermúdez et al. 2019). It will also be important to further investigate the mechanisms behind potential genetic causes for variations in susceptibility to heat stress, for example through investigation of gene expression for the de novo synthesis of specific fatty acids in the mammary gland during milk production (Knutsen et al. 2018).

Variation in productivity on organic farms
Our findings showed a high level of inter-farm differences in milk yields and their variation over time, confirming previous observations of the diversity of systems included in the 'organic' bracket. Understanding such inter-farm variation in performance is vital, as the environmental benefits of using organic farming approaches can be substantially negated by low productivity relative to conventional systems (Smith et al. 2019). Differences in productivity may suggest the existence of opportunities for improving performance through the wider adoption of best practices. However, it can be hard to tease apart differences in management that represent important adaptations to local conditions from those that represent deviations from best practice. Achieving resilience may involve a trade-off of production levels and yield stability under environmental change. The herds studied here represent particularly extensive production systems compared to organic dairy farms in other European countries (Blanco-Penedo et al. 2019), with the majority (89.6%) of the Spanish herds studied belonging to a cluster with the lowest average input of feed concentrate per cow (i.e., cows are usually fed a non-grain-based ration) and the highest level of access to grazing. For such farms, resilience to local conditions may be prioritized over maximizing yields, while variation in environmental conditions can be expected to drive variations in performance in ways that are less likely in highly controlled indoor systems. These factors make it all the more important to differentiate necessary (and potentially innovative) local adaptation to conditions from unnecessary variations in management that represent suboptimal decision-making in organic systems.

Modeling approach
The differences between our findings and those of other studies on heat stress in Spanish dairy farms may be associated with differences in the statistical models and methods used here. One limitation of studies using THI as a proxy for heat stress is that THI does not consider important factors relevant to the extensive management of the organic farm system, incorporating only environmental parameters and not considering the characteristics of the individual cow. This is problematic, as patterns of animal response may differ across climatic conditions, production systems, or both (Carabaño et al. 2016). Work is needed to develop heat stress indicators that incorporate these factors.

Next steps
Continued climatic change adds urgency to the need to study the effects of heat stress on dairy cows under less intensive production systems. Such studies can underpin strategic improvements in management and selection to enhance the ability of farm animals to cope with environmental challenges in all systems (see Supplementary Material Fig. 1). A framework for a more 'holistic' quantification of external environmental impacts on dairy cow health and productivity is key for organic production systems where many outcomes are influenced by environmental factors and are beyond direct managerial control. Future tools for measuring, assessing the risks of, and optimizing responses to heat stress will need to characterize the diversity of the agroecosystem in order to be effective. Differences between individual farms and between parities observed require further research to identify best and worst performance and to focus on ways to improve farms with the worst performance. Sustainability strategies should consider the optimization of production conditions and the potential of the animals, keeping them far from their biological limitations while avoiding conflicts between welfare and productivity (McInerney 2004). In order for organic farming to continue in its significant contribution to green growth, a comprehensive understanding of these farming systems, including the economic and environmental impacts of managerial decisions, is paramount. Tackling this challenge will enable a better understanding of the extent to which organic farming represents a resilient production strategy in the context of global warming.

Conclusions
Our results suggest that in regions such as the Mediterranean basin, where climate change is having-and will continue to have-significant impacts on weather conditions, the organic dairy sector may bear at least some costs through reductions in milk quality. However, further, larger-scale studies are needed to verify these findings and to explore their implications. In particular, researchers need to focus on the specificities of agricultural systems and on farm management decisions as important sources of innovation to increase resilience to climate change, in order to generate knowledge of real practical value to farmers.
Author contributions Alejandro Ruete (AR) and Isabel Blanco-Penedo (IB-P) developed the study design. IB-P collected the data at the farm. Secondary databases from the Spanish Milking Record scheme and the Spanish Meteorological Agency were obtained by IBP and Antonio Velarde (AV). AR performed the statistics and IBP helped with the biological interpretations. IB-P, AR, Richard Kipling (RPK), and AV were involved in writing the original draft and IB-P and RPK in manuscript review and editing under supervision of AR). All authors have read and approved the final version of the manuscript to be published.
Funding Information Open access funding provided by Swedish University of Agricultural Sciences.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.