The composition of gut microbiome is not altered by vaccination against COVID-19

To investigate the impact of the COVID-19 vaccines on the gut microbiome, shotgun metagenomic sequencing was performed on a total of 239 fecal samples from 59 patients from our three cohorts (43 healthy control, 160 cancer, and 36 primary immunodeficient patient samples) (Fig S1). Samples were taken over the course of three vaccine doses, at one of three vaccine timepoints: pre-dose (before vaccination), acute (2–3 d after vaccination) or late (16–28 d after vaccination) for each vaccine dose (Fig 1A). Pre-dose sampling provides a baseline assessment of the gut microbiome before vaccination, whereas acute samples provide the opportunity to assess the effect of vaccination on the gut microbiome at the height of the initial inflammatory response, and late samples allow a determination of any resolution or maintenance of acute effects.

When first observing the α-diversities of all samples taken from each of our cohorts, we noted significant differences between the cohorts (Fig S2A), however when assessing samples taken at different vaccine timepoints from within each cohort, we did not observe any significantly differences (Fig 1B). In agreement with this, when assessing samples from the same patient at different vaccine timepoints using a paired sample analysis (tracking individual patient samples across multiple timepoints, e.g., pre-dose and acute), we also did not see any significant differences in the α-diversities of our patient samples (Fig S2B and Table S1). This indicates that the COVID-19 vaccine is not affecting the diversity of the gut microbiome, despite the distinct microbial diversity between the cohorts and individuals.

We next applied principal component analysis to visualize the β-diversity of our microbiome composition data using the abundance of all detected operational taxonomic units (Fig 1C). The principal components (PCs) describe the largest variation components in the dataset, representing shifts in microbiome composition and potentially reflecting to the abundance changes of bacterial species between the samples. The first five principal components were responsible for the 2.9%, 2.4%, 2.1%, 2.0%, and 1.9% of variation in the data, respectively, and were further analysed using mixed-effect linear models with multiple input variables from our available metadata. In our linear models, we asked whether the vaccination timepoint of the samples in each cohort could improve the explained variance of the PCs when compared with a baseline model describing the explained variance using patient samples as the grouping variable. We found that there was no significant improvement on the baseline model (Fig S2C). Furthermore, a PERMANOVA analysis did not find the covariate of the timepoints from which the samples were taken to be significantly affecting the microbial composition (Table S2). These analyses together support that the variance we see in our samples is not a signature of the COVID-19 vaccines, rather of the patients presenting with different microbiome compositions themselves.

Next, we asked whether the COVID-19 vaccines induce any changes in the phylum-level composition of the gut microbiome and profiled the relative abundance of phyla across all samples (Fig 1D), observing variation in our patient samples. Moreover, when comparing the top six most prevalent phyla, no significant differences were observed between vaccine timepoints (Fig 1E) despite significant differences in these phyla between cohorts and individuals (Fig S2D). When using a paired analysis approach, we observed only marginal, small effect size, differences (after multiple testing correction) in seven (out of 231) timepoint comparisons in the cancer cohort (Fig S2E) (Table S3). Of those seven, however, only one had a log₂FC greater than 1 in either direction (+1 or −1) (Fig S2F), indicative of notable change in median relative abundance of Lentisphaerae.

This demonstrates that the COVID-19 vaccines do not appreciably alter the composition of the gut microbiome, irrespective of the unique compositions found in our cohort samples.

COVID-19 vaccination does not induce species level changes in the gut microbiome

We next sought to analyse differentially abundant microbial species between vaccine timepoints, that is, pre-dose, acute, and late, using DESeq2 (15). All cohorts were analysed independently for the abundance changes in samples taken at each timepoint with the most differentially abundant species in a representative heatmap. For the cancer cohort, when assessing the abundance of these top differential-responding bacterial species between samples taken pre-dose and acutely, unsupervised clustering does not demonstrate evident grouping of timepoints (Fig 2A). Among all the species, only two were significantly increased in acute samples compared with the pre-dose samples, Klebsiella pneumoniae and Butyrivibrio crossotus found in 15 (P = 1.01 × 10⁻²⁴) and seven samples (P = 7.63 × 10⁻¹²) out of the 97 cancer patient samples, respectively (Fig 2B). The former is only representative in a quarter of the cohort (11 patients), only melanoma patients, and within those has an average relative abundance of 0.7% (Fig 2C); the latter in two renal cancer patients, representing on average 3% of the relative abundance (Fig 2D). Considering ∼2,500 species are represented across all patient samples, change in two low-abundant and sparsely represented species signifies negligible changes. Similar findings were seen for our other two cohorts, healthy controls and primary immunodeficient patients (Fig S3A–D, respectively). When performing paired sample differential abundance analysis in all three cohorts using DESeq2, we find that no significantly altered species between pre-dose and acute, or pre-dose and late sample in the healthy control and cancer cohorts. Samples from only one primary immunodeficient patient showed a significant reduction in Enterobacter sp. in a late sample compared with pre-dose (Fig S3E). Thus, our findings demonstrate that on a species level there is no unified, biologically relevant change in abundance of microbial species induced by the COVID-19 vaccines.

As the differential abundance analysis considers the change in abundance of all species irrespective of their relative abundance within each sample, we were curious whether there were any noticeable changes in the most abundant species found within each patient cohort that could be attributed to the vaccine timepoints. There was no significant difference in any of the most abundant 15 species in samples taken at any of the three vaccine timepoints (Fig 2E), representing, on average, 47% of the relative abundance of the species within patient samples in the cancer cohort, 50% within the healthy controls, and 53% within the primary immunodeficient patient cohort. This indicates that we see no effect of the vaccine on the species occupying the highest proportion of the microbial niche.

There was considerable concern both at the time of the initial vaccine programme, and to this day, on the safety of COVID-19 vaccines. Given that, we sought to observe if there was any differential outgrowth of bacterial species that have been associated to various immune-related diseases, including gastric cancer and autoimmunity, and metabolic and neurological diseases. Although our study cannot address long-term outgrowth of bacterial associated to these diseases, we are able to highlight if there is any temporal, acute changes in these species which is still of physiological relevance. Faecalibacterium prausnitzii, which is reported to be reduced in both gastric cancers, autoimmunity, and Crohn’s disease (16), showed no significant abundance changes at the vaccine timepoints in our cohorts (Fig 2F). We also found no presence of Helicobacter pylori, which is strongly associated with the initiation and development of gastric cancers (17). Akkermansia muciniphila, found to be increased in obesity (18) and correlated with response rates to immune checkpoint blockade therapies in various malignancies (19), was not altered by vaccination in all three cohorts (Fig 2G). In studies of Alzheimer’s, Escherichia coli has been demonstrated to promote neurodegeneration (20); in our samples we did not see significant difference induced by the vaccine in any of our cohorts (Fig 2H). This supports that the COVID-19 vaccine does not promote the change in abundance of microbes that are associated with various immune-related diseases within any of the three cohorts and is indicative of no greater risk of the aforementioned diseases as a result of COVID-19 vaccination.

Gut microbiome diversity is not correlated with the magnitude of the response to the COVID-19 vaccines

Within the current literature, a few studies have reported changes in the gut microbiome that correlate with vaccine efficacy (12, 13, 14), so we sought to determine whether the gut microbiome composition was related to the magnitude of the COVID-19 vaccine response. We performed an assessment of vaccine efficacy using a live-virus neutralisation assay, as a predictive measure of vaccine protection (21), to assess whether vaccine efficacy had correlation with microbial diversity.

When taking the neutralising capacity of patient serum at both the second dose (Fig 3A) and third dose (Fig 3B), we asked whether the Shannon diversity of gut microbiome at different vaccine timepoints was affected by or correlated with neutralisation. We did not see any correlation between diversity and vaccine efficacy in any of our patient cohorts; the same is true for the quantity of anti-spike IgG antibodies (Fig S4). This indicates that the gut microbiome diversity was not correlated with the magnitude of the immune response in our patient cohorts, thus suggesting that improved efficacy of the vaccine does not come at a cost of microbial disturbance.

The gut microbiome functional capacity was not affected by COVID-19 vaccines

Having investigated the composition and relative abundance of the microbial species that constitute the gut microbiome, we next sought to investigate whether the functional capacity of the microbial species was altered by the COVID-19 vaccines. Using the EggNOG database (22), we assigned functional annotations to the sequenced metagenomes. The highest level of functional annotation depicts three functional groups, cellular processes, and signalling, information storage and processing, and metabolism. In these, we did not see any significant differences between the vaccine timepoints within our cohorts (Fig 4A); similar to taxonomic data presented earlier, when assessing cohort samples separately, there are significant changes (Fig S5A).

We next observed the abundance of the 22 defined functional groups in the next functional annotation level down in the separate vaccine timepoints within each patient cohort (Fig 4B). Representative graphs of the most abundant functional annotations within each of the highest three level functional levels remained unchanged after vaccination within our patient cohorts (Fig 4C), the same is true for the remaining 19 (Fig S5B).

At the lowest functional annotation level, we interrogated the abundance of cluster of orthologous genes at separate vaccine timepoints within our cohorts. Remarkably, no cluster of orthologous genes out of a possible 2,142 genes presented in our patient samples, was significantly altered (when corrected for multiple testing) as a result of the COVID-19 vaccines in any of the three cohorts. This demonstrates that the functional annotations of the gut microbiome are not altered by the administration of the COVID-19 vaccines.

Study recruitment and ethics

Participants volunteered for an observational study and were enrolled to one of 3 cohorts, healthy controls, cancer patients which presented with either melanoma or renal malignancies, or primary immunodeficient patients with defined mutations in key immunoregulatory genes, and patients with clinical presentation aligning to that of defined primary immunodeficiency, including antibody deficiency (Table 1). Almost all participants in each cohort received BNT162b2 Pfizer vaccine, aside from one patient at second dose receiving AstraZeneca, and six patients receiving Moderna at the third (Table 1). Patients were excluded if presenting with positive COVID-19 serology or if presenting in hospital with clinal symptoms/features related to their disease which may influence the physiological response to the COVID-19 vaccination. The research was conducted in accordance with the principles of Good Clinical Practice and following approved protocols of the NIHR National Bioresource. Samples were collected at NIHR Cambridge Clinical Research Facility, Addenbrookes Hospital, Cambridge, UK. Samples were provided with the written informed consent of all study participants under the NIHR National BioResource - Research Tissue Bank (NBR-RTB) Ethics (Research Ethics Committee [REC]:17/EE/0025). All participants were consented under the East of England Cambridge South national consent (REC:13/EE/0325).

Participant sampling

Participant samples were anonymised by clinical staff before sample delivery to the research laboratory. Peripheral blood and fecal samples were collected longitudinally over the course of up to three doses of the vaccines against COVID-19 (sample coverage varied across doses [Fig S1]). Samples were taken from 9^th February, 2021, to 22^nd February, 2022. PBMCs were extracted from blood samples using density gradient centrifugation, stored temporarily at −80°C, before being transferred to long-term storage in liquid nitrogen. Serum was isolated from peripheral blood via centrifugation and stored at −80°C until required. Fecal samples were collected at three timepoints around each vaccine dose: pre-dose (94% of samples taken within 3 d before vaccination, the remaining three samples taken up to 14 d prior), acute (day 2 or 3 after vaccination) or late (day 16–28 after vaccination). Fecal matter was collected in OMNIgene•GUT kits (DNA Genotek) whereby samples are stored in a stabilizing, inactivating solution. Samples were transported to the laboratory and homogenized upon arrival before being stored at −80°C until required.

Fecal DNA extraction and sequencing

DNA was extracted from fecal samples using QIAamp PowerFecal Pro DNA kits (QIAGEN). Samples were thawed and ∼250 mg of fecal sample was lysed via bead beating. According to the kit protocol, the sample was then cleaned of non-DNA organic and inorganic material, and then washed using ethanol. DNA was eluted into 10 mM Tris and quantified using the Qubit 1X dsDNA High-Sensitivity (HS) Assay Kit (Thermo Fisher Scientific) using the Qubit fluorometer (Thermo Fisher Scientific). DNA at a concentration of 10 mg/µl was sent for sequencing. Samples were sent in bulk for sequencing blinded to sample timepoints of each patient. Shotgun metagenomic sequencing was performed with Illumina NextSeq 2000 sequencing platform using paired-end reads of 150 bp in length, resulting in a median of ∼28 million reads per sample.

Shotgun metagenomic analysis

Raw sequencing data were preprocessed with PRINSEQ++ (39 Preprint) (v.1.2.4) in paired read mode, quality trimming to a minimal Phred score of 30 in a window of 15 bases and removing reads of less than 75 bp length after trimming. In addition, host contamination was removed by mapping against the GRCh38 reference human genome using Bowtie2 (40) (v.2.4.5) and removing any mapped reads from the dataset.

Raw, trimmed, and filtered reads were checked for quality using FastQC (v.0.11.9) (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). From the remaining read pairs, taxonomic profiling was determined using mOTUs3 (41) (v.3.0.1) profiler. For functional profiling, the remaining read pairs after filtering were assembled using metaSPAdes (42) (v.3.15.4) with a k-mer size of 55. The resulting scaffolds were filtered for at least 200-bp length and weighted by their average coverage of the filtered reads to enable quantitative analysis. The remaining scaffolds were aligned to the EggNOG database (22) (downloaded on 2022/04/08) using DIAMOND (43) (v.2.0.13). Microbiome analysis was performed in R (v.4.2.2) using relative abundances, within the phyloseq (44) (v.1.42.0), microbiome (45) (v.1.20.0) and vegan (46) (v.2.6-4) packages. Differential abundance analysis was performed using DESeq2 (15) (v.1.38.3).

Serological assessment of immune response to COVID-19 vaccines

Serum samples were thawed, heat-inactivated at 56°C for 30 min. Measurements for the dilution of serum that reduces viral activity by 50% (NT₅₀) against WT SARS-CoV-2 were as previously reported (47, 48). For anti–SARS-CoV-2–specific IgG antibodies, we used a previously described method (49, 50), in which Luminex bead sets are covalently coupled to the recombinant SARS-CoV-2 proteins nucleocapsid protein, spike, and receptor-binding domain to quantify antibody levels.

Statistical analysis

Wilcoxon tests with multiple testing corrections using the Bonferroni test were deployed throughout, using pairwise comparison where appropriate. Mixed-effect linear modelling was performed using lmer4 (51) in R. Correlation was determined using Spearman’s Rank coefficient.