Spatial-temporal pattern of tuberculosis mortality in China and its relationship with long-term PM2.5 exposure based on a causal inference approach

Evidence on the spatial-temporal distribution of tuberculosis (TB) mortality across China and its relationship with long-term particulate matter (PM2.5) exposure is limited. We aimed to address significant gaps in our understanding of the spatial-temporal clustering patterns of TB mortality in China and provide evidence for its causal links with long-term PM2.5 exposure. Annual pulmonary TB mortality, PM2.5 concentrations, and socioeconomic factors for provinces in mainland China between 2004 and 2017 were obtained. Turning points in the temporal trend and spatial clustering patterns of pulmonary TB mortality were identified. A difference-in-differences causal inference approach was applied to estimate the long-term effect of PM2.5 exposure on the mortality. The average annual percent change of pulmonary TB mortality in China was −2.5% (95% CI: −5.6%, 0.7%), with an 11.1% annual increase in the Northwest since 2012 (P= 0.029). The hot and cold spots, determined by the local Moran’s I index, were all located in northern China, where Xinjiang in the Northwest had the highest mortality across the study period. We found a significant association between long-term PM2.5 exposure and pulmonary TB mortality, with percent increase risk of mortality (IR%) being 0.74% (95 CI%, 0.04%, 1.45%) for 1 µg/m3 increase in PM2.5 concentration. This association varied across multiple socioeconomic groups, with the highest IR% in provinces with lower level of latitude (IR% = 0.83%, 95% CI: 0.01%, 1.65%), lower quartile of gross domestic product (IR% = 1.01%, 95% CI: 0.23%, 1.80%) or higher proportion (⩾14%) of people >65 years of age (IR% = 1.24%, 95% CI: 0.44%, 2.04%). Comprehensive sensitivity analyses showed a robust adverse effect of long-term PM2.5 exposure on pulmonary TB mortality. Attention needs to be paid to the rising trend of pulmonary TB mortality in Northwest China. Our study provides the stable evidence to date of the causal association between long-term PM2.5 exposure and the risk of death from pulmonary TB, especially in low-altitude, underdeveloped, and aged provinces.


Introduction
Tuberculosis (TB) is a major infectious disease and one of the leading causes of deaths worldwide, with 10 million infections and 1.5 million deaths in 2020 [1]. Although extensive efforts have been made in China to reduce the burden of TB, there were still 0.67 million infected people and almost 2000 deaths from TB in 2020 [2]. China is still among the ten countries with the highest TB burden [1].
The spatial and temporal characteristics of TB incidence in China have been well documented [3][4][5][6]. The decreasing trends of TB incidence in mainland China in the last two decades were reported in many studies [3,4], and some information on the spatialtemporal distribution of TB mortality is also available [7]. However, gaps remain in our understanding of the spatial clustering patterns of TB mortality across the country. Changes in the spatial clustering over the past two decades that are critical for decision-making and optimizing health resources allocation, are also unclear.
Fine particulate matter (PM 2.5 ) exposure is considered an important environmental predictor driving TB morbidity [8,9]. PM 2.5 has been found to be associated with increased expression of inflammatory cytokine receptors to activate IL-1R and IL-6R-mediated signaling pathways in lung cells [10]. It may also dysregulate a variety of metabolic enzyme activities, including superoxide dismutase and nitric oxide synthase [11]. Even limited exposures to PM 2.5 were found to be associated with increased lung oxidative stress and inflammation [12]. These mechanisms may result in a higher vulnerability of people to M. tuberculosis infection, related disease progression, and death.
Although the micro-mechanisms of the damage caused by PM 2.5 are well understood, the epidemiological evidence of increased TB mortality associated with long-term PM 2.5 exposure is limited and conflicting. Only two studies reported long-term PM 2.5 -TB mortality associations in China. A cohort study in Shanghai suggested that per interquartile range (IQR) increase in annual PM 2.5 concentrations (2.06 µg/m 3 ) was significantly associated with mortality from TB (HR = 1.46, 95% CI: 1.15, 1.85) [13]. In contrast, a time-series analysis in Shandong did not observe any significant association between the long-term PM 2. 5 and TB mortality (IR% = 0.00%, 95% CI: −0.02, 0.02) [14]. This is existing evidence contributing to our understanding of the potential health impact of long-term PM 2.5 exposure. However, previous studies generally were based on the association assessment by conventional statistical models. Existing effects estimates may be biased by the remaining confounding issue and result in contradictory conclusions [15]. Recently, some causal inference approaches, which can better control confounding by comparing observed and counterfactual outcomes, are available and have been well used in estimating the effect of PM 2.5 on other outcomes [16][17][18]. In addition, the aforementioned studies were based on local data, limiting researchers' capability to extrapolate their findings to other parts of China, or capture the nationwide scenarios of the long-term PM 2.5 -TB mortality relationship.
To fill these knowledge gaps, we first explored the temporal and spatial distribution and clustering pattern of TB mortality in mainland China from 2004 through 2017. Then, we assessed the association between long-term exposure to PM 2.5 and TB mortality using a causal inference approach named difference-in-differences (DID) design. We also assessed whether any associations were modified by gross domestic product (GDP), population aging degree, and other socioeconomic factors.

Mortality data
Because PM 2.5 damages the respiratory system first with most TB patients dying of pulmonary TB, we used the pulmonary TB mortality among the entire population as the outcome in this study. The province-specific pulmonary TB mortality between 2004 and 2017 were collected from the Public Health Science Data Center operated by the Chinese Center for Diseases Control and Prevention [19]. TB case and death were determined according to standard clinical guidelines and all patients died from pulmonary TB were coded in the International Classification of Disease, Ten Revision (ICD-10: A15-A19). The data in this platform were extracted from the surveillance system for infectious diseases, which regularly integrates and manages the pulmonary TB data separately reported by health care facilities across the country.

Exposure assessment and socioeconomic factors
Annual PM 2.5 concentrations at the 1 km × 1 km grid scale over China during 2004-2017 were obtained from the NASA Socioeconomic Data and Applications Center [20]. The dataset consists of annual concentrations (µg/m 3 ) of ground-level PM 2.5 derived from a complex estimation procedure. First, the GEOS-Chem chemical transport model was used to relate aerosol optical depth retrievals from multiple satellite algorithms to near-surface PM 2.5 concentration. Then, geographically weighted regression was utilized with ground-based measurements to adjust the initial satellite-derived values [21]. We computed the annual province-specific mean PM 2.5 concentrations by averaging all grid data in the province. Compared with ground site data, these mean concentrations from gridded data can better represent the average population exposure levels for different provinces. Following the definition in numerous studies in this field [13,14,17,18], we also define the annual average PM 2.5 concentrations as 'long-term' exposures.
Daily temperature gridded data (1 × 1 km) were obtained from the Berkeley Earth Foundation dataset (BERKEARTH) in the Climate China Store maintained by the European Center for Medium-range Weather Forecasts [22]. BERKEARTH merges temperature records from 16 archives into a single coherent dataset, and we obtained the daily mean temperature for each province by averaging over all grids in the province. Then, the mean and standard deviation of summer and winter temperatures in each province were estimated using daily province-specific temperatures. Specifically, we defined summer as June, July and August, and winter as January, February and December.
We also collected socioeconomic factors from the Statistical Yearbook for each province integrated by the China Economic and Social Big Data Research Platform [23], including GDP per capita, gender ratio (the ratio of male over female population), population density (the number of people per square meter of area), elderly proportion (the proportion of people >65 years), education rate (the proportion of people with elementary education and beyond), college rate (the proportion of people with a college degree and beyond), passenger turnover (the summation of passengers travelling between provinces multiplied by the average traveling distance), health personnel allocation (the number of health technicians per thousand people) and urbanization rate (the proportion of urban population).

Statistical analysis
We performed a joinpoint regression analysis [3] to explore the temporal trend of pulmonary TB mortality and identify the points in time when a significant change in the slope occurred for mainland China overall as well as for each of the seven geographic regions [24] in the country (figure S1). The annual percent change (APC) for each identified time period and weighted average annual percent change (AAPC) for the entire study period were estimated to represent the rate of change in slope.
The global spatial autocorrelation of pulmonary TB mortality in each year was investigated using Global Moran's I statistics. Spatial clustering analysis was performed via Anselin's Local Moran's I to identify high-risk areas. Both global and local Moran's I range from −1 to 1 and could be interpreted similarly. An index close to −1 (or 1) indicates a high dissimilarity (or similarity) across the spatial area in the mortality of pulmonary TB, while an index around 0 represents no spatial autocorrelation. We further classified the areas into four types according to local Moran's I: (1) high-high cluster or hot spot (i.e. aggregated areas with the highest mortality), indicating the high cluster areas being surrounded by other high cluster areas; (2) high-low cluster, indicating the high cluster areas being surrounded by low cluster areas; (3) low-high cluster, and (4) lowlow cluster or cold spot (i.e. aggregated areas with the lowest mortality).
We applied a modified DID approach to estimate the effect of long-term PM 2.5 exposure on pulmonary TB mortality, which has been successfully applied in estimating the association between PM 2.5 and other outcomes [17,25,26]. The essence of this design is that temporally stable confounders that have not been measured can be controlled by comparing the same population to itself at different time points [18,27]. Typically, the DID model is written as follows: where Y a s,t is the outcome in area s and year t for exposure equal to a; Z s denotes spatial confounders that have very small (no) variations over time; U t reflects temporal confounders with very small variations among areas; and W s,t represents confounders that vary across areas and time. Further, we can get the difference of differences in outcomes over time between location s and s ′ as follows: where b and b ′ are the exposures in area s ′ at times t and t − 1, respectively. If the changes in W s,t between time among areas are the same, the second term in the equation will be near to zero. Therefore, the known and unknown confounders have been controlled under proper assumptions in this approach. In this study, we assessed this assumption by calculating the relative change rate differences of annual concentrations of PM 2.5 , temperature, and GDP per capita [17], and estimated the Pearson correlations of the relative change differences between PM 2.5 and three potential confounders.
In this study, we fit a model which has been generalized to multiple areas and time periods by Wang and colleagues [26]. The function is written as follows: where Y s,t is the number of pulmonary TB death in province s and year t; I s and I t are the indicator variables for provinces and years, respectively; PM s,t is the annual PM 2.5 concentration in province s and year t; Tsum s,t , Twin s,t and SD T sum s,t , SD T win s,t represent the mean and standard deviation of summer and winter temperatures, respectively. Due to the rapid economic development in China over the past two decades, we further included GDP per capita to reflect the economic development that might vary across both provinces and years. ln(Pop s,t ) is an offset term of background population size for the Quasi-Poisson distributed outcome variables. Through comparisons between provinces and years, unmeasured spatial and temporal confounders can be controlled within this DID design. To properly estimate the coefficient of PM 2.5 , the conditional quasi-Poisson model was used to perform the DID design. Compared with the unconditional Poisson model, conditional models are able to account for autocorrelation and overdispersion existing in time-series data [28].
To evaluate the disparity of PM 2.5 -mortality association across different groups, we converted the continuous environmental and socioeconomic variables into categorical variables. Because the distribution of population significantly varies across different altitude zones in China and human activities have a great impact on the emission of air pollution, including PM 2.5 , it is very likely that the compositions as well as the toxicity of PM 2.5 vary with the altitude. Thus, we did a modification analysis to investigate the potential moderation of altitude. The provinces were divided into three groups according to the terrain and altitude (table S1). The first altitudinal gradient group covers the Qinghai-Tibet Plateau, with altitudes higher than 4000 m. The second group includes most basins and mountains in mainland China and altitudes usually ranges from 1000 to 2000 m. The third group is located in eastern China and has an average altitude lower than 500 m. According to the World Health Organization, we classified provinces as not aging (<7% of people >65 years of age), aging (7%-14% of people >65 years of age), and aged (⩾14% of people >65 years of age) society [29]. Other factors, including GDP per capita, population density, education rate, college rate, passenger turnover, health personnel allocation and urbanization rate, were divided into four groups by quartile (Q1, Q2, Q3, and Q4) in each year. We used two approaches to examine effect modification by these factors: (1) stratification analysis for altitudinal gradient groups, which would not change over time across provinces; and (2) added product terms of other time-varying factors and PM 2.5 (e.g. Q4 of education rate × PM 2.5 ) into the model.
We performed sensitivity analyses to test the robustness of the results. We modeled summer and winter temperatures using natural splines with 3 and 4 degrees of freedom, respectively [17], to assess the nonlinear effect of temperature on the association between PM 2.5 and pulmonary TB mortality. To examine potential lagged effects of PM 2.5 concentrations, we used the moving average of PM 2.5 in the current year and 1-2 years before (lag 0-1 and lag 0-2) to refit the main model. We applied a meta-regression to check the statistical differences among the effect estimates of the primary model and sensitivity models [30,31]. The effect estimates were modeled against the meta-predictor (e.g. a categorical variable with five levels, which can indicate different models), then the meta-regression would provide P-values for difference using a likelihood ratio test [32].
All results were expressed as percent increase risk of mortality (IR%) with 95% confidence intervals (95% CI) per 1 µg/m 3 increase in annual PM 2.5 concentration. The joinpoint analysis was done with Joinpoint software (version 4.9.0.1). The maps and results of the spatial analysis were generated using ArcGIS (version 10.5). The statistical models were developed in R (version 4.1.0) with the package 'gnm' for the main analysis and 'mvmeta' for the meta-analysis. A P-value < 0.10 was considered as broader significant and the P-value < 0.05 was statistically significant.  Table 2 shows the temporal trends of pulmonary TB mortality in different areas across mainland China. We did not find any change point over the study period in the country (overall AAPC: −2.5%, 95% CI: −5.6%, 0.7%). After 2006, a significant decrease in mortality in Central China occurred (APC: −6.9%, P = 0.001). In North China, the mortality continuously decreased with an AAPC of −5.8% (95% CI: −9.0%, −2.4%). In contrast, the mortality in Northwest China significantly increased after 2012 with an APC of 11.1% (P = 0.029).

Temporal trend and spatial distribution of mortality
The results of global and local spatial autocorrelation analyses are shown in table 3 and figure 1. Except     ). An adverse effect of long-term exposure to PM 2.5 in the third gradient group was observed, with IR% of 0.83% (95% CI: 0.01%, 1.65%). We found that there was a monotonic trend in the effect of long-term PM 2.5 exposure under different levels of GDP per capita, aging degree, and passenger turnover (figure 2). We found a higher effect of PM 2.5 in provinces with lower GDP per capita, with IR% decreasing from 1.01% (0.23%, 1.80%) in Q1 to 0.44% (−0.35%, 1.22%) in Q4. Moreover, we observed a higher and significant IR% of TB mortality following long-term PM 2.5 exposure in aged provinces (1.24%; 95% CI: 0.44%, 2.04%). Besides, the effect in Q1 passenger turnover group was 0.40% (−0.54%, 1.36%) with no significance, but the effect in Q4 passenger turnover group was significant with an IR% of 0.95% (0.21%, 1.68%).
The results of the sensitivity analyses can be found in supplemental table S2. Refitting the confounding effect of temperatures nonlinearly did not substantially alter the effect estimates with very similar IR% of 0.78% (95% CI: 0.06%, 1.50%) for 3 degrees of freedom and 0.74% (95% CI: 0.02%, 1.46%) for 4 degrees of freedom, respectively. Although the effect estimates for moving averages of PM 2.5 were generally larger than the main results, the differences were not statistically significant (both P > 0.05).

Discussion
This study systematically explored the spatial and temporal distribution of pulmonary TB mortality in mainland China during 2004-2017. This is the first study to estimate the risk of pulmonary TB mortality associated with long-term PM 2.5 exposure in China using a causal inference method.
We found that the pulmonary TB mortality in China decreased from 2004 to 2017, which is consistent with previous studies. Wang et al [3] suggested that TB mortality in China significantly decreased from 2004 to 2019. Zhu et al [4] also found a decreasing trend in TB mortality from 1990 to 2015. However, the similar trend in our study was not significant. A possible reason could be the difference in the definition of mortality between this and previous studies. We only accounted for pulmonary TB, which could result in a lower statistical power due to a smaller number of deaths. Previous studies focused on all TB cases, in which the declining trend of nonpulmonary TB mortality may be more pronounced. More data would be needed to clarify this issue.
Although no significant break point was observed during 2015-2017 in mainland China, the mortality in Northwest China had significantly increased since 2012. When focusing on the temporal trend of pulmonary TB mortality among provinces in Northwest China (figure S2), we found that the mortality mainly increased in Xinjiang Province. The increase might be related to the increased incidence of TB, which would increase the risk of death from pulmonary TB due to more prevalent patients. The TB incidence in Xinjiang Province continued to rise after 2012 [7], which could partially explain the increasing trend of pulmonary TB mortality.
We also found that the pulmonary TB mortality in Xinjiang remained high during the study period. Besides the high pulmonary TB incidence and prevalence in this area which resulted in a large number of TB patients [5,7,33], the high proportion of ethnical minorities, mainly Uyghurs, may be another important reason [34]. It was reported that minorities in Xinjiang had a higher prevalence of drug-resistant TB than the Han Chinese [35], where unsuccessfully treated cases may be suffering from a higher mortality risk. Moreover, it was reported that rural residents in the Southern region of Xinjiang were more likely to have poor access to health care [7], and there were more than 70% Uyghurs living in these areas, which could result in a low successful rate of Directly Observed Treatment, Short course among this group [36]. Previous studies also suggested that cultural norms, religious practices, and social determinants of health, such as education, poverty, also contributed to the limited development and high mortality burden in Xinjiang [5,7]. However, more data are needed to further understand the causes of high pulmonary TB mortality in this area.
We found the hot and cold spots were all located in northern China. Although existing studies have not reported any spatial pattern of TB mortality, a previous study reported similar spatial patterns of TB incidence [37]. The similarity in the spatial clustering characteristics of (pulmonary) TB incidence and mortality suggested that some common factors underpin the spatial distribution, such as the distribution of poverty [38,39]. As mentioned above, the health system and preventive care management may be underdeveloped in provinces with low economic levels. In our study, all provinces in hot spots had lower GDP than provinces with cold spots, indicating that increasing investment in TB prevention and control in underdeveloped provinces is still very important.
Consistent with previous cohort studies [13], we found an adverse effect of long-term PM 2.5 exposure on pulmonary TB mortality. A Chinese TB cohort study estimated an adjusted HR of 1.20 (95% CI: 1.07, 1.35) for TB mortality per 1 µg/m 3 long-term exposure to PM 2.5 [13]. We found that each 1 µg/m 3 increase in annual PM 2.5 concentration was associated with a 0.74% increase in the risk of pulmonary TB death. And our assessment of the assumption for the modified DID suggested that this effect could be considered as a causal effect (table S3). Both evidence indicated the adverse effect of long-term PM 2.5 exposure on TB mortality. Despite typical inflammatory and oxidative stress responses to long-term exposure, there is specific PM 2.5 -induced pathogenesis suggested for pulmonary TB [40]. Long-term PM 2.5 exposure could decrease the expression of hemoglobin subunit d2 and 3 in type II alveolar epithelial cells, which could increase the risk of M. TB proliferation [41]. This would produce more severe pulmonary TB cases after proliferation and therefore more deaths. It was worth noting that the P value for the effect estimate of long-term PM 2.5 exposure in mainland China was 0.038 in our study, which was relatively close to 0.05. Thus, the interpretation of the results should be cautious, and more studies were needed to confirmed our findings.
We found the effect of long-term PM 2.5 exposure was significant in provinces with low altitudes. The floating population, migrants who do not stay in an area permanently and always move from place to place [42], may be an explanation. Provinces in the third gradient shared more than 40% of the floating population in the whole China [43]. There are enormous challenges in managing floating TB patients. Floating patients were more likely to delay treatment and abandon medical treatment, so the cure rate of TB among them was much lower than the permanent population [44,45]. Further, the compositions of PM 2.5 among different altitude gradients would be different. Most people in China reside in the provinces of the third gradient, so the traffic volume in this area is much higher than that of the other gradient provinces [2]. Traffic related emissions can substantially increase the proportion of elemental and black carbon in PM 2.5 , which has a stronger association with mortality [46,47].
We found that pulmonary TB patients living in provinces with lower per capita GDP levels and higher aging degrees were more sensitive to long-term PM 2.5 exposure. As mentioned above, poverty is an important factor influencing the susceptibility of TB patients [38,39], and provinces with lower per capita GDP levels would have poorer nutrition, poorer health services and medical resources. In our study, the overall Spearman correlation between per capita GDP and health personnel allocation was 0.86 (P < 0.05), and the slight downward trend in the modification effect of health personnel allocation could also support this view. Many studies have ascertained that PM 2.5 does more harm to the elderly [48,49]. Frailty is common among older people and can make them more vulnerable to PM 2.5 exposure [50]. In 2020, the proportion of people older than 65 in China had reached 13.5%, while the birth rate has fallen below 10% [2]. Thus, population aging will continue to increase the burden of TB deaths from PM 2.5 .
This study systematically explored the spatial and temporal distribution of pulmonary TB mortality in mainland China during 2004-2017. We used a causal inference method to estimate the risk of pulmonary TB mortality associated with long-term PM 2.5 exposure in China. There were several limitations in this study. We could not adjust individual-level factors (e.g. job, marital and smoking status), and potential ecological bias might exist because the study was done on the population level. Further, the statistical power might be limited because all the exposure and factors in this study were averaged values at the provincial level. Therefore, non-differential exposure misclassification would make the results conservative. However, we still found a significant adverse effect of long-term PM 2.5 exposure on pulmonary TB mortality. In addition to mortality, case fatality ration is also a meaningful statistic, representing the mortality specifically among the TB patients. However, data on the number of TB patients is not available. Further studies on this new topic would be helpful.

Conclusions
We found that the pulmonary TB mortality in China slowly decreased during 2004-2017. High-risk areas for pulmonary TB mortality were all in northern China, especially Xinjiang. Our study suggests that long-term PM 2.5 exposure increases the risk of death from pulmonary TB. Further, this association may vary among environmental and socioeconomic subgroups, including altitude, GDP per capita, the degree of aging in a provincial population. These findings suggest that strengthening economic investment and actively achieving healthy aging may help reduce pulmonary TB mortality following long-term PM 2.5 exposure.

Data availability statements
The data cannot be made publicly available upon publication due to legal restrictions preventing unrestricted public distribution. The data that support the findings of this study are available upon reasonable request from the authors.