Temperature-mortality relationship in dairy cattle in France based on an iso-hygro-thermal partition of the territory

The issue of global warming and more specifically its health impact on populations is increasingly concerning. The aim of our study was to evaluate the impact of temperature on dairy cattle mortality in France during the warm season (April–August). We therefore devised and implemented a spatial partitioning method to divide France into areas in which weather conditions were homogeneous, combining a multiple factor analysis with a clustering method using both weather and spatial data. We then used time-series regressions (2001–2008) to model the relationship between temperature humidity index (an index representing the temperature corrected by the relative humidity) and dairy cattle mortality within these areas. We found a significant effect of heat on dairy cattle mortality, but also an effect of cooler temperatures (to a lesser extent in some areas), which leads to a U-shaped relationship in the studied areas. Our partitioning approach based on weather criteria, associated with classic clustering methods, may contribute to better estimating temperature effects, a critical issue for animal health and welfare. Beyond the interest of its use in animal health, this approach can also be of interest in several situations in the frame of human health.


Introduction
The effect of temperature on human health has been evidenced for decades, but only recently has the development of distributed lag non-linear models (DLNM) [1] allowed a better estimation of the temperaturemortality relationship by considering both the delayed effect of the exposure and the non-linear structure of the relationship. Typically, the relationship between an environmental factor (temperature, pollution, etc) and mortality is first assessed in several sites and a metaanalysis is then performed to pool the model parameters obtained for the different sites [2][3][4][5][6]. The first step is usually performed at the city level, which has the double advantage of being a spatial scale where the temperature is uniform and the population density is high, which ensures power and precision in the models.
The influence of temperature on cattle mortality has been evidenced by different methods [7][8][9]. Two recent studies, based on DLNM, described the pattern of the relationship between temperature and cattle mortality [10,11]. However, for animals, the choice of the spatial unit to assess the relationship in different sites, with both a homogeneous distribution in temperature and a high population density, remains problematic.
In this study, we used a partitioning procedure that addresses this issue to evaluate the impact of temperature on dairy cattle mortality in France during the warm season (April-August). We proceeded in two stages: we first performed a partitioning of the national territory into clusters of homogeneous temperature and relative humidity; we then evaluated the relationship between temperature and dairy cattle mortality within each cluster.

Methods
Spatial partitioning: Weather data were provided by Meteo-France (https://donneespubliques. meteofrance.fr/). The minimum daily temperature, maximum daily temperature, average daily temperature and average daily humidity were available on a square lattice with 8 km sides for the period 2001-2008.
We partitioned France on a grid of approximately 1300 regular hexagons (side = 13.2 km, area = 453.6 km 2 ) and computed average daily values of each weather variable for each hexagon. The clustering of the territory was carried out using weather data for 2002, because no major climate event was recorded for that year and we assume that the climatic partitioning did not evolve between 2001 and 2008. The iso-hygro-thermal clusters were built by grouping hexagons similar in temperature and humidity during this year that were geographically close. The hexagons were grouped in two steps: We first performed a multiple factorial analysis (MFA) on five blocks of variables: three blocks of 365 daily temperature values (minimum, maximum and average temperatures), one block of 365 daily relative humidity values and one block of geographical coordinates (latitude and longitude) of the centroids of the hexagons.
We then performed a hierarchical ascendant clustering (HAC) on the factorial coordinates of the hexagons derived from the MFA. The generalized Ward's Criterion was used as the aggregation criterion for HAC. This consists of aggregating two clusters in a way that minimizes intra-cluster variance and maximizes inter-cluster variance [12].
We arbitrarily defined two trim levels in the hierarchical clustering tree to get two levels of hexagon aggregation. We chose the first trim so that the first partition was composed of 100 sites. This choice ensures that the sites obtained are sufficiently small to be homogeneous in terms of weather conditions and, at the same time, large enough to contain a sufficiently high population, two conditions for obtaining accurate sitespecific estimates of the temperature humidity index (THI)-mortality relationship [13]. We excluded the sites with an average of fewer than ten dead dairy cattle per day to fulfil the second condition. The second trim defined nine large climatic areas used to estimate the pooled relationships. Given the properties of hierarchical clustering, the 100 sites were nested in the nine study areas.
Estimation of the THI-mortality relationship: Data on dairy cattle were extracted from the French National Cattle Register managed by the French Ministry of Agriculture. We selected dairy cattle present in mainland France during the warm season (April-August) from 1 January 2001 to 31 December 2008. We computed daily numbers of dead and live dairy cattle per site. We focused the analysis on unplanned deaths (on-farm or euthanasia) and thus did not take slaughtered dairy cattle into account. THI was computed as follows: THI = air temperature [ For each site, we fitted a DLNM to assess the association between the daily number of dead dairy cattle and THI. A Poisson generalized linear model allowing for over-dispersion in daily counts of deaths was constructed: = DOW t + PH t + PH t+1 + PH t+2 + NS (THI, df.THI, df .lag) where E(Y t ) was the expected number of daily dead dairy cattle at day t; NS(.) denoted a natural cubic spline; NS(t, df.t), the spline function for time, controlled for a smooth long-term trend; NS(t_in_Year, df.t_in_Year), with t_in_Year the day of the year, controlled for a regular seasonal trend that was forced to be identical each year using df.t = 3 and df.t_in_Year = 5 degrees of freedom, respectively. Dichotomous variables indicating the day of the week (DOW t ), public holidays (PH t ) and the first and second day after a public holiday (PH t+1 and PH t+2 ) were also included in the model; it has been shown that mortality rates are lower on public holidays and more variable on the following two days, because dairy cattle that die on public holidays are generally collected and declared on the two days following the holiday. The spline function for THI was bi-dimensional with five degrees of freedom for the lag (df.lag) and four for THI (df.THI). Spline knots were placed at equally-spaced quantiles along the national-level average temperature range, and knots in the lag space were set at equally spaced values on the log scale of lags to allow more flexible lag effects at shorter delays [14]. The log of the number of dairy cattle at risk on day t (N t ) was included in the model as an offset. The maximum lag was set at 21 days to explore the lag structure of the THI effect.
We then performed multivariate meta-analyses to combine results across the different sites at the level of the study areas [15,16]. The THI-lag-mortality association was first reduced to the overall THI-mortality association, cumulating the risk during the lag period. This step reduced the number of parameters to be pooled in the second-stage meta-analysis, while still preserving the complexity of the estimated dependency between THI and dairy cattle mortality [15]. The regression coefficients of the relationship estimated in the first stage were used as outcomes for the second stage and multivariate meta-analyses were defined by intercept-only models. In the same way, lag-specific effects at the 5th and 95th percentiles of THI were computed for every site and the relationship between relative risk (RR) of death and lags was then pooled in each study area.
In each area, the THI corresponding to the level of minimum mortality (MMTHI) was used as the reference to compute the RR. We computed RR and the associated 95% confidence intervals for extreme and moderate cold values, corresponding respectively to  Statistical analyses were performed with R software [17] using the 'FactoMineR' [18], 'dlnm' [14] and 'meta' [15] packages.

Results
Spatial partitioning: The HAC on the coordinates of the hexagons in the factorial space constructed by the MFA resulted in the characterization of nine study areas (figure 1). Each study area was formed by gathering together 5 to 16 sites, which were themselves formed by the grouping of hexagons. The distribution of THI in each study area (median and percentiles) is shown in table 1.
Estimation of the THI-mortality relationship: The 100 sites defined by the clustering counted on average between 0 and 152 dead dairy cattle per day for the 8 years of study (table 2). Based on Bhaskaran's recommendations on time-series regression studies in environmental epidemiology [13], we did not estimate the temperature-mortality relationship in the clusters with fewer than 10 deaths per day on average, because this would not have provided enough power and precision to accurately estimate the relationship in these clusters. This concerned all of the study areas 1 and 9, which are known to have a low cattle population, and 36 sites (46%) distributed across the remaining seven study areas. The low mortality figures in these clusters resulted from a low number of dairy cattle, and were not related to higher survival rates. The pooled relationship between THI and dairy cattle mortality had a U-shape, whose flatness or depth depended on the study areas (figure 2). It generally matched the form of the within-sites relationships, although the effect of cool temperatures was more heterogeneous among sites in study areas 4, 5 and 7. MMTHI ranged between 13 • C THI and 15 • C THI ( extreme cold, but was not significant in half of the areas. The lag-response association for moderate heatestimated from the models in the different sites and study areas-indicated that a THI increase from MMTHI to the 95th percentile was associated with a steep increase in mortality risk on the day of exposure and a decreasing risk during the 2-5 following days (figure 3). Area 6 was characterized by an increase in the risk after exposure, followed by a decrease and a subsequent protective effect for a lag between 4 and 9 days, but the confidence intervals were very large. The pattern for the effect of cold was less pronounced (figure 4). Unlike the heat effect, the exposure to low THI values had no immediate effect on mortality risk and was even protective in some study areas. Figure 2. THI-mortality relationship for the sites (dotted lines) and for the study areas (solid red lines). Shaded areas are 95% confidence intervals for study areas. RR of mortality for the THI relative to the THI with the lowest mortality.
An increase in the mortality risk appeared a few days after exposure and lasted only a few days.

Discussion
In this study, we divided the territory into homogeneous temperature and humidity clusters to estimate the relationship between THI, a temperature-related index, and dairy cattle mortality more accurately than by using pre-defined geographical units, e.g. administrative units, which do not have climate consistency. We focused the partitioning of the territory primarily on the temperature pattern. In practice, we gave more weight to temperatures by including the three temperature variables (minimum, average and maximum) in three different blocks of the MFA. Although strongly correlated, the inclusion of the three variables in the MFA also enabled us to take into account the average temperature, the extreme temperatures and the range of daily temperature changes. The MFA also considered the geographical coordinates of the hexagons' centres to integrate a spatial proximity dimension, and limited, but did not exclude, the creation of completely spatially-fragmented clusters. Accordingly, our final partitioning was formed by six compact study areas and three partially fragmented areas, which corresponded to mountainous areas (where altitude-related climatic features were more important than the compactness of the area). The first-level partitioning into 100 sites also displayed partially fragmented sites in the eastern and central parts of France (supplementary file 1 available at stacks.iop.org/ERL/12/114022/mmedia); however, the sites were more compact than when the partitioning was performed without accounting for the geographical coordinates (results not shown). . RR of mortality estimated by the multivariate meta-analytical model for moderate heat (95th percentile of the study areaspecific THI distribution), relative to the minimum mortality THI for lags between 0 and 20 d. Dotted lines stand for the sites and solid red lines for the pooled relationship in study areas. Shaded areas are 95% confidence intervals around the pooled relationship.
A recent study has already shown the effect of temperature on dairy cattle mortality in France [10]. The main limitation of that study was that weather data were available from 50 meteorological stations, among which only 12 stations were retained on the basis of having a minimum dairy cattle population on which to run the models. Despite this selection, models for half of the 12 stations suffered from a lack of precision and power. In the current study, we overcame this issue by using the partitioning approach. Because we used weather data on a fine and exhaustive grid, we could aggregate hexagons in order to have a satisfactory trade-off between intra-cluster homogeneity of the weather conditions and the population size inside the clusters. The advantage of the method is that the choice of the trim level enables fine-tuning between the two properties.
Our study showed that both high and low THI were associated with an increased risk of mortality among dairy cows in France during the warm season. Although not expected, the predicted effect of the low THI could be attributed to the association of cold temperatures with other weather parameters, such as precipitation and wind speed, in particular in freerange cows, as suggested in a recent study in Belgium [11]. The authors of the Belgian study found an RR of 1.09 (95% CI: 1.02-1.17) for moderate heat (95th percentile, 19.7 • C) and 1.26 (1.08-1.48) for extreme heat (99th percentile, 22.6 • C). These estimates tended to be lower than the RR predicted in our seven study areas . RR of mortality estimated by the multivariate meta-analytical model for moderate cold (5th percentile of the study areaspecific THI distribution), relative to the minimum mortality THI for lags between 0 and 20 d. Dotted lines stand for the sites and solid red lines for the pooled relationship in study areas. Shaded areas are 95% confidence intervals around the pooled relationship.
(although not statistically different, since CIs were not exclusive). Differences between the two studies in values considered as moderate and extreme heat, in study populations (restricted to dairy cattle over two years old in the Belgium study) as well as in farming practices (i.e. animal management, animal housing, use of heat abatement technology to reduce heat stress) may explain the observed difference in mortality risk. Furthermore, as in the Belgian study [11], we showed an acute effect of heat during the first days, followed by a progressive return of the risk to near zero. We predicted a similar pattern in most areas for the coldest temperatures, with an effect that appeared after a delay of a few days and also lasted a few days. In contrast, this effect lasted for more than two weeks in the Belgian setting, probably because it was measured during winter.
A knowledge of temperature-related mortality in dairy cattle is important for animal health and welfare, as well as for economic reasons. The proposed partitioning method, based on weather criteria, is innovative and solves the issue of the definition of accurate climate areas in a simple way that is useful for studies in animals. We believe this approach will be particularly useful in further studies, especially to quantify the effect of heat waves. Such analyses indeed require a very accurate estimation of the temperature-mortality relationship to discriminate between the added effect of extreme temperature (heat wave) and the 'normal' effect of temperature.
Beyond its relevance in animal health, our geographic aggregation approach can also be of interest in human health, in particular to examine the impact of environmental stressors in rural areas. While many health studies have been conducted on climate stress or air pollution at the level of big cities, where exposure is relatively homogeneous and populations are large, only a limited number of studies have explored these effects in low-populated areas. If it is a matter of methodological issue, we believe our approach can help overcoming these limitations. It can also be useful within the framework of environmental topics that are specific to rural populations, such as exposure to pesticides.