Vulnerability and Burden of All-Cause Mortality Associated with Particulate Air Pollution during COVID-19 Pandemic: A Nationwide Observed Study in Italy

Background: Limited evidence is available on the health effects of particulate matter (PM including PM2.5 with an aerodynamic diameter ≤ 2.5 μm; PM10, ≤ 10 μm; PM2.5–10, 2.5–10 μm) during the pandemic of COVID-19 in Italy. The aims of the study were to examine the associations between all-cause mortality and PM in the pandemic period and compare them to the normal periods (2015–2019). Methods: We collected daily data regarding all-cause mortality (stratified by age and gender), and PM concentrations for 107 Italian provinces from 1 January 2015 to 31 May 2020. A time-stratified case-cross design with the distributed lag non-linear model was used to examine the association between PM and all-cause mortality. We also compared the counts and fractions of death attributable to PM in two periods. Results: Italy saw an increase in daily death counts while slight decreases in PM concentrations in pandemic period. Each 10 µg/m3 increase in PM was associated with much higher increase in daily all-cause mortality during the pandemic period compared to the same months during 2015–2019 (increased mortality rate: 7.24% (95%CI: 4.84%, 9.70%) versus 1.69% (95%CI: 1.12%, 2.25%) for PM2.5; 3.45% (95%CI: 2.58%, 4.34%) versus 1.11% (95%CI: 0.79%, 1.42%) for PM10; 4.25% (95%CI: 2.99%, 5.52%) versus 1.76% (95%CI: 1.14%, 2.38%) for PM2.5–10). The counts and fractions of deaths attributable to PM were higher in 2020 for PM2.5 (attributable death counts: 20,062 versus 3927 per year in 2015–2019; attributable fractions: 10.2% versus 2.4%), PM10 (15,112 versus 3999; 7.7% versus 2.5%), and PM2.5–10 (7193 versus 2303; 3.7% versus 1.4%). Conclusion: COVID-19 pandemic increased the vulnerability and excess cases of all-cause mortality associated with short-term exposure to PM2.5, PM2.5–10, and PM10 in Italy, despite a decline in air pollution level.


Data Collection
We collected the ERA5 hourly surface (at 2 meters above the land surface) ambient temperature and ambient dew point temperature at 0.1° × 0.1° spatial resolution from January 1 st , 2015 to May 31 st , 2020 for Italy and its surrounding areas. This data was downloaded from the ERA5-Land hourly data (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=form).
To valid and calibrate the ERA5 data, we downloaded hourly ambient temperature and ambient dew point temperature observed by weather stations from the Integrated Surface Database (https://www.ncdc.noaa.gov/isd) hosted by National Oceanic and Atmospheric Administration of US. A total of 519 weather stations with hourly weather observations during 2015-2020 for Italy and its surrounding areas were used to validate the EAR5 hourly dataset ( Figure S1). Hourly observations were transformed into daily observations by averaging all hourly observations within each day. Validation and calibration were based on a total of 692,002 daily observations. We did not use the daily temperature observations from weather stations for our analyses mainly for two reasons: 1) the spatial coverage of those stations is not satisfactory, with many provinces not covered; 2) the daily temperature observations from many weather stations were often incomplete. In contrast, the ERA5 data has perfect spatial and temporal coverage, despite its less accuracy.
To improve the performance of the calibration model, we also collected data on elevation and leaf area index (LAI) that may affect temperatures. LAI is an indicator for the density of vegetation cover. The monthly LAIs for both high and low vegetation at 0.1° × 0.1° spatial resolution were downloaded from ERA5-Land monthly averaged data (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5land-monthly-means?tab=form). The global elevation data at 7.5 arc second (i.e., 7.5/3600 degree × 7.5/3600 degree, about 250 m × 250 m) spatial resolution was downloaded from the NASA Shuttle Radar Topographic Mission (SRTM http://srtm.csi.cgiar.org/srtmdata/). We aggregated the elevation data to 0.1° × 0.1° spatial resolution by averaging elevation values of the 48 × 48 blocks of 7.5/3600 degree × 7.5/3600 degree pixels. ERA5 daily temperature, ERA5 LAI, and elevation data were linked to all weather stations according to date, longitude and latitude. Figure S1. Distribution of weather stations that used to calibrate the ERA5 daily temperature and dew point temperature.

Calibration Model
We used random forest model to calibrate the ERA5 daily mean temperature and dew point temperature against the observations from weather stations. In the random forest models, the dependent variable was either the daily mean temperature or daily mean dew point temperature, and the independent variables (predictors) included ERA5 daily mean temperature, ERA5 daily mean dew point temperature, monthly LAI, elevation, year, day of year (1 to 366), latitude, longitude. The random forest models were fitted by the "ranger" (version 0.10.1) R package with default settings. The models built based on all the 692,002 daily observations from 519 weather stations were then used to calibrate the ERA5 daily temperature observations. The calibrated ERA5 daily mean temperature and ERA5 daily mean dew point temperature were then used to calculate daily mean relative humidity, using the algorithm provided by the "humidity" R package.
To evaluate the accuracy of the calibration model, we performed a 10-fold cross-validation by randomly splitting the 519 stations into 10 folds. As showed in Figure S2, the calibrated ERA5 temperatures (i.e., predicted values by the calibration model) were highly correlated with the observed values from weather stations, suggesting a high accuracy of the calibration model. Notes: other covariates included elevation, year, day of year (1 to 366), latitude, longitude, and leaf area index. We performed a 10-fold cross-validation by randomly splitting the 519 stations into 10 folds. RMSE, root mean square error. MAE, mean absolute error.    Note: p-value was estimated by fixed effect meta-regression. * p-value for differences in subgroups; # p-value for differences in different time periods where the normal period of 2015-2019 was taken as the reference group.