COVID-19 incidence and mortality in Lombardy, Italy: An ecological study on the role of air pollution, meteorological factors, demographic and socioeconomic variables.

Lombardy, the most populated and industrialized Italian region, was the epicentre of the first wave (March and April 2020) of COVID-19 in Italy and it is among the most air polluted areas of Europe. We carried out an ecological study to assess the association between long-term exposure to particulate matter (PM) and nitrogen dioxide (NO2) on COVID-19 incidence and all-cause mortality after accounting for demographic, socioeconomic and meteorological variables. The study was based on publicly available data. Multivariable negative binomial mixed regression models were fitted, and results were reported in terms of incidence rate ratios (IRRs) and standardized mortality ratios (SMR). The effect of winter temperature and humidity was modelled through restricted cubic spline. Data from 1439 municipalities out of 1507 (95%) were included in the analyses, leading to a total of 61,377 COVID-19 cases and 40,401 deaths from all-causes collected from February 20th to April 16th and from March 1st to April 30th, 2020, respectively. Several demographic and socioeconomic variables resulted significantly associated with COVID-19 incidence and all-cause mortality in a multivariable fashion. An increase in average winter temperature was associated with a nonlinear decrease in COVID-19 incidence and all-cause mortality, while an opposite trend emerged for the absolute humidity. An increase of 10 μg/m3 in the mean annual concentrations of PM2.5 and PM10 over the previous years was associated with a 58% and 34% increase in COVID-19 incidence rate, respectively. Similarly, a 10 μg/m3 increase of annual mean PM2.5 concentration was associated with a 23% increase in all-cause mortality. An inverse association was found between NO2 levels and COVID-19 incidence and all-cause mortality. Our ecological study showed that exposure to PM was significantly associated with the COVID-19 incidence and excess mortality during the first wave of the outbreak in Lombardy, Italy.

Introduction

The Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is responsible for the Coronavirus disease 2019 (COVID-19), firstly identified in Wuhan (China) in December 2019 and declared pandemic by the World Health Organization (WHO) on March 11th, 2020 (World Health Organization, 2020). COVID-19 incidence and mortality rates are still rising in the world, with more than 60 million cases and about 1.5 million deaths as of November 30th, 2020 (Johns Hopkins University, 2020).

Italy was dramatically hit by the COVID-19 pandemic, with 1.6 million cases and nearly 55.000 deaths as of November 30th, 2020 (Johns Hopkins University, 2020). Lombardy, the most populated and industrialized Italian region, was the epicentre of the outbreak. It experienced one of the highest case-fatality rate (14%) in the world (Istituto Superiore di Sanità, 2020a; Johns Hopkins University, 2020), with more than 400.000 cases (about one fourth of those registered in Italy) and nearly 55.000 deaths out of 10 million inhabitants (Istituto Superiore di Sanità, 2020a). In March and April 2020, Lombardy registered a nearly three-fold excess in all-cause mortality, which was much higher than that observed on a national level (about 50%), as compared to the same period in the quinquennia 2015–2019 (Istituto Nazionale di Statistica, 2020a; Odone et al., 2020; Tosi et al., 2020).

Little is known about the risk factors related to the widespread outbreak and the unexpectedly high lethality of COVID-19 in Lombardy. A large retrospective study conducted on nearly 4000 critically ill COVID-19 patients hospitalized in a Lombardy network of Intensive Care Units (ICUs) showed a significant association between mortality and age, male sex, and pre-existing comorbidities, including chronic obstructive pulmonary disease, hypercholesterolemia and diabetes (Grasselli et al., 2020a, 2020b). The same figure emerged from a study by the Italian National Institute of Health (“Istituto Superiore di Sanità (ISS)”), which showed that male sex, advanced age and presence of chronic diseases were predictors of death (Istituto Superiore di Sanità, 2020b).

In addition to medical history and lifestyle factors, the scientific community has hypothesised a role of environmental factors in the diffusion of SARS-CoV-2, individual susceptibility to infection, and progression and severity of COVID-19 disease (Copat et al., 2020; Domingo et al., 2020). Lombardy and other northern Italian regions were severely hit by the COVID-19 pandemic, and also ranks among the most air polluted areas of Europe, being located in the Po valley (European Environment Agency, 2020). For these reasons, several studies were recently carried out to investigate the association between air pollution and COVID-19 outbreak spread and severity in Lombardy and other northern Italian regions (Bontempi, 2020a; Coker et al., 2020; Comunian et al., 2020; Conticini et al., 2020; Filippini et al., 2020, Ogen, 2020; Zoran et al., 2020a, 2020b).

In addition to air pollution, other population characteristics such as demographic, socioeconomic and community variables may contribute to COVID-19 spread and severity, but little is still known on their role (Bontempi et al., 2020; Bontempi, 2020b; Liang et al., 2020; Wu et al., 2020a). Also, the role of meteorological factors on COVID-19 incidence and mortality have been investigated in several studies (Guo et al., 2020; Islam et al., 2020; Wu et al., 2020b; Yuan et al., 2020). In fact, low temperature and humidity are considered to favour the spread of respiratory virus infections, but current evidence on their role in COVID-19 spread and severity is still unclear (Harmooshi et al., 2020; Mecenas et al., 2020).

Our study aims to investigate the impact of long-term exposure to PM2.5, PM10 and NO2 on COVID-19 incidence and excess in all-cause mortality in Lombardy through an ecological approach, accounting for demographic, socioeconomic and meteorological variables.

Methods

Study design and source of data

Study design

This is an ecological study based on publicly available data (Agenzia Regionale Protezione Ambiente Lombardia, 2020a; Borruso, 2020; Istituto Nazionale di Statistica, 2020b, 2011a, 2011b; Regione Lombardia, 2020, 2018). All the variables included in the analyses, as well data sources, are listed in Supplementary Table 1. Since all the data were publicly available and aggregated on a municipality level, ethical approval by the Ethical Committee of Brescia was not required.

Outcomes

The two outcome variables were: i) the daily number of COVID-19 cases diagnosed through SARS-CoV-2 swab tests from February 20th to April 16th, 2020. Most COVID-19 cases were likely severe, because SARS-CoV-2 swab tests were administered only to hospitalized patients in that period (Borruso, 2020); ii) the excess in all-cause mortality in March and April 2020, defined as the ratio between the observed and expected number of deaths. The former was centrally collected by the Italian National Institute of Statistics (“Istituto Nazionale di Statistica” (ISTAT)”) (Istituto Nazionale di Statistica, 2020b) from municipality-based register offices, while the latter was estimated from the average daily all-cause mortality registered in Lombardy in the quinquennia 2015–2019. Both observed and expected all-cause mortality data were retrieved from the ISTAT datawarehouse (Istituto Nazionale di Statistica, 2020b) (Supplementary Table 1). The excess in all-cause mortality was a more appropriate outcome measure as compared to official COVID-19 mortality data, because the former included both direct and indirect effects of SARS-CoV-2 infection, whereas the latter did not include people who died at home or in nursing homes without being diagnosed for COVID-19 (Michelozzi et al., 2020).

Demographic, socioeconomic and community variables

We considered the following demographic and socioeconomic variables: population size (the number of residents); population density (inhabitants per Km2); proportion of population over 75 years old; sex ratio (the number of males for 100 females); average family size; house crowding index (defined as the ratio between the number of occupied houses with less than 40 m2 and over 4 occupants, or with 40–59 m2 and over 5 occupants, or with 60–79 m2 and over 6 occupants, and the total number of occupied houses). All these variables were extracted from the 2011 national population Census data (Istituto Nazionale di Statistica, 2011a, b).

In addition, we considered two socioeconomic indicators: i) the high to low education ratio, defined as the ratio between the number of residents aged 25–64 years with high school and/or university degree and the total number of residents, and ii) the “Imposta Regionale Persone Fisiche” (IRPEF), a regional tax collected by the “Ministero dell’Economia e delle Finanze” (MEF), as a proxy of the household income (Istituto Nazionale di Statistica, 2020b).

As a community variable, we considered the percentage of private mobility use, defined as a ratio between the residing population that commutes daily for work or study using a private motor vehicle (car or motorbike) and the total residing population that commutes daily for work or study in each municipality (Istituto Nazionale di Statistica, 2011a, b).

The following proxy variables describing the health services available in each municipality were considered: i) the number of beds in nursing homes and ii) the distance of each municipality from the closest hospital. Data on the number of beds in nursing homes was retrieved from the regional welfare department datawarehouse (Regione Lombardia, 2020), while hospital's geographical coordinates (longitude and latitude) were retrieved from data made freely available by the Lombardy regional epidemiological observatory (Regione Lombardia, 2018). The distance of a municipality from the closest hospital, identified through the polygon centroid, was computed through a distance matrix (QGIS Project, 2020) and identified as the minimum of the Euclidean distances between each municipality and each hospital in the region. When a municipality had one or more healthcare facilities within the municipal boundaries, the default distance was assumed to be 0.

ISTAT published yearly the business register “Archivio Statistico delle Imprese Attive” (ASIA), which includes employment data of all enterprises carrying on economic activities in the fields of industry, trade and services (Istituto Nazionale di Statistica, 2011a, b). As proxy variables of social activities, we considered for each municipality the number of employees in bars, restaurants and mobile catering activities (ATECO code 56) and the number of employees in sports, entertainment and recreational activities (ATECO code 93). The number of employees in health and social assistance activities (ATECO code 86-87-88) for each municipality was also considered since health professionals have been on the frontline fighting the epidemic within all settings of the healthcare system.

Meteorological variables and ambient air pollutants concentrations

We considered two meteorological variables: the average winter temperature, as measured in Celsius degrees (°C), and the average winter absolute humidity, measuring grams of water vapor (moisture) per cubic meter of air (g/m3), regardless of temperature. These data are regularly monitored, processed and published by the Lombardy “Agenzia Regionale Protezione Ambiente” (ARPA). The Lombardy ARPA meteorological service computes hourly meteorological fields on a 1.5 × 1.5 km2 resolution by applying optimal interpolation techniques to measured data (Uboldi et al., 2008). We computed average temperature and average absolute humidity values from February 15th to April 15th, 2020, and aggregated them on a municipality level using the GIS raster zonal statistics algorithm (QGIS Project, 2020). This time window included a 15 days lagged effect to take into account the COVID-19 incubation period.

PM10, PM2.5 and NO2 spatial annual average concentrations were estimated over the Lombardy region on a 4 × 4 km2 spatial resolution grid by means of the Chemical Transport Model ARIA Regional, and integrated, through assimilation techniques (Bocquet et al., 2015), with the air quality data monitored by the regional network (Agenzia Regionale Protezione Ambiente Lombardia, 2020b). Model runs were validated comparing measurements and model estimates. ARPA Lombardia tests for the accuracy of the air quality network measurements in two steps: daily and yearly. PM concentrations model estimates were released for the previous years (2016–2019) assimilating yearly tested measurements.

Modelled PM and NO2 concentrations were processed to estimate area-weighted mean concentration for each municipality (Agenzia Regionale Protezione Ambiente Lombardia, 2020a). Since PM10 and NO2 yearly tested modelled data was available from 2011 to 2019 and PM2.5 data was available from 2016 to 2019, we computed the long-term PM and NO2 exposure as the mean of annual PM concentrations over the four years period 2016–2019.

Statistical analysis

The statistical analysis aimed to assess the associations between demographic, socioeconomic and community, meteorological variables and environmental air pollutants concentrations and the two outcomes: i) the incidence of COVID-19 (model 1) from February 20th to April 16th, 2020, and ii) the excess in all-cause mortality (model 2) from March 1st to April 30th, 2020.

We conducted a descriptive analysis by summarizing the continuous variables as mean and standard deviation (SD), or median and interquartile range (IQR). For each municipality, we plotted the number of COVID-19 cases and the standardized mortality ratio (SMR) in the study period, as well as the concentrations of PM2.5 and PM10, on a georeferenced map of Lombardy, to provide a visual representation of the spatial distribution of the COVID-19 cases and all-cause deaths and the air pollutant levels.

We fitted model 1 and model 2 through a multivariable negative binomial (NB) random-intercept mixed regression model. We chose the NB model based on the nature of the two outcomes variables, which include count data, and because of the presence of overdispersion, which make the use of a Poisson regression model inappropriate. We used the population size and the expected all-cause deaths registered in March–April in the quinquennia 2015–2019 as offsets for model 1 and model 2, respectively. A mixed model with a random intercept was fitted to take into account the longitudinal component of the model, as the number of COVID-19 cases (model 1) and the number of all-cause deaths (model 2) were grouped on a 15-days’ time window. The model statistical unit was the municipality-specific biweekly COVID-19 incidence and excess in all-cause mortality.

For both model 1 and model 2, we fitted separate multivariable models investigating the role of long-term exposure to PM2.5, PM10 and NO2 concentrations. These models included the time since the first COVID-19 case report (a proxy for the epidemic stage) and the following independent variables (reported in Supplementary Table 1): population size (log2-transformed and standardized); population density (log2-transformed and standardized); sex ratio (standardized); percent of population over 75 years old (standardized); average family size (standardized); house crowding index (categorized in tertiles); high to low education ratio (log2-transformed and standardized); median municipality personal IRPEF tax (standardized); percentage of private mobility use (standardized); number of licensed beds in nursing homes (categorized in no facilities, 0–80 and > 80 number of beds); distance to the closest hospital; number of employees in bars, restaurants and mobile catering activities per capita (categorized in tertiles); number of employees in health and social assistance activities per capita (categorized in tertiles); number of employees in sports, entertainment and recreational activities per capita (categorized in tertiles); winter average temperature (°C) and absolute humidity (g/m3) in February and April 2020. The effect of the latter meteorological variables was lagged by 15 days to account for the COVID-19 incubation period and modelled through a 3-knots restricted cubic spline to allow for nonlinear relationships with the two investigated outcomes.

The results of model 1 are presented in terms of Incidence Rate Ratios (IRRs) and those of model 2 in terms of Standardize Mortality Ratios (SMRs), together with their 95% CIs. The IRR or the SMR can be interpreted as the relative decrease (if lower than 1) or increase (if greater than 1) of the COVID-19 incidence or excess in all-cause mortality associated with a one unit increase for each independent variable. The model fitting was assessed through the R2 for generalized linear mixed model (Nakagawa and Schielzeth, 2013).

A type I alpha level of 0.05 was considered as the threshold for statistical significance. We carried out all analyses with the R version 4.0.0 using the “lme4” package (Bates et al., 2015).

Results

Complete data were available for 1439 out of 1507 (95%) Lombardy municipalities. A total of 61,377 COVID-19 cases and 40,401 all-causes deaths were recorded in Lombardy from February 20th to April 16th and from March 1st to April 30th, 2020, respectively. The incidence rate of COVID-19, the observed and the observed to expected all-cause mortality rates over a 15-days’ time window are shown in Fig. 1 . Incidence and mortality rates trends follow the same pattern and show a peak in the two last weeks of March (16th-31st). The effect of regional and national restrictive measures for contrasting human-to-human transmission, particularly the lockdown, established on March 9th, determined a sharp decline of the epidemic curve starting from April 1st.

Fig. 1

COVID-19 incidence and observed and expected all-cause mortality rates in Lombardy in March–April 2020 (15-days’ time windows).

Table 1 reports municipality-level descriptive statistics of study outcomes and investigated variables. Lombardy has more than 10 million inhabitants, a mean municipality population size of 6495.7 inhabitants and a mean population density of 564.7 inhabitants per Km2. The mean percentage of population over 75 years old is 9.8%. The number of COVID-19 cases and the number of observed all-cause deaths varied substantially among municipalities, with median values of 16 COVID-19 cases and 11 deaths per municipality. In March and April 2020, the mean SMR was equal to 2.9, meaning that there was a nearly three-fold excess in all-cause mortality in Lombardy municipalities as compared to the same period of the quinquennia 2015–2019.

Table 1

Municipality-level descriptive statistics of study outcomes, demographic, socioeconomic and community, meteorological variables and air pollutants concentrations in 1459 Lombardy municipalities.

	Mean (SD)	Median (Q1, Q3)
Study outcomes
Number of observed deaths (Mar–Apr, 2020)	27.3 (123.6)	11.0 (4.0, 25.0)
Number of expected deaths (Mar–Apr, 2020)	11.0 (63.7)	4.4 (2.0, 9.8)
Standardized mortality ratio (SMR)	2.9 (2.2)	2.4 (1.5, 3.8)
Number of COVID-19 cases	41.5 (184.4)	16.0 (5.5, 40.0)
Demographic, socioeconomic and community variables
Population size (inhabitants)	6495.7 (33826.1)	2842.0 (1227.5, 6033.5)
Population density (inhabitants per Km2)	564.7 (787.9)	269.8 (101.0, 745.6)
Sex ratio	98.2 (5.4)	97.9 (95.0, 101.0)
Proportion of population over 75 years old (%)	9.8 (3.3)	9.3 (7.7, 11.3)
Average family size	2.4 (0.2)	2.4 (2.3, 2.5)
House crowding index	0.4 (0.3)	0.3 (0.2, 0.5)
High to low education ratio	133.4 (48.8)	126.4 (100.8, 157.3)
IRPEF per capita (€)	14449.6 (2770.6)	14311.0 (12824.2, 16017.2)
Percentage of private mobility use (%)	68.0 (7.3)	68.5 (63.6, 73.1)
Number of beds in nursing homes	43.0 (222.2)	0.0 (0.0, 60.0)
Distance to the closest hospital (meters)	5976.6 (3856.1)	5443.1 (3426.8, 8090.2)
Number of employees in bars, restaurants and mobile catering activities per capita (per 1000 inhabitants)	16.9 (30.3)	11.2 (7.7, 16.9)
Number of employees in health and social assistance activities per capita (per 1000 inhabitants)	0.5 (0.6)	0.3 (0.0, 0.7)
Number of employees in sports, entertainment and recreational activities per capita (per 1000 inhabitants)	0.8 (0.0)	0.6 (0.0, 1.1)
Meteorological variables
Winter (Feb–Apr, 2020) average temperature (°C)	8.5 (2.3)	9.4 (8.2, 9.9)
Winter (Feb–Apr, 2020) average absolute humidity (g/m3)	5.5 (0.9)	5.7 (5.1, 6.3)
Air pollutants concentrations
PM2.5 (μg/m3)	18.6 (4.8)	19.8 (15.6, 22.1)
PM10 (μg/m3)	21.8 (5.4)	22.7 (18.2, 25.8)
NO2 (μg/m3)	33.0 (13.0)	33.0 (25.5, 40.0)

Abbreviations: IRPEF, Imposta Regionale Persone Fisiche; NO2, nitrogen dioxide; PM, particulate matter; Q1, first quartile; Q3, third quartile; SD, standard deviation; SMR, standardized mortality ratio.

The spatial plot in Fig. 2 shows that the highest incidence rates (panel A) and SMRs (panel B) were observed in the central and eastern part of the Lombardy region. The first Italian COVID-19 case was officially diagnosed on February 18th, 2020 in the town of Codogno, Lodi province, located in the central part of the region. In the next weeks the epidemic evolved quickly, and the provinces of Cremona, Bergamo and Brescia, in the central-eastern part of the region, were the most severely hit.

Fig. 2

Municipality-level COVID-19 incidence rates (panel A) and standardized mortality ratios (SMRs, panel B) in Lombardy in March–April 2020.

The pattern of air pollution in Lombardy showed in Fig. 3 indicates that the highest concentrations of PM2.5 (panel A) and PM10 (panel B) were distributed homogeneously in the central part of the region, along a horizontal axis, corresponding to the Po Valley, the core of industrial activities in Lombardy.

Fig. 3

Municipality-level particulate matter (PM) concentrations in Lombardy. Panel A shows estimates of long-term (2016–2019) PM2.5 concentrations, panel B estimates of long-term (2016–2019) PM10 concentrations.

COVID-19 incidence as outcome

PM concentrations were associated with COVID-19 incidence at univariate analysis, while no association emerged for NO2 (Supplementary Table 2). The results of the multivariable model 1 are shown in Table 2 . An inverse relationship resulted for the high to low education ratio (IRR 0.76), meaning that municipalities with a higher rate of secondary school and university graduates experienced a less severe spread of the disease as compared to those with lower rates of highly educated residents. The percentage of private mobility use was inversely related (IRR 0.85) to COVID-19 incidence, as well as the distance from the municipality to the closest hospital (IRR 0.92). A 49% increase in the IRR in municipalities with 80 or more licensed beds in nursing homes emerged as compared to those with none. The number of employees in bars, restaurants and mobile catering activities per capita (IRR 1.11) and the number of employees in sports, entertainment and recreational activities per capita (IRR 1.18) showed a positive association with COVID-19 incidence. Temperature and absolute humidity were associated in a nonlinear J-shaped fashion with COVID-19 incidence (Fig. 4 ). An increase of 10 μg/m3 in the mean annual PM2.5 and PM10 concentrations in 2016–2019 was associated with a statistically significant 58% and 34% increase in COVID-19 incidence, respectively. Mean annual NO2 levels were inversely associated with COVID-19 incidence, although the result was significant only within the model including PM2.5. The pseudo R2 of the model was 0.71, showing a good fit of the model to the data.

Table 2

Multivariable estimates of incidence rate ratios (IRRs) for COVID-19 in Lombardy municipalities between February 20th and April 16th, 2020, according to demographic, socioeconomic and community variables and air pollutants concentrations.

	Model with PM2.5			Model with PM10
	IRRa	95% CI	p-value	IRRa	95% CI	p-value
Demographic, socioeconomic and community variables
Population sizeb	0.90	0.83–0.98	0.01	0.90	0.83–0.98	0.01
Population densityb	0.92	0.83–1.01	0.07	0.92	0.84–1.01	0.10
Sex ratioc	0.99	0.94–1.05	0.83	1.00	0.94–1.05	0.87
Proportion of population over 75 years oldc	1.08	1.01–1.15	0.02	1.09	1.02–1.16	0.01
Average family sizec	1.01	0.95–1.08	0.79	1.03	0.97–1.10	0.33
House crowding index (high vs low)	0.93	0.84–1.03	0.19	0.92	0.83–1.02	0.12
High to low education ratiob	0.76	0.72–0.81	<0.01	0.76	0.72–0.81	<0.01
IRPEF per capitac	1.19	1.12–1.27	<0.01	1.18	1.11–1.26	<0.01
Percentage of private mobility usec	0.85	0.81–0.88	<0.01	0.85	0.81–0.89	<0.01
Number of beds in nursing homes (>80 vs 0)	1.49	1.32–1.69	<0.01	1.49	1.32–1.68	<0.01
Distance to the closest hospitalc	0.92	0.87–0.97	<0.01	0.92	0.88–0.97	<0.01
Number of employees in bars, restaurants and mobile catering activities per capita (>14 vs ≤ 9 per 1000 inhabitants)	1.11	1.00–1.23	0.04	1.12	1.02–1.24	0.02
Number of employees in health and social assistance activities per capita (>0.6 vs 0 per 1000 inhabitants)	1.02	0.92–1.14	0.70	1.02	0.92–1.14	0.66
Number of employees in sports, entertainment and recreational activities per capita (>0.9 vs ≤ 0.4 per 1000 inhabitants)	1.18	1.07–1.30	<0.01	1.18	1.07–1.31	<0.01
Air pollutants concentrations
PM (per 10 μg/m3)	1.58	1.31–1.90	<0.01	1.34	1.16–1.55	<0.01
NO2 (per 10 μg/m3)	0.93	0.88–0.99	0.02	0.95	0.89–1.00	0.07

Abbreviations: CI, confidence interval; IRPEF, Imposta Regionale Persone Fisiche; IRR, incidence rate ratio; NO2, nitrogen dioxide; PM, particulate matter.

Adjusted by weeks since first COVID-19 case report (a proxy for epidemic stage), winter average temperature and winter average absolute humidity as modelled through a 3-knots restrcited cubic spline (Fig. 4).

Variable was log2- transformed and standardized.

Variable was standardized.

Fig. 4

Effects of winter average temperature and winter absolute humidity on COVID-19 incidence as modelled through a 3-knots restricted cubic spline in a multivariable model including long-term exposure to PM2.5 (Panel A and B, respectively) and PM10 (Panel C and D, respectively).

Excess in all-cause mortality as outcome

The univariable models showed significant associations between PM concentrations and excess in all-cause mortality (Supplementary Table 2).

The results of model 2 are shown in Table 3 . An inverse relationship resulted for the high to low education ratio (SMR 0.75), while a 27% increase in the SMR resulted in municipalities with 80 or more licensed beds in nursing homes compared to those with none. The number of employees in bars, restaurants and mobile catering activities per capita (SMR 1.12) and the number of employees in sports, entertainment and recreational activities per capita (SMR 1.15), were positively associated with all-cause deaths. An increase in average winter temperature was associated with a nonlinear decrease in all-cause mortality (Fig. 5 Panel A and C), while an opposite trend emerged for the absolute humidity (Fig. 5 Panel B and D). A 10 μg/m3 increase in mean annual PM2.5 levels in 2016–2019 was associated with a 23% excess in all-cause mortality, while no association emerged for exposure to PM10. NO2 was inversely related to all-cause mortality, although the result was significant only within the model including PM2.5. Model fit as measured by the pseudo R2 was 0.43, showing a lower fit to the data as compared to the model fitted with COVID-19 as outcome.

Table 3

Multivariate estimates of all-cause standardized mortality ratios (SMRs) in Lombardy municipalities between March 1st and April 30th, 2020, according to demographic, socioeconomic and community variables and air pollutants concentrations.

	Model with PM2.5			Model with PM10
	SMRa	95% CI	p-value	SMRa	95% CI	p-value
Demographic, socioeconomic and community variables
Population sizeb	0.92	0.86–0.98	0.01	0.92	0.86–0.98	0.01
Population densityb	1.12	1.04–1.22	<0.01	1.13	1.04–1.22	<0.01
Sex ratioc	1.00	0.96–1.04	0.97	1.00	0.96–1.05	0.99
Proportion of population over 75 years oldc	1.01	0.96–1.06	0.76	1.01	0.96–1.07	0.75
Average family sizec	1.00	0.95–1.06	0.98	1.01	0.96–1.07	0.63
House crowding index (high vs low)	0.95	0.87–1.03	0.22	0.94	0.87–1.03	0.18
High to low education ratiob	0.75	0.72–0.79	<0.01	0.75	0.72–0.89	<0.01
IRPEF per capitac	1.16	1.10–1.22	<0.01	1.15	1.09–1.21	<0.01
Percentage of private mobility usec	0.93	0.90–0.96	<0.01	0.93	0.90–0.96	<0.01
Number of beds in nursing homes (>80 vs 0)	1.27	1.16–1.40	<0.01	1.27	1.16–1.40	<0.01
Distance to the closest hospitalc	0.98	0.95–1.02	0.44	0.99	0.95–1.03	0.49
Number of employees in bars, restaurants and mobile catering activities per capita (>14 vs ≤ 9 per 1000 inhabitants)	1.12	1.03–1.22	<0.01	1.13	1.04–1.23	<0.01
Number of employees in health and social assistance activities per capita (>0.6 vs 0 per 1000 inhabitants)	0.96	0.88–1.06	0.43	0.96	0.88–1.05	0.42
Number of employees in sports, entertainment and recreational activities per capita (>0.9 vs ≤ 0.4 per 1000 inhabitants)	1.15	1.06–1.25	<0.01	1.15	1.06–1.25	<0.01
Air pollutants concentrations
PM (per 10 μg/m3)	1.23	1.05–1.44	<0.01	1.06	0.94–1.20	0.31
NO2 (per 10 μg/m3)	0.94	0.90–0.99	0.01	0.96	0.92–1.01	0.08

Abbreviations: ARPA, Agenzia Regionale Protezione Ambiente; CI, confidence interval; IRPEF, Imposta Regionale Persone Fisiche; NO2, nitrogen dioxide; PM, particulate matter; SMR, standardized mortality ratio.

Adjusted by weeks since first COVID-19 case report (a proxy for epidemic stage), winter average temperature and winter average absolute humidity as modelled through a 3-knots restrcited cubic spline (Fig. 5).

Variable was log2- transformed and standardized.

Variable was standardized.

Fig. 5

Effects of winter average temperature and winter absolute humidity on excess in all-cause mortality as modelled through a 3-knots restricted cubic spline in a multivariable model including long-term exposure to PM2.5 (Panel A and B, respectively) and PM10 (Panel C and D, respectively).

Discussion

This study found that incidence of COVID-19 and all-cause mortality were influenced by several factors, with a high goodness of fit of the models. The frequency of COVID-19 cases decreased for increasing educational level, private mobility use, distance to the closest hospital and temperature, and increased for increasing humidity, beds in nursing homes, number of employees in bars, restaurants and mobile catering activities and employees in sports, entertainment and recreational activities, proxy variables of social activity levels. Most of these variables were also associated with excess in all-cause mortality, although the model fit was lower than that achieved by the model on COVID-19 incidence.

These findings confirm the role of demographic, socio-economic, community, and meteorological variables found in other ecologic studies on COVID-19 mortality and fatality rate (Guo et al., 2020, Wu et al., 2020a, Yuan et al., 2020).

The effects of meteorological factors, particularly temperature and humidity, on COVID-19 spread have been widely investigated since the beginning of the pandemic. In fact, agents that cause respiratory tract infections usually have their maximum infection potential in cold seasons, and experimental studies have shown that increasing temperature and humidity decreases SARS-CoV-2 survival on surfaces (Biryukov et al., 2020; Magurano et al., 2020). However, although the results of published studies on this issue are rather homogeneous, the overall evidence for a role of these factors on COVID-19 spread is still weak, due to the low quality of most studies (Mecenas et al., 2020).

In our study we found that a decrease in temperature and an increase in absolute humidity was associated with increasing COVID-19 incidence and all-cause mortality in a non-linear fashion, when also controlling for demographic, socio-economic and community variables. These findings are substantially in agreement with those of vast country-levels ecological studies, some of which also found a non-linear relationship between meteorological factors and COVID-19 incidence (Guo et al., 2020; Islam et al., 2020; Wu et al., 2020b; Yuan et al., 2020). However, although temperature and humidity seem to influence incidence and severity of COVID-19, they cannot explain much of the variability in SARS-CoV-2 spread in the world (Mecenas et al., 2020). In fact, the virus can also be transmitted with warm and moist weather, as demonstrated by the high number of cases registered in Europe and North America in July and August 2020 (Johns Hopkins University, 2020).

Air pollution has been considered a potential risk factor for COVID-19 incidence and mortality since the initial outbreak of COVID-19 in some areas of China, Iran and Northern Italy characterized by low ambient air quality. The detection of SARS-CoV-2 RNA on ambient air PM in one of the Italian areas with the highest COVID-19 incidence has also suggested a role of PM in sustaining virus survival in the environment (Setti et al., 2020). Furthermore, the high COVID-19 case-fatality rate observed in some highly industrialized areas of Northern Italy, mostly due to circulatory and respiratory system diseases, suggested that subjects with SARS-CoV-2 infection could experience a more severe disease if also exposed to high levels of air pollutants (Domingo et al., 2020).

We found that a 10 μg/m3 increase in PM2.5 and PM10 concentrations was associated with a 58% and 34% increase in COVID-19 incidence. Further, a 23% increase in all-cause mortality was observed for a 10 μg/m3 increase in PM2.5 levels. This last finding is in agreement with the 11% increase in the county's COVID-19 mortality rate for 1 μg/m3 increase in PM2.5 observed in the US ecological regression analysis by Wu et al. (2020a), who also controlled for 20 county-level covariates (Wu et al., 2020a). The hypothesis that PM may increase the risk of COVID-19 occurrence, severity and lethality is biologically plausible. Indeed, the effects of PM exposure on circulatory and respiratory systems have long been established, and PM exposure has been associated with the occurrence and severity of other respiratory tract infections (Domingo and Rovira, 2020; Mehta et al., 2013). PM exposure has been found to cause chronic inflammation and cellular damage through oxidative stress, and to dysregulate the human immune response making people more susceptible to infections (Sood et al., 2018, Wang et al., 2020). Furthermore, air pollutants can also increase the angiotensin-converting enzyme 2 (ACE-2), which is a receptor for SARS-CoV-2, thus increasing the risk of developing severe COVID-19 (Wang et al., 2020).

Our study showed an inverse association between long-term exposure to NO2 and COVID-19 incidence and all-cause mortality, in line with the findings of a study conducted in the metropolitan area of Milan (Zoran et al., 2020b).

Nevertheless, in spite of a plea of scientific publications and reprints on the air-pollution and COVID-19 association, this hypothesis still needs confirmation. Although the majority of studies found associations between either recent or historical levels of air pollutants and COVID-19 morbidity or mortality (Copat et al., 2020; Domingo et al., 2020), some discrepancies deserve attention. A European study showing a relationship between NO2 levels and COVID-19 fatality rate (Ogen, 2020) has also raised criticism due to methodologic flaws (Chudnovsky, 2020; Pisoni and Van Dingenen, 2020). An US county-based study found an association between NO2, but not PM2.5, levels and COVID-19 mortality and fatality (Liang et al., 2020).

As pointed out by Heederick et al. (Heederik et al., 2020) and Villeneuve and Goldberg (2020), recent published studies investigating the impact of air pollution on COVID-19 had some methodologic flaws, including, among the others, a lack of consideration of the role of confounders and the use of inappropriate statistical methods, as well as an excess of emphasis on statistical significance rather than on the magnitude of the effect (Boutron and Ravaud, 2018; Pearce et al., 2020b; Villeneuve and Goldberg, 2020).

In our study, we aimed to overcome all these potential flaws. First, by limiting place and time of the analysis. Italian regions have a wide autonomy in health management and adopted different policies in managing the COVID-19 emergency. In order to avoid confounding or effect modification by health policy, we restricted our analysis to Lombardy, the most severely affected European region. We analysed only the data of the first months of the epidemic, March and April 2020, because subsequent measures established for contrasting the epidemic, i.e. the total lockdown imposed in Italy, determined a substantial reduction of COVID-19 incidence and lethality, making the assessment of risk factors on a community basis more difficult. Second, we analysed both incidence and excess in mortality by fitting adequate models and taking into account the population size and the expected deaths as model offsets. We used official data for COVID-19 cases to avoid the risk of misclassification of cases, which is a common limit of ecological studies (Villeneuve and Goldberg, 2020). Furthermore, the use of all-cause mortality allowed to overcome the possible misclassification of causes of death in subjects with COVID-19. Finally, we controlled for a comprehensive set of potential confounders, including the epidemic stage, several demographic, socioeconomic and community variables, proxy variables of the intensity of inter-human contacts, and meteorological factors.

Nevertheless, our study has some limitations, too. First, the ecological study design, due to lack of individual data, is prone to ecologic bias (Villeneuve and Goldberg, 2020). However, this study design allows to estimate the effect of exposures with a homogeneous diffusion on a community basis, like air pollutants (Pearce et al., 2020a). Moreover, although we considered several proxy variables of population characteristics and activities, no data were available on other factors influencing COVID-19 transmission and severity, from individual habits (tobacco smoking, physical activity, obesity, etc.) to comorbidities. We could not also consider the occupational exposure to airborne particulates, which is historically relevant for most provinces of the Lombardy region, because no aggregate data are publicly available on this potential additional source of exposure. For these reasons, high-level evidence from longitudinal prospective studies is necessary to provide a comprehensive evaluation of risk factors for COVID-19 incidence and mortality (Heederik et al., 2020; Villeneuve and Goldberg, 2020).

In conclusion, our study suggests that PM exposure, meteorological factors and several socioeconomic and community variables may contribute to the incidence and mortality of COVID-19. Public health interventions for better air pollution control might contribute to reduce the global burden of related diseases, and possibly of severe COVID-19 cases, especially in highly polluted areas such as Lombardy, according to the WHO air quality guidelines (World Health Organization, 2017). Further studies are needed at the individual level to assess more representative exposure metrics of environmental and occupational exposure to airborne particulates, and to assess the role of comorbidities.

Author contribution

Elena De Angelis: Conceptualization; Data curation; Writing – review & editing. Stefano Renzetti: Conceptualization; Methodology; Software; Formal analysis; Data curation; Writing – review & editing. Marialuisa Volta: Conceptualization; Writing – original draft; Writing – review & editing; Supervision; Visualization; Project administration. Francesco Donato: Conceptualization; Writing – original draft; Writing – review & editing; Supervision; Visualization; Project administration. Stefano Calza: Conceptualization; Methodology; Supervision; Visualization; Project administration. Donatella Placidi: Conceptualization; Supervision; Project administration. Roberto Lucchini: Conceptualization; Writing – review & editing; Supervision; Visualization; Project administration. Matteo Rota: Conceptualization; Writing – original draft; Writing – review & editing; Visualization; Project administration

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.