A Batch-Dependent Safety Signal Related to All-Cause Mortlity Associated with COVID-19 Vaccination

1. Introduction

In response to COVID-19, the rapid development and global deployment of vaccines that express the prefusion conformation of the SARS-CoV-2 spike protein in the recipient [1] have been hailed as a pivotal intervention in the pandemic, substantially reducing severe disease and mortality [1,2]. However, these vaccines based on novel platforms, such as modified RNA (modRNA) and adenoviral vector-mediated expression, are associated with an unusually broad spectrum of adverse events affecting multiple organ systems. [3,4,5,6], as well as increased mortality in young people in the weeks following vaccination [7]. A recent study of the Japanese population has identified a significant increase in excess mortality after repeated COVID-19 vaccination [8].

Robust and timely pharmacovigilance is essential for any medical product, and traditional post-marketing safety surveillance for vaccines has largely relied on spontaneous reporting systems for adverse events (AEs) and large-scale epidemiological studies that monitor pre-defined outcomes [9]. Reported suspected adverse events (SAEs) for COVID-19 vaccines account for one quarter of all SAEs recorded in the European Union’s pharmacovigilance database, EudraVigilance [10,11]. Spontaneous reporting systems like EudraVigilance are open to public scrutiny, but are inherently susceptible to under-reporting, reporting biases, and difficulties in establishing causality [12,13]. Conversely, epidemiological studies using healthcare databases provide stronger evidence for associations, but require ethical approval and typically investigate hypotheses based on specific biological mechanisms or signals obtained from spontaneous reports [14].

A relatively underexplored dimension in vaccine safety is the potential for heterogeneity be-tween individual production batches of the same product [15]. According to Good Manufacturing Practice (GMP), batches must meet stringent quality specifications; nevertheless, subtle variations in components, potency, or purity could, in principle, influence reactogenicity or safety profiles [16,17].

Emerging analyses from several countries, including Denmark [18], Sweden [17], and the United States [19], have reported non-random, batch-associated patterns in the reporting rates of serious SAEs for COVID-19 vaccines, suggesting possible batch-dependent heterogeneity. Similar signals were identified in the Czech Republic using spontaneous SAE data [20]. While provocative, these findings are constrained by the limitations of voluntary reporting systems.

To investigate this phenomenon using a more robust and definitive endpoint, we turned to all-cause mortality (ACM). ACM is a hard endpoint, comprehensively and reliably captured in national registries, and less susceptible to the ascertainment biases that affect SAE reporting. This study aims to determine whether specific batches of COVID-19 vaccines administered in the Czech Republic were associated with heterogeneous short-term ACM rates. Further-more, we sought to correlate any observed batch-ACM patterns with previously reported batch-SAE patterns to assess the concordance of signals from these orthogonal data sources.

2. Materials and Methods

Ethical Statement

The study relied solely on anonymized secondary data and was therefore exempt from research ethics board review.

Data Set

Data was obtained on all citizens in the Czech Republic from 1/1/2020 to 13/03/2024 via a freedom of information request (FOIR) to the Institute of Health Information and Statistics (IHIS). Registry data was provided including a line listing for each individual included variables for gender sex (male/female), birth year (year), vaccination date (day/month/year), COVID-19 vaccine product (name), manufacturer (name), batch identification (ID) and date of death (DD-MM-YYYY), if applicable (only available for the period up to and including 31/12/2022). Similarly, a list of all COVID-19 vaccine batches administered from 27/12/2020 to 13/03/2024 was obtained via FOIR to the Czech Republic State Institute for Drug Control (SÚKL), which included the product name and batch ID [20].

Data Selection

A flow chart of data sources and curation is shown in Figure 1. Data for individuals with a missing or incomplete batch number were excluded. The number of deaths per batch was calculated for 10-year functional age-groups for each gender and for deaths occurring within 3 months of the last vaccination date. The number of administered doses per batch was calculated per month. Only batches with a maximum monthly usage during the vaccination peri-od from 27/12/2020 until 31/9/2022 were selected for further analysis. Thus, mortality rates per batch were calculated over a 3-month period from the date of the last vaccination, using the available mortality data (1/1/2020 to 31/12/2022). Age was estimated from year of birth and date of death, where applicable. Crude mortality rates and gender- and age-standardized mortality rates to the total distribution for the dataset were calculated for all batches. The number of administered doses per batch was calculated directly from the number of occurrences in the IHIS registry data.

Data Curation

Pfizer-BioNTech

In total, 68 batches of Comirnaty were identified during the study period (doses administered: mean n = 170,937; mean σ = 178,412; max n = 823,343; min n = 1,245; Table S1). Batch FR8477 was administered only twice and was omitted from the analysis. Batch PCB0017 was administered mainly to younger (n = 1086 administrations; y = 29.4 years of age) and a few elderly individuals. This batch was associated with a few deaths in older age groups (1 death in each age group 70-79, 80-89, 90-99), leading to a high age-standardized rate and was also excluded from further analysis (not shown). A total of n = 66 batches remained for analysis.

Moderna

A total of 40 batches of mRNA-1273 (SPIKEVAX) were administered during the study period (doses administered: mean n = 39,483; mean σ = 26,085; max n = 82,070; min n = 502; Table S1). Batches 000256A, 016G21A, 000162A, 000202A, and 000384A were administered n = 42, 111, 144, 67, and 3 times, respectively, and were not associated with any deaths and were excluded. Thus, n = 35 batches were included for further analysis (Table S1).

Astra Zeneca

A total of 20 batches of ChAdOx1 nCoV-19 (Vaxzevria) were administered. (doses administered: mean n = 42,185; mean σ = 38,840; max n = 144,504; min n = 4,968; Table S1). All batches showed a peak of administration between February and July 2021. Batch ABY4055 showed a peak administration in October 2021, and, as only 10 doses were administered and no associated deaths occurred, it was omitted, leaving n = 19 batches for further analysis (Table S1).

Jansen

In total, 15 batches of Ad26.COV2.S (Jcovden) were administered, with peak usage from May to November 2021 (doses administered: mean n = 26,780; mean σ = 35,242; max n = 146,120; min n = 9,371; Table S1). Batch ACA4541 was administered only once and was omitted, leaving n = 14 batches for further analysis (Table S1).

Data Analysis

Data curation and statistical analyses were performed in IBM SPSS Statistics (Version 26) [21]. Charts and figures were produced in Microsoft Excel.

To assess heterogeneity in the data, hierarchical cluster analysis (HCA) was performed on 3-month ACM rates to determine the appropriate number of clusters, followed by non-hierarchical cluster analysis (N-HCA) to generate segments, and data for a generalized linear model (GLM) to test for cluster separation.

Statistical significance for ACM rates per batch was tested using a Z-test. Since the Z-test assumes a normal distribution, this test is appropriate when testing rates far from zero. For rates near zero, testing was performed using the Poisson distribution. Correlation between ACM and SAEs was performed using bivariate Pearson correlation combined with recursive removal of outliers outside 3 standard deviations (number of removed outliers/batches: Pfizer: 3, Spikevax: 4, AstraZeneca: 0, Jansen: 0). Since SAE data was not adjusted for age and sex, the correlation was calculated/reported using the crude ACM rates, although using the age and sex adjusted ACM showed similar results.

3. Results

Following data curation (Figure 1), we analyzed age- and sex-standardized mortality rates (deaths per 100,000) for the four vaccine products administered in the Czech Republic over the study period. Rates for Comirnaty varied significantly across batches (range = 0 to 1124.2, mean = 260.4, σ = 236.7; Table S1). Initial HCA identified 3-5 distinct clusters, along with several individual batches that did not cluster. A final solution with 3 clusters was chosen (Figure 2A), as it was the simplest option with significant separation between clusters using GLM (p < 0.0001) and also effectively captured the heterogeneity in the sample. Batches within the blue cluster (n = 8) accounted for 4.6% of all doses and were linked to 11.3% of all deaths (R2 = 0.99, β = 0.00898, 95% CI = 0.00799–0.00996; Table S1). Batches within the green cluster (n = 38) made up the majority of doses (80.3%) and the majority of deaths (79.6%): R2 = 0.99, β = 0.00353, 95% CI 0.00340–0.00367; Table S1), while batches in the yellow cluster (n = 20) accounted for 15.1% of doses and were associated with 9.1% of deaths (R2 = 0.98, β = 0.0023834, 95% CI 0.002219–0.002548; Table S1).

We previously described a temporal association between the severity of SAEs and different vaccine-batch clusters [17,18,20]. Therefore, we plotted batch number against mortality rate based on the peak month of administration and observed a temporal pattern where mortality rates and clustered batches aligned (Figure 2B). Batches in the blue cluster were administered early in the vaccination campaign, followed by most batches in the green cluster, and later batches in the yellow cluster, with some overlap between green and yellow during the later study period. Additionally, all batches except one in the blue cluster showed a highly significant positive deviation from the mean mortality rate (n = 954.3, σ = ±101.1; p ≤ 0.0364; Table S1). Notably, several batches in the yellow cluster showed a significant negative deviation from the mean mortality rate, including EX8680, which was specifically used for Pfizer employees (mean mortality rate = 0.0; p = 0.000; Table S1), and three pediatric vaccine batches (FN4072, FN4071, and FN6988; mean mortality rate = ±0.1949; p ≤ 0.000; Table S1).

For SPIKEVAX, the mortality rate also varied considerably (range = 0 to 2385.1, mean = 419.9, σ = 402.3; Table 1). Initial hierarchical clustering analysis (HCA) identified one large cluster (n = 33 batches) and one small cluster (n = 2 batches). Then, a non-HCA method with four clusters was used, resulting in two clusters with single batches each, which were then combined into a common blue cluster (n = 2). The remaining batches were shared between two additional clusters (green, n = 22, and yellow, n = 11), and the three clusters (Figure 2C) were tested for significant differences using GLM (p ≤ 0.0001). Batches within the blue cluster (n = 2 batches) accounted for 2.0% of all doses and were associated with 7.5% of all deaths (R2 = 0.97, β = 0.01406 (95% CI -0.01361–0.04173); Table S1). Batches within the green cluster (n = 18 batches) made up the majority of doses (70.7%) and related deaths (74.7%) (R2 = 0.98, β = 0.00431 (95% CI 0.00401–0.00462); Table S1), while batches within the yellow cluster (n = 15 batches) represented 27.4% of doses and were linked to 17.8% of deaths (R2 = 0.99, β = 0.00271 (95% CI 0.00256–0.00286); Table S1).

As observed with Comirnaty, SPIKEVAX batches that clustered together also exhibited a tem-poral pattern, with a clear separation between batches in the blue and green clusters and some overlap between batches in the green and yellow clusters (Fig. 2F). A significant positive deviation from the mean mortality rate was found for the two batches in the blue cluster, followed by a decrease over time.

For Vaxzevria, the mortality rate varied significantly (range = 324.6 to 1072.2, mean = 538.1, σ= 161.3; Table S1). HCA suggested solutions of 2 or 3 clusters. Non-HCA analysis produced a 3-cluster solution, each with significant separation, using GLM (p < 0.0001). Similar to Comirnaty and SPIKEVAX, heterogeneity in Vaxzevria’s mortality rate showed a comparable temporal pattern, with batches separated in the blue (n = 1) and green (n = 9) clusters and overlap between the green and yellow (n = 9) clusters. A significant positive deviation from the mean mortality rate was observed for the single batch in the blue cluster (ABV2856; Table S1).

Batches within the blue cluster (n = 1 batch) accounted for 2.3% of all doses and were associated with 4.7% of all deaths (R2 = 1, β = 0.01072 (95% CI not available (n = 1)); Table S1). Batches within the green cluster (n = 9 batches) made up approximately half of the doses 42.5% and related deaths 48.3% (R2 = 0.99, β = 0.00575 (95% CI 0.00340–0.00613); Table S1), while batches within the yellow cluster (n = 9 batches) represented 55.3% of doses and were linked to 47.0% of deaths (R2 = 0.99, β = 0.00457 (95% CI 0.00420–0.00494); Table S1).

For Jcovden, the mortality rate also varied significantly (range = 304.4 to 756.3, mean = 496.0, σ = 118.0; Table S1). HCA indicated a solution consisting of 2-3 clusters, and the final solution for non-HCA consisted of 2 clusters with a significant separation using GLM (p ≤ 0.0001). In contrast to Comirnaty, SPIKEVAX, and AZ1222, Jcovden batches within the same cluster did not show a clear temporal relationship with the mortality rate. A single batch (21C17-05) showed a significant positive deviation from the mean (rate = 765.3, p ≤ 0.0112).

Batches within the blue cluster (n = 8 batches) accounted for 75.0% of all doses and were associated with 81.2% of all deaths (R2 = 0.999, β = 0.0055 (95% CI -0.00537–0.00564); Table S1). Batches within the green cluster (n = 6 batches) made up the minority of doses (25.0%) and related deaths (18.8%) (R2 = 0.985, β = 0.00383 (95% CI 0.003295–0.004373); Table S1). We did not see a seasonal pattern in mortality rates for any of the products (Table S1).

To examine the relationship between mortality rates and reported SAEs, we utilized SAE data from the Czech Republic [20] to calculate crude SAE rates, as this data could not be adjusted for age or gender, and compared these with the crude mortality rates for the corresponding batch (Table S1). After removing outliers (see Method; Comirnaty n = 3 and SPIKEVAX n = 4 batches), a strong correlation was observed between SAE and ACM rates for Comirnaty (r = 0.82; Fig. 3A). SPIKEVAX (r = 0.69; Fig. 3B) and Vaxzevria (r = 0.82; Fig. 3B) also showed high correlation, but not for Jcovden (r = 0.27; Fig. 3D).

Figure 3. Correlation of batch-related mortality and suspected adverse event report data. Scatter plots showing the relationship between three-month all-cause mortality (ACM) and suspected adverse event (SAE) report data for (A) Comirnaty, (B) SPIKEVAX, (C) Vaxzevria, and (D) Jcovden. Each point represents a single vaccine batch, color-coded according to that presented in Fig. 2. Dashed trendlines are linear regression lines. The coefficient of determination (R²) is shown in the bottom right of each plot.

Figure 3. Correlation of batch-related mortality and suspected adverse event report data. Scatter plots showing the relationship between three-month all-cause mortality (ACM) and suspected adverse event (SAE) report data for (A) Comirnaty, (B) SPIKEVAX, (C) Vaxzevria, and (D) Jcovden. Each point represents a single vaccine batch, color-coded according to that presented in Fig. 2. Dashed trendlines are linear regression lines. The coefficient of determination (R²) is shown in the bottom right of each plot.

4. Discussion

We examined the relationship between specific COVID-19 vaccine batches administered over a 21-month period, from the start of the vaccination rollout through September 2022, and the temporal ACM data up to the end of that year, covering the entire population of the Czech Republic. Using descriptive statistics and cluster analysis on this comprehensive, high-quality national dataset (which includes approximately 19 million doses given to 11 million people), we identified detailed three-month ACM rates linked to specific batches. Building on our previous analyses of batch-related differences in severe adverse event (SAE) reports [20], we identified batches with significantly higher or lower ACM rates than other batches of the same product. This relationship persisted after adjusting for age and sex and demonstrated a clear temporal trend. For three of the four products studied (Comirnaty, SPIKEVAX, and Vaxzevria), notably higher ACM was associated with batches administered early in the vaccination campaign. Similarly, strong correlations with SAE results were found for these three products, but not for Jcovden.

Remarkable platform-dependent differences in the safety profiles of COVID-19 vaccines have been reported [22]. Whilst our findings do not establish causality at the individual level between COVID-19 vaccine administration and death, they highlight important questions for the safety monitoring of COVID-19 vaccines, especially those using the modified RNA (modRNA) platform. The strong correlation between batch-associated ACM and batch-associated SAE rates suggests a genuine safety signal rather than residual confounding. This correlation is significant because the ACM and SAE data are independent; SAEs come from a spontaneous reporting system susceptible to under-reporting and bias [9], while ACM is a definitive endpoint from reliable national mortality records. The convergence of signals from these separate sources greatly strengthens the evidence for a non-random phenomenon. The bias present in results based on SAE data is considerably weakened, if not eliminated, by the findings from the ACM data, and the ACM results are further supported by the SAE data regarding potential causality [23].

The results of this study, when viewed in the context of previous SAE analyses, support earlier hypotheses about SAE-related safety signals. Our findings align with studies reporting batch-associated SAE patterns in the Czech Republic [20], Denmark [17,18], Sweden [17], and the USA [19]. They differ from the findings of Hviid et al., who later reported no significant rise in ACM for Danish SAE-associated batches [24]. However, the methodological limitations of this study, which might have diluted a potential ACM-associated signal concentrated in high-risk, early-vaccinated individuals as well as the uncertain relationship between SAEs, short-term hospitalization, and long-term prognostic consequences [23] preclude drawing a definitive conclusion. Our data confirm not only a temporal decline in batch-associated ACM but also a correlation between ACM and SAE patterns within the same population. Both modRNA products showed several batches (EL1484, EJ6796, and EX8680 for Comirnaty and 300493, 300496, and 3001531 for SPIKEVAX; Fig. 3A and B, respectively) that deviated substantially and were removed as outliers in the correlation analyses. These batches all showed high SAE rates relative to ACM rates and were administered early on, except for Comirnaty EX8680 (a batch reserved for Pfizer employees), which had peak utilization in August 2021. We were unable to identify any defining features of these outlier batches that might explain their lack of correlation with closely matched batches in series or peak usage. However, the strong correlations between ACM and SAE rates for modRNA vaccine batches reinforce the possibility of significant batch-to-batch product differences. Potential factors include modRNA integrity [25,26], Spike protein expression dynamics [27,28], residual DNA [29], and the reformulation of the buffer component of the drug product for Comirnaty [30].

The main strength of this study is the analysis of a comprehensive population-level ACM dataset that complements spontaneous SAE data. Key limitations must be acknowledged. First, although we adjusted for age and sex, population-level data cannot account for all individual-level confounders, such as specific comorbidities. Second, we only assessed mortality within a three-month period after vaccination; longer-term effects were not evaluated due to the limited duration of mortality data. Lastly, while ACM exhibits seasonal trends, we did not identify a pattern linking batch-related ACM to seasonality, indicating that it was not a primary factor in our results. It is also important to highlight that vaccine-induced mortality is likely only a small part of total ACM. However, even minor, systematic differences in ACM across batches deserve careful examination as possible safety signals.

5. Conclusions

We conclude that clusters of COVID-19 vaccine batches administered in the Czech Republic were associated with significantly different rates of short-term ACM. For Comirnaty, SPIKEVAX, and Vaxzevria, batches administered earlier were linked to higher ACM. The correlation between ACM rates and reported SAE rates for these batches strengthens the evidence for a batch-dependent safety signal. These signals, unreported during the initial rollout, warrant further investigation through detailed medical record reviews, analysis of pathological specimens where available, and quality-control re-evaluation of retained batch samples, especially concerning mRNA integrity and residual DNA content.