An external exposome-wide association study of COVID-19 mortality in the United States.

The risk factors for severe COVID-19 beyond older age and certain underlying health conditions are largely unknown. Recent studies suggested that long-term environmental exposures may be important determinants of severe COVID-19. However, very few environmental factors have been studied, often separately, without considering the totality of the external environment (i.e., the external exposome). We conducted an external exposome-wide association study (ExWAS) using the nationwide county-level COVID-19 mortality data in the contiguous US. A total of 337 variables characterizing the external exposome from 8 data sources were integrated, harmonized, and spatiotemporally linked to each county. A two-phase procedure was used: (1) in Phase 1, a random 50:50 split divided the data into a discovery set and a replication set, and associations between COVID-19 mortality and individual factors were examined using mixed-effect negative binomial regression models, with multiple comparisons addressed, and (2) in Phase 2, a multivariable regression model including all variables that are significant from both the discovery and replication sets in Phase 1 was fitted. A total of 13 and 22 variables were significant in the discovery and replication sets in Phase 1, respectively. All the 4 variables that were significant in both sets in Phase 1 remained statistically significant in Phase 2, including two air toxicants (i.e., nitrogen dioxide or NO2, and benzidine), one vacant land measure, and one food environment measure. This is the first external exposome study of COVID-19 mortality. It confirmed some of the previously reported environmental factors associated with COVID-19 mortality, but also generated unexpected predictors that may warrant more focused evaluation.

Introduction

The 2019 novel coronavirus disease (COVID-19) is a now global pandemic with severe health, social, and economic consequences. There are over 200,000 deaths in the United States (US) as of October 2020 (CDC, 2020b); however, little is known about the risk factors of severe COVID-19 beyond older age and a few known comorbidities such as hypertension, diabetes mellitus (DM), cardiovascular diseases (CVD), asthma, and chronic lung diseases (Guan et al., 2020; Wu et al., 2020). In addition, large geographic disparities in COVID-19 mortality exists in the US, which cannot be totally explained by the established risk factors alone (Stephanie et al., 2020). A recent study showed that 29% and 22% of the COVID-19 patients in the US who were hospitalized or admitted to the intensive care unit (ICU), respectively, did not have any of the comorbidity included in the current guideline (Chow et al., 2020). It is thus critical to identify additional risk factors contributing to COVID-19 severity and mortality.

There is a large overlap between the currently known risk factors (i.e., comorbidities such as hypertension, DM, CVD) of severe COVID-19 and the diseases that are affected by long-term exposure to environmental factors. Multiple studies have shown the adverse effects of exposures from the natural (e.g., air pollution (An et al., 2018; Brook et al., 2010; Eze et al., 2015; Guarnieri and Balmes, 2014)), built (e.g., walkability (Creatore et al., 2016; Lovasi et al., 2011; Simons et al., 2018)), and social (e.g., neighborhood deprivation (Gale et al., 2011; Laraia et al., 2012)) environments on hypertension, DM, CVD, asthma, and chronic lung diseases. Recent studies found that long-term environmental exposures such as air pollution (i.e., fine particulate matter with diameters that are 2.5 μm and smaller or PM2.5, and nitrogen dioxide or NO2) are associated with COVID-19 mortality (Coker et al., 2020; Conticini et al., 2020; Fattorini and Regoli, 2020; Konstantinoudis et al., 2021; Liang et al., 2020; Ogen, 2020; Wu et al., 2020b; Yongjian et al., 2020), suggesting that environmental exposures may be important determinants of COVID-19 severity. Further, environmental exposures can induce inflammation (Halonen et al., 2010; Keita et al., 2014; Ostro et al., 2014; Shen et al., 2014), which is believed to be an important contributing factor to rapid clinical deterioration of COVID-19 patients (Huang et al., 2020). It is plausible that environmental exposures may be important determinants of severe COVID-19.

Traditional environmental health studies usually examine environmental factors separately without considering the totality of the external environment (i.e., the external exposome) (Niedzwiecki et al., 2019; Wild, 2012). Such approach is limited to generate timely findings for COVID-19 because of (1) the time-consuming process to conduct separate studies on individual exposures, and (2) the difficulties to control for potential confounding by co-exposures. The external exposome is defined as all exposures from the external environment that an individual experiences across the lifetime (Wild, 2012), which is an ideal framework to identify novel environmental factors associated with COVID-19 mortality as it can systematically and efficiently screen hundreds of environmental exposures (Hu et al., 2020; Niedzwiecki et al., 2019; Wild, 2012; Zheng et al., 2020). In this study, we analyzed the county-level COVID-19 mortality data in the US to examine the association between the external exposome and COVID-19 mortality by integrating data on a wide range of environmental factors. Using the agnostic and hypothesis-free external exposome-wide association study (ExWAS) approach with integration of multi-source environmental big data, we aimed to identify novel external exposome factors associated with COVID-19 mortality.

Materials and methods

COVID-19 mortality

Nationwide county-level COVID-19 mortality data in the US were obtained from Johns Hopkins University, Center for Systems Science and Engineering Coronavirus Resource Center (JHU CSSE) (JHU CSSE, 2020). JHU CSSE aggregates daily county-level COVID-19 data from multiple sources such as the Centers for Disease Control and Prevention (CDC) and state health departments. We obtained the cumulative number of deaths for each county (or county equivalent) in the contiguous US (i.e., 48 adjoining US states plus the District of Columbia) up to and including October 31, 2020, which include both confirmed and probable deaths. In the US, a confirmed death is defined by meeting confirmatory laboratory evidence for COVID-19, and a probable death is defined by one of the following: (1) meeting clinical criteria and epidemiologic evidence with no confirmatory laboratory testing performed for COVID-19, (2) meeting presumptive laboratory evidence and either clinical criteria or epidemiologic evidence, and (3) meeting vital records criteria with no confirmatory laboratory testing performed for COVID-19 (CDC, 2020a). We also obtained the 2019 estimate of county population totals from the US Census, which is the latest data available for county-level population in the US. County-level COVID-19 mortality rates were calculated as the ratio of COVID-19 deaths to county level population size.

The external exposome

Data on a variety of natural, built, and social environment measures were obtained from 8 well-validated sources and spatiotemporally linked to each county to determine its residents' long-term exposures to the external exposome before the COVID-19 pandemic. Table 1 shows a summary of the external exposome data sources. A total of 337 external exposome factors covering 9 categories were included.

Table 1

Summary of external exposome measures.

Category	Data source	Time period	Spatial scale	Temporal scale	Number of variables
PM2.5	Atmospheric Composition Analysis Group, WUSTL	2006–2018	0.01 degree in lon/lat	1-year	1
PM2.5 compositions	Atmospheric Composition Analysis Group, WUSTL	2006–2017	0.01 degree in lon/lat	1-year	7
PM10/O3/NO2/CO/SO2	The Center for Air, Climate, & Energy Solutions	2006–2015	Census block group	1-year	5
Air toxicants	National Air Toxic Assessment, EPA	2005, 2011, 2014	County	1-year	164
Walkability	Walkability Index, EPA	2006–2013	Census block group	Cross-sectional	1
Food environment	Food Environment Atlas	2007–2018	County	Cross-sectional	98
Vacant land	Aggregated USPS Administrative Data on Address Vacancies, HUD	2006–2019	Census tract	3-month	19
Social capital	Census Business Pattern	2006–2018	County	1-year	10
Crime and safety	Uniform Crime Reporting Program, FBI	2006–2016	County	1-year	32

Natural environment

Data on PM2.5 and its compositions (i.e., sulfate or SO4 2−, ammonium or NH4 +, nitrate or NO3 −, organic matter or OM, black carbon or BC, mineral dust or DUST, and sea-salt or SS) were obtained from the University of Washington at St. Louis Atmospheric Composition Analysis Group (ACAG) (Van Donkelaar et al., 2019), which estimated annual PM2.5 and its compositions at a spatial resolution of 0.01 degree in longitude and latitude. The estimates were derived using data from a chemical transport model (GEOS-Chem) and satellite observations of aerosol optical depth statistically fused by geographically-weighted models that have been extensively cross-validated (Van Donkelaar et al., 2019). Area-weighted averages were calculated to aggregate the exposure estimates to the county-level, and temporally averaged PM2.5 levels (2006–2018) and PM2.5 compositions (2006–2017) were generated. We also obtained county-level annual estimates of other criteria air pollutants (i.e., particulate matter with diameters that are 10 μm and smaller or PM10, ozone or O3, carbon monoxide or CO, sulfur dioxide or SO2, and NO2) regulated by the US Environmental Protection Agency (US EPA) from the Center for Air, Climate, & Energy Solutions (CACES), which were derived using data from the US Environmental Protection Agency (EPA) regulatory monitors, land use, and satellite-derived estimates of air pollution with well-validated land use regression models (Kim et al., 2020). In addition to the criteria air pollutants, other air toxicant measures were generated using data from the National Air Toxics Assessment (NATA) developed based on a national emissions inventory of outdoor air toxics sources (Logue et al., 2011). Exposure estimates of 164 air toxicants are available at the county level for 2005, 2011, and 2014. Averages were calculated to indicate each county's long-term exposure to the toxicants.

Built environment

Food environment measures were obtained from the US Department of Agriculture (USDA)’s Food Environment Atlas (USDA, 2019), which included county-level food environment indicators from 2007 to 2018. Among the 112 available indicators, we excluded 13 indicators with missing values >10% and 1 indicator with all the counties having the same value. A total of 98 variables were included. Walkability were assessed using the National Walkability Index developed by the US EPA (US EPA, 2015), which measures walkability on a scale from 1 to 20 for each census block group, with 1 indicating the least walkable block group and 20 indicating the most walkable block group. We spatially aggregated the census block group level measures to the county level based on area-weighted averages.

Social environment

Vacant land measures at the census-tract level in 2006–2019 were obtained from the US Department of Housing and Urban Development (HUD) aggregated US Postal Service (USPS) administrative data (Garvin et al., 2013). A total of 19 measures that are available across all years were included. We aggregated the data to the county-level. In addition, ten social capital measures were constructed using the Census Business Pattern data based on the North American Industry Classification System (NACIS) codes (Rupasingha et al., 2006) at the county level. Furthermore, 32 county-level annual crime measures were obtained from the Uniform Crime Reporting Program in 2006–2016 (Barnett-Ryan, 2007).

Potential confounders

A total of 19 county-level variables and 2 state-level variables were considered as potential confounders (Supplemental Table 1): (1) 11 county-level sociodemographic factors obtained from the 2014–2018 American Community Survey (Mather et al., 2005), including population density, percent of population ≥ 65 years, percent of the population 45–64 years, percent of the population 15–44 years, percent living in poverty, median household income, percent black, percent Hispanic, percent of the adult population with less than a high school education, median house value, percent of owner-occupied housing, (2) 2 county-level health factors obtained from the 2020 County Health Rankings (Remington et al., 2015), including percent current smokers and percent obese, (3) 4 county-level meteorological factors generated using data from the gridMET via Google Earth Engine (Blankenau et al., 2020), including average daily temperature and relative humidity for summer (June–September) and winter (December–February), (4) average number of hospital beds per unit population for each county since the first COVID-19 case reported, using daily hospital bed capacity data from the Definitive Healthcare via Amazon Web Services (AWS) Marketplace (Wang et al., 2018), (5) number of days since first COVID-19 case reported for each county obtained from JHU CSSE (JHU CSSE, 2020), and (6) days since state COVID-19 emergency orders and state reopening, obtained from the Boston University COVID-19 United States state policy database (Boston University, 2020).

Statistical analysis

Normalization transformations of all continuous external exposome variables and potential confounders were performed using the bestNormalize package in R (Peterson, 2018), which implements several transformation methods including the log, square root, exponential, arcsinh, Box Cox, and Yeo-Johnson transformations. The best transformation was determined based on the Pearson P statistics. Supplemental Table 2 shows the chosen transformations and parameters for all the external exposome variables. All continuous variables were also z-score standardized (mean = 0 and standard deviation = 1). Missing data for all external exposome factors and potential confounders were imputed using the chained equations method by the mice package in R. A variable was considered as a predictor in the imputation model if its proportion of nonmissing values among counties with missing values in the variable to be imputed was larger than 40% and they were correlated (i.e. with the absolute correlation value>0.4) with the variable to be imputed or the probability of the variable being missing. We imputed a single dataset given the minimal impacts of the imputation procedure due to the large sample size and small fractions of missing data.

In this external ExWAS, we used the standard two-phase EWAS-MLR approach (environment-wide association study followed by a multivariable regression step including the identified hits) (2016; Patel et al., 2013). In Phase 1, we randomly split the data into a 50% discovery set and a 50% replication set. We considered all the 337 external exposome variables for associations with COVID-19 deaths after accounting for multiple comparisons. Mixed effect negative binomial regression models were fitted for each external exposome factor after adjusting for all the potential confounders with population set as the offset and a random intercept by state. The probability function is given bywhere y is the number of COVID-19 death for the ith county in state j, μ is the mean parameter, and α is the over dispersion parameter. The mean parameter is then modelled as:where X and C are the external exposome variable and confounders for the ith county in state j, λ is a random intercept for state, and ε is the population size for the ith county in state j, which is included as an offset. To account for the multiple testing, the Benjamin-Hochberg procedure was used to control the false discovery rate (FDR) at 5% (Benjamini and Hochberg, 1995). A variable is regarded as significant if it has an FDR-adjusted p-value (or q-value) < 0.05 in both the discovery and replication sets. We also generated a correlation heatmap showing the pairwise Pearson correlations of the variables retained from Phase 1. In Phase 2, we used a multivariable mixed-effect negative binomial regression model including all significant external exposome variables from Phase 1as well as all the potential confounders to estimate the effect sizes. Variables remained significant in Phase 2 are retained. The mortality rate ratios (MRRs) and 95% confidence intervals (CIs) were reported. To examine the potential spatial confounding for counties within the same state, we further obtained the residuals and calculated Moran's I for each state (except for states with ≤5 counties/county-equivalents, i.e., Delaware, District of Columbia, Hawaii, and Rhode Island). To evaluate how the aggregated COVID-19 mortality to different transmission stages across counties could affect the results, we performed sensitivity analyses using COVID-19 mortality data reported by September 30, 2020. Fig. 1 shows the flow chart summarizing the external ExWAS methods employed in this study. All analyses were performed using the R statistical software (version 3.6; R Development Core Team). The study has been approved by the Institutional Review Board at University of Florida (IRB202001156).

Fig. 1

Flowchart.

Results

All the 3108 counties and county equivalents in the contiguous US were included in this study, of which 2795 (89.9%) have reported at least one COVID-19 death by October 31, 2020. Supplemental Fig. 1 shows the county-level COVID-19 mortality, and Table 2 shows the COVID-19 mortality as well as county- and state-level characteristics of the included counties.

Supplemental Fig. 1

County-level COVID-19 mortality (till October 31, 2020).

Table 2

Characteristics of 3108 counties and county equivalents in the contiguous US.

	Mean (SD)
COVID-19 deaths (till October 31, 2020)
Number of COVID-19 deaths	73.4 (341.7)
COVID-19 mortality (per 100,000 population)	57.6 (61.6)
County-level characteristics
Population density (per km2)	106.3 (697.7)
Percent of the population ≥ 65 years	18.4 (4.5)
Percent of the population 45–64 years	27.1 (2.9)
Percent of the population 15–44 years	36.1 (5.2)
Percent living in poverty	15.6 (6.5)
Median household income (US dollar)	38,303.0 (116,188.9)
Percent Black	10.0 (14.8)
Percent Hispanic	9.3 (13.9)
Percent of the adult population with less than a high school education	13.7 (6.1)
Median house value (US dollar)	146,109.1 (89,065.4)
Percent of owner-occupied housing	71.5 (8.1)
Percent of current smokers	17.5 (3.6)
Percent of obese	32.9 (5.4)
Average daily temperature for summer (K)	303.1 (3.2)
Average daily temperature for winter (K)	280.4 (6.6)
Average daily relative humidity for summer (%)	89.0 (9.7)
Average daily relative humidity for winter (%)	87.5 (4.8)
Average number of hospital beds since the first reported COVID-19 case (per 100,000 population)	234.3 (333.2)
Number of days since first reported COVID-19 case	156.6 (64.9)
State-level characteristics
Number of days since state COVID-19 emergency orders	233.3 (3.8)
Number of days since state COVID-19 reopening	176.2 (10.9)

Fig. 2 shows the volcano plot summarizing the results from Phase 1. After accounting for multiple comparisons using the Benjamin-Hochberg procedure, a total of 13 and 22 variables were significantly associated with COVID-19 mortality rate in the discovery and replication sets, respectively. Among them, 4 variables were significant in both the discovery and replication sets. Supplemental Table 2 shows the MRRs, 95% CIs, p-values, and q-values for each of the 337 external exposome variables from Phase 1.

Fig. 2

Volcano plot showing the results from Phase 1 of the external ExWAS of county-level COVID-19 mortality in the contiguous US.

In Phase 2, all the 4 significant variables that were significant in both the discovery and replication sets from Phase 1 were simultaneously included in a multivariable mixed-effect negative binomial regression model after adjusting for the potential confounders, and all the 4 variables remained statistically significant. Supplemental Fig. 2 shows the spatial distribution of the 4 variables. Fig. 3 shows the pairwise correlations of the 4 variables retained from Phase 1. All correlation coefficients have absolute values below 0.3.

Supplemental Fig. 2

Spatial distributions of significant external exposome variables from Phase 1 of the external ExWAS.

Fig. 3

Correlation heatmap of significant external exposome variables from Phase 1 of the external ExWAS.

Table 3 shows the adjusted MRRs and 95% CIs (for each standard deviation increase in continuous variables) for these variables. Two of the 4 variables were air toxicants, including exposure to NO2 (MRR: 1.19, 95% CI: 1.13, 1.26), and benzidine (MRR: 0.92, 95% CI: 0.88, 0.95). One was a food environment measure: county-level percent of students eligible for reduced-price lunch in 2015 (MRR: 0.90, 95% CI: 0.87, 0.93). The other one was a vacant land measure: percent of addresses with no-stat in the previous quarter but currently in service (MRR: 0.89, 95% CI: 0.85, 0.92). Supplemental Table 3 shows the Moran's I statistic of residuals for each state from Phase 2, and none of the Moran's I statistic is statistically significant.

Table 3

Results from the external ExWAS of county-level COVID-19 mortality (till October 31, 2020) in the contiguous US.

Exposure		Transformation	Standard deviation	Phase 1						Phase 2
				Discovery set			Replication set
Variable	Category			MRR(95% CI)	p-Value	q-Value	MRR(95% CI)	p-Value	q-Value	MRR(95% CI)	p-Value
Benzidine (ug/m3)	Air toxicant	Box Cox(λ = −0.74)	4.11 × 105	0.90(0.85, 0.95)	2.70 × 10−4	0.011	0.88(0.83, 0.94)	6.17 × 10−5	0.003	0.92(0.88, 0.95)	3.23 × 10−5
NO2 (ppb)	Criteria air pollutant	Yeo-Johnson(λ = −0.49)	1.25 × 10−1	1.22(1.13, 1.32)	3.98 × 10−7	<0.001	1.25(1.15, 1.35)	2.38 × 10−8	<0.001	1.19(1.13, 1.26)	1.41 × 10−10
Percent previous quarter no-stat currently in service	Vacant land	log10(x + 0.001)	2.26 × 10−1	0.86(0.81, 0.91)	2.90 × 10−7	<0.001	0.85(0.80, 0.90)	2.56 × 10−8	<0.001	0.89(0.85, 0.92)	1.04 × 10−8
Percent students eligible for reduced-price lunch, 2015	Food environment	No Transformation	3.71	0.88(0.84, 0.93)	5.29 × 10−6	<0.001	0.90(0.85, 0.95)	5.14 × 10−5	0.003	0.90(0.87, 0.93)	6.22 × 10−9

Mortality rate ratio (MRR) and 95% confidence interval (CI) for each standard deviation increase.

Supplemental Table 4, Supplemental Table 5 shows the results from the sensitivity analyses using COVID-19 mortality data reported by September 30, 2020. The 4 significant variables found in the main analyses were also identified from the sensitivity analyses with similar effect sizes. In addition, the sensitivity analyses additionally identified three significant variables, including chloroform (MRR: 1.05, 95% CIL 1.01, 1.10), methyl bromide (MRR: 1.07, 95% CI: 1.02, 1.12), and antimony compounds (MRR: 1.05, 95% CI: 1.01, 1.08).

Discussions

This is the first external ExWAS to estimate the associations between long-term exposure to the external exposome and COVID-19 mortality. Using the standard two-phase EWAS-MLR approach, we assessed the association of 337 external exposome variables with COVID-19 mortality (till October 31, 2020) in the contiguous US. After accounting for multiple testing and high correlations among the exposures, 4 county-level external exposome variables characterizing the natural (i.e., criteria air pollutants and air toxicants), built (i.e., food environment), and social environment (i.e., vacant land) were identified to be significantly associated with county-level COVID-19 mortality.

Recent studies in the US and the Europe have suggested that long-term exposures to criteria air pollutants regulated by the US EPA such as PM2.5, NO2, SO2, and O3 are associated with higher COVID-19 mortality (Coker et al., 2020; Cole et al., 2020; Conticini et al., 2020; Fattorini and Regoli, 2020; Hendryx and Luo, 2020; Konstantinoudis et al., 2021; Liang et al., 2020; Ogen, 2020; Wu et al., 2020b; Yongjian et al., 2020). Long-term air pollution exposures can induce inflammation (Chuang et al., 2007), which may lead to severe symptoms among COVID-19 patients. In Phase 1 of our study, O3 was significant in the discovery set only, and PM10 and SS were significant in the replication set only. Only NO2 remained significantly associated with county-level COVID-19 mortality after addressing multiple testing and account for the potential confounding due to co-exposures to other external exposome factors. This finding is consistent with several other studies conducted in the US and Europe (Conticini et al., 2020; Liang et al., 2020; Ogen, 2020). To our knowledge, no study has examined air toxicants beyond the criteria air pollutants (i.e., carbon monoxide, lead, ground-level ozone, nitrogen dioxide, particulate matter, and sulfur dioxide). Using the county-level population exposure estimates from NATA, we found that benzidine was associated with county-level COVID-19 mortality. Benzidine is widely used for the detection of blood and as a reagent in the production of dyes (EPA, 2016a). In addition, the sensitivity analyses using COVID-19 mortality data reported by September 30, 2020 additionally found methyl bromide, chloroform, and antimony compounds to be associated with COVID-19 mortality, although chloroform and antimony compounds were not significant in either the discovery or replication sets and methyl bromide was significant only in the discovery set in Phase 1 of the main analyses. Methyl bromide is extensively used as a pesticide (Committee on Acute Exposure Guideline Levels - National Research Council, 2012), and antimony is primarily used in grid metal for lead acid storage batteries (Winship, 1987). Both of them have been known to cause inflammation and fibrosis (Boorman et al., 1986; Committee on Acute Exposure Guideline Levels - National Research Council, 2012; Winship, 1987), which are believed to be important contributing factors to rapid clinical deterioration of COVID-19 patients (Huang et al., 2020). Chloroform is primarily used in industry as a solvent and can also be present in chlorinated water as a byproduct of water chlorination (McCulloch, 2003). While long-term exposure to chloroform and benzidine is known to cause adverse effects on the liver and central nervous system (EPA, 2016a; EPA, 2016b), it is largely unknown how they may contribute to COVID-19 mortality, and future studies are needed to confirm the associations and understand the underlying mechanisms.

To date, no study has examined whether long-term exposures to the built and social environments are associated with COVID-19 mortality. In this study, we found that increased long-term access to reduced-price lunch was protective of COVID-19 mortality, suggesting that long-term food security may play an important role in COVID-19 outcomes. While the underlying mechanisms are largely unknown, previous studies showed that food environment is a major contributing factor for obesity and many other health outcomes (Dubowitz et al., 2012; Morland and Evenson, 2009). In addition, we also found that a higher percent of addresses currently in service but were classified as “No-Stat” (i.e., vacant, under construction, or inactive (USPS, 2020)) in the previous quarter is associated with lower COVID-19 mortality, suggesting that living in counties with decreasing vacant lands are protective of COVID-19 mortality. Previous studies have shown that vacant land may lead to negative health outcomes due to their association with higher levels of chronic stress and less social interactions (Garvin et al., 2013).

Our study has several strengths. By using the agnostic external ExWAS approach, we examined the associations between COVID-19 mortality and long-term exposures to a variety of external exposome factors, which addressed limitations of previous studies that assessed only a small fraction of these factors separately (Coker et al., 2020; Cole et al., 2020; Hendryx and Luo, 2020; Konstantinoudis et al., 2021; Liang et al., 2020; Wu et al., 2020b). In addition, we spatiotemporally linked data characterizing the natural, built, and social environments of each county in the contiguous US to account for the spatiotemporal dynamic nature of external exposome factors. Furthermore, a wide range of potential confounders have been adjusted.

Several limitations also need to be acknowledged. First, this is an ecological study with COVID-19 mortality data only available at the county level. It is important to note that the results of this study should not be used to make individual-level inferential statements due to the potential ecological bias. However, as risk factors beyond older age and comorbidities for severe COVID-19 are largely unknown, the novel external exposome factors identified from this study may generate new hypotheses to be tested in future investigations when nationwide individual-level COVID-19 data are available. Second, although we obtained county-level COVID-19 death data from JHU CSSE, which is considered as the most comprehensive data source for COVID-19 in the US (JHU CSSE, 2020), the quality of the COVID-19 data may change substantially over time across different counties within the US. However, the consistent findings observed in the sensitivity analyses suggest that the results we observed are robust. Third, the standard two-phase EWAS-MLR approach we used did not consider nonlinear association and potential interactions, and many of the exposures considered may be subject to potential measurement error. In addition, although many environment factors have been included to characterize the external exposome, this list is not exhaustive, and continuing efforts are needed to further improve the measurement of the external exposome.

Conclusions

This external ExWAS of county-level COVID-19 mortality in the contiguous US provides new insights into the role of long-term exposures to the external exposome in COVID-19 mortality. We confirmed a previously reported association (i.e. NO2), and identified novel environmental factors associated with COVID-19 mortality, including air toxicants (i.e., methyl bromide, benzidine, and chloroform), vacant land, and food environment measures. Although this is an ecological study with only county-level data available, the novel external exposome factors identified from this study may inform future investigations with individual-level data to confirm and understand the mechanisms underlying these associations.

The following are the supplementary data related to this article.

Supplemental Table 1

Data sources for potential confounders.

Supplemental Table 2

Results from Phase 1 of the external ExWAS of county-level COVID-19 mortality (till October 31, 2020) in the contiguous US.

Supplemental Table 3

Moran's I statistic of residuals for each state from Phase 2 of the external ExWAS of county-level COVID-19 mortality (till October 31, 2020) in the contiguous US.

Supplemental Table 4

Results from Phase 1 of the sensitivity analyses of the external ExWAS of county-level COVID-19 mortality (till September 30, 2020) in the contiguous US.

Supplemental Table 5

Results from the sensitivity analyses of the external ExWAS of county-level COVID-19 mortality (till September 30, 2020) in the contiguous US.

Funding

This work was supported by the National Institute of Environmental Health Sciences of the under award number 1R21ES032762. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

CRediT authorship contribution statement

Hui Hu: Data curation, Formal analysis, Funding acquisition, Methodology, Supervision, Writing – original draft. Yi Zheng: Data curation, Writing – review & editing. Xiaoxiao Wen: Data curation, Writing – review & editing. Sabrina S. Smith: Data curation, Writing – review & editing. Javlon Nizomov: Writing – review & editing. Jennifer Fishe: Writing – review & editing. William R. Hogan: Writing – review & editing. Elizabeth A. Shenkman: Writing – review & editing. Jiang Bian: Funding acquisition, Writing – review & editing.

Declaration of competing interest

The authors disclose that they have no actual or potential competing interest.