Correlation of PM_2.5 pollution and daily mortality rate of cardiovascular and cerebrovascular diseases in Northeast China through PM_2.5 sources analysis

2.1 Study Area

This study focused on Harbin, a city in northeastern China located between 44°04’ to 46°40’N and 125°42’ to 130°10’E, covering an area of 53 076 square kilometers. Harbin experiences a temperate continental monsoon climate, characterized by long, cold winters and short, cool summers, with distinct seasonal changes. The study specifically examined air pollution in the main urban area, establishing two PM_2.5 component monitoring sites based on the annual wind direction and the dispersion of air pollutant.

According to the annual wind frequency rose chart, Daoli District is situated upwind of Harbin, while Xiangfang District is located downwind, where numerous heavy industry factories are concentrated. Consequently, Daoli District was selected as a light pollution area, and Xiangfang District as a heavy pollution area.

2.2 Meteorological factors and ambient air pollutants

The daily meteorological data for 2016 to 2018 were obtained from the Harbin Meteorological Bureau. These data includes average temperature (°C), relative humidity (%), average atmospheric pressure (hpa), precipitation (mm), average wind speed (m/s), and sunshine hours (h). The air pollution data were sourced from the National Urban Air Quality Real-Time Release Platform. Harbin has a total of 11 air pollution monitoring stations covering the entire main urban area. The primary environmental air pollutants measured include PM_2.5(μg/m³), NO₂(μg/m³) and SO₂(μg/m³), O₃(μg/m³), and CO (mg/m³).

2.3 Composition of PM_2.5

Based on wind direction and population density factors in Harbin, Daoli District and Xiangfang District were selected as sampling sites for PM_2.5 components. A particle sampler with a flow rate of 100 L/min, equipped with quartz fiber and glass fiber filter membranes, was used to collect PM_2.5 samples. The concentrations of water-soluble ions, heavy metals, and PAHs in the PM_2.5 samples were determined.

2.4 Outcome measure

The daily mortality data from 2016 to 2018 were obtained from the Harbin Center for Disease Control and Prevention. All deaths related to cardiovascular and cerebrovascular diseases (CCVD) were coded according to the International Classification of Diseases, 10th Edition (ICD 10). This study was approved by the Ethics Committee of the Harbin Center for Disease Control and Prevention (Approval No. 2023015).

2.5 Statistics analysis

The mean and standard deviation were used to describe quantitative data with a normal distribution, while the median and interquartile range were used for data with a skewed distribution. The correlation between air pollution, meteorological data, and mortality data was assessed using the Spearman correlation coefficient.

To investigate the short-term relationship between PM_2.5 concentration and mortality from cardio-cerebrovascular diseases, a time series generalized linear model (GLM) based on Poisson regression was established using R3.0.1 to calculate the expose-response relationship. PM_2.5 monitoring data were used as the exposure concentration, while mortality data from CCVD in Harbin from 2016 to 2018 were used as the health outcome. Assuming a Poisson distribution, a log-link function was applied. The model controlled for long-term seasonal trends, meteorological factors, and day-of-the-week effects to establish the expose-response relationship between PM_2.5 pollution and CCVD mortality. The model is expressed as follows:

where E (Yt) represents the expected number of CCVD deaths on day t. Zt is the concentration level of air pollutants on day t. β is the exposure-response coefficient, indicating the increase in daily CCVD mortality associated with a unit increase in pollutant concentration. Ns represents the natural spline function, with df indicating the degree of freedom. Time is a date variable, and the appropriate df value helps control for long-term and seasonal fluctuations of mortality data. DOW is an indicator variable for ‘day of the week’. Xt represents the meteorological factor on day t, including average temperature and relative humidity.

To test the robustness of the model, a sensitivity analysis was conducted as follows: (1) The degrees of freedom (df) for the time smoothing function were sequentially varied from 6 to 8 to assess the impact of controlling for time trends on effect estimation. (2) The impact of different lag times (0-7 days) on the estimation was analyzed to determine the lag period during which PM_2.5 affects CCVD mortality. (3) The health effects were reanalyzed after incorporating coexisting pollutants, and a multi-pollutant model was developed to evaluate the impact of other atmospheric pollutants on the estimation of PM_2.5 health effects.

The PMF model was used to analyze the characteristics of pollutants and the contribution of various pollution sources. The dataset was represented as a data matrix X (i × j), where i is the number of samples measured and j is the number of chemical components measured. The matrix X (i × j) can be decomposed into two matrices: G (i × p) and F (p × j), where p represents the number of main pollution sources (factors), G is the pollution source contribution matrix, and F is the source component spectrum matrix. E represents the residual matrix. The PMF model can be represented as:

In this study, the EPA PMF5.0 model analysis tool was applied, using a concentration uncertainty of 10% as the initial estimate for each species. The criteria for evaluating the model included the proportion of data below the method detection limit (MDL), data loss, signal-to-noise (S/N) ratio, proportion of residuals, and r² value. After multiple iterations, the Q value stabilized, and the species were classified as “strong”, “weak”, or “bad” based on these criteria. Following this, a reasonable number of factors were selected for modeling and analysis.

3.1 Characteristics of meteorological factors, air pollutants and cardio-cerebrovascular system mortality

As shown in Table 1, the annual average values of meteorological factors in Harbin from 2016 to 2018 showed no significant changes, except for a slight increase in temperature. The annual average concentrations of PM_2.5 were 54.12 μg/m³, 57.16 μg/m³, and 37.37 μg/m³, respectively. The daily average number of deaths from cardio-cerebrovascular diseases was 92, 93, and 89, respectively. Between 2016 and 2018, the maximum daily average concentration of PM_2.5 reached 460.45 μg/m³, exceeding the Class II standard of the Environmental Air Quality Standard (GB 3095, 2012) (75 μg/m³) by more than six times. In 2016, the number of days with daily average concentrations exceeding the standard was 96 days, in 2017, it was 75 days, and in 2018, it was 43 days. T These exceedances were primarily concentrated during the heating season. Additionally, NO₂ and O₃ concentrations also exceeded the standard. The maximum daily average concentration of NO₂ was 134.36 μg/m³, exceeding the standard (80 μg/m³) by 1.6 times, with 25 days surpassing the standard. The maximum daily average concentration of O₃ was 229.73 μg/m³, exceedingthe standard (160 μg/m³) by more than 1.4 times, with 19 days exceeding the standard. However, SO₂ and CO concentrations did not exceed their respective standards.

3.2 Correlation between air pollutants and meteorological factors

Fig. 1 shows the Spearman rank correlation coefficient analysis between atmospheric pollutants and meteorological factors. PM_2.5 was found to be negatively correlated with temperature, relative humidity, wind speed, and O₃, and positively correlated with SO₂, NO₂, CO, and atmospheric pressure.

Fig. 1

Correlation analysis between air pollutants and meteorological factors

3.3 The single factor exposure-lag-response associations of PM_2.5

As shown in Table 2, there is a significant correlation between PM_2.5 exposure and the number of deaths from CCVD on the day of exposure, as well as with delays of 1 day, 2 days, and 5 days. The exposure-response (ER) value with a 2-day delay reached its maximum at 0.35% (95% CI: 0.15%, 0.55%).

After stratifying by gender, the results revealed a significant correlation between PM_2.5 exposure and the number of deaths from CCVD in males, with delays of 1 day, 2 days, 4 days, and 5 days. The ER value with a 2-day delay was at its maximum, reaching 0.39% (95% CI: 0.14%, 0.63%). In the female group, a significant correlation was found on the day of exposure, as well as with delays of 2 days and 5 days. The ER value was highest on the day of exposure, reaching 0.37% (95% CI: 0.09%, 0.65%).

After stratifying by age, the results indicated a significant correlation between PM_2.5 exposure and the number of deaths from CCVD in individuals under 65 years old on the day of exposure, with an ER value of 0.37% (95% CI: 0.02%, 0.71%). For individuals aged 65 and older, significant correlations were observed on the day of exposure, as well as with delays of 2 days, 4 days, 5 days, and 7 days. The maximum ER value was observed on the day of exposure, reaching 0.37% (95% CI: 0.15%, 0.60%).

3.4 The multi-factor exposure-lag-response associations of PM_2.5

PM_2.5 lag days data with statistically significant results in the single-pollutant model were selected for further analysis in the multi-pollutant model, as shown in Table 3.

On the day of PM_2.5 exposure, after correcting for one or more pollutants, the impact of PM_2.5 on the mortality rate from CCVD increased compared to the single-pollutant model.

For a one-day lag in PM_2.5 exposure, after adjusting for SO₂, the impact of PM_2.5 on CCVD mortality increased compared to the single-pollutant model, with an ER value of 0.27% (0.03%-0.50%). No significant correlation with CCVD mortality was found after adjusting for other pollutants or combinations of pollutants.

For a two-day lag in PM_2.5 exposure, the impact of PM_2.5 on CCVD mortality was reduced compared to the single-pollutant model when adjusted for SO2 or O3, with ER values of 0.29% (95% CI: 0.06%, 0.53%) and 0.30% (95% CI: 0.10%, 0.50%), respectively. No significant correlation with CCVD mortality was found after adjusting for other pollutants or combinations of pollutants.

After five days of delayed exposure to PM_2.5, the impact of PM_2.5 on CCVD mortality was reduced compared to the single-pollutant model when adjusted for NO₂, SO₂, or O₃, with ER values of 0.29% (0.02%-0.55%), 0.27% (0.03%-0.50%), and 0.29% (0.08%-0.49%), respectively. When adjusting simultaneously for SO₂ and NO₂, SO₂ and O₃, the impact of PM_2.5 on CCVD mortality was further reduced compared to the single-pollutant model, with ER values of 0.28% (0.01%-0.55%) and 0.25% (0.01%- 0.49%), respectively. No significant correlation with CCVD mortality was found after adjusting for other pollutants or combinations of pollutants.

3.5 Source Apportionment of PM_2.5

The classification of each species was determined based on the proportion of data below the method detection limit (MDL), as well as the criteria including missing data, signal-to-noise (S/N) values, proportion residuals, and r² values. Species classified as “strong” include SO₄ ²⁻, NO₃ ⁻, NH₄ ⁺ , As, Cd, Cr, Pb, Ni, TI, Ace, Ant, Flua, Pyr, Chr, BaA, BbF, BkF, BaP, DahA, BghiP, and IcdP. Species classified as “weak” include Cl-, Sb, Al, Be, Hg, Mn, Se, Nap, Flu, and Phe. Acy was classified as “bad”. By multiple adjustments to the number of factors through the model, the optimal number of 7 factors was determined. These 7 main sources of atmospheric PM_2.5 pollution in Harbin are illustrated in Fig. 2 and Fig. 3.

Fig. 2

Component spectra of each factor obtained by positive matrix factorization

Fig. 3

Contribution rates of each factor obtained by positive matrix factorization

In Factor 1, the characteristic elements Pyr, Flua, and BaA are identified. Pyr and BaA are both 4-ring PAHs, which generally originate from coal-fired emissions. Flua, a 3-ring PAHs, is typically associated with the industrial combustion of wood or oil^[19]. Moreover, Pyr and Flua are known biomarkers of coal-fired emis- sions^[20]. Thus, Factor 1 is identified as coal smoke dust, with a contribution rate of 6.7%.

In Factor 2, the characteristic elements Phe, Ace, and Flu areidentified. Phe and Flu are both 3-ring PAHs and generally originate from biomass combustion sources^[21]. In the winter and spring, the burning of farmland around Harbin, particularly straw burning, can generate a large amount of Ace. Therefore, Factor 2 is identified as biomass combustion dust, with a contribution rate of 6.9%.

In Factor 3, the characteristic elements NO₃ ^-, NH₄ ⁺ , and SO₄ ^2- are identified. These ions are all formed through the secondary conversion of their precursors. NO₃ ^- typically originates from the conversion of nitrogen oxides emitted by motor vehicle and other combustion processes. SO₄^2- is a secondary pollutant formed by SO₂ through gas-solid conversion mechanism. NH₄⁺ is mainly generated through the neutralization reaction of atmospheric NH₃^[22]. Therefore, Factor 3 is identified as secondary particles, with a contribution rate of 29.7%.

In Factor 4, the characteristic elements Cr and Ni are identified. Both Cr and Ni are heavy metals produced during metallurgical process^{[23, 24]}. Therefore, Factor 4 is identified as smelting dust, with a contribution rate of 8.4%.

In Factor 5, the characteristic elements TI, Pb, and As are identified. TI mainly originates from industrial emissions. Pb is a metal component used in various factory processes, characterized by its low melting point, which makes it prone to dust generation during smelting^[25]. As is mainly derived from decolorizing agents used in mining, metallurgy, pharmaceuticals, glass industries, as well as from various insecticides, coal burning, fertilizers, leather, and pesticides^[26]. Therefore, Factor 5 is identified as industrial dust, with a contribution rate of 6.3%.

In Factor 6, the characteristic elements Mn, Se, Al, and Hg are identified. Mn, Se, and Al are representative elements of soil sources^{[27, 28]}. Hg may originate from emissions related to activities such as cement production, mining, and the processing of mineral materials. Therefore, Factor 6 is identified as surface dust and construction dust, with a contribution rate of 31.0%.

In Factor 7, the characteristic elements BghiP, DahA, IcdP, BkF, BbF, and BaP are identified s. Research indicates that high ring PAHs (5 ring and 6 ring) are typically emitted by motor vehi- cles^{[29, 30]}. Therefore, Factor 7 is classified as motor vehicle exhaust dust, with a contribution rate of 10.9%.