Fluctuating temperature modifies heat-mortality association around the globe

Studies have investigated the effects of heat and temperature variability (TV) on mortality. However, few assessed whether TV modifies the heat-mortality association. Data on daily temperature and mortality in the warm season were collected from 717 locations across 36 countries. TV was calculated as the standard deviation of the average of the same and previous days’ minimum and maximum temperatures. We used location-specific quasi-Poisson regression models with an interaction term between the cross-basis term for mean temperature and quartiles of TV to obtain heat-mortality associations under each quartile of TV, and then pooled estimates at the country, regional, and global levels. Results show the increased risk in heat-related mortality with increments in TV, accounting for 0.70% (95% confidence interval [CI]: −0.33 to 1.69), 1.34% (95% CI: −0.14 to 2.73), 1.99% (95% CI: 0.29–3.57), and 2.73% (95% CI: 0.76–4.50) of total deaths for Q1–Q4 (first quartile–fourth quartile) of TV. The modification effects of TV varied geographically. Central Europe had the highest attributable fractions (AFs), corresponding to 7.68% (95% CI: 5.25–9.89) of total deaths for Q4 of TV, while the lowest AFs were observed in North America, with the values for Q4 of 1.74% (95% CI: −0.09 to 3.39). TV had a significant modification effect on the heat-mortality association, causing a higher heat-related mortality burden with increments of TV. Implementing targeted strategies against heat exposure and fluctuant temperatures simultaneously would benefit public health.


INTRODUCTION
Global warming pervasively affects human life and undermines the years of gains in public health. 1,,2, Owing to the increasing rate of 0.2 C in temperature per decade, human-induced warming has been associated with an increase in frequency and intensity of hot days, reaching 2.9 billion additional person-days of exposure to heatwave events of vulnerable populations older than 65 years in 2019. [3][4][5][6] Increasing heat exposure, in turn, results in excess morbidity or mortality. 7 37.0% (95% CI: 20.5-76.3) of heat-related deaths in the warm season can be attributed to human-induced heating. 10 Another challenging issue driven by climate change is temperature variability (TV), an indicator of short-term temperature fluctuations or stability. Previous studies have projected an increasing trend of TV in some regions, in particular in tropical countries. 11,,12, Extensive studies have established evidence of the health effects of TV, showing a significant association between TV and mortality in many parts of the world and substantial public health burden stemming from TV. 13,,14, Traditionally, studies of the health effects of TV have controlled for mean temperature as a confounder. 13,,15,,16, Based on similar biological mechanisms underlying the health effects of TV and heat exposure, it is worth investigating whether there exists a synergistic effect between them. 17, Previous studies observed season-differentiated effects of TV on mortality, which suggests potential effect modification of mean temperature on TV-related mortality. 13,,18 However, to the best of our knowledge, few studies have examined whether TV modifies the heat-mortality association. For example, the temperatures fluctuations from 5 C to 25 C and from 13 C to 17 C represent the same mean temperature of 15 C, but their effects on mortality are very likely to be different. Assessing heat-mortality association without considering the modification of TV may fail to recognize the most severe weather in relation to heat and to implement an effective early warning system.
In this study, using data from the Multi-Country Multi-City (MCC) Collaborative Research Network, we systematically evaluated the contribution of TV to the heat-mortality association in 717 locations across 36 countries over the period 1972-2018. Through this study, we aimed to provide a more complete picture of the TV-differentiated heat-related mortality burden and to provide scientific evidence that could improve the sensitivity of current heat-health warning systems in hot temperatures with dramatic temperature fluctuations.

RESULTS
The descriptive statistics from each country are shown in Table 1. A total of 36.42 million deaths from total or non-external causes were identified during the warm season from 1972 to 2018 (Table 1). On average, the median daily TV across 717 locations was 5.8 C (25th-75th percentile: 4.9-6.7). The average daily mean temperature under each quartile of TV was 20.6 C (Q1), 21.7 C (Q2), 22.3 C (Q3), and 23.0 C (Q4) ( Table 1). The overall correlation coefficient between daily mean temperature and TV was 0.10 (Table S2). The summary descriptions of the daily mean temperature and TV for each location are shown in Tables S3-S5. Figure 1 shows exposure-response curves between daily mean temperature and mortality in warm season. In general, J-shaped associations between daily mean temperature and mortality were found under different groups of TV, with the risks increasing dramatically at extreme hot temperatures ( Figure 1A). The MMTs for different quartiles of TV were 17.62 C (Q1), 17.62 C (Q2), 18.39 C (Q3), and 19.80 C (Q4). The difference in four curves was tested as statistically significant. Regional exposure-response curves indicate potential geographical patterns in the modification effect of TV on the heat-mortality association (Figure 1B). In general, higher mortality risks were observed for higher TV levels across different regions. Southern Europe and central Europe generated a greater difference in heat-related mortality risks between Q1 and Q4 of TV. For most countries, the country-specific curve demonstrated increased mortality risks as TV rose ( Figure S1).
In the sensitivity analyses, our results were robust. The same patterns were observed after changing the length of exposure to TV (Figures S2-S5; Tables   Table 1. Mortality data and description of daily TV and mean temperature in each stratum (from Q1 to Q4) of TV in 717 locations from 36 countries during the warm season S9-S12) and after adding separate predictors in meta-regression with BLUP (Table S13). Using incremental lag periods, shortening the duration of the warm season, and adding relative humidity to the model, the modification effect of TV still existed (Tables S14-S16). The AFs changed slightly after using different methods to handle missing values (Table S17).

DISCUSSION
Our study showed that heat exposure together with high TV could significantly increase the mortality risk in the warm season. We saw an upward trend in premature death due to heat exposure with the increase in TV. For temperatures higher than TVST (96.34th percentile on average), AF showed a greater difference across 4 TV groups. The TV-modified heat-mortality burden showed disparate geographical variations.
The physiological mechanisms that explain the synergistic effects of TV and heat exposure on mortality are not yet clearly defined. However, there are several ways in which the two exposures may interact. When exposed to heat, people expend more of their reserves on thermoregulation to respond to heat. 19, This process involves elevations in heart rate and blood pressure, vasodilatation to transfer heat to the skin, and respiration to lose heat with the expired air. [20][21][22] Physiological adaptation to higher temperatures takes time. If the temperature suddenly changes in a short period of time, then people may have difficulty with internal thermoregulation, resulting in inflammatory responses and coagulation abnormalities induced by heat stress. 23, In addition, sudden temperature changes may also activate bronchopulmonary vagal afferent nerves and the inflammatory response. 24, For people with underlying conditions (e.g., preexisting illness, chronic diseases, poor fitness level), heat exposure may place extra pressure on the cardiovascular and respiratory systems, 25 especially when temperature changes dramatically in a short period. Consequently, heat-related deaths may occur.
Our findings are generally consistent with previous studies focusing on TV or heat, indicating that both heat exposure and TV were positively associated with mortality. 13,,26,,27, Some studies have explored the modification effect of temperature on the association between mortality and diurnal temperature range (DTR) or TV. 17,,28-30 A study in Japan showed a higher risk of cardiovascular mortality associated with TV during extremely hot days in comparison with extremely cold days. 17, Another study in England and Wales showed a J-shaped curve of the relationship between percentiles of DTR and mortality in the warm season, while an inverted-V-shaped association was observed during the cold season. 28 Although many studies have quantified the mortality burden associated with heat exposure, to the best of our knowledge, none of them explored the potential heterogeneity attributable to TV. [31][32][33] Without taking TV into account when assessing the association between heat exposure and mortality, the heat-related mortality burden may be underestimated.
Some heterogeneity across regions was also found. For countries such as Guatemala and Colombia, temperatures above TVST together with low TV showed a higher mortality burden. The potential reasons for these results may derive from the low variation in TV and thus a small difference between Q1 and Q4 of TV. TV can be affected by many factors, such as greenness, soil moisture, and precipitation. 34,,35, For example, vegetation removal and soil aridation would act to increase daily temperature fluctuation with more rainfall. 34, Temperate desert steppe generally has comparable warming effects on T min and T max , while temperate meadow and temperate steppe may have larger cooling effects on T max than T min . 35, In addition, population aging may mediate the association. The variation in heat exposure is caused by the differences in both vulnerable populations (in particular, the elderly) and temperature across regions. 36, For the countries such as South Africa (7.7%) and Cambodia (6.8%), the vulnerable population aged 60 or older accounted for only a small proportion of total population, lower than the global average level of 12.3% in 2015. 37, The lower the fraction of older persons, the less sensitivity to heat exposure for the whole country, and the lower the mortality burden attributable to heat exposure. With the deepening of the aging of society, vulnerable populations are expected to increase. 37 Further research is warranted to explore the geographic variation in TV-differentiated heat-related mortality burden and call for targeted strategies in mitigation and adaptation against climate change.
As a response to the increasing heat conditions, heat-health warning systems are developed in some parts of the world. For example, the WHO Regional Office for Europe developed the heat-health action plans (HHAPs) in 2008, covering 35 out of 53 member states of the WHO European Region by the end of 2018. 38, Although it is hard to assess how much heat-health prevention was associated with HHAPs, a substantial reduction in heat-related deaths was observed since the implementation of preventive measures (either HHAPs or other types of intervention). [38][39][40][41] However, the majority of warning systems use mean temperature or T max as indices to trigger the warnings and inform the general public through the mass media without user-oriented attractive notification. 42,,43, As suggested by our findings, the health effects of heat exposure could be magnified if they are accompanied by a higher TV and the findings were consistent across countries in Europe where HHAPs are implemented, 44 implying that it is difficult for current warning systems to obtain effective heat prevention when the TV is high.

Strengths and limitations
This study has several strengths. First, to the best of our knowledge, this is the first study to systematically explore the modifying effects of TV on heat-related Article mortality risks. The assessment of TV-differentiated mortality risks associated with heat exposure provides a better understanding of heat vulnerability. When quantifying the heat-related mortality burden, the effect of TV should not be ignored. Second, this study benefits from the large-scale investigation across multiple countries and a long period of time. It enables us to provide a global vision of the heat-related mortality burden under different TV exposure and to explore potential variation in estimations in terms of socioeconomic status and climatic and geographic features. Finally, our study targeted the heat-related mortality burden more precisely by highlighting the modification effect of TV when temperatures were higher than TVST.
Several limitations should be acknowledged. As the time series design and temperature data from fixed monitoring stations were used, ecological fallacy and measurement errors in exposure seem to be inevitable. Due to the lack of data for some areas of the world, some regions contain only one country. The extrapolation and interpretation of the findings are restricted to the generalization of our results. In addition, we are unable to investigate the modification effect of TV on the association between heat exposure and age-or cause-specific mortality. As the biological mechanisms are sensitive to causes of death, further research is warranted to investigate the disease-differentiated modification effect of TV on heat-related mortality burden.

CONCLUSIONS
Our findings, which are based on multi-country data, revealed that higher TV over a short period of time increased the mortality risks associated with heat exposure. It is imperative to raise public awareness of the potential health risks of TV. Targeted adaptation strategies against heat-related mortality burden should be implemented after taking into account the fluctuation of temperatures and geographical patterns.

MATERIALS AND METHODS Data collection
The MCC Collaborative Research Network database (http://mccstudy.lshtm.ac.uk/) was used. Daily death counts and meteorological data, including mean temperature, maximum temperature (T max ), minimum temperature (T min ), and relative humidity were extracted. The International Classification of Diseases, 9th and 10th revisions (ICD-9 and ICD-10) codes were used to identify the causes of death. We extracted the data series on non-external causes of death (ICD-9: 0-799; ICD-10: A00-R99), or, if not available, all-cause mortality.
Our analyses were restricted to locations with complete weather data in the warm season (the warmest 4 consecutive months) for at least 2 consecutive years. Detailed information on data cleaning is described in Text S1. Finally, 717 locations across 36 countries were included. Daily mortality, mean temperature, T min , and T max data had overall missing rates of 0.14%, 0.84%, 0.86%, and 0.73%, respectively (Table S1).

Calculation of temperature variability
TV was calculated as the standard deviation (SD) of the average of T min and T max for the current day and 1 day before (T max-lag0 , T max-lag1 , T min-lag0 , T min-lag1 ). 13 In the sensitivity analysis, higher lengths of exposure (0-2 days and 0-3 days) were applied.

Statistical analysis
We used a two-stage time series design to assess the modification effect of TV on the heat-related mortality burden in the warm season. In the first stage, a generalized linear regression with the quasi-Poisson family allowing overdispersion in the death counts was applied for each location to obtain location-specific estimates for the heat-mortality association. To capture the modification effect of TV on the heat-mortality association, we introduced an interaction term between a cross-basis function of daily mean temperature and quartiles (Qs) of TV. We used relative rather than absolute levels of TV to accommodate different levels of adaptive capacity to TV across locations. The equation was as follows: where Y it is daily deaths counts in location i on day t; TV it is the linear function of TV; and cbðTemp it ; lag = 10Þ, built by distributed-lag nonlinear models (DLNMs), is a cross-basis function of daily mean temperature featuring the nonlinear and delayed association over 10 days of lag. 45, It incorporates a natural cubic spline function with two internal knots placed at the 50th and 90th percentiles of the location-specific temperature distributions during the warm season and a natural cubic spline function of lag days with 2 internal knots placed at equally spaced values in the log scale to capture the delayed effect of temperature. 10 Quartile TV stands for the dummy coded categorical variable of TV groups (from Q1 to Q4). To control for unmeasured temporal trends such as seasonality and long-term trend, nsðTime it ;df = 4 =yearÞ, a natural cubic function of calendar days with 4 degrees of freedom (df) per year was included. In addition, an indicator for the day of the week ðDOW it Þ was adjusted in the model to control for weekly variations in risk.
In the second stage, multivariate random-effects meta-analysis without predictors was used to pool the location-specific estimates at the global, regional, and national levels. Precise estimates for each location were obtained using the best linear unbiased prediction (BLUP) estimations. BLUPs can provide more accurate estimates in locations with small daily death counts or short time series by borrowing information across locations. 46 Quantification of temperature-specific risks and attributable fraction For each quartile of TV, we showed the risk over a 10-day lag period at each temperature value compared with the risk at minimum mortality temperature (MMT), at which the risk of mortality was the lowest. The statistical significance of the difference in heat-mortality risks across quartiles of TV was assessed using repeated-measures multivariate meta-analysis. Briefly, based on location-specific estimates of four quartiles of TV, we performed a random effects meta-regression with TV quartiles as the only meta-predictor to test the difference in overall exposure-response curves of heat-mortality risks across quartiles of TV. Furthermore, in each country, to identify the potential targeted percentile ranges of temperature modified by TV, we tested the statistical significance of the difference in heat-mortality risks between Q1 and Q4 of TV using a fixed-effects meta-regression model at each temperature value. The modification effect of TV was identified as being statistically significant if there was a significant difference in heat-mortality risks between Q1 and Q4 of TV. The fixed-effects model was used because these country-specific estimates were based on the same samples. The country-specific temperature percentiles above which statistically significant modification effects of TV were observed, were called TV sensitive heat thresholds (TVSTs) and used to separate components of mortality burden attributable to heat exposure.
We compared the attributable deaths and attributable fractions (AFs) associated with heat exposure above MMT for each quartile of TV. Two components of heat exposure were used to separate the overall AF: from the quantile of MMT to TVST and above TVST. Daily attributable deaths due to each component of heat were calculated using BLUP location-specific association and then summed to obtain the total attributable deaths during the study period. 47 AF was computed by dividing the attributable deaths by the total death counts. Monte Carlo simulation (n = 1,000) was used to derive 95% CIs.

Sensitivity analysis
We conducted several sensitivity analyses to check the robustness of our results: (1) using different lengths of exposure to TV (TV 0-2 and TV 0-3); (2) adding separately location-specific predictors (region, average mean temperature, range of mean temperature, indicators for Köppen-Geiger climatic zones, gross domestic product [GDP] per capita, latitude, and longitude) to the meta-analytical model with BLUP in the second stage to check whether the pooled estimates would change while adjusting for the effects of each predictor on the location-specific estimates; (3) choosing alternative 2 and 3 warmest consecutive months to define the warm season; (4) changing the lag days of heat from 7 to 13 days; and (5) controlling the potential effect of relative humidity using a natural cubic spline with 3 df. Finally, although the overall missing rates for mortality and temperature data were generally small (Table S1), we also performed sensitivity analyses by restricting our analyses to locations with complete data and by using complete data after imputation. Missing values in time-series data were imputed by the spline interpolation method. R software (version 3.6.2) with packages "dlnm" (for the construction of the cross-basis functions), "mvmeta" (for meta-regression), and "imputeTS" (for spline interpolation of the time-series data) was used to perform all of the analyses. A two-sided p < 0.05 was set as statistically significant. The code is available at the personal website of the first author (Github: https://github.com/yaowu-ops/Modification-effect-of-TV-on-heat-mortalityassociation.git).

Data sharing
Data were collected within the MCC Collaborative Research Network under a data-sharing agreement and cannot be made publicly available. Researchers can refer to MCC participants, who are listed as coauthors of this article, for information on accessing the data for each country.  Figure 2 Fractions of all-cause mortality attributable to heat exposure by temperature variability (A) TVST identified for each country for TV. *For countries without identifiable TVST, the quantile threshold of the 96.34th percentile (the average of all identifiable country-specific TVSTs) in the temperature distribution for each country was used. (B) Comparison of AFs of mortality due to heat exposure for Q1 and Q4 of TV in each country, stratified by TVST. Yellow bars: the AFs for the Q1 of TV; purple bars: the AFs of mortality due to heat exposure above MMT for the Q4 of TV; dark blue bars: the AFs of mortality due to heat exposure from MMT to TVST for the Q4 of TV; light blue bars: the AFs of mortality due to heat exposure above TVST for the Q4 of TV. AF = attributable fraction; MMT = minimum mortality temperature; Q1 = the 1st quartile; Q4 = the 4th quartile; TV = temperature variability; TVST = temperature variability sensitive heat threshold.

DECLARATION OF INTERESTS
The authors declare no competing interests.

Table of contents
Text S1. Information about Data collection.      Table S1. Descriptive information on missing rates of mortality and temperature data at the country level. Table S2. Correlation coefficients between daily mean temperature and temperature variability. Table S3. Description of TV and daily mean temperature in each quartile of TV for each location included in our study. Table S4. Description of TV 0-2 and daily mean temperature in each quartile of TV 0-2 for each location included in our study. Table S5. Description of TV 0-3 and daily mean temperature in each quartile of TV 0-3 for each location included in our study. Table S6. Attributable deaths of mortality due to heat exposure, stratified by quantiles of TV in each region. Table S7. Attributable fractions of mortality due to heat exposure stratified by quantiles of TV in each country. Table S8. Attributable fractions of mortality due to heat exposure stratified by TVST for each quartile of TV in each country. Table S9. Attributable fractions of mortality due to heat exposure stratified by quantiles of TV 0-2 in each country. Table S10. Attributable fractions of mortality due heat exposure stratified by quantiles of TV 0-3 in each country. Table S11. Attributable fractions of mortality due to heat exposure stratified by TVST for each quartile of TV 0-2 in each country. Table S12. Attributable fractions of mortality due to heat exposure stratified by TVST for each quartile of TV 0-3 in each country. Table S13. Overall attributable fractions after adding separate predictor in meta-regression with BLUP, stratified by quartiles of TV. Table S14. Sensitivity analyses of overall attributable fractions after changing the lag of mean temperature in the main model, stratified by quartiles of TV. Table S15. Sensitivity analyses of overall attributable fractions using different definitions of warm season, stratified by quartiles of TV. Table S16. Sensitivity analyses of overall attributable fractions after adjusting relative humidity in the main model using 556 cities with relative humidity data, stratified by quartiles of TV. Table S17. Sensitivity analyses of overall attributable fractions using different complete data after removing cities with missing vlaues or after imputaion, stratified by quartiles of TV.
Text S1. Information on data collection and data cleansing We collected data on daily mortality and mean temperature (Tmean) in warm season for 750 cities across 43 countries or regions from the database of the MCC Collaborative Research Network (http://mccstudy.lshtm.ac.uk/). Detailed information on the process of data collection has been described in previous publications 1 . As daily minimum temperature (Tmin) and maximum temperature (Tmax), the primary exposures of interest in our study, were not available for some countries, we downloaded data from official websites instead. For European countries including U.K., Finland, Germany, Spain, we use the daily Tmean, Tmin, and Tmax data (10 km × 10 km resolution) extracted from the European daily high- For each location, the following data cleansing procedures were used. First, we collected locations that have daily measurements of Tmean, Tmin, Tmax, and mortality for more than 99 days of warm season per year (corresponding to 80% of the complete data in warm season), and with a coverage over 2 years. For locations with several discontinuous periods over 2 years, we separated data into different datasets according to each period. Second, we identified the type of missing values for each location. As timeseries data was used in this study, two types of missing manners were identified: the lack of measurements for certain variables and a small amount of missing values scattered during the study time period of each location. If it's unavailable of measurements for certain variable in one location, we excluded this location in the first step. For the second type of missing data, there is no suggestion that item being missing is correlated with death counts or any other predictors of interest. Previous studies using the same dataset also suggested that measurements being missing were mainly caused by data logging errors or abnormal operations of the monitoring equipment 3 . Thus, they are more likely to go missing at random. A summary of missing rates of temperature data and mortality data for each country (Table S1). Total mortality, mean temperature, minimum temperature and maximum temperature data have an overall missing rate of 0.13%, 0.82%, 0.84%, and 0.71%, respectively. The total missing rates were relatively small, which has nonsubstantial influence on our estimates. To provide more convincing evidence, we also perform sensitivity analyses by restricting our analyses to locations with complete data and by using complete data after imputation.  Figure S1. Country-specific exposure-response curve between mean temperature and mortality in warm season, stratified by quantiles of TV. Definition of abbreviations: TV = temperature variability; Q1 = the first quartile; Q2 = the second quartile; Q3 = the third quartile; Q4 = the fourth quartile.                   Definition of abbreviations: TV = temperature variability; Q1 = the first quartile; Q2 = the second quartile; Q3 = the third quartile; Q4 = the fourth quartile. a warm season used in the primary model. Definition of abbreviations: TV = temperature variability; Q1 = the first quartile; Q2 = the second quartile; Q3 = the third quartile; Q4 = the fourth quartile. Definition of abbreviations: TV = temperature variability; Q1 = the first quartile; Q2 = the second quartile; Q3 = the third quartile; Q4 = the fourth quartile.