An ecological analysis of long-term exposure to PM2.5 and incidence of COVID-19 in Canadian health regions.

BACKGROUND: Ambient fine particulate matter (PM2.5) is associated with a wide range of acute and chronic health effects, including increased risk of respiratory infection. However, evidence specifically related to novel coronavirus disease (COVID-19) is limited.
METHODS: COVID-19 case counts for 111 Canadian health regions were obtained from the COVID-19 Canada Open Data portal. Annual PM2.5 data for 2000-2016 were estimated from a national exposure surface based on remote sensing, chemical transport modelling and ground observations, and minimum and maximum temperature data for 2000-2015 were based on a national interpolated surface derived from thin-plate smoothing splines. Population counts and sociodemographic data by health region were obtained from the 2016 census, and health data (self-rated health and prevalence of smoking, obesity, and selected chronic diseases) by health region, were obtained from the Canadian Community Health Survey. Data on total number of COVID-19 tests and changes in mobility comparing post-vs. pre-introduction of social distancing measures were available by province. Data were analyzed using negative binomial regression models.
RESULTS: After controlling for province, temperature, demographic and health characteristics and days since peak incidence by health region, long-term PM2.5 exposure exhibited a positive association with COVID-19 incidence (incidence rate ratio 1.07, 95% confidence interval 0.97-1.18 per μg/m3). This association was larger in magnitude and statistically significant in analyses excluding provinces that reported cases only for aggregated health regions, excluding health regions with less than median population density, and restricted to the most highly affected provinces (Quebec and Ontario).
CONCLUSIONS: We observed a positive association between COVID-19 incidence and long-term PM2.5 exposure in Canadian health regions. The association was larger in magnitude and statistically significant in more highly affected health regions and those with potentially less exposure measurement error. While our results generate hypotheses for further testing, they should be interpreted with caution and require further examination using study designs less prone to bias. Crown

Introduction

There is considerable evidence from multiple lines of research including toxicology, human clinical studies and epidemiological studies that air pollution in general and fine particulate matter (PM2.5) more specifically are associated with a wide range of acute and chronic health effects (Thurston et al., 2017). Based on estimates from the Global Burden of Disease initiative, PM2.5 accounts for the greatest burden of mortality of any environmental exposure, accounting for approximately 3 million worldwide deaths annually (GBD, 2017 Risk Factor Collaborators, 2018). It is well established that acute exposure increases the risk of emergency visits and hospital admissions for respiratory infections including pneumonia (Atkinson et al., 2014; Domingo and Rovira, 2020; Nhung et al., 2017). There is also a growing body of evidence that long term exposure increases the risk of morbidity and mortality from respiratory infection (Mehta et al., 2013; Neupane et al., 2010).

Evidence specifically related to PM2.5 and novel coronaviruses such as severe acute respiratory syndrome coronavirus (SARS-CoV-1) and middle east respiratory syndrome (MERS) is limited. Studies based on the SARS-CoV-1 outbreak suggest that meteorology and exposure to air pollution increased transmission (Cai et al., 2007) and worsened patient prognosis (Kan et al., 2005). Notably SARS-CoV-1 patients from more polluted regions were twice as likely to die as those in less polluted places (Cui et al., 2003). There is also evidence that air pollution exposure more generally adversely affects respiratory immune defences (Domingo and Rovira, 2020; Yang et al., 2020), and emerging evidence suggesting that novel coronavirus disease (COVID-19) incidence and mortality may be increased in relation to both acute (Zhu et al., 2020) and chronic exposure (Andree, 2020; Liang et al., 2020; Ogen, 2020; Wu et al., 2020). This evidence implies that deterioration in air quality over short time periods (e.g. from wildfire smoke, other local burning, specific meteorological events such as temperature inversions) may lead to more cases of severe COVID-19 infections, adding further demand to the healthcare system. Conversely, improving air quality by reducing both the occurrence of acute events and long term average concentrations, may help to protect communities from COVID-19 and reduce the burden on hospitals.

In this study, we conduct an ecological analysis of COVID-19 cases and 17 year average PM2.5 concentrations among Canadian health regions. While ecological analyses have many limitations which preclude attribution of cause and effect, they can be readily conducted once data are available and permit the generation of hypotheses to be more rigorously examined in subsequent studies.

Materials and methods

COVID-19 case counts compiled from publicly available reports for 111 health regions were obtained from the COVID-19 Canada Open Data portal (Berry et al., 2020). Health regions are defined by provincial ministries of health; in some jurisdictions, they correspond to areas served by local public health departments or authorities (Statistics Canada, 2018a). The provinces of British Columbia and Saskatchewan reported counts for groups of health regions (up to four per group), which we distributed to individual health regions in proportion to population. Similarly, health region was not reported for many cases in Nova Scotia. These were distributed among four health regions in proportion to population. In an attempt to account for stage of outbreak, response to distancing measures, and correlation between case counts and deaths, we also obtained data on days elapsed since the first case, days since peak daily incidence of new cases, and deaths at the health region level (Berry et al., 2020), as well as date of declaration of public health emergency or state of emergency at the provincial level (Boire-Schwab et al., 2020). Total number of COVID-19 tests was only available by province. Annual PM2.5 data for 2000–2016 were derived from a surface combining a 0.01° × 0.01° (approximately 1 km × 1 km) resolution Aerosol Optical Depth (AOD) retrieval from the National Aeronautics and Space Administration (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS), Multi-angle Imaging SpectroRadiometer (MISR), and Sea-Viewing Wide Field-of-View Sensor (SeaWiFS) instruments, Goddard Earth Observing System chemical transport model (GEOS-Chem) simulations, and ground observations (van Donkelaar et al., 2019). Values prior to 2004 were temporally adjusted (Boys et al., 2014). These data have been used extensively in air pollution epidemiology studies in Canada (Pappin et al., 2019). Greenness data for 2000–2019 were based on growing season maximum normalized difference vegetation index (NDVI) data from the MODIS onboard the Terra satellites (Didan, 2015; Gorelick et al., 2017). NDVI values vary from −1 to 1, negative values indicating features such as water and positive values indicating vegetation. We employed positive values only in calculating average greenness by health region. Temperature data for 2000–2015 (annual minimum of lowest monthly maximum temperature – henceforth referred to as “minimum temperature,” and annual maximum of highest monthly minimum temperature – henceforth referred to as “maximum temperature”) were based on a national interpolated surface of available observations derived using thin-plate smoothing splines, implemented in Australian National University Spline (ANUSPLIN) climate modeling software (Wang et al., 2018). All exposure data were available by 6 character postal codes (CANMAP 2015), which were mapped to 2018 health region boundaries in R (GISTools (Brunsdon and Chen, 2014), rgdal (Bivand et al., 2019) and raster (Hijmans, 2020) packages). Postal codes are used by Canada Post corporation for mail delivery, and are analogous to American zip codes – there are currently approximately 875,000 postal codes in Canada; as they are point locations, they do not have a specified surface area (CANUE, 2018). Maps were generated using the ggplot2 (Wickham, 2016) and broom (Robinson and Hayes, 2020) packages. Population counts and sociodemographic data by health region on percent of population 65 or older, with income less than the low income cutoff (LICO), and Black (Adams et al., 2020; Dyer, 2020; Yancy, 2020), were obtained from the 2016 census (Statistics Canada, 2018b). LICOs are defined as income levels below which families spend a disproportionate share of their income on necessities, and are family-size and community-size specific (Statistics Canada, 2016). Health data by health region on factors thought to increase susceptibility to COVID-19 (Adams et al., 2020) (percent of population who rate health as fair or poor, are daily or occasional smokers, overweight, obese, have asthma, chronic obstructive pulmonary disease (COPD), hypertension or diabetes), were based on data from the 2017 and 2018 Canadian Community Health Survey (CCHS) (Statistics Canada, 2020), which is an annual national cross-sectional survey of individuals 12 years of age and over. Changes in mobility by province comparing post- vs. pre-introduction of social distancing measures were based on aggregated data from Google Account users who opted-in to location history for their account (Google LLC, 2020). Correlations among variables were examined using a correlogram based on Spearman correlations (Harrell et al., 2020, Wei and Simko, 2017).

Data were analyzed using negative binomial regression models, specifying PM2.5 and covariates (including province) as fixed effects. Log population was included as an offset. PM2.5 and covariates were first regressed individually vs. case counts, then those exhibiting statistically significant associations were included in multivariate models. Prevalence of asthma, COPD, hypertension and diabetes were excluded from multivariate models with PM2.5, since they could be intermediate in a putative causal pathway with COVID-19 incidence. The most parsimonious model was selected based on the Akaike Information Criterion (AIC) (Akaike, 1974). Presence of residual spatial autocorrelation was examined by mapping model residuals and computing Moran's I (Bivand and Wong, 2018). Sensitivity analyses were conducted by excluding Montréal (which accounted for 27.5% of cases but only 5.5% of population); excluding British Columbia, Saskatchewan and Nova Scotia (which reported cases counts only for aggregated health regions); excluding health regions with less than median population density (with presumably greater exposure measurement error over larger, more sparsely populated health regions); restricting the analysis to Ontario and Quebec, the two provinces with the highest incidence, attributed in part to provincial level policies on testing (Weeks, 2020), and timing of school vacation periods (Perreaux, 2020) respectively; and specifying province as a random rather than fixed effect. Analysis was conducted in R (R Core Team, 2019) using the lme4 package (Bates et al., 2015). Research ethics board approval was not required because all data were publicly available and aggregated at health region level. The work described here has been carried out in accordance with the Uniform Requirements for manuscripts submitted to Biomedical journals.

Results

There were 73,390 cases up to May 13, 2020 and overall incidence was 208.8 cases/100,000 (2016 population). Incidence was highest in Quebec (489.0/100,000) and Ontario (166.5/100,000) and lowest in Nunavut (0/100,000) and the Northwest Territories (12.0/100,000) (Fig. 1 ). The ten health regions with the highest incidence were distributed among multiple provinces, and included three large cities (Montréal, Toronto, Calgary), regions surrounding Montréal (Laval, Lanaudière, Montérégie, Laurentides, Mauricie/Centre-du-Québec), and one remote northern area (Far North, Saskatchewan) (Supplementary Table S1). Deaths were very strongly correlated with case counts by health region (R2 = 0.95). PM2.5 concentrations averaged 6.1 (standard deviation 2.1) and were highest in urban areas as well as more generally in southern Ontario and Quebec (Fig. 2 ). Variability in other exposures, sociodemographic and health characteristics by health region are summarized in Table 1 . There was little variability in days since declaration of emergency. A correlogram is provided in Fig. S2. PM2.5 was strongly positively correlated with percent Black, minimum and maximum temperature. Population density exhibited strong positive correlations with percent Black and low income, PM2.5, minimum and maximum temperature, as well as strong negative correlations with prevalence of obesity and smoking. Prevalence of obesity was strongly positively correlated with prevalence of diabetes, hypertension and smoking, while percent age 65 and older was strongly positively correlated with prevalence of hypertension and COPD. Changes in mobility during the period March 13 (the date after which social distancing recommendations were introduced in most provinces)-April 26 vs. January 3-February 6, were generally comparable among provinces (Supplementary Table S3). Notably, changes in visits/length of stay were somewhat more modest in Nunavut, which had no cases, and visits/length of stay in parks increased substantially in British Columbia, the most populous regions of which typically experience milder weather during the March/April period compared to elsewhere in Canada.

Fig. 1

COVID-19 incidence by health region. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 2

Long term average PM2.5 concentrations by health region. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Table 1

Summary of sociodemographic characteristics, health measures and exposures by health region.

	Minimum	25th %ile	Median	Mean	75th %ile	Maximum
Population (n)	2632	76,626	169,244	316,682	421,538	2,731,571
Population density per km2	0.0	2.2	10.6	212.4	39.0	4848.3
Population density per km2a	10.6	23.2	39.0	418.4	227.9	4848.3
Incidence/100,000	0.0	20.4	54.5	114.0	118.8	1041.8
Days since first case	24	54	59	61	63	109
Days since emergency declared	52	56	57	57.	57	61
Days since peak incidence	4	23	39	35	44	109
Age ≥65 (%)	3.8	15.2	18.0	17.8	21.6	26.4
Black (%)	0.0	0.5	0.7	1.5	1.7	9.5
< Low Income Cut-off (%)	2.5	4.9	6.3	6.9	7.6	18.1
Poor or fair self-rated health (%)	7.0	11.1	12.6	12.9	14.3	21.9
Overweight (%)	28.5	33.7	35.5	35.8	38.0	46.1
Obese (%)	12.1	28.1	32.3	31.9	36.7	47.2
Asthma (%)	4.5	7.2	8.5	8.4	9.5	14.1
COPD (%)	0.0	3.4	4.8	4.8	6.3	9.1
Hypertension (%)	12.1	17.1	19.0	19.4	21.6	29.8
Diabetes (%)	2.7	6.4	7.8	8.0	9.4	14.4
Daily or occasional smoker (%)	8.8	15.7	18.3	18.9	20.6	63.1
Physically active (%)	38.0	50.1	55.0	54.7	58.6	71.9
PM2.5 (μg/m3)	2.0	4.2	6.0	6.1	7.5	11.5
Minimum temperature	−34.9	−24.6	−18.0	−17.8	−12.4	−0.9
Maximum temperature	10.3	17.5	19.4	19.2	21.4	23.8
NDVIb>0a	0.54	0.59	0.62	0.62	0.65	0.66

For health regions with population density ≥ median.

Normalized Difference Vegetation Index.

Results of regression models are summarized in Table 2 . In bivariate models, PM2.5, minimum and maximum temperature, percent Black population, percent of population < LICO, population density, and days since first case were significantly positively associated with COVID-19 counts, while days since peak incidence, percent of population 65 or older, and prevalence of asthma, COPD, hypertension, overweight and obesity exhibited significant negative associations. COVID-19 counts were not significantly associated with total number of tests or percent change in mobility to work locations by province (not shown). In the best fitting multivariate model, minimum temperature exhibited a significant positive association, and percent age 65 and older and days since peak incidence exhibited a significant negative association with COVID-19 counts. PM2.5 and percent low income remained positively associated with COVID-19 counts, and prevalence of overweight and obesity remained negatively associated, but the associations were no longer significant. Additional multivariate models are summarized in Table S4. The association of PM2.5 with COVID-19 counts was not sensitive to the sequential addition to the best-fitting model of population density, percent Black population, or days since first case; the PM2.5 coefficient changed by <10% and additional covariates were not significantly associated with COVID-19 counts. However, the association of PM2.5 with COVID-19 counts was sensitive to the removal of minimum temperature from the best-fitting model (IRR increased to 1.10, 95% CI 1.00–1.21, p < 0.05). Mapping of model residuals and Moran's I did not indicate presence of residual spatial autocorrelation.

Table 2

Summary of regression model results.

Variable	Bivariatea		Multivariatea
	IRR	95% CI	IRR	95% CI
PM2.5 (per μg/m3)	1.25	1.13–1.38	1.07	0.97–1.18
Minimum temperatureb	1.64	1.23–2.19	1.42	1.05–1.93
Maximum temperatureb	1.42	1.08–1.86
Population densityb	1.40	1.19–1.63
Percent age 65+b	0.57	0.48–0.69	0.72	0.59–0.88
Percent < LICOb	1.57	1.33–1.84	1.12	0.91–1.38
Percent Blackb	1.49	1.27–1.75
Percent asthmab	0.81	0.68–0.97
Percent COPDb	0.61	0.52–0.72
Percent hypertensionb	0.64	0.52–0.79
Percent diabetesb	1.03	0.84–1.27
Percent physically activeb	0.92	0.70–1.19
Percent overweightb	0.78	0.65–0.93	0.86	0.71–1.03
Percent obeseb	0.71	0.57–0.87	0.90	0.68–1.20
Percent smokersb	0.97	0.72–1.30
Days since first caseb	1.52	1.24–1.86
Days since peak incidenceb	0.60	0.50–0.72	0.71	0.59–0.84
NDVIc	0.58	0.44–0.77

Includes log(population) as offset and province as factor.

scaled.

only health regions with population density greater ≥ median.

In sensitivity analyses (Table 3 ), the association of PM2.5 with COVID-19 incidence was unchanged when Montréal was excluded or when province was included as a random effect rather than a fixed effect. In models excluding British Columbia, Saskatchewan and Nova Scotia, excluding health regions with less than median population density, and restricted to Ontario and Quebec, the association of PM2.5 with COVID-19 incidence was larger in magnitude than the national analysis, and statistically significant. These analyses accounted for 85–95% of cases and 62–82% of population. Notably, in the subset of health regions with population density greater than the median, greenness exhibited a significant negative association with case counts on its own, but was no longer significant in a multivariate model which had a larger AIC than the final model. Results were also similar when based on data from earlier in the course of the pandemic when there were 35,986 cases (April 19).

Table 3

Summary of sensitivity analyses of association of PM2.5 and COVID-19 incidence.

Model	IRR	95% CI	% cases	% population
Best fitting multivariate model (province as fixed effect)	1.07	0.97–1.18	99.9a	99.6a
Exclude Montréal	1.07	0.96–1.18	72.4	94.5
Exclude 3 provinces with aggregated health regions	1.15	1.00–1.32b	94.6	81
Exclude health regions with population density < median	1.16	1.00–1.34b	94.8	82
Ontario and Quebec only	1.21	1.03–1.40	84.9	61.5
April 19 data (35,986 cases)	1.10	1.01–1.19	99.9	99.6
Best fitting multivariate model (province as random effect)	1.09	0.98–1.21c	99.9a	99.6a

Five health regions excluded due to missing data.

Lower bound >1 with additional decimal places.

Corrected for overdispersion by multiplying standard error by square root of dispersion factor (Bolker, 2020).

Discussion

We found that after controlling for province, temperature, health and demographic characteristics, and time since peak incidence by health region, long-term PM2.5 exposure exhibited a non-statistically significant positive association with COVID-19 incidence among 111 Canadian health regions. This association was larger in magnitude and statistically significant in models excluding provinces which reported cases only for aggregated health regions, excluding health regions with less than median population density, and restricted to Ontario and Quebec, in each case accounting for a majority of COVID-19 cases. The association we observed parallels that reported in an analysis of PM2.5 and COVID-19 mortality in approximately 3000 US counties (Wu et al., 2020). While this is striking, especially given the lower COVID-19 incidence and narrower range of PM2.5 exposure in Canada compared to the US, findings should be interpreted with caution; the analyses examined different outcomes (incidence vs. mortality) and associations with PM2.5 in both studies could result from similar biases related to ecologic analyses. PM2.5 exposures and other risk factors for COVID-19 incidence can vary over a smaller scale than health region, particularly for large health regions. Greenness in particular varies at a much smaller scale than health region, even in predominantly urban health regions. Our analysis entailed coarser spatial resolution and a much smaller number of health regions than US counties. However, analyses excluding provinces that reported cases only for aggregated health regions, excluding health regions with less than median population density, and restricted to Quebec and Ontario resulted in a larger magnitude and statistically significant association, which could reflect reduced exposure measurement error for these generally smaller, more densely populated health regions. The US study also found that the mortality risk ratio was elevated in an analysis restricted to urban areas, but it was no longer significant (Wu et al., 2020). Owing to the smaller sample, we were also constrained in the number of covariates that could be included in multi-variate models. Still, we did account for unmeasured confounders by province by specifying it as a fixed or random effect, and found that the magnitude of the association was larger in more highly urban and more highly affected areas. Mobility data indicated that effects of social distancing restrictions were generally similar among provinces. We employed publicly available laboratory confirmed case counts, thus our results are not generalizable to milder or asymptomatic cases which did not undergo laboratory confirmation, and case criteria may differ by province and change over time. Data on COVID-19 testing were only available at the provincial level. Although criteria for testing are standardized, access may differ by health region, which could increase the observed incidence in more highly urban areas with potentially better access and generally higher PM2.5 concentrations, potentially resulting in a spurious association between PM2.5 exposure and COVID-19 incidence. While COVID-19-related hospitalizations and deaths may be less affected by artefacts introduced by differential eligibility for and/or access to testing by health region, the number of hospitalizations and deaths in Canada is relatively small. At the health region level, a large number of small or zero counts would be expected, which introduces other limitations in the analysis. In any case, hospitalization data were not available nationally at the health region (or smaller) level. Our observation that deaths were very strongly correlated with case counts by health region provides evidence that differences in case numbers by health region were not driven by differences in testing, which would be expected to introduce more scatter in the relationship between deaths and cases. Finally, associations observed based solely on area level measures may not exist at the individual level (ecological fallacy).

The magnitude of the association both in our study and the US study was several fold larger per unit PM2.5 than hazard ratios typically observed in cohort studies of mortality (Crouse et al., 2015). It is premature to speculate on pathophysiological mechanistic explanations for this observation. While there is existing evidence that PM2.5 increases the risk of respiratory infections, including the SARS virus, the most likely explanation for the large magnitude association is residual confounding by unmeasured factors.

In addition to the Wu et al. (2020) study, another American ecological study based on approximately 3000 counties found that nitrogen dioxide (NO2) was significantly associated with COVID-19 mortality and case-fatality rates, while PM2.5 exhibited a marginally significant association with mortality (Liang et al., 2020). A significant positive association of PM2.5 with COVID-19 incidence and hospital admissions was also reported in the Netherlands (Andree, 2020). Ogen (2020) reported co-location of high NO2 concentrations estimated from remote sensing and high COVID-19 mortality counts in northern Italy and Madrid. However, the analysis did not account for the underlying population at risk, population density, timing of onset of cases or introduction of control measures, or sociodemographic or health characteristics. In a time-series study in 120 Chinese cities, significant positive associations were observed between short term (two week) exposure to PM2.5, NO2 and ozone and COVID-19 incidence (Zhu et al., 2020).

We also found that percent Black population was positively associated with COVID-19 incidence (although it was not included in the final multivariate model). This is consistent with reports from the US of disproportionately high COVID-19 incidence and mortality among Black Americans (Dyer, 2020; Yancy, 2020). Percent of population less than LICO was also positively associated with COVID-19 incidence. Those with lower income have more barriers to self isolating and social distancing. Inequality in access to care affecting access to COVID-19 testing and outcome following infection would be expected to be less in Canada with universal healthcare. Negative associations of percent of population age 65 and older, prevalence of asthma, COPD, hypertension, overweight and obesity with COVID-19 incidence may reflect a higher prevalence in predominantly rural health regions with few cases.

Subsequent studies should examine effects of both acute and chronic exposure to PM2.5, and air pollution more generally, on COVID-19 morbidity and mortality using more highly spatially resolved area level data, as well as individual level data to determine whether results from ecological analysis are borne out.

Conclusions

We report an analysis of COVID-19 incidence by Canadian health region. After controlling for temperature, demographic and health characteristics and days from peak incidence by health region, long-term PM2.5 exposure exhibited a positive association with COVID-19 incidence, paralleling results of a recent American analysis of mortality. The association was larger in magnitude and statistically significant in more highly affected health regions and those with potentially less exposure measurement error. While parallels with US findings are noteworthy, and generate hypotheses for further testing, results should be interpreted with caution and require further examination using study designs less prone to bias.

Credit author statement

David M. Stieb: Conceptualization, Methodology, Formal analysis, Data curation, Writing - original draft, Writing - review & editing; Greg J. Evans: Conceptualization, Methodology, Writing - original draft, Writing - review & editing, Teresa M. To: Conceptualization, Methodology, Writing - original draft, Writing - review & editing, Jeffrey R. Brook: Conceptualization, Methodology, Data curation, Writing - original draft, Writing - review & editing, Richard T. Burnet: Conceptualization, Methodology, Writing - original draft, Writing - review & editing

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.