The association of wildfire air pollution with COVID-19 incidence in New South Wales, Australia

The 2020 COVID-19 outbreak in New South Wales (NSW), Australia, followed an unprecedented wildfire season that exposed large populations to wildfire smoke. Wildfires release particulate matter (PM), toxic gases and organic and non-organic chemicals that may be associated with increased incidence of COVID-19. This study estimated the association of wildfire smoke exposure with the incidence of COVID-19 in NSW. A Bayesian mixed-effect regression was used to estimate the association of either the average PM10 level or the proportion of wildfire burned area as proxies of wildfire smoke exposure with COVID-19 incidence in NSW, adjusting for sociodemographic risk factors. The analysis followed an ecological design using the 129 NSW Local Government Areas (LGA) as the ecological units. A random effects model and a model including the LGA spatial distribution (spatial model) were compared. A higher proportional wildfire burned area was associated with higher COVID-19 incidence in both the random effects and spatial models after adjustment for sociodemographic factors (posterior mean = 1.32 (99% credible interval: 1.05–1.67) and 1.31 (99% credible interval: 1.03–1.65), respectively). No evidence of an association between the average PM10 level and the COVID-19 incidence was found. LGAs in the greater Sydney and Hunter regions had the highest increase in the risk of COVID-19. This study identified wildfire smoke exposures were associated with increased risk of COVID-19 in NSW. Research on individual responses to specific wildfire airborne particles and pollutants needs to be conducted to further identify the causal links between SARS-Cov-2 infection and wildfire smoke. The identification of LGAs with the highest risk of COVID-19 associated with wildfire smoke exposure can be useful for public health prevention and or mitigation strategies.


H I G H L I G H T S
• Estimated the association of  with PM 10 level and wildfire burned area extent • Wildfire burned area was used as an alternative indicator of wildfire smoke exposure. • Bayesian and spatial regression models used to estimate the associations and map risk • An association of wildfire burned area with COVID-19 incidence in NSW was found. • Identification of higher risk areas can be used for public health strategies.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
The 2020 COVID-19 outbreak affected all Australian states and territories, with New South Wales (NSW) having a notification rate of 57.3 per 100,000 persons in the context of the national rate of 110.4 per 100,000 persons in 2020 (Australian Department of Health, 2021). The NSW COVID-19 outbreak began in March 2020, after an unprecedented wildfire season where more than 18 million ha of land were burned between September 2019 and March 2020. Air pollution has been associated with a higher risk of COVID-19 in other wildfire affected regions such as California with a 57% increase in cases after about 4 million ha of forest were lost in the 2020 wildfires (Meo et al., 2021). Previous research has identified the transmission of coronaviruses via dust particles and particulate matter (Sedlmaier et al., 2009;Zhao et al., 2019), and a recent study evidenced the presence of SARS-CoV-2 RNA on outdoor PM, suggesting that it can enhance the persistence of the virus in the atmosphere (Setti et al., 2020). Wildfire smoke typically contains a complex mixture of particulate matter (PM), toxic gases, inorganic elements and ionic constituents, and organic compounds such as polycyclic aromatic hydrocarbons, methoxyphenol, alkanes and levoglucosan generated from the burning of biomass fuels (Rager et al., 2021). Many of these agents have been associated with altered immune responses that can play a role in increasing the risk of coronavirus infections (Velazquez-Salinas et al., 2019).
An increased incidence of respiratory infections has been associated with exposures to pollutants occurring in wildfire smoke, including particulate matter (PM) and nitrogen dioxide (NO 2 ) (Cohen et al., 2017;Forouzanfar et al., 2016). Wildfire pollutants also include gases such as methane (CH 4 ) and carbon dioxide (CO 2 ), and precursors of ozone (O 3 ), secondary inorganic aerosols and carbonaceous particles that are potentially associated with respiratory infections (Yang et al., 2021b;Zhang et al., 2019;Zhang et al., 2013). Recent studies have found a higher risk of SARS-Cov-2 transmission in areas with exposure to high levels of CO, PM and O 3 in the United States (Meo et al., 2020;Meo et al., 2021). These exposures are associated with increased number of COVID-19 cases during and after the wildfires with a higher risk of COVID-19 infection after two weeks of exposure to PM 2.5 , and a cumulative lag-effect of PM 10 , PM 2.5 , carbon monoxide (CO) and NO 2 (Wang et al., 2020;Zhu et al., 2020). The mechanisms to explain the higher incidence of COVID-19 after weeks of exposure are not yet known although multiple studies have found altered levels of leukocytes and immunoregulators such as interleukin 6 and interferon beta in people and animals weeks and or months after exposure to wildfire smoke and PM Ma et al., 2017;Zhu et al., 2021).
Whereas the exposure to wildfire smoke can determine an increased risk of viral infections, the association of wildfire air pollution with the outbreak of COVID-19 after the 2019-2020 wildfires in NSW has not yet been investigated. Previous studies in Australia have estimated an increased rate of respiratory emergency attendances and hospital admissions during and after wildfires (Chen et al., 2006;Crabbe, 2012;Ryan et al., 2021;Tham et al., 2009). These studies used data from air quality monitoring stations that is limited in terms of the range of air pollutants measured and the few areas covered. NSW has only 57 high standard air monitoring stations for PM 10 and gases occurring in wildfire smoke, covering an area of 801,150 km 2 . Despite the geographical extension of the state, air monitoring stations are located mostly in highly populated areas. Although air monitoring measures in NSW are limited in terms of the relatively small number of monitoring stations, the geographic extent of the wildfires can be used as an alternative indicator of wildfire smoke exposure. Satellite images can be used to predict PM concentrations using nonlinear models (Yang et al., 2021a), and have been incorporated to develop satellitederived air pollution and climate indicators by efforts that include the Lancet Commission on Pollution and Health and the NASA Health and Air Quality Applied Science Team (Anenberg et al., 2020;Reid et al., 2015). This study aims to investigate whether there was an association between either PM 10 levels or the geographic extent of the 2019-2020 wildfires with the incidence rate of COVID-19 infection in NSW.

Methods
The study followed an ecological design, using the NSW Local Government Areas (LGA) as the ecological units. A cross-sectional analysis was used to assess the association of the LGA-level COVID-19 incidence rate (ir-COVID-19) with either the percentage of wildfire burned area (WBA) or the average concentration of PM 10 (avg-PM 10 ) in the 2019-2020 bushfire season, as two different proxies of wildfire exposure.

Data
The ir-COVID-19 was calculated as the number of new locally acquired COVID-19 cases in the first 5 months of the outbreak (March 2-August 4, 2020) per 100,000 people per LGA of residence, using data from the NSW Department of Health. Data on PM 10 from monitoring stations using the National Environment Protection Measure standard, for the period Sept-2019-March-2020 was obtained from the Department of Planning Industry and Environment, DPIE (DPIE, 2021). The levels of PM 10 associated with wildfires in Australia are highly correlated with other PM, especially PM 2.5 , and have been identified as a consistent indicator of the effects of wildfires . In addition, PM 10 is the most measured air pollution indicator in NSW with significantly more monitoring stations compared to those measuring PM 2.5 . The avg-PM 10 was calculated with an inverse distance weighting interpolation model in ArcMap v10.6.
The WBA was calculated in ArcMap using high resolution raster data (satellite images, pixel size: 10 m) from the fire extent and severity mapping data for the same period (DPIE, 2020). The COVID-19 records did not include individual characteristics therefore the ir-COVID-19 could not be age standardised. To adjust for age and sex, the percentage of population by sex, and age groups 0-15 years and 65+ years was calculated per LGA using the Australian Bureau of Statistics (ABS) census data (ABS, 2021). Population density was estimated as the proportion of total population per km 2 per 100,000 people per LGA. Other sociodemographic data included the Index of Relative Socioeconomic Disadvantage (IRSD) per LGA (ABS, 2021) to adjust for socio-economic status that is considered a significant confounder in analyses of air pollution associated with COVID-19 (Chakrabarty et al., 2021). There were 129 LGAs with an average size of 6208 km 2 (quartile-1: 334 km 2 , quartile-3: 7141 km 2 ).

Statistical analysis
A mixed effects Bayesian regression model was used to measure the association of the ir-COVID-19 with either the avg-PM 10 or WBA, adjusting for sociodemographic and environmental covariates. Bayesian models are a robust approach to produce reliable estimates in epidemiological studies that consider health outcomes distributed at geographical area levels and can provide unobserved parameters of interest using the posterior probability distribution (Kang et al., 2016). The models were run in R using the R-INLA package that implements a computationally efficient approach comparable to Markov Chain Monte Carlo methods (Rue et al., 2009). Collinearity between the predictors can affect the prediction of the ir-COVID-19 therefore a variance inflation factor test was used to verify there was no multicollinearity (Menard, 1995). For the i-th LGA, the number of new COVID-19 cases was modelled as with the linear predictor defined on the logarithmic scale for each exposure (avg-PM 10 and WBA): where α is the intercept; exp is either the avg-PM 10 or WBA percentage; pop is the population, included as the offset; υ i~N (0, σ 2 ) is the random effect (i-th LGA); and X represents the list of independent variables (population density, IRSD, and percentages of females, population 0-15 years and 65+ years) with their respective regression coefficients β x . Population density was included on the logarithmic scale to facilitate the interpretation of the regression coefficients.
As there is important variability in the distribution of the LGAs with a higher concentration of smaller and densely populated LGAs in the greater Sydney region (Fig. 1a), a model accounting for the spatial structure of the LGAs was defined as where s represents the spatial structure and u represents the unstructured (i.e. random effects) component according to the Besag-York-Mollie specification (Besag et al., 1991). The spatial structure was defined as an adjacency matrix using a queen specification (i.e., including all neighbours sharing a border with each LGA). This spatial specification is considered optimal (Duncan et al., 2017) and has been tested for Australian LGAs (Cortes-Ramirez et al., 2021). All models implemented the default minimally informative priors (Blangiardo et al., 2013;Rue et al., 2009). Associations were identified if the regression coefficient credible intervals (CI) did not cross the null value 1. In a preliminary analysis, the Poisson regression was compared against a negative binomial and a zero inflated Poisson regressions to assess overdispersion and zero inflation. All predictors were scaled and the models were compared using predictive information criteria to identify the best fit (Millar, 2009; Morales-Otero and Núñez-Antón, 2021). There were no differences in the associations found in all models compared. Both Poisson models (Random Effects and Spatial models) had the best fit compared to all other models therefore the Poisson regression was identified as the best model to estimate the association of ir-COVID-19 with either avg-PM 10 or WBA, adjusting for sociodemographic covariates (Supplementary material). Finally, the LGA-specific relative risk compared to the whole of NSW and the probability of the risk to be greater than 1 were mapped using the R-package T-map (Tennekes, 2018) since area-specific estimates are often mapped to assess the spatial distribution of risk and the probability of increased risk in specific geographical areas (Richardson et al., 2004). This study was approved by the Queensland University of Technology Human Research Ethics Committee (11-11-2020).

Results
The ir-COVID-19 was higher in most LGAs within the metropolitan Sydney region and eastern LGAs alongside the coast, compared to the rest of the state (Fig. 1a). Fig. 1b and c shows the distribution of the 2019-2020 wildfires, the WBA, and the avg-PM 10 , per LGA. The descriptive statistics of the ir-COVID-19 and the predictors included in the analysis are shown in Table 1. Table 2 shows the fixed effects estimated in the regression models. An association between WBA and the ir-COVID-19 was found in both the Random Effects and Spatial models. No association of avg-PM 10 with the ir-COVID-19 was identified. Using the unscaled mean of WBA, the regression coefficients are interpreted as that each 1.8% increase in the WBA was associated with 32% and 31% increased risk of COVID-19 in the Random Effects and Spatial models respectively, after adjusting for sociodemographic factors. There were no important differences in the goodness-of-fit between the models with WBA and avg-PM 10 , and very small goodness-of-fit differences between the Random Effects and Spatial models for either WBA or avg-PM 10 (DIC values were close). A higher population density was associated with increased ir-COVID-19 in all models. Fig. 2 shows the spatial distribution of the LGA-specific relative risk of COVID-19 compared to the whole of NSW after adjustment for avg-PM 10 and WBA and sociodemographic risk factors respectively (a and c); and the probability of the risk to be greater than 1, in the Spatial Models with PM 10 and WBA (b and d). A greater ir-COVID-19 increase (LGA-specific risk >2.5 and CI not crossing 1) associated with avg-PM 10 was found in LGAs in the greater Sydney region (Blue Mountains, Penrith, Waverley) and the Hunter region (Mid-Coast). An increased probability of a positive association of avg-PM 10 with ir-COVID-19 was identified in LGAs close to coastal areas, with the highest probability (>98%) in Shoalhaven, Liverpool, Penrith, Waverley, Blue Mountains and Mid-Coast. A greater association of WBA with ir-COVID-19 (LGA-specific risk >2.5 and CI not crossing 1) was found in LGAs in the greater Sydney region (Penrith, Waverley) and the Hunter region (Mid-Coast). A higher probability of a positive association between WBA and ir-COVID-19 was identified in LGAs in the Southeast and Tablelands, greater Sydney, Hunter and New England Norwest regions with the highest probability (>98%) in Liverpool, Penrith, Waverley and Mid-Coast.

Discussion
We estimated an increase of the COVID-19 incidence rate (ir-COVID-19) associated with a higher percentage of wildfire burned area (WBA) in NSW, after adjusting for sociodemographic factors. We did not find evidence of an association between higher average PM 10 levels and the ir-COVID-19. Using a spatial regression, we identified the local government areas (LGAs) in the greater Sydney and Hunter regions showing a stronger association between WBA and ir-COVID-19 compared with the whole of NSW. Our findings may indicate that exposure to pollutants in wildfire smoke increases the risk of COVID-19 infection. Toxic gases and airborne particles released in wildfires have been associated with higher incidence of COVID-19 in other regions severely affected such as California (Meo et al., 2021). The increased risk of COVID-19 in areas affected by the unprecedented 2019-2020 wildfire season in NSW identified in our study concur with previous research that estimates the association of wildfire pollutants with SARS-Cov-2 infection (Meo et al., 2020;Meo et al., 2021).
We used high resolution images (pixel size: 10 m) to calculate the proportional burned area as an indicator of wildfire smoke, a known approach to estimate wildfire air pollution, different to other indicators such as PM measures (Anenberg et al., 2020;Reid et al., 2015). Using PM measures in NSW is challenging given the larger and uneven geographical extension of the LGAs and the scant and highly concentrated distribution of monitoring stations around the most populated areas. A measure such as the WBA is an alternative indicator of the mix of pollutants in wildfire smoke (Reid et al., 2015). Our findings identified an increased incidence of COVID-19 in NSW areas where polluting agents occurring in wildfires would have higher concentrations. Wildfires air pollutants include gases and chemical compounds such as polycyclic aromatic hydrocarbons, in addition to airborne particles, that can increase the risk of respiratory infections (Rager et al., 2021). Epidemiological studies have consistently estimated the association of COVID-19 with wildfires air contaminants, especially PM (Chakrabarty et al., 2021;Navarro et al., 2021;Wang et al., 2020). Although we did not find a statistical association of PM 10 concentration with the COVID-19 incidence, this might be because the PM 10 levels were calculated with a spatial interpolation model that could not provide highly accurate estimates, while the WBA is an indicator of PM amid other air pollutants.  The higher incidence of COVID-19 weeks after the occurrence of the wildfires shows a lagged association between the wildfire smoke exposure and the higher risk of SARS-Cov-2 infection. Previous research has identified that multiple agents released in wildfires are linked to altered biological mechanisms with potential effect on later viral immune responses (Rager et al., 2021;Velazquez-Salinas et al., 2019). For example, higher levels of PM are significantly associated with increased blood levels of interleukin-6 (IL-6) in people exposed to ambient air pollution over short and long-term (Zhu et al., 2021). This cytokine is produced by immune system cells such as macrophages, and B and T cells to start a cascade of events promoting the transcription of genes linked to cellular signalling processes in the inflammatory response to viral infections (Velazquez-Salinas et al., 2019). In-vitro studies have identified that PM also induces the production of other pro-inflammatory cytokines such as IL-8 and Tumour Necrosis Factor alpha involved in signalling events leading to necrosis or apoptosis and the initiation of chemotaxis and phagocytosis (Shahbaz et al., 2021). If the exposure to wildfires agents in NSW was associated with biological mechanisms linked to the cytokine cascade and altered blood interleukin levels, the immune response of people living in areas affected by the wildfires could have been compromised (Wu et al., 2015).
Multiple epidemiological studies have shown that these immune responses are dose-related, with higher levels of pro-inflammatory cytokines present after PM exposures ranging from 1 day to several years (Zhu et al., 2021); and animal studies have identified increasing levels of pro-inflammatory cytokines and cumulative lung damage after 15 to 90 days of exposure to PM and carbon black particles Ma et al., 2017). The severe and long-term exposure of the NSW population to airborne pollutants in the 2019-2020 wildfires could have played a role in the later spread of COVID-19 infections, linked to immune mechanisms. Other chemicals occurring in wildfires such as polycyclic aromatic hydrocarbons (PAH), methoxyphenol and alkanes are also associated with abnormal immunological responses. Higher PHA levels have been associated with incongruent immune response against viruses such as measles and hepatitis B in studies of vaccine induced immunity (Andrews et al., 2021;Nilamsari et al., 2020). Some species of methoxyphenol and alkanes have been found to induce immunosuppression and enhance pathogenicity of some gram-positive bacteria and to mediate oxidative stress during influenza A virus infection in insects and in-vitro studies respectively (Mollah et al., 2021;Schivo et al., 2014). The above studies support the possibility of an increased risk of viral infections after lagged exposures to airborne particles and chemicals occurring in wildfires via impaired immune responses. However other potential factors mediating the SARS-Cov-2 infection such as social interaction and or mobility between at risk population groups, comorbidities and individual susceptibility can play a determinant role in a COVID-19 outbreak. Although the association between wildfires smoke and COVID-19 incidence found in this study might be partially explained by alterations in the immune response of infected individuals, further research is required to assess individual biomarkers associated with wildfires air pollutants to determine their causal links with COVID-19 infection.
Our analysis also estimated the LGA-specific effects, a valuable addition to the fixed effects that can be used to identify spatial patterns of risk of COVID-19 across NSW. A higher risk was found in LGAs in north regions, beside LGAs with larger proportional area burned and highly populated in the greater Sydney and Hunter regions. This represents the variability of the effect of the wildfire exposures between LGAs, after considering the LGA spatial structure and adjustment for sociodemographic factors. Since the whole posterior distribution can be used to estimate the uncertainty associated with the LGA-specific regression coefficients, we identified the LGAs in these regions with the highest probability of excess risk (relative risk >1). These estimates can be used as proxies of higher risk of COVID-19 incidence per LGA as done in previous research to rank geographical areas and set disease mapping decision rules to prioritize screening and detection (Courtemanche et al., 2015). Although less conservative rules can be set (Richardson et al., 2004), we highlighted LGAs with an excess risk over 98% that match LGAs with a posterior mean greater than 2.5. Another important finding was the association of population density with the COVID-19 incidence. This concurs with previous research that identifies population density as a relevant factor to be included in analyses to determine the risk of COVID-19 infection (Bhadra et al., 2020;Coskun et al., 2021;Kadi and Khelfaoui, 2020;Rashed et al., 2020). These studies found a significant association between population density and SARS-Cov-2 infection although did not adjust for sociodemographic covariates. Our models were adjusted for important factors such as relative socioeconomic disadvantage, sex, and younger and older population groups (<15 and >65 years) to produce an estimate of the COVID-19 incidence associated with WBA. Given the predicted increase in severity and frequency of wildfires in Australia and globally (Beggs et al., 2019), our findings can be used to inform policy making and support public health strategies aimed at curbing COVID-19 and other respiratory infections associated with wildfire pollution.

Limitations
A limitation of this study was the lack of covariates at the individual level within the data available. The inclusion of these variables would have allowed the identification of comorbidities and individual patterns of exposure that have an important effect on the incidence of COVID-19 infection. An ecological study was therefore used to address the absence of individual level variables, using LGAs as the ecological units. This allowed the estimation of spatial trends using spatial regressions for some of the models. However, using LGA-aggregated data carries a risk of ecological bias (i.e., the potential to miss confounding factors at the group level that may produce spurious associations). To reduce this risk of bias, we adjusted for known sociodemographic factors including socioeconomic indicators and population percentages for sex and relevant age groups. We also implemented a mixed effects model that can reduce the risk of ecological bias (Wakefield, 2009) and included a Bayesian spatial regression that incorporates a spatial component to smooth estimates between geographical areas and increase the robustness of the regression coefficients (Kang et al., 2016).
Nevertheless, our findings provide statistical estimates of association but are not measures of causality between wildfire smoke exposures and COVID-19 infection.

Conclusions
A proportional increase of the area affected by wildfires, used as a proxy of the exposure to wildfire smoke, was associated with a higher COVID-19 incidence in NSW. Local Government Areas in the greater Sydney and Hunter regions were shown to have a higher increase in COVID-19 incidence compared to the whole of NSW. The use of wildfire burned area as an alternative indicator of wildfire smoke exposure can be helpful to identify patterns of higher risk in spatial epidemiological analyses to support decision making in the prevention and or mitigation of future health impacts of wildfires. Further research on individual responses to specific wildfires chemicals is necessary to identify the causal links between SARS-Cov-2 infection and wildfires air pollution.