Association between live childhood vaccines and COVID-19 outcomes: a national-level analysis.

We investigated whether countries with higher coverage of childhood live vaccines [BCG or measles-containing-vaccine (MCV)] have reduced risk of coronavirus disease 2019 (COVID-19)-related mortality, while accounting for known systems differences between countries. In this ecological study of 140 countries using publicly available national-level data, higher vaccine coverage, representing estimated proportion of people vaccinated during the last 14 years, was associated with lower COVID-19 deaths. The associations attenuated for both vaccine variables, and MCV coverage became no longer significant once adjusted for published estimates of the Healthcare access and quality index (HAQI), a validated summary score of healthcare quality indicators. The magnitude of association between BCG coverage and COVID-19 death rate varied according to HAQI, and MCV coverage had little effect on the association between BCG and COVID-19 deaths. While there are associations between live vaccine coverage and COVID-19 outcomes, the vaccine coverage variables themselves were strongly correlated with COVID-19 testing rate, HAQI and life expectancy. This suggests that the population-level associations may be further confounded by differences in structural health systems and policies. Cluster randomised studies of booster vaccines would be ideal to evaluate the efficacy of trained immunity in preventing COVID-19 infections and mortality in vaccinated populations and on community transmission.

Introduction

Trained immunity, or long-term boosting of innate immune responses, by live vaccines [Bacillus Calmette–Guérin (BCG), measles-containing-vaccine (MCV)] can induce heterologous protection against other pathogens including RNA viruses [1]. The concept of trained immunity has been established by animal/human experimental studies as well as epidemiological studies. It has been speculated that widely administered vaccines could be an important tool for reducing susceptibility to and for decreasing severity of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection [1]. One hypothesis is that heterogeneity of mortality rates among different age groups and different countries may be explained to some extent by differing degrees of trained immunity from vaccines that could provide some basal level of protection from this novel pathogen.

Several ecological studies have evaluated the impact of BCG on Coronavirus disease 2019 (COVID-19) outcomes (case incidence, death and case fatality ratio) [2, 3]. The results are conflicting, likely due to insufficient control of important factors, including differential entrance of the epidemic to each country, vaccine coverage and other health metrics. If live pathogen vaccinations are associated with decreased risk of poor COVID-19 outcomes on a country basis, one might postulate that vaccine coverage could be beneficial both through direct effects on individuals as well as indirectly by decreasing community viral load and hence transmission to others [1, 3]. We used publicly available national-level data from 140 countries to investigate whether countries with higher proportion of people who received childhood live vaccines (BCG or MCV) have reduced COVID-19 mortality or fewer cases, adding new considerations to previously published analyses.

Methods

Data source

All of the data used in this analysis consists of national-level estimates. We acquired data on COVID-19 through 13 July 2020 from the publicly available databases, Our World In Data (ourworldindata.org) for COVID-19 statistics, and Foundation for Innovative New Diagnostics (https://www.finddx.org/at-a-glance/) for COVID-19 testing rates (cumulative number of tests per million people), which we standardised as below. Data on BCG and MCV coverage from 1980 to 2018 were extracted from WHO/UNICEF website (https://apps.who.int/immunization_monitoring/globalsummary/timeseries/tswucoveragebcg.html). We used published estimates of the Healthcare Access and Quality Index (HAQI) 2016 which were obtained from the Global Burden of Disease (http://ghdx.healthdata.org/record/ihme-data/gbd-2016-healthcare-access-and-quality-index-1990-2016). Other health metrics, including national life expectancy, hospital beds per capita were obtained separately from Our World In Data (ourworldindata.org), publicly available databases.

Outcomes and exposures of interest

We defined the start of the epidemic in each country as the date that country reported a total of 100 cases, and evaluated cumulative deaths and cases reported during the first 60 days of each respective epidemic. A proportion of vaccinated people in the entire country during each year was calculated by multiplying the proportion of the eligible population actually vaccinated by the estimated proportion of vaccine eligible people (proportion of population under age 5 annualised, given these are both childhood vaccines). With vaccine coverage from each year, we calculated a cumulative vaccine coverage index in each country from each specified year through 2018 (BCG and MCV coverage indices) where we prespecified several possible durations of trained immunity effects, and therefore selected the following years to begin the vaccine exposure: 2005 (cumulative vaccinations in the last 14 years since 2005 inclusive), 1995 (in the last 24 years since 1995 inclusive), 1985 (in the last 34 years since 1985 inclusive) and 1980 (in the last 39 years since 1980 inclusive). We also assigned zero BCG vaccine coverage to countries for each year for which universal vaccination was not implemented as per national policy records (http://www.bcgatlas.org/).

To calculate a vaccine coverage index, formally, let V denote the vaccine of interest (BCG or MCV), i denote the country of interest, t denote the calendar year and VC denote the coverage of vaccine V for country i during year t. For each country i, Eq. (1) shows the vaccination coverage index spanning M prior to (and excluding) 2019.where is the estimated proportion of the national population of country i who were eligible to receive vaccine V in year t, and is the estimated proportion of the eligible population of country i who received the vaccine V in year t. Prespecified values of M equal to 14, 24, 34 and 39 were chosen to explore different postulated durations of trained immunity. Hence, each assumed value of M resulted in distinct vaccine indices with different cutoff years, representing different possible durations of trained immunity. For example, the 2005 vaccine indices represent the proportion of people who could directly benefit from augmented immunity to other pathogens (including SARS-CoV-2) in the entire country if the effect of trained immunity were to last 14 years. As expected, the number of countries with complete vaccine data in more distant years (e.g., 1980 vaccine indices) is smaller than that available for more recent years (since 2005). We chose the 2005 BCG/MCV indices (14-year coverage) as our vaccine variables based on their best prediction of the primary outcome (log-transformed COVID-19-related deaths per million) in bivariate models compared to the more expansive indices (1980−1995) and the more conservative assumption that trained immunity may not last as long as vaccine pathogen-specific immunity. However, all BCG and MCV indices were highly correlated with the indices included that data from 1980, showing relative stability of vaccination rates in each country over time (Supplementary Fig. S1). While the publicly available country-level data that are used to produce these vaccine variables are typically reliable, it is worth noting that these are all estimates and the data quality varies across years and countries.

Analysis

Our analysis sought to infer (i) if COVID-19-related death was associated with vaccine exposure variables after controlling for known national characteristics of healthcare systems and population health such as life expectancy, number of hospital beds, physicians per population, gross domestic product and composite measures of public health policies and (ii) whether the association between vaccination exposure and COVID-19 deaths differed as a function of these known national-level variables. Due to our small sample size (n = 140 countries), we used the HAQI, as a single metric of population health in our models to handle and include the high-dimensional correlated predictors, while isolating the effects of our exposures. HAQI 2016, a publicly available, validated summary score of healthcare systems developed by the University of Washington's Institute for Health Metrics and Evaluation, has gone through rigorous external validation, and encompasses all variables under consideration [4]. While the variables used to construct the HAQI are not all publicly available, published estimates of the HAQI are publicly available.

We included all countries with data on vaccine coverage, COVID-19 outcomes, and HAQI available to conduct multiple linear regression models. We standardised our vaccine exposure variables as well as the HAQI to have a mean of zero and standard deviation of 1. We used nested linear models to explore:

Whether each of the BCG or MCV index is associated with COVID-19-related deaths conditional on HAQI.

Whether BCG and MCV indices are simultaneously associated with COVID-19-related deaths conditional on HAQI.

Whether HAQI modifies the association between both BCG and MCV indices with COVID-19related deaths.

Whether HAQI modifies the association between BCG index and COVID-19-related deaths.

Whether MCV index modifies the association between BCG index and COVID-19-related deaths.

Results

One hundred and forty countries had complete publicly available data with a wide range in COVID-19 cumulative cases (3.3−9862.2 per million people) and mortality (0−721.2 per million people) by 60 days after epidemic start (Supplementary Table 1). Variation of BCG coverage (BCG index 2005 with 136 countries: median 0.24, range 0−0.5) was wider than that of MCV coverage (MCV index 2005 with 138 countries: median 0.23, range 0.11−0.47). Wide variation of HAQI among 140 countries was also observed (median 68.1, range 18.6−97.1).

BCG was marginally associated with fewer reported COVID-19 death rates, with the association remaining after adjusting for HAQI (Table 1). However, the magnitude of the association between BCG and COVID-19 death rates increased as a function of HAQI, as reflected by a statistically significant interaction between the BCG vaccine variable and HAQI. Similar patterns were observed for BCG in multivariable models that included both BCG and MCV vaccine exposure variables. MCV was also significantly associated with fewer reported COVID-19 death rates, but this association was no longer significant after adjusting for HAQI. MCV coverage did not notably change the association between BCG and COVID-19 death rates. Neither BCG nor MCV coverage were associated with reported COVID-19 case rates after adjusting for HAQI (not shown).

Table 1.

Associations between vaccine coverage and COVID-19 mortality (total COVID-19-related deaths per million at day 60 following epidemic start)

	Estimate (s.e.)	P value
Models for BCG
Exposure only
Intercept	2.357 (0.114)	<0.001
BCG2005a	−1.013 (0.114)	<0.001
Additive with HAQI
Intercept	2.370 (0.113)	<0.001
BCG2005a	−0.664 (0.245)	0.008
HAQI	0.396 (0.246)	0.110
HAQI interaction
Intercept	2.137 (0.158)	<0.001
BCG2005a	−0.573(0.246)	0.021
HAQI	0.484 (0.246)	0.052
BCG2005*HAQIa	−0.270 (0.129)	0.038
Models for MCV
Exposure only
Intercept	2.362 (0.121)	<0.001
MCV2005b	−0.859 (0.121)	<0.001
Additive with HAQI
Intercept	2.373 (0.114)	<0.001
MCV2005b	−0.266 (0.182)	0.147
HAQI	0.764 (0.182)	<0.001
HAQI interaction
Intercept	2.185 (0.158)	<0.001
MCV2005b	−0.389 (0.195)	0.048
HAQI	0.698 (0.185)	<0.001
MCV2005*HAQIb	−0.243 (0.142)	0.091
Models for BCG and MCVc
Exposures only
Intercept	2.354 (0.114)	<0.001
BCG2005a	−1.142 (0.250)	<0.001
MCV2005b	0.146 (0.251)	0.561
Additive with HAQI
Intercept	2.367 (0.114)	<0.001
BCG2005a	−0.777 (0.340)	0.072
MCV2005b	0.120 (0.250)	0.633
HAQI	0.388 (0.247)	0.238
HAQI interactions
Intercept	2.089 (0.171)	<0.001
BCG2005a	−0.349 (0.389)	0.744
MCV2005b	−0.216 (0.291)	0.744
HAQI	0.518 (0.251)	0.128
BCG2005*HAQIa	−0.330 (0.152)	0.128

The results from linear regression models using national-level aggregate data for both predictors and outcomes are presented. Three sets of nested models are examined for the BCG index 2005 as exposure, the MCV index 2005 as exposure, and both BCG index 2005 and MCV index 2005 as exposures. The models are shown in ascending order of complexity, starting from unadjusted models with only the exposures to the most complex models including interactions between each exposure and HAQI.

MCV, measles-containing-vaccine; BCG, bacillus Calmette–Guérin; HAQI, Healthcare Access and Quality Index (2016).

Estimated proportion of vaccinated people with BCG in the country during the last 14 years.

Estimated proportion of vaccinated people with MCV in the country during the last 14 years.

Reported P-values for the models using both BCG index 2005 and MCV index 2005 as exposures account for multiple testing using Holm's method.

We performed additional univariate analyses with COVID-19 testing rate and other outcomes and exposures to elucidate relationships between COVID-19 metrics and other country-level health and development metrics for better understanding of the above observations (Supplementary Fig. S2). Because HAQI represents a composite of many relevant health metrics, we sought to evaluate the relationships between several individual variables with publicly available data to assess for possible confounders or effect modifiers for the relationship between vaccine coverage and COVID-19 outcomes. COVID-19 testing rate was associated with higher HAQI. COVID-19 testing rate and other markers of better health infrastructure (life expectancy, number of hospitals per population) were negatively associated with BCG index 2005 and MCV index 2005. Furthermore, both COVID-19 testing rate and life expectancy were associated with higher COVID-19-related mortality rates. These relationships between healthcare metrics, vaccine coverage and COVID-19-related mortality rates suggest that several healthcare system metrics (frequency of COVID-19 testing, life expectancy) are inversely associated with vaccine coverage but positively associated with reported COVID-19 death rates.

Discussion

MCV and BCG vaccination have demonstrated benefits on reducing non-infection specific childhood mortality in multiple randomised-controlled trials (RCTs), the mechanism of which is thought to be non-specific protection against heterogeneous pathogens (trained immunity) [5]. The emergence of SARS-CoV-2 heightened the interest in trained immunity, and recently BCG-induced cross-reactive T cell response to SARS-CoV-2 has been reported [6], although prior ecological studies to investigate the association between BCG vaccines and COVID-19 outcomes have had conflicting results. Much of the heterogeneity in prior studies may have been contributed to by simplistic analyses that did not account for timing of epidemic onset or different healthcare systems in each country [2]. Although it remains unclear how long the trained immunity effect could last: other studies have either unrealistic estimations for duration of effect (e.g., up to 75 years) or simply categorised a country with a binary designation of being a BCG-promoting country, without regard to yearly vaccine coverage or coverage across time [1, 7].

The strengths of this analysis include use of the same relative ascertainment window for COVID-19 outcomes for each country considering differential timing of epidemic arrival. Limiting outcome measures to the first 60 days of the epidemic may also mitigate the effects of the highly variable political and public health responses following the ‘first wave’. We selected deaths as a primary endpoint as a downstream parameter to account for all possible vaccine effects (preventing acquisition of the virus, symptomatic episodes, severe illness or onward transmission) [1, 3, 7]. Our models included BCG, MCV and HAQI considering both effect modification and confounding [4]. Although we selected 2005 BCG/MCV indices since they had the strongest association in bivariate analysis and quantified the cumulative vaccine exposure for each vaccine over the last 14 years, vaccine coverage dating back to 1980 for each vaccine were also highly correlated (Supplementary Fig. S1). Ranging vaccine coverage back four decades allows for extrapolation to vaccine coverage for a larger proportion of the middle-age and older population at risk for COVID-19 mortality, and allows for a sensitivity analyses to evaluate longer trained immunity effects. The observed weakening of effects when longer duration vaccine coverage is included could have multiple interpretations. One possibility is that there is less protection seen with inclusion of more distant vaccination. It is possible that simply using older vaccine coverage data biases results toward the null due to higher variability and less completeness in older data, or less strong associations with current health metrics. Regardless of the metrics used in our analyses, the actual vaccine exposure in the population remains the same, but inclusion of vaccine coverage back to four decades may simply be less specific for both possible vaccine effects or as a surrogate for other unmeasured health statistics.

The current study demonstrated marginal associations between coverage of both childhood vaccines with COVID-19-related mortality. The significance of BCG coverage remained in the conditional model after controlling for HAQI. The association of MCV coverage with COVID-19 death rates was weaker than that of BCG coverage in all models. This may reflect lower variability of MCV coverage between countries or indicate this is a surrogate marker for other factors associated with COVID-19 outcomes, if not a biological difference. Somewhat counter-intuitively, COVID-19-related mortality rates were higher in countries with higher HAQI. Since countries with higher HAQI tended to have higher COVID-19 testing rate early in the epidemic (Supplementary Fig. S2), this likely describes a phenomenon of under ascertainment of both cases and deaths in countries with less robust health infrastructure, or else employed COVID-19 containment strategies less reliant on testing (masking, social distancing). Similarly, countries with lower BCG vaccine coverage were more likely to have higher rates of COVID-19 testing and reported COVID-19 associated mortality. Due to small sample size (n = 140 countries) and compounding of uncertainty within each country-level metric, we elected to not undertake a full multivariable analysis with a more inclusive list of available health metrics. Interestingly, the magnitude of the association between BCG coverage and COVID-19 death rates increased as a function of HAQI. This suggests either that there is a true relationship between BCG coverage and COVID-19 outcomes which are masked by health-system heterogeneity, or else the association could be attributed to an unmeasured confounder associated with BCG coverage but not HAQI, such as confidence in government or low parity in access to care despite adequate resources [8]. Because the majority of countries with high HAQI no longer use BCG, sensitivity analyses focusing on high-HAQI countries were not suitable.

We explicitly acknowledge the principal limitation that country-level data does not represent individual-level exposures or outcomes [9]. There are other potential confounders, such as variability in other biological exposures including host-related, microbial and environmental exposures; national COVID-19-specific public health responses (e.g., per cent mask use, policy and adherence to movement restrictions) and heterogeneity in population density and disparities in access to care; the majority of which are difficult to accurately measure and quantify [2, 8]. Although studies have indicated certain BCG strains might induce more effective trained immunity than others, there was insufficient data on vaccine type to evaluate these hypotheses [10]. Moreover, our analysis treated the national-level summaries as fixed variables, ignoring the uncertainty in these estimates and possibly their biases, and our models only examined linear associations.

In conclusion, we found an association between higher cumulative BCG coverage and COVID-19-related death rates but not cases, and only a marginal effect of MCV coverage on either. We cannot rule out that these observations are attributable to differential healthcare infrastructure, including COVID-19 testing rate, population age distribution or other unmeasured confounders. Several RCTs (primarily BCG) are currently underway (https://clinicaltrials.gov/), and one compiled (https://www.umcutrecht.nl/en/about-us/news/article/jan-18-tuberculosis-vaccine-does-not-protect-vulnerable-elderly-against-covid-19), to evaluate the impact of heterologous live pathogen vaccines on COVID-19-related outcomes. The first of these RCTs was in >6,000 elderly persons randomised to BCG or placebo in the Netherlands, which showed neither any protection of vaccinees from COVID-19 nor reduction in disease severity. These ongoing studies may however not be able to adequately assess full vaccine effects, which depend on characteristics of the population vaccinated, including whether children, middle-age or older adults are vaccinated, not be able to adequately assess the presence and heterogeneity of preceding BCG coverage [1, 7, 10]. Cluster randomised trials of booster BCG or MCV vaccines could best evaluate both individual and community-level effects of trained immunity.