key: cord-0778701-3em9vq1c
authors: de Angelis, Elena; Renzetti, Stefano; Volta, Marialuisa; Donato, Francesco; Calza, Stefano; Placidi, Donatella; Lucchini, Roberto G.; Rota, Matteo
title: COVID-19 incidence and mortality in Lombardy, Italy: an ecological study on the role of air pollution, meteorological factors, demographic and socioeconomic variables
date: 2021-01-22
journal: Environ Res
DOI: 10.1016/j.envres.2021.110777
sha: 36b480d5f8d317ba4e1da72be1d42b3416bcbdef
doc_id: 778701
cord_uid: 3em9vq1c

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 (NO(2)) 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 1,439 municipalities out of 1,507 (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 20(th) to April 16(th) and from March 1(st) to April 30(th), 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/m(3) 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/m(3) increase of annual mean PM2.5 concentration was associated with a 23% increase in all-cause mortality. An inverse association was found between NO(2) 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.

The Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) is responsible for the Coronavirus disease 2019 , firstly identified in Wuhan (China) in December 2019 and declared pandemic by the World Health Organization (WHO) on March 11 th , 2020 (World Health Organization, 2020). 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 30 th , 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 30 th , 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 4,000 critically ill 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 recent 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 (E. Bontempi, 2020 (E. Bontempi, , 2020 Coker et al., 2020; Comunian et al., 2020; Conticini et al., 2020; Filippini et al., 2020; Zoran et al., 2020a Zoran et al., , 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; E. Bontempi, 2020; Liang et al., 2020; X. Wu et al., 2020) . 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; 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 NO 2 on COVID-19 incidence and excess in all-cause mortality in Lombardy through an ecological approach, accounting for demographic, socioeconomic and meteorological variables.

This is an ecological study based on publicly available data (Agenzia Regionale Protezione Ambiente Lombardia, 2020a; Borruso, 2020; Istituto Nazionale di Statistica, 2020b 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.

J o u r n a l P r e -p r o o f 5 The two outcome variables were: i) the daily number of COVID-19 cases diagnosed through SARS-CoV-2 swab tests from February 20 th to April 16 th , 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) .

We considered the following demographic and socioeconomic variables: population size (the number of residents); population density (inhabitants per Km 2 ); 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 m 2 and over 4 occupants, or with 40-59 m 2 and over 5 occupants, or with 60-79 m 2 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, 2011).

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). J o u r n a l P r e -p r o o f 6 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, 2011).

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, 2011). 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.

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 7 cubic meter of air (g/m 3 ), 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.5x1.5 km 2 resolution by applying optimal interpolation techniques to measured data (Uboldi et al., 2008) . We computed average temperature and average absolute humidity values from February 15 th to April 15 th , 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 NO 2 spatial annual average concentrations were estimated over the Lombardy region on a 4x4 km 2 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 NO 2 concentrations were processed to estimate area-weighted mean concentration for each municipality (Agenzia Regionale Protezione Ambiente Lombardia, 2020a). Since PM10 and NO 2 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 NO 2 exposure as the mean of annual PM concentrations over the four years period 2016-2019.

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 20 th to April 16 th , 2020, and ii) the excess in all-cause mortality (model 2) from March 1 st to April 30 th , 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 NO 2 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 (log2transformed 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 J o u r n a l P r e -p r o o f employees in sports, entertainment and recreational activities per capita (categorized in tertiles); winter average temperature (°C) and absolute humidity (g/m 3 ) 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 R 2 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) .

Complete data were available for 1,439 out of 1,507 (95%) Lombardy municipalities. A total of 61,377 COVID-19 cases and 40,401 all-causes deaths were recorded in Lombardy from February 20 th to April 16 th and from March 1 st to April 30 th , 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 Figure 1 .

Incidence and mortality rates trends follow the same pattern and show a peak in the two last weeks of March (16 th -31 th ). The effect of regional and national restrictive measures for contrasting human-to-human transmission, particularly the lockdown, established on March 9 th , determined a sharp decline of the epidemic curve starting from April 1 st . 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 Km 2 . 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 J o u r n a l P r e -p r o o f 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.

The spatial plot in Figure 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 18 th 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.

The pattern of air pollution in Lombardy showed in Figure 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.

PM concentrations were associated with COVID-19 incidence at univariate analysis, while no association emerged for NO 2 (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 (Figure 4 ). An increase of 10 μg/m 3 in the mean annual PM2.5 and PM10 concentrations in 2016-2019 was associated with a statistically significant 58% and 34% increase in COVID-although the result was significant only within the model including PM2.5. The pseudo R 2 of the model was 0.71, showing a good fit of the model to the data.

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 ( Figure 5 Panel A and C), while an opposite trend emerged for the absolute humidity ( Figure 5 Panel B and D). A 10 ug/m 3 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. NO 2 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 R 2 was 0.43, showing a lower fit to the data as compared to the model fitted with COVID-19 as outcome.

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 (X. .

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 decreasing 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; 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 J o u r n a l P r e -p r o o f 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/m 3 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/m 3 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/m 3 increase in PM2.5 observed in the US ecological regression analysis by , who also controlled for 20 county-level covariates (X. . 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) . 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 .

Our study showed an inverse association between long-term exposure to NO 2 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 NO 2 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 NO 2 ,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 (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.

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. 

Agenzia Regionale Protezione Ambiente Lombardia, 2020a. Data repository

Agenzia Regionale Protezione Ambiente Lombardia, 2020b. I sistemi modellistici in ARPA -Aria / Qualità dell'Aria

ARPA.aspx?firstlevel=Modellistica (accessed 7.15.20)

Fitting Linear Mixed-Effects Models Using lme4

Increasing Temperature and Relative Humidity Accelerates Inactivation of SARS-CoV-2 on Surfaces

Data assimilation in atmospheric chemistry models: current status and future prospects for coupled chemistry meteorology models. Atmospheric Chemistry and Physics

Air quality in Europe -2020 report. Luxembourg: Publications Office of the European Union

Satellite-detected tropospheric nitrogen dioxide and spread of SARS-CoV-2 infection in Northern Italy

COVID-19 Lombardy ICU Network, 2020a. Risk Factors Associated With Mortality Among Patients With COVID-19 in Intensive Care Units in

COVID-19 Lombardy ICU Network, 2020b. Baseline Characteristics and Outcomes of 1591 Patients Infected With SARS-CoV-2 Admitted to ICUs of the Lombardy Region

Go slow to go fast: a plea for sustained scientific rigour in air pollution research during the COVID-19 pandemic

COVID-19 and climatic factors: A global analysis

Impact of the COVID-19 epidemic on the total mortality of the resident population in the first quarter of 2020 [WWW Document

Istat.it | Registro statistico delle imprese attive (ASIA -Imprese) [WWW Document

8mila Census -una selezione di indicatori per ogni comune di Italia

COVID-19 integrated surveillance data in Italy

Characteristics of COVID-19 patients dying in Italy

SARS-CoV-2 infection: the environmental endurance of the virus can be influenced by the increase of temperature

Effects of temperature and humidity on the spread of COVID-19: A systematic review

Ambient particulate air pollution and acute lower respiratory infections: a systematic review and implications for estimating the global burden of disease

Mortality impacts of the coronavirus disease (COVID-19) outbreak by sex and age: rapid mortality surveillance system

A general and simple method for obtaining R2 from generalized linear mixed-effects models

COVID-19 deaths in Lombardy, Italy: data in context

Assessing nitrogen dioxide (NO2) levels as a contributing factor to coronavirus (COVID-19) fatality

Comparisons between countries are essential for the control of COVID-19

Are Essential for Policy Guidance and Decisions

Assessing nitrogen dioxide (NO2) levels as a contributing factor to coronavirus (COVID-19) fatality

Elenco RSA Accreditate | Open Data Regione Lombardia [WWW Document

Georeferenziazione strutture | Open Data Regione Lombardia [WWW Document

SARS-Cov-2RNA found on particulate matter of Bergamo in Northern Italy: First evidence

ERS/ATS workshop report on respiratory health effects of household air pollution

Clarification of Misleading Perceptions of COVID-19 Fatality and Testing Rates in Italy: Data Analysis

Three-dimensional spatial interpolation of surface meteorological observations from high-resolution local networks

Methodological Considerations for Epidemiological Studies of Air Pollution and the SARS and COVID-19 Coronavirus Outbreaks

Is there an association between the level of ambient air pollution and COVID-19?

speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing-on-covid

Evolution of WHO air quality guidelines: past, present and future

Air pollution and COVID-19 mortality in the United States: Strengths and limitations of an ecological regression analysis

Effects of temperature and humidity on the daily new cases and new deaths of COVID-19 in 166 countries

Non-linear correlation between daily new cases of COVID-19 and meteorological factors in 127 countries

Assessing the relationship between surface levels of PM2.5 and PM10 particulate matter impact on COVID-19 in

Assessing the relationship between ground levels of ozone (O3) and nitrogen dioxide (NO2) with coronavirus (COVID-19