Are environmental pollution and biodiversity levels associated to the spread and mortality of COVID-19? A four-month global analysis.

On March 12th, 2020, the WHO declared COVID-19 as a pandemic. The collective impact of environmental and ecosystem factors, as well as biodiversity, on the spread of COVID-19 and its mortality evolution remain empirically unknown, particularly in regions with a wide ecosystem range. The aim of our study is to assess how those factors impact on the COVID-19 spread and mortality by country. This study compiled a global database merging WHO daily case reports with other publicly available measures from January 21st to May 18th, 2020. We applied spatio-temporal models to identify the influence of biodiversity, temperature, and precipitation and fitted generalized linear mixed models to identify the effects of environmental variables. Additionally, we used count time series to characterize the association between COVID-19 spread and air quality factors. All analyses were adjusted by social demographic, country-income level, and government policy intervention confounders, among 160 countries, globally. Our results reveal a statistically meaningful association between COVID-19 infection and several factors of interest at country and city levels such as the national biodiversity index, air quality, and pollutants elements (PM10, PM2.5, and O3). Particularly, there is a significant relationship of loss of biodiversity, high level of air pollutants, and diminished air quality with COVID-19 infection spread and mortality. Our findings provide an empirical foundation for future studies on the relationship between air quality variables, a country's biodiversity, and COVID-19 transmission and mortality. The relationships measured in this study can be valuable when governments plan environmental and health policies, as alternative strategy to respond to new COVID-19 outbreaks and prevent future crises.

Introduction

Outbreaks of emerging infectious diseases, such as the 1918 influenza pandemic (Taubenberger and Morens, 2006), the 2014 Ebola (Rewar and Mirdha, 2014) virus disease, the white-nose-syndrome in bats (Blehert et al., 2009), the ash dieback (Pautasso et al., 2013) fungal disease in ash trees, and the pandemic chytridiomycosis, which killed amphibians worldwide (Fisher et al., 2009; Fisher et al., 2012), are occurring with an increasing frequency and terrible consequences. One of the main causes of pandemic events and epidemic diseases is the close interaction between human populations and both domesticated and wildlife animals that carry pathogens (Woolhouse and Gowtage-Sequeria, 2005). Most pathogens pass from their wildlife reservoirs onto human populations through hunting, the consumption of wild species, wild animal trade, and other contacts with the wildlife. Additionally, changes in the Earth’s climate and weather continue to impact the planet’s ecosystems, which include the environmental communities with infectious disease agents and hosts that acts as vectors (Yeh K et al., 2020). Therefore, the intensified emergence of infectious pathogens can also be attributed to climate change, biodiversity loss, habitat degradation, and a rate increase of human-wildlife interactions (HWI) (Schmeller et al., 2020).

A reduction of biodiversity richness and evenness causes the disappearance of a key part of the ecosystem that serves as a buffer to the spread of infectious diseases onto humans, animals, and plants (Peixoto and Abramson, 2006; Pongsiri et al., 2009; Ostfeld and Keesing, 2017). In that sense, several studies have suggested that the transmission of diseases increases with the loss of biodiversity (Keesing et al., 2010; Wood et al., 2014; Lacroix et al., 2014; Johnson et al., 2015). Additionally, the fifth edition of UN’s Global Biodiversity Outlook report (GBO-5) published by the Convention of Biological Diversity (CBD), remarks the importance of biodiversity when addressing climate change and long-term food security. GBO-5 concludes that action to protect biodiversity is essential to prevent future pandemics. Moreover, the report notes that biodiversity loss might also lead to a faster rate of emergence and re-emergence of infectious diseases.

There have been six large-scale epidemics in the 21st century (i.e., SARS, swine flu, MERS, Ebola, Zika, and Avian). Here, we briefly present their characteristics:

The Severe Acute Respiratory Syndrome (SARS) occurred in 2003, which led to more than 8000 infections with a mortality rate of approximately 10% and an impact limited only to local and regional economies (LeDuc and Barry, 2004). This epidemic ended abruptly in July 2003 and no human cases of the SARS coronavirus have been detected since then.

The 2009 H1N1 influenza virus, which causes swine flu, was a pandemic that first appeared in Mexico in March 2009 and then in April in the United States. The 2009 swine flu became a pandemic as a result of global mobility and airline travel and led to an estimated 0.4% case fatality (Al Hajjar and McIntosh, 2010).

The Middle East respiratory syndrome (MERS) was first identified in humans in Saudi Arabia and Jordan in 2012 (Memish et al., 2020). MERS is considered a zoonotic pathogen that jumps from infected dromedary camels to humans (El-Kafrawy et al., 2019; Gardner et al., 2019). By contrast to SARS, which was contained within a year after emerging, MERS continues to have a limited circulation in the Middle East region and causes intermittent sporadic human infection cases, infected community clusters, and nosocomial outbreaks, all of which hold a high risk of spreading globally (Zumla et al., 2015).

The Ebola virus was first detected in 1976 in Zaire (presently known as the Democratic Republic of Congo). Since the virus was first detected, over 20 known outbreaks of Ebola have been identified in sub-Saharan Africa, mostly in Sudan, Uganda, Democratic Republic of Congo, and Gabon (Malvy et al., 2019). At present, no vaccine or efficient antiviral management strategy exists for Ebola (Hasan et al., 2019). Although the Ebola virus has a substantial epidemic and pandemic potential due to the ease of international travel, as demonstrated by the 2013–2016 West-African Ebola virus epidemic with approximately 28,000 confirmed cases and 11,000 deaths (Garske et al., 2017).

The Zika fever (2015–2016) was first isolated in 1947 from a febrile rhesus macaque monkey in the Zika Forest of Uganda. Since 1954, when the first cases in humans were reported, the Zika virus caused only limited sporadic infections in Africa and Asia. However, a large outbreak with approximately 440,000 to 1,300,000 cases spread from Brazil to 29 countries in the Americas in 2015 (Plourde and Bloch, 2016). In November 2016, WHO announced the end of the Zika outbreak.

Avian flu (or bird flu) was first reported in 1997 in Hong Kong with only 18 infections and 6 human deaths. However, more than 700 cases of the avian flu have been reported from over 60 countries (Alexander and Brown, 2009), which include the 2016 outbreaks that occurred in China (Chatziprodromidou et al., 2018).

Some of those epidemics have been studied to anticipate their societal and environmental impacts (Crameri et al., 2015; Kim et al., 2018; Qiu et al., 2018). The current novel coronavirus disease (COVID-19) dwarfs those six large-scale epidemics of the 21st century in terms of spatial extent and societal consequences (Ali and Alharbi, 2020), and COVID-19 is the only pandemic with widespread and complex environmental impact (Lal et al., 2020; Nakada and Urban, 2020).

On March 12, 2020, the World Health Organization (WHO) declared the COVID-19 as a pandemic (WHO, 2020). By May 18th, 2020, more than 210 countries reported confirmed cases of COVID-19 (Cucinotta and Vanelli, 2020). Extended virus transmission outside China was reported among various European Union (EU) countries, the United States of America (USA), Latin American countries (e.g. Brazil and Peru) and African countries (e.g. South Africa) (Davis, 2020).

Several social distancing measures were implemented to intervene and contain the alarming spread of COVID-19 (Pan et al., 2020). Various clinical trials to develop a vaccine or a pharmaceutical treatment are being directed to fight against the virus. Some initial data on the effect of environmental factors (i.e. temperature and humidity) on virus spread and mortality have been presented (Prata et al., 2020; Tobías and Molina, 2020; Triplett, 2020; Wang et al., 2020a; Wang et al., 2020b). However, none of these studies included as confounders factors such as the mutual impact of both social distancing and government movement restriction policies on virus spread. Additionally, while various researchers pointed out the role of biodiversity on COVID-19 spread (Corlett et al., 2020; Grandcolas and Justine, 2020; Lorentzen et al., 2020; Outlook et al., 2010), until now, there is no information on that role at the local or international level. Although still a debate in ecological forums, theoretically, a higher biodiversity acts as an increased protection factor that enhances the human immune system against unknown viruses (Maas et al., 2006; Rook, 2013).

The collective impact of environmental and ecosystem factors, as well as biodiversity, on COVID-19 spread remain empirically unknown, particularly in regions with a wide ecosystem range. To the best of our knowledge, there is no information on the relationship of biodiversity with COVID-19 at a worldwide level. Estimates, over time, of COVID-19 spread and mortality that consider ecosystem and biodiversity determinants could help identify which level of these factors can be beneficial to slow down spread and mortality, and with how much impact. In our view, an estimate of the effect of environmental and biodiversity parameters, jointly along with other studied factors, can give insights that may help guide authorities when establishing an early decision for the containment of the future outbreaks. Thus, the aim of this study was to assess the relationship between biodiversity, environmental, and other ecosystem factors with COVID-19 spread and mortality at global levels.

Methods

Study design

We conducted a retrospective, observational, longitudinal study. We obtained data on COVID-19 spread and mortality, and related risk factors from 218 countries. We compiled a dataset of COVID-19 daily cases and deaths spanning January 21st to May 18th, 2020, based on the most recent publicly available population-level information (per country), as reported by WHO (https://www.who.int/emergencies/diseases/novel-coronavirus-2019/situation-reports/). The current study was approved by Parc Sanitari’s Sant Joan de Déu, Ethics Committee (PIC-67-20, Barcelona, Spain) and complies with the ethical guidelines of the 1975 Declaration of Helsinki.

COVID-19 international data and other baseline measures

The WHO daily situation reports were used to compile data between January 21st to May 18th, 2020, on daily confirmed cases, total confirmed cases, daily confirmed deaths, the total amount of confirmed deaths, and time since the last reported case for each of the 218 countries/regions. Cases identified in cruise ships were excluded from the analysis. Cases among all China’s provinces were grouped all together. COVID-19 cases were classified separately, in particular, administrative regions of China such as Hong Kong, Macao, and Taiwan since they applied different government interventions and policy measures than mainland China. Based on the WHO database, Puerto Rico, Northern Mariana Islands, Guam, and United States Virgin Islands were classified separately from the US.

The effective date of each social distancing intervention, per country, was initially extracted from online policy databases (World Health Organization, 2020) to create a 4-level government policy intervention score that varied from levels 0 to 3, which represented “low”, “intermediate”, “high”, and “very high” intervention levels. Level 0 (“low”) is defined as those countries with no restrictions at entry points or temperature check or additional medical screening (at entry) or completion of travel health questionnaires at entry points or quarantine in suspected cases imported from affected areas in each time; Level 1 (“intermediate”) is defined as those countries that announced the “low” measures plus visa suspension or suspension of entry to specific passengers or flow suspension to and from COVID-19 affected areas; Level 2 (“high”) is defined as those countries that applied all the aforementioned and proceeded to annulations and suspensions of large gatherings, events or educational activities (schools and/or universities closures) or isolation of specific areas; and Level 3 (“very high”) is defined as those countries that applied all the aforementioned and proceeded to quarantine the entire country (defines as full or mandatory lockdown) or announced a stay-at-home-order or applied nationwide curfew (defined as the specific government order referred to the period of time when the population was required stay at home or was allowed to move outside their homes - night or day curfews were considered). More detail on the abovementioned methodology could be found in Tyrolovas et al. (2020).

The World Bank classification system was used to classify each of the 218 countries in a distinct country income level: High (HICs), Upper (UMICs), Lower Middle (LMICs), and Low (LICs) income. Also, to define different geographical regions, the WHO classification was used as follows: European Region, Western Pacific Region; Region of the Americas; African Region, Eastern Mediterranean Region, and South-East Asia Region. Population density per square kilometer was also assessed by country from the World Bank data (Bank, 2016).

Finally, we used two more variables: “days since the first case” and “days since the last case”. The former variable measures the number of days since the first COVID-19 case was reported in each country. This variable allows us to compare all countries even if they have different starting times of the disease. Thus, it describes the disease arrival into a region. Additionally, “days since the last case” is a variable created by WHO and measures the count of days in which a country has not reported a new COVID-19 case. Therefore, it describes the acceleration or deceleration of the disease spread in a region.

Biodiversity, environmental, and ecosystem assessments by country

The National Biodiversity Index (NBI) reported by the Convention on Biological Diversity (https://www.cbd.int/gbo1/annex.shtml) is a measure of the variation of genetic, species, and ecosystem levels for each country and is based on estimates of a country’s richness and endemism of four terrestrial vertebrate classes and vascular plants; vertebrates and plants are ranked equally. The NBI includes an adjustment for country size and values range 0–1. Countries with a land area less than 5000 sq km, overseas territories, and dependencies were excluded.

Daily temperature and precipitation measures were downloaded from NASA’s Goddard Earth Sciences Data and Information Services Center (GES DISC). Specifically, we obtained data from the algorithm called Integrated Multi-satellite Retrievals for Global Precipitation Measurement (https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGDE_06/summary) (IMERG) (Nasa ED, 2020a). The measurement of precipitationCa measures the precipitation of the combined microwave-IR spectrum. Additionally, the maximum, minimum, and average, daily temperatures from 2m above the ground, for each country, were obtained from the MERRA-2 (a Modern-Era Retrospective analysis for Research and Applications version 2) (Nasa ED, 2020b) (https://disc.gsfc.nasa.gov/datasets/M2SDNXSLV_5.12.4/summary). We used only the maximum temperature in our models to avoid multicolinearity. We note that we obtained the 14-day lagged date of the temperature and precipitation measures to account for 14 days between those measures and case and mortality confirmation. This lag approach has been similarly applied before in other studies (Shi et al., 2020; Runkle et al., 2020).

The factors associated with the ecosystem vitality and environmental health, per country, are listed in Table 1 and were obtained from the 2020 Environmental Performance Index (EPI) report (Index EP, 2018) (https://epi.yale.edu/downloads). A data-driven summary of the state of sustainability around the world is provided by EPI, which uses 32 performance factors across 11 issue categories. The data comes from trusted third-party sources, such as international governing bodies, nongovernmental organizations, and academic research centers. Credible datasets rely on established collection methods that have been peer-reviewed by the scientific community or endorsed by international authorities.

Table 1

The list below shows the environmental and ecosystem vitality factors from the 2020 Environmental Performance Index Framework. The framework organizes 32 factors into 11 issue categories and two policy objectives. The code designates each factor variable. Source. (Index EP, 2018).

Policy objective	Issue category	Factor	Code
Environmental Health	Air Quality	Ambient particulate matter pollution	PMD
		Household air pollution from solid fuels	HAD
		Ozone	OZD
	Sanitation & Drinking Water	Unsafe drinking water	UWD
		Unsafe sanitation	USD
	Heavy Metals	Lead Exposure	PBD
	Waste Management	Solid Waste	MSW
Ecosystem Vitality	Biodiversity & Habitat	Terrestrial Biome Protection – National weights	TBN
		Terrestrial Biome Protection – Global weights	TBG
		Marine protection	MPA
		Protected Areas Representativeness Index	PAR
		Species Habitat Index	SHI
		Species Protection Index	SPI
		Biodiversity habitat Index – Vascular Plants	BHV
	Ecosystem Services	Tree cover loss, %	TLC
		Grassland Loss	GRL
		Wetland Loss	WTL
	Fisheries	Fish Stock Status	FSS
		Regional Marine Trophic Index	RMS
		Fish caught by Trawling	FGT
	Climate Change	CO2 intensity trend	CDA
		Methane intensity trend	CHA
		F-gases intensity trend	FGA
		N2O intensity trend	NDA
		Black Carbon intensity trend	BCA
		GHG emission intensity growth rate	GIB
		GHG emission per capita	GHP
		CO2 from Land Cover, trend	LCB
	Pollution Emissions	SO2 intensity trend	SDA
		NOX intensity trend	NXA
	Agriculture	Sustainable Nitrogen Management Index	SNM
	Water Resources	Wastewater treatment level	WWT

This study emphasizes air deficiency 2020 EPIs, which are broken down into five factors. Household solid fuels (HAD) measures the actual outcomes from exposure to indoor air pollution from household use of solid fuels. Ambient particulate matter pollution (PMD) measures the average annual concentration of PM2.5 to which the typical citizen of each country is exposed. Ozone (OZD) measures the intensity of ground-level ozone to which the typical citizen of each country is exposed. The units of the former three measures are the number of age-standardized disability-adjusted life-years lost per 100,000 persons. And SO2 Emissions (SDA) and NOx Emissions (NXA) measure the intensity of SO2 and NOx emissions respectively, from the entire economy, as a blend of current-year intensity and a 10-year trend. Technical details of how these measures have been calculated can be found in the EPI technical appendix (https://epi.yale.edu/downloads/epi2020technicalappendix20200604.pdf). For all variables from the EPI report, we have used the 2020 EPI values and also the 10-year change rate to include increments or decrements of the measures. The 2020 EPI report calculates the air quality and pollutant measures from the previous year.

Atmospheric measurements by cities

We obtained information on the daily average of temperature, level of humidity, ground-level ozone, atmospheric particulate matter of 10 μm or less in diameter (PM10), and 2.5 μm or less in diameter (PM2.5) from the World Air Quality Index project (https://aqicn.org/data-platform/COVID-19). These measures are city-based, and we used three cities with different levels of COVID-19 spread: Denver (with a medium case rate), and Barcelona and Milan (with high case rates). We used data that spanned from March 8th until May 18th, 2020.

Statistical analysis

Bayesian space-time analysis

We aimed to identify the influence of biodiversity, temperature, and precipitation factors related to the number of positive COVID-19 cases and deaths. As we utilize data collected across different countries in time, we employ a Bayesian spatio-temporal approach to capture those effects between the outcomes and the covariates. In particular, this approach accounts for temporal correlations as well as spatio-temporal interactions, which have proven to be important due to the nature of the spread and mortality cases of COVID-19. Additionally, we assumed a Zero-Inflated Negative Binomial (ZINB) distribution in the outcomes of interest that exhibited overdispersion and excess zeros to lead to statistically-sound inferences.

The spatio-temporal models for assessing the evolution of COVID-19 spread and mortality were applied with the adjustment of the following confounders: government intervention level (4 stringency levels), the number of days since the last COVID-19 new case (days since last case), precipitation and temperature measurements, country income level, count of days since the first COVID-19 case is reported in each country (days since first case), and population density (per sq. km). The natural logarithm of the total population was added to the linear predictor function (as an offset) to account for the infection and mortality rate per country, as is more relevant to model spread or mortality rates (scaled based on population) than counts of cases. We selected the best-fitting model based on both the Deviance Information Criterion (DIC) (Spiegelhalter et al., 2002) and the Watanabe-Akaike Information Criterion (WAIC) (Watanabe, 2010). All model details are presented in the Supplementary Appendix S1. The Bayesian spatio-temporal models were used as a first step in the current analysis to capture spatial and in-time effects. All computations were carried out using the R package R-INLA (Rue et al., 2009) in R Version 4.0.2.

Generalized linear mixed models

To identify the effects of ecosystem vitality and environmental health variables related to COVID-19 spread and mortality, we also fitted generalized linear mixed models with Template Model Builder (glmmTMB) and assumed a ZINB (ZINB) distribution in the outcomes of interest. The analysis period was from January 21st to May 18th.

The factors used in this analysis were those listed in Table 1. As a preliminary analysis, given that all variable types are numeric, we checked correlations among them to avoid multicollinearity issues (see Supplementary Figure S1 for the correlation matrix with the final set of variables used). The glmmTMB models were applied with the adjustment of the following confounders: government intervention level (low, intermediate, high, and very high), country income level, count of days since the first COVID-19 case in each country, the World Bank geographical region categorization, and the interaction between the last two confounders. We selected the best-fitting model based on the Akaike Information Criterion (AIC) (Akaike, 1974). Additionally, because our data comes from a longitudinal study with measurements over time per country, we included random effects for the country grouping variable. However, we also assumed that the correlations within a country over time are not constant. Therefore, we included random intercepts and random slopes model which implies that correlations between the observations within a country are functions of time. The confounders stated above were incorporated as fixed effects of the model. As we applied in the Bayesian space-time analysis, the population offset was taken into account in the models. Similar methodological approaches have been applied by previous studies in the field (Wu et al., 2020; Travaglio et al., 2020). All model details are presented in the Supplementary Appendix S2.

We fitted the same models without the assumption of random slopes with the aim of fitting parsimonious models. We calculated the likelihood ratio test to compare the fit of the two models (with and without random slopes). The results show evidence that a model including random slopes is significantly better than the simple models (p-value < 0.001). Finally, we also fitted similar regression models stratified by country-income level to control the virus’ spread and mortality variability of surveillance infrastructures and monitoring systems among regions. Thus, we created two groups: the Low and Lower-middle income countries (94 countries – 58.7%) and the High and Upper-middle income countries (66 countries – 41.3%), in which each income level group reflects healthcare quality in a homogenous way. All modelling was carried out using the R package glmmTMB (Brooks et al., 2017) in R Version 4.0.2.

Analysis of count time series following generalized linear models

We fitted a count time series and followed Negative Binomial models with an autoregressive term to characterize the association between COVID-19 spread and air quality factors. Particularly, we were interested in finding if the following air quality factors: atmospheric particulate matter 10 μm or less in diameter (PM10), 2.5 μm or less in diameter (PM2.5), and ground-level ozone (O3), have a significant effect related to the number of COVID-19 cases. We assessed all count time series models adjusted by humidity and temperature confounders and the analysis was carried out using the R package tscount (Liboschik et al., 2015) in R Version 3.6.3. Besides, we calculated Spearman correlations between the COVID-19 spread and all air quality factors. We note that the date of the air quality factors was that of two weeks before respective COVID-19 case dates to account for 14 days between transmission and case confirmation.

Results

We obtained data on the social demographic, country-income level, and government policy intervention factor from all 208 countries. The presence of missing values reduced our sample of countries to 192 (90.8%). Out of these 192 countries, 174 (90.6%) and 180 (93.8%) had the National Biodiversity Index and the EPI report informed, respectively. Due to differential patterns of missing data, our final data pertained to 160 countries with all variables relevant to the analyses reported here. The distribution of country-income level is quite similar among all groups with exception of lower income countries (HICs: 29.4%; UMICs: 28.8%, LMICs: 25.6%; LICs: 16.3%). The region distribution is in descending order: Europe & Central Asia (29.5%), Sub-Saharan Africa (27.5%), Latin America & Caribbean (16.4%), East Asia & Pacific (11.9%), Middle East & North Africa (9.4%), South Asia (4.4%), and North America (1.3%).

Association of COVID-19 spread and mortality with national biodiversity index; a global and spatio-temporal analysis

The raw data for COVID-19 spread (case rates by 100,000) and NBI are displayed in the world maps in Fig. 1 . The color scale was set to range between 0 and 1 in both maps for the COVID-19 case rates and NBI to allow comparison. We see that color distribution between the maps are reversed, i.e. large values in one map (more reddish tones) correspond to smaller values in the other map (more bluish tones) and vice versa. These maps are visual indication of a negative correlation between NBI and COVID-19 case rate.

Fig. 1

World map of the NBI values (map on the left) and the COVID-19 spread (map on the right) from January 21st to May 18th.

The summary table given in Fig. 1 confirms this univariate association. The table shows the average of COVID-19 case rates for different groups of NBI (categorized using quartiles). For instance, the COVID-19 case rate corresponding to the lower values of NBI is 1.56, which drastically increased from a case rate of 0.134 for high values of NBI.

The spatio-temporal regression analysis, which is summarized in Table 2 , assesses the COVID-19 spread along with government interventions, level of income, biodiversity, environmental, and other factors, among 160 countries.

Table 2

Bayesian spatio-temporal regression analysis to evaluate the COVID-19 spread.

Items	Estimated coefficient	95% HPDI
NBI	−0.606	−0.946, -0.268
precipitationCa	−0.001	−0.003, 0.001
Temperature (max.)	−0.010	−0.017, -0.002
Population Density (sq/km)	0.0003	−0.0002, 0.0006
Days since last case	−0.015	−0.041, 0.010
Days since first case	−0.009	−0.016, -0.002
HICs	Reference Category
LICs	−0.334	−0.801, 0.132
LMICs	−0.275	−0.636, 0.086
UMICs	−0.305	−0.612, 0.000
NBI: Temperature (max.)	−0.011	−0.023, 0.002

Significant effects where the 95% HPDI does not include the zero-value are shown in boldface. HPDI: Highest Posterior Density Interval, is the equivalent CI in a Bayesian framework. LMICs: Lower Middle-income countries; UMICs: Upper Middle-income countries; LICs: Low-income countries. NBI: National Biodiversity Index as reported by the Convention on Biological Diversity. precipitationCa: measures the precipitation of the combined microwave-IR spectrum. Days since last case: the number of days since the last COVID-19 new case. Days since first case: count of days since the first COVID-19 case is reported in each country.

Spatio-temporal models were also adjusted for government policy interventions.

The biodiversity degree of a country was reversely associated with to COVID-19 spread globally [NBI: −0.61, 95%CI (−0.95, −0.27)]. NBI impact was statistically strongest in this model: the lower the country’s variations in genetic, species, and ecosystem levels, the higher was the impact on the levels of COVID-19 spread. Although with a lower significance, there was also an inverse effect of the maximum 2-m air temperature (see Table 2), suggesting the spread of the virus is lower in countries with higher temperatures. There were no significant interaction effects between biodiversity and temperature related to COVID-19 spread. The same observation was made with precipitation measures. However, number of days since the first reported COVID-19 case was inversely related to the number of COVID-19 cases [days since first case: −0.01, 95%CI (−0.02, −0.002)].

We also utilized the Bayesian spatio-temporal, and the generalized linear mixed models to analyze the evolution of mortality (see Table S1 of the Supplementary Appendix S3 for the detailed results). The results showed that biodiversity and environmental factors do not have a direct influence on mortality.

Association of COVID-19 spread and mortality with ecosystem vitality and environmental health variables

Fig. 2 shows the set of factors having a significant effect on the COVID-19 spread, which are estimated using a generalized mixed ZINB model for ecosystem measures as presented in Table 3 .

Fig. 2

The plots show the effect on the spread of COVID-19 when associated to (a) Air quality deficiency, (b) country income group, days since the first COVID-19 case by (c) region and by (d) country-income group. The COVID-19 spread units are based on rates via the offset as produced by the related models.

Table 3

Generalized mixed ZINB model regression analysis to evaluate the COVID-19 spread.

Items	Estimated coef	95% CI
HICs	Reference Category
LICs	−2.080	−3.842, -0.311
LMICs	−1.880	−3.425, -0.332
UMICs	−1.310	−2.515, -0.114
Air Deficiency	0.028	0.004, 0.053
Air Deficiency 10-year change	−0.020	−0.142, 0.101
Sanitation & Drinking Water 10-year change	0.001	−0.140, 0.143
Heavy Metals 10-year change	0.030	−0.116, 0.177
Biodiversity & Habitat	−0.001	−0.019, 0.017
Biodiversity & Habitat 10-year change	0.005	−0.028, 0.038
Ecosystem Services	0.006	−0.009, 0.020
Climate Change 10-year change	0.0002	−0.021, 0.021
Pollution Emissions	−0.014	−0.033, 0.006
Pollution Emissions 10-year change	0.009	−0.004, 0.023
Agriculture	−0.010	−0.031, 0.011
Agriculture 10-year change	−0.004	−0.036, 0.029
Days since first case:East Asia & Pacific	Reference Category
Days since first case:Europe & Central Asia	0.050	0.024, 0.075
Days since first case:Latin America & Caribbean	0.052	0.023, 0.080
Days since first case:Middle East & North Africa	0.053	0.023, 0.084
Days since first case:North America	0.111	0.050, 0.172
Days since first case:South Asia	0.077	0.036, 0.117
Days since first case:Sub-Saharan Africa	0.061	0.032, 0.090
Days since first case:HICs	Reference Category
Days since first case:LICs	0.024	−0.006, 0.054
Days since first case:LMICs	0.032	0.009, 0.054
Days since first case:UMICs	0.027	0.007, 0.046

Significant effects are shown in boldface. LMICs: Lower Middle-income countries; UMICs: Upper Middle-income countries; LICs: Low-income countries.

Generalized mixed ZINB models were also adjusted for days since the first case, World Bank region and government policy interventions.

After we adjusted for various confounders, we found only a single, but relevant, 2020 EPI variable: the level of air deficiency. This variable is constructed using weighted average of exposure to household air pollution (55%), fine air particulate matter smaller than 2.5 μm (PM2.5, 40%), and ground-level ozone pollution (5%).

Fig. 2a depicts the increasing effect on the spread of COVID-19 when associated with air quality deficiency. Out of the four factors, the one that impacts health the worst was the quality of the air, the greater the value of air deficiency, the larger was the effect on increasing COVID-19 spread [air deficiency: 0.028, 95%CI (0.004, 0.053), p-value = 0.021]. There were additional covariates significantly associated with the spread of COVID-19. In terms of the spread of the virus, HICs had the worst effect, as depicted in Fig. 2. Moreover, the plot shows a consistent gradient on COVID-19 spread from highest in HIC to lowest in LIC [LICs: −2.080, 95%CI (−3.842, −0.311), p-value = 0.021], when compared to HICs. However, it should be noted that this relationship might be driven by some other underlying variables such as the nature of global travel, which possibly skewed towards HIC in the period studied (Menkir et al., 2020; Schellekens and Sourrouille, 2020). Number of days since the first reported case of COVID-19 was positively associated with the increase of COVID-19 cases across all World Bank regions, except for the East Asian & Pacific region (see Fig. 2c).

We also observed some interesting interactions. In the case of the number of days since the first COVID-19 case was reported across country income levels (Fig. 2d), LMICs and the UMICs had a significant impact on virus’ spread compared to HICs [Days since first case: LMICs 0.032, 95%CI (0.009, 0.054), p-value = 0.005; Days since first case: UMICs: 0.027, 95%CI (0.007, 0.046), p-value = 0.006].

Analysis of mortality with the same set of covariates under the same modeling strategy (results shown in Table S2 and Figure S2 in Supplementary Appendix) shows that the lack of air quality had a small impact on mortality evolution. Air quality was less significant in the rate of mortality than in the spread of the disease analysis. Further, number of days since the first reported case across the World Bank regions is the only significant predictor in the COVID-19 mortality analysis with similar conclusions as in the analysis of spread of COVID-19.

The analysis of the random intercepts for both models shows country-to-country variability even after accounting for all the differences in the underlying covariates reported in our analyses [spread model: 3.205, 95% CI (2.444, 4.204); mortality model: 4.052, 95% CI (2.939, 5.586)]. This indicates that random variability when all covariates are set to zero is significant.

Association of COVID-19 spread and mortality with air quality measures

Given the observed impact of air deficiency on COVID-19 spread and mortality, we fitted similar mixed-effects models as described in the previous section to further investigate this impact. Air deficiency variable is broken down to three specific air quality measures: Air pollution (HAD), PM2.5 exposure (PMD) and Ground-level ozone exposure (OZD). Our analysis also used the other ecosystem vitality and environmental health variables described in Table 1 and adjusts them by the same aforementioned set of confounders (model estimates for COVID-19 spread and mortality are given in Tables S3 and S4 in Supplementary Appendix S4). Fig. 3 presents the effects of the significant factors for both outcomes: spread (see Fig. 3a) and mortality (see Fig. 3b).

Fig. 3

Air pollution (HAD) effects for COVID-19 spread in (a) and mortality in (b). In (a) the blue line represents the effect of air pollution on COVID-19 spread and the upper and lower bands represent the 95%CI. Equivalently in (b) for COVID-19 mortality evolution, the effect and bands appear in color orange. The COVID-19 spread and mortality units are based on rates via the offset as produced by the related models. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 3a and b shows a similar tendency in exposure to air pollution, where the effect is more significant on mortality than on COVID-19 spread [Air pollution (spread): 0.021, 95%CI (0.001, 0.043), p = 0.049; Air pollution (mortality): 0.029, 95%CI (0.005, 0.054), p = 0.019]. The larger is the level of air pollution, the greater is the influence of that factor on COVID-19 spread and mortality. No other analyzed air quality factor seems to have an effect on spread or mortality.

Association of COVID-19 spread and mortality with the level of air deficiency. Stratification by country-income level

Some countries possibly do not capture at an equal level, like other countries, the spread of COVID-19 and mortality, because there might be a certain variability of surveillance infrastructures and monitoring systems among the regions. Based on that assumption, we stratified the analysis of the mixed models presented in the previous two sections by country-income level.

The analyses of results are given in Table S5 through Table S8 in the Supplementary Appendix S4 and show that the lack of air quality consistently made an impact on the spread of the disease and its mortality across the High and Upper-middle income countries [Air deficiency (spread): 0.037, 95%CI (0.009, 0.064), p-value = 0.008; Air deficiency (mortality): 0.032, 95%CI (0.003, 0.062), p-value = 0.031]. The plot, in Figure S3, in Supplementary Appendix S3 illustrates the effects of air deficiency on the spread of the disease and its mortality. However, the impact of air deficiency was not observed when the factor of Low, and Lower-middle, income countries were considered.

COVID-19 spread association with specific air pollutants; analysis of three cities

Our time series modeling, for each city, shows a significant effect for ground-level ozone in Barcelona [O3: 0.04, 95%CI = (0.02, 0.08), p-value<0.001] and in Milan [O3: 0.02, 95% CI = (0.01,0.03), p-value = 0.003], which implies that high values of ground-level ozone in the previous 14 days increases COVID-19 spread. Similar patterns in terms of atmospheric particulate matter of 10 μm or less (PM10) were observed in the case of Denver. Thus, PM10 was positively correlated to COVID-19 daily cases for the period analyzed [PM10: 0.02, 95%CI = (0.06, 0.26), p-value = 0.002] even after the underlying models controlled for the set of covariates mentioned before. Details of the model estimates are shown in Supplementary Table S9 of the Supplementary Appendix S4.

Fig. 4 shows the COVID-19 time series overlapped with the 14-day lag daily PM10 or O3 levels for the three cities. The grey shaded areas depict the period when each city was under a strict intervention level (“very high”) following a similar methodology as used at the country level. This point only concerns the time series of COVID-19 daily cases as the time series of PM10 or O3 are 14-day lagged. Thus, the time series of those pollutant measures were in the period when a “high” level of intervention was implemented and, therefore, when there still was a business and industrial activity that may generate air pollutants. Only the factor with a significant effect on COVID-19 daily cases is shown in each graph. We observe similar trend patterns within each city and some peaks occur concurrently. For instance, in plot (c) for Denver, there are simultaneous peaks of PM10, and COVID-19 daily cases, in the period between the ends of March and April. We then calculated Spearman correlations between the COVID-19 cases and the air pollutant measures. The results supported what we had observed with the time series [Barcelona COVID-19 spread and O3: ρ = 0.38, p-value = 0.001; Milan COVID-19 spread and O3: ρ = 0.55, p-value<0.001; Denver COVID-19 spread and PM10: ρ = 0.48, p-value<0.001].

Fig. 4

Time series of COVID-19 daily cases with ground-level ozone (O3) and atmospheric particulate matter PM10 and PM2.5 for Barcelona (a), Milan (b), and Denver (c). The grey shaded area highlights the period when the city was under a strict intervention level.

Discussion

Like other major known epidemics, such as Ebola or SARS, the emergence of the COVID-19 is not unrelated to the climate and biodiversity crises we are experiencing. There is evidence that human health is intimately connected to the intervention humans have in the natural world (Seymour, 2016; Thompson Coon et al., 2011). Our research focused on the association of environmental, biodiversity, and ecosystem factors with COVID-19 spread and mortality between January 21st and May 18th, 2020, where we adjusted for several sociodemographic, social distancing, and county-income level factors. There are several aspects of the results worthy of more discussion than presented in our manuscript. First, a country’s level of biodiversity was of moderate impact on the spread of the disease. Secondly, the COVID-19 spread appeared to be smaller in countries with higher temperatures. Thirdly, there was a direct impact between the level of air deficiency and the spread of COVID-19 and mortality evolution. Particularly, the higher the exposure is to indoor air pollution from household use of solid fuels (see the HAD variable in Table 1), the higher is the impact in both the COVID-19 spread and mortality. Fourthly, the lack of air quality consistently made an impact on the spread of the disease and mortality in HICs and UMICs, rather than LICs and LMICs. Fifthly, as presented in the case studies on cities, COVID-19 transmission is associated is with the 14-day lag ground-level ozone and atmospheric particulate matter of 10 μm, or less, levels.

The appearance of COVID-19 has highlighted the extreme importance of the combat against the loss of biodiversity (WHO, 2017). Mounting evidence of the relationship between this new disease and the reduction in biodiversity requires urgent attention (Lorentzen et al., 2020). In an ecosystem with more biodiversity, a quick spread of the pathogen is harder. Loss of biodiversity can affect the transmission of infectious diseases (Keesing et al., 2006) and provides an opportunity for viruses to pass between animals and people (Keesing et al., 2010). Researchers have reported that the disturbance of natural ecosystems increases wildlife-to-human disease jumps, which has been suggested as the principal cause of neglected, forgotten, and the recent occurrence of unknown human diseases (Epstein et al., 2003). For example, three studies detected a strong association between low bird diversity and an increased human risk or incidence of West Nile encephalitis in the United States (Allan et al., 2009; Ezenwa et al., 2006; Swaddle and Calos, 2008). The results of our study based on the NBI are aligned with these assessments. Thus, preservation and sustainable management of biodiversity might be necessary to mitigate climate disruption and prevent pandemics in order to protect health and wellbeing for the generations to come (United Nations Secretary-General, 2020). Our findings suggest there is a relationship between biodiversity factors and the spread of COVID-19. Our study aims to highlight this possible relationship with the underlying data collection conditions and over the period from January 21 to May 18. However, we are aware that this relationship could be driven by some other underlying variables, such as socio-economic development (Travaglio et al., 2020) and, therefore, further analysis is required.

A new study showed that wild animals that are known to be host to pathogens and parasites, which can jump to humans, become a greater share of the local animals on those sites where humans have turned natural habitats into secondary, agricultural, or urban ecosystems (Gibb et al., 2020). The study found that the latter effect was strongest for zoonotic diseases with host species such as rodents, bats, and passerine birds. Some of those species are the ones that survive after humans diminish biodiversity, and as a consequence, there is a higher risk of dangerous pathogens that can make the leap to humans (Tollefson, 2020) However, the latter connection between the loss of biodiversity due to human development and disease outbreaks does not predict the next pandemic or the current COVID-19 spread, which prompts to keep investigating to establish a mechanism with the current pandemic.

The WHO estimates that around 7 million people die every year from exposure to fine particles in polluted air that lead to diseases such as stroke, heart disease, lung cancer, chronic obstructive pulmonary diseases, and respiratory infections, which include pneumonia (WHO, 2016). Air pollution is also known to weaken the immune system because and to compromise a person’s ability to fight off infection according to the European Public Health Alliance (Vettore, 2020). In particular, household air pollution has contributed to 3.8M deaths (WHO, 2016). A 2003 study found that patients with SARS, a respiratory virus disease closely related to COVID-19, were 84% more likely to die if they lived in areas with high levels of pollution (Cui et al., 2003). Our study showed that there is also an association between air pollution and COVID-19. However, the impact is not the same among countries, because in HICs and UMICs the impact of air pollution was higher. One possible reason for this discrepancy between countries with different income levels is that HICs and UMICs have higher industrialization levels, which in turn causes higher levels of air deficiency (Vigo et al., 2020).

Industrialized countries have driven 5.3 million of the confirmed cases of COVID-19 worldwide and 350,000 global deaths to date (WHO, 2020b). Those countries tend to have higher levels of atmospheric particulate matters (PM10 and PM2.5) and ground-level ozone (O3) pollutants, which originate from man-made sources, such as vehicles and industrial emissions. PM10 and PM2.5 are fine particles, which tend to stay longer in the air, and ground-level ozone (O3) is an irritant gas. All those pollutants trigger or worsen respiratory chronic diseases (Vigo et al., 2020). This is a serious health issue with scientific evidence. For instance, one recent study found that an increase of just 1 μg per cubic meter of PM2.5 corresponded to a 15% increase in COVID-19 deaths (Wu et al., 2020). Researchers have also observed that the consequence of an increased level of O3 is worse than the impact of PM10 level on COVID-19 spread (Cui et al., 2003; Australian Governament, 2020; Staehelin et al., 2001; United States Environmental Protection Agency, 2020). We studied the levels of PM10, PM2.5, and O3 in three industrialized cities and observed an association of those pollutants with COVID-19 spread. The number of COVID-19 cases correlated with levels of O3 in the two European cities (Barcelona and Milan) and an association was observed with the levels of PM10 for Denver. In support of our findings, recent studies marked the relation of PM levels and other pollutants with COVID-19 spread (Travaglio et al., 2020; Zoran et al., 2020) mentioning as a potential mechanism the possible virus attachment to large pollutants (Reche et al., 2018).

Apart from the mentioned findings related to environmental and biodiversity factors, there were other significant effects. When our study was stratified by World Bank regions, then as the number of days became larger, since the first COVID-19 case was reported in a particular World Bank region, the effect on COVID-19 spread was greater. The only exception was the East Asian & Pacific region. One plausible explanation is that we know the disease struck first in Asia, then Europe, North America, and so on. Therefore, the region of Asia has had a longer period to reduce the COVID-19 transmission than the other regions where the disease arrived later. Moreover, the disease has made an impact in regions where more HICs countries exist. Our analysis also determined that the wealthier a country is, the lesser the effect is on disease transmission, as the number of days from the first case increases (Clouston et al., 2020). A possible argument is that HICs perform a higher number of COVID-19 test, have better prepared national health systems and organizations, on the other hand, Governments of LICs lack the resources to implement mitigation measures.

Finally, we conducted a sensitivity analysis to assess our inferences for large countries in terms of area extension. For the sensitivity analysis, we fitted the models again for data on COVID-19 spread and mortality removing the five top countries with the largest area (i.e., Russia, Canada, China, USA, and Brazil). Taking away those large countries, we were able to check if they influenced the inferential analyses. In both cases, mortality and spread measures were not only in line with the definitive model estimates, but also the significant environmental covariates were estimated with similar effects to those when the five countries were included. We only observed changes in the income covariate, which is possibly due to the exclusion of three UMIC countries (Brazil, Russian, and China) with high level of cases. The results of the spatio-temporal analysis are shown in Tables S10 and S11 in Supplementary Appendix S4 for spread and mortality, respectively.

Limitations

To the best of our knowledge, this is the first study that uses publicly available COVID-19 data to analyze the spread of the disease and mortality related to environmental factors and biodiversity levels, adjusted by social demographic, country-income level, and government policy intervention confounders, among 160 countries, globally. However, there are several limitations which must be remarked:

Our inferences are drawn using observational data. To the extent data available, our inferences adjust for the differential covariates across observational units. Ideally, we should have fitted the models in a randomized design, but such design is impossible to pursue in the current settings. However, an extension of this work could use post-randomization techniques based on matching or weighting-based random sampling methods that specifically target potentially varying background characteristics. Unfortunately, there is a lack of readily-available methods to address the intricate nature of correlations and the nature of spatio-temporal data designs. Based on our analytical findings of the correlation structure, one should anticipate that the role of the study design would be significant in such considerations.

This is an observational study, which includes a potential risk of bias that these kinds of studies carry as reported by previous similar studies (Wu et al., 2020; Travaglio et al., 2020). Thus, the results of this study should not be used to make individual-level inferential conclusions. Moreover, this study could also be proned to unmeasured confounding bias. However, we tried to adjust for the most important confounding factors such as population density, time since the beginning of the epidemic, government interventions, weather, and socio-economic factors. Additionally, due to the observational nature of the current study, reported results could be sensitive to specific modeling choices. To assess the sensitivity of such modeling approaches, we also fitted models stratified by country income level (see Tables S5-S8 in Supplementary Appendix S4).

A conclusive capture of temporal COVID-19 spread and mortality trends may not exist due to several factors such as completeness of WHO COVID-19 datasets and government interventions that may be announced one particular day and applied effectively after several days.

Some of the countries have a reliable reporting system, while others do not. Thus, there is variability in the monitor and surveillance of the number of COVID-19 cases and deaths per country and region. For instance, COVID-19 cases and deaths could be underreported (Piovani et al., 2020; Islam et al., 2020). To this extent, this study was unable to take into account asymptomatic cases (Oran and Topol, 2020). Nonetheless, the countries still must report spread and mortality figures at the national level to the WHO with the use of a certain criteria. In that sense, the association of COVID-19 spread and mortality with ecosystem vitality and environmental health factors was also analyzed stratified by country-income level with the aim to capture diversities among surveillance and healthcare systems (Barber et al., 2017). The results showed a consistency with the global analysis.

Daily screening COVID-19 tests per country (Udugama et al., 2020) was not used, because only a limited number of countries had reported that information. Other studies have also noted that the COVID-19 testing rate by country for the same period was challenging (Islam et al., 2020). As we also had a reduced number of countries with environmental and NBI variables, we preferred not to use the COVID-19 test data as it could alter the parameter estimations due to the use of a small sample of countries.

We used covariates in the spatio-temporal model related to precipitation and temperature. Those covariates only summarize the precipitation and temperature levels for the centroid of each country (e.g., we only have the values of Rivas-Vaciamadrid (near Madrid) in Spain), which is not representative of the entire country. This is only used as a proxy and any error associated is incorporated in the zero-inflated negative binomial probability distribution. However, an extension of this work could use a more precise geographical proxy focused on subnational data of a set of countries.

Our analysis is based on the National Biodiversity Index (NBI) by country. We understand that the country’s measure of biodiversity has a great variation between areas of a country (urban or not), but probably the transmission of COVID-19 occurred mainly on urban areas. A correlation analysis of our confounded models with city-level data would be interesting. To the best of our knowledge, there are no public city-level data sets available with biodiversity indexes. Additionally, we applied sensitivity analysis in our spatio-temporal models, excluding the top 5 countries with the largest area (i.e. Russia, Canada, China, USA, and Brazil). The results remained in the same direction.

The 2020 Environmental Performance Index is the most updated data set at the country-level. However, the air quality and pollutant measures were calculated from the EPI report of the previous year. Therefore, they were not measured at the same period as the time series of COVID-19 daily cases. Therefore, our study tries to assess an association between environmental trends and COVID-19 transmission and mortality, during a four-month period.

A limitation of our mortality models is that they did not consider a measure of the age-density of the population or an indirect measure of life expectancy per country. Additionally, the factors related to an individual’s micro-environment (such as the workplace, schools, etc) and behavior (i.e., outdoor activities, smoking habits, etc) are highly related to people’s exposure to air pollution (Travaglio et al., 2020) and could not be assessed in the current study.

A limitation of our analysis is we did not use a higher spatial resolution in each country to perform an analysis at a local scale. To have the same data set but for the local scale is almost impossible. Even if such data were available, significant parameter correspondence and reliability issues would lead to substandard quality in the statistical inference. In addition, these issues would require statistical interventions to compensate them and, unfortunately though, they are beyond the scope of this study.

Our investigation and analysis focused on data variations in the spread and mortality of COVID-19 from January 21st to May 18th, 2020. For country subsets LIC, LMIC, UMIC, and HIC, we ran a sensitivity analysis with the use of a Granger test that compared trends until May 18th and until June 30th of the median of number of cases. The results showed that the trends for HIC and LIC countries until June 30th were similar to the trends until May 18th. However, the trends before and after May 18th for LMIC and UMIC were significantly different. Therefore, the results of this article should be interpreted with caution as they only relate to the underlying data collection conditions and period. As COVID-19 is an infection with a dynamic transmission and all the covariates we used might change, we do not think it would be appropriate to make conclusions beyond May 18th as further data and analysis would be required.

Conclusions

The COVID-19 pandemic is probably a consequence related to the global crises of biodiversity loss and environmental health. The extent and significance in the association between air quality variables, a country’s measure of biodiversity, and COVID-19 transmission and mortality were measured in this study, which also prompts to continue with further investigations to reveal the mechanism of this relationship. Air pollution reduction and biodiversity preservation are complex problems that depend on government actions and financial resources. Our findings give insights that may help governments plan environmental and health policies, as alternative strategy to respond to new COVID-19 outbreaks and prevent future crises.

Author contribution

Contributors: Daniel Fernández, Stefanos Tyrovolas, Marianthi Morena, and William Pan designed the study. Daniel Fernández, Stefanos Tyrovolas, Marianthi Morena, Iago Giné-Vázquez, and Marta Nai Ruscone obtained the data from public available sources. Daniel Fernández, Stefanos Tyrovolas, Iago Giné-Vázquez, Recai Yucel, Ivy Liu, and Marta Nai Ruscone analyzed the data and drafted the manuscript. All authors contributed to the interpretation of the data and revision of the manuscript. All authors had primary responsibility for final content and act as guarantors. All authors read and approved the final manuscript.

Funding sources

Daniel Fernández supported by grant 2017 SGR 622 (GRBIO) administrated by the Departament d’Economia i Coneixement de la Generalitat de Catalunya (Spain), by grant RTI2018-100927-J-I00 administrated by Ministerio de Ciencia e Innovación (MCI, Spain), by the Agencia Estatal de Investigación (AEI, Spain), and by the European Regional Development Fund (FEDER, UE). Stefanos Tyrovolas was supported by the Foundation for Education and European Culture, the Miguel Servet programme (reference CP18/00006), and the Fondos Europeos de Desarrollo Regional. William Pan is supported by NASA-ROSES Grant NNX15AP74G. Ivy Liu was supported by Marsden grant E2987-3648 administrated by the Royal Society of New Zealand.

Data availability statement

All data used for this study are publicly available (information on how to access the data is provided in text). The compiled COVID-19 dataset and algorithms of the current study are available upon request to the corresponding author. The corresponding author (Daniel Fernández) affirms that the manuscript is an honest, accurate, and transparent account of the study being reported; that no important aspects of the study have been omitted; and that any discrepancies from the study as planned (and, if relevant, registered) have been explained.

Ethical approval

The study’s protocol has been approved by the research ethics board at the Parc Sanitari Sant Joan de Déu (PIC-67-20, Barcelona, Spain).

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.