All-cause mortality attributable to long-term changes in mean temperature and diurnal temperature variation in China: a nationwide quasi-experimental study

Previous studies have demonstrated an association between short-term exposure to ambient temperature and mortality. However, the long-term effects of elevated temperature and temperature variability on mortality have remained somewhat elusive in epidemiological studies. We conducted a comprehensive epidemiological study utilizing Chinese population census data from 2000 and 2010. Census-derived demographic and socioeconomic factors were paired with temperature data from the European Re-Analysis Land Dataset across 2823 counties. We employed a difference-in-difference approach to quantitatively examine the relationship between all-cause mortality and annual exposure to mean temperature and diurnal temperature range (DTR). Additionally, we evaluated the potential effects of socioeconomic and environmental covariate modifications on this relationship and calculated the attributable mortality. Lastly, we projected excess deaths attributable to annual temperature exposure under various shared socioeconomic pathways (SSPs, e.g. SSP126, SSP370, and SSP585). For each 1 °C rise in annual mean temperature and DTR, the mortality risk could increase by 6.12% (95% CI: 0.84%, 11.69%) and 7.72% (95% CI: 3.75%, 11.84%), respectively. Counties with high labor-force ratios and high NO2 and O3 concentrations appeared to be sensitive to the annual mean temperature and DTR. Climate warming from 2000 to 2010 may have resulted in 5.85 and 14.46 additional deaths per 10 000 people attributable to changes in annual mean temperature and DTR, respectively. The excess mortality related to changes in annual mean temperature and DTR is expected to increase in the future, with special attention warranted for long-term temperature changes in Southwest China. Our findings indicate that long-term mean temperature and DTR could significantly impact mortality rates. Given the spatial heterogeneity of increased mortality risk, the formulation of region-specific strategies to tackle climate change is crucial.


Introduction
As the repercussions of climate change become increasingly apparent, the health effects associated with annual mean temperature and temperature variability have become critical areas of public concern.
Numerous studies have revealed non-linear 'U-' and 'J-' shaped associations between short-term exposure to ambient temperature and a range of health outcomes, but time series models and case-crossover designs were employed in many of these investigations to control for potential confounding influences of long-term trends and air pollution [1][2][3].These approaches limit the ability to infer causation or adequately control for confounding factors on an annual scale [4].
To overcome these limitations, methodologies designed to control for unmeasured confounders, such as the difference-in-difference (DID) approach, have been proposed and applied in several recent environmental epidemiology studies [4][5][6].A classic DID design is applied using a quasi-experimental approach, with outcome measurement and comparison between pre-and post-intervention time points.This method is beneficial for the inference of causality via adjustment for unmeasured and omitted confounding factors [7,8].A DID or its variant design has been used in several studies to explore the health effects of long-term exposure to air pollution, with daily/yearly pollution concentrations across spatial units treated as measurements at different continuous time points [4,5,9].However, research in which a causal model/design has been applied to analyze the long-term health effects of ambient temperature remains scarce.To our knowledge, the DID approach has been used to examine the relationship between long-term mean temperature and mortality across multiple communities in China in only one study, which was confined to eastern China and thus did not yield findings representative of the entire population [10].
Numerous studies have demonstrated that shortterm exposure to varying temperatures can be associated with low temperature-mediated sympathetic nervous system excitation, dehydration induced by cold and heat, and systemic inflammatory responses triggered by heat stimulation.These factors contribute to the development of chronic cardiorespiratory diseases [11].However, epidemiological evidence for the effects of long-term temperature changes on all-cause mortality and the evolution of such risks over large temporal and spatial scales remains insufficient.Given the observed and expected increases in long-term temperature trends due to climate change, assessment of the long-term effects of high temperature and temperature variability is crucial to understand the public health impacts of climate change.
Beyond the annual mean temperature, temperature variability may contribute significantly to increased mortality risk.Evidence suggests that residents accustomed to typical temperatures are more vulnerable to harm from abnormal temperature fluctuations [12].The diurnal temperature range (DTR), defined as the difference between the maximum and minimum temperatures in a day, has been identified as a key metric of temperature variation (TV) [13].Although numerous studies have revealed positive associations of short-term increases in TV with morbidity and mortality [2,14,15], studies of the effects of long-term TV increments, particularly in developing countries such as China, are scarce.
As part of the latest framework for climate change scenarios, shared socioeconomic pathways (SSPs) were developed based on the traditional representative concentration pathways (RCPs); they incorporate a variety of socioeconomic assumptions such as population, economic growth, and education [16].Due to China's vast territory and large population, significant socioeconomic and demographic differences may exist among cities and even counties [5].Thus, the forecasting of future temperature-related mortality burdens based on the exposure-response relationship and population growth in Chinese counties is of paramount importance for the devising of climate change adaptation strategies and response measures.Some studies have involved the assessment of future risks related to the mean temperature or DTR, but many are limited by factors such as restricted study locations [17,18], the use of emission scenarios without socioeconomic assumptions (e.g.RCP scenarios) [19,20], disregard of future population changes [17,21], and the use of single scenarios/general circulation models (GCMs) for analyses [22,23], which complicate the generation of results that are accurate and comparable across the entire country.
In this study, we used a DID design and census data from 2000 and 2010 to estimate the associations of changes in annual exposure to ambient temperature and DTR with changes in all-cause mortality across 2823 counties.Additionally, we aimed to (1) evaluate the modifying effects of socioeconomic and environmental factors, (2) calculate changes in attributable deaths under a counterfactual framework, and (3) project future excess mortality attributable to changes in annual mean temperature and DTR.

Data collection 2.1.1. Population and mortality data
This study was conducted with publicly available Chinese census data from 2000 and 2010, sourced from the website of the National Bureau of Statistics (www.stats.gov.cn/).These datasets encompass various indicators, such as population size, mortality rates, and household living conditions, for countylevel units across the country.The Regulations on the National Population Census stipulate that Chinese citizens must cooperate with the population census and provide the required information accurately.Post-statistical surveys revealed underreporting rates for the 2000 and 2010 censuses of 0.12% and 1.81%, respectively, significantly below the internationally accepted standard of 3% [24].
Our study area encompassed 2823 county-level administrative units in mainland China (excluding Hong Kong, Macao, and Taiwan).Countylevel census data were geocoded to correspond with geographic information system shapefiles.As some administrative regions merged or split between 2000 and 2010, we used the national county-level administrative divisions in 2000 as the reference to calibrate the vector data on administrative divisions in 2010.In the matching of demographic data to county-level administrative units, the population of split administrative areas was presumed to be distributed evenly.The population was then estimated proportionately to the area of the region.
To facilitate comparability, we calculated the sex-age standardized mortality rate (SMR) for each county to standardize mortality counts, using the equation (equation ( 1)): where i, j, and k represent the administrative code, census year, and specific gender and age group, respectively, and D i,j and E i,j denote the actual and expected numbers of deaths in a county.SMR i,j = 1 indicates that the mortality risk in county i matches the national average.P i,j,k represents the population of a specific sex-age group in a county, and R k,j is the sex-age group mortality rate, denoted by the national average from the two census datasets.
To control potential confounding effects of demographic factors, we extracted and derived the following variables from the two censuses: (1) the urbanization rate (%), defined as the proportion of urban residents and calculated by dividing the number of urban residents by the total population; (2) the college education rate (%), represented as the percentage of the population that has completed a college education or higher; (3) the immigration rate (%), determined as the percentage of immigrants (population migrating from other places to the local county) and calculated by dividing the number of immigrant residents by the total population; (4) the mean age (years), defined as the population-weighted average age in a county; and (5) the dependency ratio (%), defined as the ratio of the population aged <15 and >65 years (non-labor force) to the number of people aged 15-65 years (labor force).Given that high ambient temperature often coincides with high electricity demand (e.g. for air conditioning use) [25], we also incorporated 1 × 1 km gridded electricity consumption data [electricity (10 8 , kW•h)] for both census years from a public database (https://figshare.com/) [26].Changes in these variables between 2000 and 2010 were integrated into the statistical model as covariates for confounding control.

Environmental variables
We procured hourly data for air and dew point temperatures 2 m above land surfaces in 2000 and 2010 from the European Re-Analysis Land Dataset (ERA5), which provides global coverage with a spatial resolution of 0.1 • .Stemming from high-resolution numerical integrations from the European Centre for Medium-Range Weather Forecasts land surface model, ERA5 generates variables describing the evolution of water and energy cycles over land [27].The extracted temperature data serve as an effective substitute for data from observation stations, adeptly addressing issues of spatial heterogeneity [28].After mapping the hourly data to mainland China, we calculated daily mean temperature and DTR, followed by annual average values.We then used an area-weighted average method to spatially align the annual mean temperature and DTR with county-level mortality data [29].Relative humidity was determined based on the 2 m air and dew point temperatures, using transformation formulas described elsewhere [30].
Fine particulate matter (PM 2.5 ) data with a resolution of 0.1 • were extracted from a published product offering annual global estimates for 1998-2016.With the use of satellite remote sensing and geographically weighted region evaluation, these estimates were found to be highly consistent with the cross-validation of monitor samples (R 2 = 0.81) [31].Such advanced PM 2.5 estimates have been utilized extensively in assessments of the impact of air pollution on health, and have been proven to be particularly suitable for areas or periods lacking regular monitoring networks [29,32,33].The area-weighted average approach was employed to estimate annual average concentrations for each county, as described previously.

Data for future projections
To facilitate future risk estimation, we collected daily simulated temperature data from 2015 to 2050 from the Inter-Sectoral Impact Model Intercomparison Project (www.isimip.org/).We considered three typical scenarios: SSP126 (radiation forcing stabilized at 2.6 W m −2 in 2100), SSP370 (radiation forcing stabilized at 7.0 W m −2 in 2100), and SSP585 (radiation forcing stabilized at 8.5 W m −2 in 2100).Each scenario encompasses five GCMs (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0, and UKES1-0-LL), including daily maximum, minimum, and mean temperatures at a global scale of 0.5 • grid resolution, validated through deviation correction and statistical downscaling processing [34].We extracted the simulated daily maximum, minimum, and mean temperatures from the grid cells corresponding to each county from 2015 to 2050, and then calculated the annual mean temperature and DTR for each county.
Gridded 1 km downscaled population data for China from 2015 to 2050 under SSP1 (sustainability), SSP3 (regional rivalry), and SSP5 (fossil-fueled development) were retrieved from the Socioeconomic Data and Applications Center (https://sedac.ciesin.columbia.edu/data/set/popdynamics-1-kmdownscaled-pop-base-year-projection-ssp-2000-2100-rev01).Population projections for each county were computed by summing the populations of grid cells within their boundaries.To account for the discrepancy between projected population and census data, we calculated a county-specific correction factor by dividing the 2000 SSP population data by the 2000 census data.This correction factor was assumed to be constant and was used to adjust the predicted population of a specific county under the three SSPs [35].

Study design and statistical model
A DID design with a log-linear regression model was adopted to explore the associations of changes in annual temperature (mean of daily means and DTRs) between 2000 and 2010 with changes in mortality over the same period among 2823 spatial units.This design enabled the filtering out of effects of unmeasured confounding (e.g.lifestyle or cultural) factors that rarely changed during the study period through the comparison of each sample at different year, thereby allowing the estimation of the true effects of ambient temperature.In this study, we compared the annual average variations of mean temperature/DTR and changes in mortality risk from 2000 to 2010 to obtain the long-term effect of temperature, as it not only includes the annual change in temperature, but also to some extent implies a decadal trend of TV.
We initiated our analysis by graphically examining the relationship between annual average temperature/DTR and mortality risk.A thin plate spline function with four degrees of freedom was fitted for the annual mean temperature and DTR separately, and relative humidity was included in the basic log-linear model (figures 1(B) and (D)).The statistical significance of the smoothing term of temperature/DTR was tested in the basic model, and the Bayesian information criterion (BIC) was used to test the goodness of fit of linear versus non-linear models.
A multivariate model containing all potential covariates was subsequently built.Specifically, we specified the response distribution as 'Gaussian' and the link function as 'log' , which describes the distribution of SMR in each county and how its expected value (after a log transformation for SMR) relates to corresponding predictors (mean temperature or DTR) in a linear way.The fully adjusted model for annual average temperature is represented as equation (2).For the DTR-mortality model, an additional thin plate spline term of annual average temperature was incorporated into equation (2), given previous research demonstrating that mean temperature may play an important confounding role in the relationship between DTR and human health [36,37]: where i denotes the county code; x symbolizes the change in annual average temperature between 2000 and 2010; β is the intercept (β 0 ) or slopes (β 1−10 ) for linear terms; s 1 () is a two-dimensional (latitude and longitude) smoothing term that controls for geographically varying unmeasured confounders [29]; s 2 () represents a smoothing function of relative humidity, with three degrees of freedom.The risk ratio percentage change [(RR-1) × 100%] for each 1 • C increase in x was quantified to display the effects' magnitudes.Diagnostic graphs (e.g.Q-Q plots, residual histograms) were used to assess the model's goodness of fit.
We also conducted stratified analyses to examine effect modifications, considering various socioeconomic and environmental factors, such as the urbanization rate, dependency ratio, education level, and levels of PM 2.5 , O 3 , and NO 2 .Equation (3) was used to test the significance of heterogeneity between strata: where β 1 and β 2 represent the coefficients from the model in each stratum, and se 1 and se 2 denote the corresponding standard errors.

Attributable risk calculation
We utilized the calculated associations to quantify the changes (∆D i ) in deaths attributable to variations in annual mean temperature and DTR between the baseline scenario (established in 2000) and a counterfactual scenario (considering temperature changes up to 2010 only): where D i,2000 denotes the number of deaths in county i in the year 2000 and ∆temp i signifies changes in annual temperature or DTR between 2000 and 2010.
We compared predicted and actual death counts in each city for 2010 to ascertain model accuracy.

Future risk projection
To predict future attributable mortality, we used the year 2000 as a baseline and estimated attributable deaths annually from 2015 to 2050: where POP s,i,year represents the adjusted population of county i in a given year under scenario s; MR i,2000 is the mortality rate in county i from the census year 2000, used as a proxy for the mortality rate in all projection years; and ∆temp s,i,year represents the temperature change from 2000 to a specific future year.
To facilitate result visualization and interpretation, attributable deaths in each city were aggregated from the corresponding county-level counts.

Sensitivity analyses
In addition to model diagnosis, sensitivity analyses were performed to test the robustness of our findings.These included: (1) including the values of various socio-economic variables in 2000 on the basis of the basic model; (2) alternating the degrees of freedom for annual mean temperature and DTR from 3 to 6; (3) alternating the degrees of freedom for relative humidity from 3 to 6; and (4) evaluating lag effects using temperature data from the previous year (1999 and 2009) and the average of current and previous year temperatures.All analyses were executed using R software (version 4.2.3;R Foundation for Statistical Computing).The 'mgcv' and 'mgcViz' packages were employed to construct the log-linear model, and the 'raster' and 'sf ' packages were used to handle spatial grid data.Two-tailed P values < 0.05 were considered to be significant.

Descriptive statistics
The study included a total of 2823 spatial units across mainland China.Statistical overviews of environmental and demographic variables across all counties are provided in tables 1 and S1.Between 2000 and 2010, the average mean temperature, DTR, and SMR rose by 0.28 • C, 0.32 • C, and 0.09, respectively.Average concentrations of PM 2.5 , O 3 , and NO 2 also increased by 1.40, 5.95, and 0.64 ug m −3 , respectively (table 1).The total population and actual mortality number increased (from 436 000 to 467 000, and 2570 to 2610, respectively).Increases were also observed in the urbanization rate, education rate, immigration rate, mean age, and electricity consumption.In contrast, the dependency ratio decreased (table S1).The spatial distributions of changes in the annual mean temperature and DTR between 2000 and 2010 are depicted in figures 1(A) and (C).

Epidemiological associations
Figures 1(B) and (D) shows that the mortality risk in all counties increased approximately linearly with increasing annual mean temperature (F for trend = 8.05, P < 0.01) and DTR (F for trend = 4.41, P < 0.01).For DTR, the BIC shows that the linear model (362.38) is slightly better than the nonlinear model (371.94).The annual mean temperature change varied between −0.3 • C and 0.7 • C and the DTR change fluctuated between −0.8 • C and 1 • C in most counties.Figure 2 shows associations of annual mean temperature and DTR with census-based mortality in China, with stepwise adjustments for confounding variables.The primary model indicates positive risk alterations for the annual mean temperature and DTR, despite CIs straddling zero.All multivariate models present significant and sturdy statistical associations of escalated mean temperature and DTR with the mortality risk.In the fully adjusted model, the risk of death increased by 6.12% (95% CI: 0.84%, 11.69%) and 7.72% (95% CI: 3.75%, 11.84%), respectively, with each 1 • C increase in the annual mean temperature and DTR.

Effect modification assessment results
Table 2 shows the mortality rate percentage changes corresponding to 1 • C increases in the annual mean temperature and DTR, stratified by socioeconomic and environmental factors.The prolonged effects of annual mean temperature and DTR were more intense in counties with smaller labor forces, and increased concentrations of NO 2 and O 3 , than in their counterparts.For the DTR, marginally greater effect modifications were also observed in low urbanization level, low education level and higher concentrations of PM 2.5 .

Counterfactual analysis results
The counterfactual analyses indicated that the increase in the annual mean temperature from 2000 to 2010 would result in 719 485 attributable deaths nationwide, equating to an increase of 5.85 deaths per 10 000 individuals.For the DTR, the attributable deaths and death change rate were 1778 090 and 14.46 per 10 000, respectively.We identified spatial heterogeneity in attributable deaths, aligning approximately with temperature changes (figure S1): annual mean temperature-related deaths were concentrated predominantly in the southwestern region and DTRrelated deaths increased primarily in the southwestern, central, and southern regions.Considering a more significant increase in excess mortality compared to the other two scenarios, we focused on the predicted results under the SSP585 (for mean temperature) and SSP126 scenario (for DTR).The estimated average deaths in cities attributable to the annual mean temperature show an increasing trend from 2015 to 2050, peaking under SSP585 (figure 3(B)).Compared with the base year (2000), the attributable deaths under SSP585 would increase by 15

Model diagnostics and sensitivity
Table S2 shows that the estimated risks do not change significantly after incorporating the initial values of some socio-economic factors.Figure S2 displays the diagnostic diagram for the multivariate model, indicating that the residuals essentially conform to a normal distribution, the error term implied in the regression conforms to the assumed Gaussian distribution.Figure S3 shows consistent alignment between county-level projected deaths from the counterfactual analysis and actual deaths reported in the 2010 census, supplemented by results from several typical cities. Figures S6-S8 show that the estimated values exhibited minimal variation following alterations in the degrees of freedom and lag settings, suggesting that the model is robust.

Discussion
This study, to our knowledge, is the first census-based research performed with a DID design to examine the causal effect of annual exposure to ambient temperature and DTR on all-cause mortality in China.
The annual mean temperature and DTR increased between 2000 and 2010, significantly increasing the death risk and number of attributable deaths nationwide.Future projection analyses suggested that the upward trend of the annual mean temperature will persist, causing a considerable number of excess deaths under the SSP585 scenario, particularly in the southwestern region.Although the future increase in the DTR is not pronounced, it will cause some excess deaths in certain areas of the southwestern, northwestern, and eastern regions under the SSP126 scenario.
Limited studies have assessed the long-term effects of temperature on mortality risk, and results were inconsistent due to differences in exposure indices and modeling strategies [38].For example, Li and Gu [39] employed a dynamic model to probe the historical  association between temperature exposure and mortality evolution in China, they observed that increases in the number of extreme hot days would significantly increase annual mortality rate: compared to days below 27 • C, one day with an average temperature above 27 • C would increase annual mortality rate by about 0.1%.Yitshak et al [40] examined the impact of annual exposure to mean temperature on the deaths of patients hospitalized with cardiorespiratory diseases, reporting 6.24% (95% CI: 5.93%, 6.54%) and 1.37% (95% CI: 1.28%, 1.47%) excess risks for respiratory admissions and 7.32% (95% CI: 6.68%, 7.96%) and 0.15% (95% CI: −0.04%, 0.34%) excess risks for ischemic stroke admissions in association with per-IQR (2.2 • C) increases in the 10th (8 • C) and 90th (10.2 • C) centiles of temperature, respectively.In our study, a positive relationship between the annual mean temperature and mortality risk was observed.However, another study revealed a negative association between the annual mean temperature and mortality risk in 364 Chinese communities [10].In Hu et al's study, the selected Chinese communities covered areas with death monitoring points (usually with good economic development), without considering the vast western and remote areas (usually with high population vulnerability), thus there was a possibility of underestimation of risk.In addition, although both studies utilized the idea of DID design for selfcomparison at different time points, there were still differences in many aspects, such as data source, mortality indicator, study period, and modeling strategy, which might be another reason for the inconsistent estimation.More county-scale, better-designed temperature-mortality studies were still needed in the future to validate our results.
While the underlying mechanism remains uncertain, a wealth of evidence from fields such as physiology, immunology, and behavioral science illustrates that abrupt changes in temperature can incite stress responses in the cardiovascular and respiratory systems, thereby increasing the risk of death [13,45].In particular, large fluctuations in the DTR may alter the heart rate, arterial blood pressure, oxygen intake, and provoke an inflammatory response, thereby intensifying the burden on the cardiovascular system and potentially triggering or worsening cardiovascular events [46].
Our stratified analyses suggested that counties with lower levels of urbanization, a higher labor force ratio, and lower education levels seem to experience an increased risk of all-cause mortality linked to annual exposure to mean temperature and temperature variability.These findings align with previous findings [44,47,48].For instance, Hu et al [47] explored the differences in temperature-mortality relationships between rural and urban regions of Zhejiang Province, China.They discovered that rural residents were more vulnerable to extreme heat and cold than were their urban counterparts, attributing this susceptibility to lower socioeconomic status and substandard living conditions, such as extended outdoor work hours, inadequate access to medical resources, and inability to afford air conditioning, in rural areas [47].
Compared with their more-educated counterparts, individuals with lower education levels are often employed in environments that offer fewer protective measures and higher toxicity levels [49].Air pollutants (PM 2.5 , NO 2 , and O 3 ) appear to contribute significantly to aggregating temperature-related health risks, having a larger estimated impact than do socioeconomic factors.This observation reinforces existing evidence suggesting that the intricate interplay between air pollution and climate change, particularly in the context of global warming, poses a greater threat to human health [49].The concept of a 'climate penalty' may also amplify health damage in individuals exposed to higher levels of air pollution [50].
In an effort to understand the number of deaths attributable solely to the temperature change between 2000 and 2010, we conducted counterfactual analyses.The results indicate that a substantial number of deaths could be ascribed to increases in the annual temperature and DTR at a national level, underlining the immediate need for the development of regionspecific climate change adaptation strategies.
We illustrated the anticipated future variations in the annual mean temperature and DTR.The mean temperature demonstrates a yearly increasing trend, with a corresponding increase in predicted deaths, particularly under the SSP585 scenario.These observations align with previous findings [35,51].Contrary to the annual mean temperature, the DTR displays a decreasing trend in all scenarios, except that it increases slightly under SSP126.These results reflect inequal increases in future annual maximum and minimum temperatures.This phenomenon, tied to a more rapid increase in nighttime than daytime temperatures, has been reported previously and is believed to be driven by factors such as precipitation, cloud cover, and net longwave radiation [52].However, it is not absolute, as predictions for some areas of Europe and Central and South America indicate future DTR increases, suggesting that future changes will depend largely on local climate characteristics [53].
For the DTR, we observed a slight increase accompanied by minor fluctuations in city-averaged excess deaths relative to 2000.To date, research on DTRrelated projections remains scarce; however, existing findings suggest a general increase in DTR-related excess mortality with varying patterns among regions and countries [22,53].
Spatially, the southwestern region (including provinces such as Xinjiang, Qinghai, Tibet, and Sichuan) seems to contribute significantly to annual mean temperature-and DTR-related excess mortality.This contribution may be attributed to factors such as high-altitude living conditions, lack of selfprotection awareness, and the inadequacy of adaptive measures for extreme weather [54,55].Population density in the eastern region (e.g.Jiangxi and Zhejiang provinces), combined with residents' adaptability to warm weather, may explain the increased vulnerability to significant DTR variations in this region [22].
Unlike this study, which included nearly all counties across China, previous large-scale research focused predominantly on surveillance sites in eastern China, while the inclusion of monitoring points in western region was insufficient [10,22].Additionally, we must acknowledge that socioeconomic structures can vary greatly among cities and even counties [5], whereas some counties may have good economic development and sufficient resources to cope with adverse weather, others, especially those in rural and remote areas, may lack effective response measures [56].Thus, further scientific evidence from western China is needed to highlight the vulnerability of smaller geographical units and to guide the implementation of adaptive strategies and measures.
This study, while providing important insights, has a few limitations that should be acknowledged.
(1) Due to the constraints of census data, we were unable to conduct a health assessment based on stratified seasons, as has been done in many previous studies [10,12].The high mortality rate in a certain county may be not only related to annual temperature changes, but also related to some extreme weather events that occur in a specific season (especially in summer and winter), such as floods, wildfires, and cold waves, which may cause a sudden increase in mortality rate in a short period of time.In addition, some infectious diseases (such as influenza) may also contribute to the increase in local mortality rate.However, due to the fact that census data can only provide variable information on an annual scale rather than a seasonal scale, we were unable to adjust for these seasonal confounding factors in the model, thus affecting the accuracy of the estimation to some extent.More season-stratified studies are warranted in the future.(2) Our future risk assessment presupposed that people's adaptability to climate change (i.e. the exposure-response relationship) remains constant.However, socioeconomic development and the widespread use of air conditioning could significantly mitigate the health damage from global warming [57], introducing another element of uncertainty.Therefore, the results should be interpreted with caution.(3) Although the SMR indicator accounts for the population's age structure, data limitations prevented us from completely discounting the impact of future population aging.(4) As the study subjects were counties, rather than individuals, unavoidable misclassification errors of environmental exposure occurred.
In summary, this study offers significant insights into the long-term effects of the mean temperature and DTR on all-cause mortality in China.Utilizing a DID approach, we determined that increases in the annual mean temperature and DTR could have significant impacts on the mortality risk, as they led to a marked increase in attributable deaths between 2000 and 2010.Certain socioeconomic and environmental factors might modify the association between annual mean temperature/DTR and mortality.Future temperature-related excess mortality may increase under unmitigated climate scenarios, particularly in the southwestern region.These findings would provide further epidemiological evidence on the health effects of long-term temperature exposure, more and better designed studies are needed in the future to provide more conclusive evidence.
) and (C) shows the trends in the alteration of the annual mean temperature and DTR under the three SSP scenarios relative to 2000.The annual mean temperature displays a notable upward trend from 2015 to 2050, with a drastic increase

Figure 1 .
Figure 1.County-level annual mean temperature and DTR change and their relationship with mortality.A and C show the spatial change of annual mean temperature and DTR (diurnal temperature range) from 2000 to 2010.B and D show the exposure-response relationship from basic model (adjusting for relative humidity).The blue dots in B and D corresponds to all counties in our study, the red curve represents the overall association.The black lines at the bottom of B and D represent the probability density of the sample distribution.

Figure 2 .
Figure 2. The estimated association between annual mean temperature, DTR and census-based all-cause mortality in China, by adjusting for multiple confounding factors.The associations are shown as percent change (%) in mortality risk corresponding to one Celsius degree increment in annual mean temperature (Tmean) and DTR change.

Figure 3 .
Figure 3.The future changes in annual mean temperature, DTR, and associated attributable deaths under three typical scenarios, from 2015 to 2050, compared with 2000.A and C show the national average changes of annual mean temperature and DTR, compared to 2000, respectively; B and D show the national averaged changes in annual mean temperature-and DTR-related deaths, respectively.

Table 1 .
Descriptive statistics of annual average of meteorological, pollution, and mortality variables.
Note: Tmean, Tmax, T min , DTR, and RH represents annual mean temperature, maximum temperature, minimum temperature, diurnal temperature range, and relative humidity, respectively.SMR represents sex-age standardized mortality rate.

Table 2 .
Percentage change (%, 95% CI) of all-cause mortality risk per 1 • C increase in annual mean temperature and DTR by socioeconomic and environmental status.