Main

The Pfizer–BioNTech COVID-19 BNT162b vaccine is highly effective at preventing disease and infection at the individual and community levels, as demonstrated in clinical trials, as well as in real-world vaccination campaigns1,2,3,4,5. Furthermore, among vaccinated individuals infected with SARS-CoV-2, a lower viral load was observed6. Reduced infection and viral load suggest that reduced transmission occurred. However, vaccination could, in principle, also increase transmission due to behavioral effects, as vaccinated individuals may not quarantine after contacting a patient with COVID-19, or may be less mindful of social-distancing measures7. It is therefore unclear whether, overall, vaccination reduces transmission at the population level, thereby conferring protection for those who are not vaccinated, such as individuals who are immunodeficient or those who are currently ineligible for vaccination4,8,9,10. Moreover, as the reproduction number varies with sociobehavioral and environmental factors even in the absence of vaccination and because global disease rates represent both vaccinated and unvaccinated individuals, it has proven challenging to determine the effect of vaccination on community-level SARS-CoV-2 transmission11,12,13.

The rapid vaccine rollout in Israel, initiated on 19 December 2020 and covering almost 50% of the population within 9 weeks, presents a unique opportunity to test these questions using real-world data. Capitalizing on differences in vaccination rates among geographically distinct communities, and on the availability of an unvaccinated bystander cohort of individuals below 16 years of age for whom the vaccine was not authorized in the first stages of vaccine rollouts, we asked whether and to what extent the fraction of patients vaccinated in each community affects the risk of infection in an unvaccinated cohort of individuals under 16 years old within this same community.

We focused our analysis on the vaccination rates and test results of 177 distinct communities with a presumed low rate of natural immunization as inferred by a low fraction of individuals infected with SARS-CoV-2. We retrieved the vaccination dates and test results, from 9 December 2020 to 9 March 2021, of members of Maccabi Healthcare Services (MHS), Israel’s second largest healthcare maintenance organization. We defined geographical communities based on residence codes, and identified 246 communities that each comprised a sufficient number of tests and people (Methods). As both vaccination and natural infection could render individuals immunized, thereby possibly conferring protection to unvaccinated individuals, high infection rates could mask the effect of vaccination-induced immunity. To minimize the confounding effect of natural immunization, we considered the test record since 1 March 2020, and included only communities in which the fraction of people tested positive by 9 March 2021 was lower than 10% (177 communities; Methods and Extended Data Fig. 1). This exclusion of communities did not introduce notable biases in the age distributions of neither the vaccinated individuals nor the patients with a confirmed SARS-CoV-2 infection (Extended Data Fig. 2).

To control for the spatiotemporally dynamic nature of the epidemic, we focused on relative changes in positive test fraction within each community between fixed time intervals, thus accounting for global temporal patterns as well as for intrinsic differences among communities. First, for each community, the mean cumulative fraction of patients inoculated with the first dose of the vaccine was calculated for two consecutive 3-week intervals (V1, V2) for individuals aged 16–50 years, who are assumed to represent the population that is likely to interact with the bystander unvaccinated cohort of under 16 years of age14 (Fig. 1a). Second, for each such vaccination time interval, a corresponding time-shifted 3-week testing interval was defined with a delay of 28 d, to allow for the putative immunization of the vaccinated individuals and subsequent cross-protection of unvaccinated individuals to take effect (Fig. 1a,b). For these two testing intervals, the positive test fraction among the unvaccinated bystander cohort (P1, P2) was calculated (Fig. 1b). Comparing periods 1 and 2, we then considered the relationship between the change in positive test fraction of the unvaccinated cohort (P2/P1) and the change in vaccination (V2 − V1) (Fig. 1c).

Fig. 1: An increase in the fraction of vaccination in a community is associated with a later reduction in the number of confirmed infections of the unvaccinated cohort in the same community.
figure 1

a,b, The accumulated fraction vaccinated among 16–50-year-old individuals (a) and the positive test fraction among the 0–15-year-old bystander unvaccinated cohort (b, 7-d moving window average; to avoid a logarithm of zero, time windows with zero positive tests are plotted as having a single positive test) are shown for each community as a function of time. The vaccination intervals (a) and their corresponding time-shifted testing intervals (b) are indicated above each panel (dark gray and light gray). The black trajectory highlights an example community, indicating its mean vaccination fraction for each of the intervals (V1, V2), and the mean positive test fractions in the testing intervals (P1, P2) (see Methods; the example community is highlighted only for demonstration). c, For each community (n = 177), the change in positive test fraction between time intervals (P2/P1) is shown as a function of the change in the vaccinated fraction in the corresponding vaccination intervals (V2 − V1). Data are mean ± 1 s.e.m. The dotted line shows linear fit (P = 0.0003) (see Methods). The community highlighted in a and b is indicated (black x, error bar and arrow). d, The association is robust to a shift in interval dates. For the same communities (n = 177), slopes and P values are shown as a function of the first vaccination interval start date (keeping the same delay of 28 d between vaccination and test intervals). Circles and error bars represent maximum likelihood estimator (MLE) ± 2 s.e. (95% confidence intervals). The slope corresponding to the intervals used in ac is indicated (filled black circle). **P < 0.005; ***P < 0.0005.

The positive test fraction in the unvaccinated cohort decreased in proportion to the rate of vaccination in each community. A linear fit showed a strong negative association between the log-transformed odds ratio of the positive test fraction and the increase in vaccination fraction (log(P2/P1) = a(V2 − V1) + b, where a = −4.1 (95% confidence interval, −6.3 to −2), P = 0.0003) (Fig. 1c; see Methods, ‘Correlation analysis’). For each 20 percentage points of vaccinated individuals, the positive test fraction for unvaccinated individuals decreased approximately twofold; such amplification could appear owing to multiple rounds of infection cycles. These results were robust to changes in interval length (a = −3.6 (95% confidence interval, −4.8 to −2.3), P < 10−4 for 4-week intervals), interval dates (Fig. 1d), and member and community exclusion thresholds (Extended Data Figs. 3, 4). Furthermore, consistent with the known time for vaccine-induced immunity1 plus one or more infection cycles, the effect appears robustly for time shifts beyond 24 d (Extended Data Fig. 5). Considering the age of the vaccinated cohort, we find that the association was stronger when the considered vaccinated cohort was younger than the reference vaccinated cohort of 16–50 years of age, and became non-significant for vaccination cohorts of individuals aged 66 years and older (Extended Data Fig. 6), consistent with an assumed lower level of social interactions between the two age groups14. We also observed a strong negative correlation when considering, for each community, the positive test fraction per capita instead of per test (a = −8.7 (95% confidence interval, −11.3 to −6.1), P < 10−4). Taken together, our analysis shows a strong and robust negative association between the vaccination of adults in the community and a later decrease in infection of the younger bystander cohort.

Our study has several important limitations. First, beyond vaccine-based immunity, infections of unvaccinated individuals could also be affected by naturally acquired immunity. While we minimize this confounding factor by focusing on communities with low cumulative prior infections, future studies could more directly control for this effect—for example, by including seroprevalence data. Second, individual behavior and public policy guidelines, including enforcement and compliance with a lockdown imposed between 8 January and 7 February 2021 (refs. 15,16,17), might correlate with vaccination rates and also affect the infection potential of the unvaccinated group. Third, the proportion between the measured positive test fraction and the actual infection rate could be different in each community and might also vary with time. Using the logarithmic derivatives log(P2/P1), our analysis accounts for the possibility of different proportions between the infection rate and the positive test fraction in each community as well as for the possibility that this proportion changes uniformly with time, but it cannot account for the possibility that these proportions change over time in different ways for different communities. Finally, our study population is limited to MHS members, who are representative of only part of the overall population for each community.

In this study, we identified a strong negative association between vaccination rate at the community level and the risk of infection for unvaccinated members of the community. We find that high vaccination rates were associated with lower infection rates at later time points among the unvaccinated cohort. Although the observed vaccine-associated protection of the unvaccinated population is encouraging, further studies are required to understand whether and how vaccination campaigns might support the prospect of herd immunity and disease eradication.

Methods

Data collection

Anonymized electronic health records were retrieved from MHS database for the study period: 6 December 2020–9 March 2021. These records include the following information. (1) Patient demographics, which included for each MHS member: a random ID used to link records, year of birth, and a coded geographical location of residence. For anonymization purposes, location of residence was given as random codes for the city or town and geographical statistical area (GSA; the country is divided by the National Bureau of Statistics into GSAs, which correspond to areas or neighborhoods). (2) The test results, which included for any SARS-CoV-2 reverse-transcription quantitative PCR test performed for MHS members: the patient random ID, the sample date and an indication of positive and negative results (total, 2,137,063 tests with 110,788 positive results). (3) Vaccination, which included the patient by random ID and date of inoculation with the first dose of the BNT162b2 mRNA COVID-19 vaccine (1.37 million vaccinated individuals).

Ethics committee approval

The study protocol was approved by the ethics committee of The MHS. The IRB-approved protocol includes an exempt from informed consent.

Vaccination and testing intervals

We defined two consecutive 3-week testing intervals between 6 January 2021 and 16 February 2021 (3-week-long intervals were long enough to allow for a sufficient number of tests; see ‘Defining communities’). For each of these two testing intervals and each community, the positive test fraction among unvaccinated individuals below 16 years old was calculated (Fig. 1b; P1, P2), and their standard errors were evaluated assuming Poisson statistics. MHS members were excluded from this calculation if in total (at any time) they were tested more than 20 times (to avoid patients participating in extensive surveys; a sensitivity analysis is included in Extended Data Fig. 4). In addition, for each MHS member, any test following a positive test was excluded from this calculation. To account for the expected time after the first dose for a vaccinated individual to acquire immunity and to potentially prevent further infections, each 3-week testing interval was matched with a corresponding 3-week vaccination interval that was time-shifted backwards by 28 d. For each vaccination interval and each community, the mean cumulative fraction of individuals aged 16–50 years vaccinated with the first vaccine dose (out of all MHS members aged 16–50 years in the community) was calculated (Fig. 1a; V1, V2).

Defining communities

Communities were defined as groups of members residing in the same GSA geographical code. For each city code, all GSAs with fewer than 2,000 members were joined into a single community. Communities were excluded if (1) they contained fewer than 2,000 members; (2) the number of tests for any 7-d window was smaller than 2 (too noisy to analyze the positive test fraction over time, Fig. 1b); (3) the number of tests in any 3-week time window was smaller than 40 (too noisy for correlation analysis, see below); (4) the number of positive tests in any 3-week time window was zero (which would cause a division by zero or a logarithm of zero when calculating y values for the correlation analysis); (5) they contained fewer than 20 members aged 16–50 years (too noisy for correlation analysis); or (6) by 9 March 2021, more than 10% of the community members tested positive (too high a presumed natural immunization confounder; a sensitivity analysis is included in Extended Data Fig. 3). This analysis resulted in 177 communities (Extended Data Fig. 7). Members not residing in one of these eligible communities were excluded from our analysis.

Statistical methods

Correlation analysis

A linear fit was performed on the log-transformed ratio of the positive test fraction y = log(P2/P1) and the change in vaccination x = V2 − V1 using the built-in function glmfit in MATLAB. The P value of the slope was calculated as the fraction of smaller or equal slopes calculated in 10,000 xy permutation bootstrapping simulations.

Sensitivity analysis

For each main parameter, a correlation analysis was performed with different values of this parameter, while all other parameters were kept fixed. The fixed parameters were those used to produce Fig. 1c. Each of these runs produced the slope and the P value of the slope (Fig. 1c; see ‘Correlation analysis’). The specific choices of parameters that were tested and their values are indicated in Extended Data Figs. 36.

Reporting Summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.