We collected fecal samples from fourteen healthy adults, who were traced since January to May 2021 at Beijing Chaoyang Hospital. Each volunteer was examined at the baseline before subjecting to injection of BBIBP-CorV vaccines (Before group, n = 14), and after they finish the completely two doses intramuscular injection of vaccines (After group, n = 14). The samples were harvested twice from each enrolled participant at two weeks prior to injection and two weeks after receiving the vaccines, when the vaccines was considered to be effectively activate the host immune system. Individuals with or with a history of cancer, heart failure, renal failure, stroke, peripheral artery disease, and chronic inflammation disorders were excluded. Moreover, pregnant or breast-feeding female subjects were excluded. Other exclusion criteria were subjects with previous SARS-CoV-2 exposure or infection, and those who received statin, aspirin, insulin, metformin, antibiotics, or probiotics etc. within the last two months. In present study, a self-controlled protocol was designed to eliminate some confounding factors including Body Mass Index (BMI), smoking status, alcohol consumption, dietary habits, and physical exercises etc. And food habits and lifestyle behaviors did not change significantly during the short survey periods. This study was approved by the Medical Ethics Committee from Beijing Chaoyang Hospital and carried out in accordance with the Helsinki declaration. Written informed consent was obtained from all the participants.

The middle section of fresh fecal samples wascollected from participants, and immediately stored in a freezer at –80 °C within half hour. Fecal microbial DNA was extracted with a TIANamp Stool DNA Kit (TIANGEN Biotech, Beijing, China) according to the manufacturer’s instructions. Isolated DNA samples were quantified using a NanoDrop® ND-2000c UV-vis spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA), and DNA integrity was estimated with agarose gel electrophoresis. Fecal microbial DNA samples were diluted to 1ng/$\mu$L and stored at –20 °C for further examination. The V3–V4 region of bacterial 16S rRNA gene was amplified through Polymerase chain reaction (PCR) with universal primer set 338F (5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’). Following library preparation and generation, then we assessed the library quality on the Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system. The libraries were subsequently sequenced on IlluminaHiSeq instrument 2500 and 250 bp paired-end reads were obtained.

According to the information of primers and unique barcode regions, sequencing data for each sample was assigned. By removing the barcode and primer, all pair-end sequences were merged with FLASH (Version 1.2.17, Baltimore, MD, USA) to obtain the raw tags. Strict quality-controlled processes including truncating tags with continuous low-quality bases and filtering out tags with 75% high quality bases were carried out with QIIME (Version 1.9.1, Denver, CO, USA). Effective clean tags were acquired when chimeras were further removed using UCHIME Algorithm (Version 7.0.1001, Tiburon, CA, USA). Utilizing UPARSE (Version 7.0.1001, Tiburon, CA, USA), effective tags from all samples were constructed into OTUs, with a cutoff for aligning identity at 0.97. The phylogenetic affiliations of clustered OTU sequences were annotated through RDP classier (Version 2.6, East Lansing, MI, USA) with a threshold at 0.8. Sequences failed to be assigned into specific taxonomy were considered as unclassified.

Rarefaction analysis and OTUs accumulation curves by randomly including certain number of individuals and counting the observed OTUs were applied to evaluate the rationality of sample size. And rarefaction curves by gradually increasing the sequencing depth were conducted to examine the rationality of sequencing quantity. Bacterial alpha-diversity based on OTU profiles was presented by the chao1 richness, shannon diversity, pielou evenness, ACE index, simpson index and invsimpson (inverse simpson) index using the R software (Version 3.3.3, Auckland, New Zealand) program package “vegan”. Significant differences of alpha-diversity indices in each individual before and after vaccination were tested by Wilcoxon matched-pairs signed rank test.

Subsequently, the relative abundance of OTUs and Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology (KOs) profiles was applied to generate the bray-curtis distance through QIIME (Version 1.9.1, Denver, CO, USA) and assess the microbial beta-diversity. Nonmetric dimensional scaling (NMDS), principal-component analysis (PCA) and principal coordinate analysis (PCoA) were conducted with the vegan package, ade4 package and WGCNA package in the R software (Version 3.3.3, Auckland, New Zealand), respectively. Significant differences between groups were tested by Analysis of similarities (Anosim).

According to the KEGG database, the potential functional capabilities of gut microbiota were predicted and identified by phylogenetic investigation of communities by reconstruction of unobserved states (PICRUSt) Version 1.0.0 (Halifax, NS, Canada). KEGG ontology and pathway profiles were generatedbased on the OTU tree and the gene information in Greengene database. Predicted functional genes were categorized into KO, and assign into KEGG pathway.

Analysis for the demographic and clinical data was conducted with the SPSS 20.0 statistical package (IBM, IL, USA), using two-tailed student’s t-test or chi-square test, as appropriate. Significant differences of alpha-diversity indices were tested by Wilcoxon matched-pairs signed rank test, and of beta-diversity were determined by Anosim. The Differential abundance of taxa and KOs for paired samples (two samples before and after treatment from the same patients) was determined by the Wilcoxon matched-pairs signed rank test. p 0.05 was set as statistically significant.

The characteristics of 14 healthy adults included in this study were summarized in Supplementary Table 1. The mean age of the study population is 38.4 years, and 71.4% were men. Meanwhile, the mean systolic blood pressure (SBP) and diastolic blood pressure (DBP) of these participants are 120.7 and 73.9 mmHg, respectively.

Twenty-eight fecal samples from all the recruited subjects were sequenced on an IlluminaHiSeq sequencer. A total of 1,970,533 qualified tags were obtained from 2,139,298 raw reads, and for each sample, 65,364 effective tags were acquired. A total of 85.41% of all qualified tags were clustered into qualified OTUs. OTUs those detected in less than 2 samples were discarded for their low frequencies in samples. Ultimately, qualified operational taxonomy units (OTUs) were clustered for downstream analysis. The details for sequencing characteristics are shown in Supplementary Table 2 of the supporting information.

Rarefaction analysis was carried out and OTUs accumulation box plots were performed to describe increased diversity and richness following expanding the sample size. The rate at which novel OTUs emerge by continuous sampling was estimated based on the box plots in Supplementary Fig. 1A. It was observed to be tending flat, but not sharply increased, indicating that the OTUs would not significantly augment with larger sample size. On the other hand, by randomly extracting a certain amount of sequencing data from the samples, the number of OTUs they represent was depicted with rarefaction curves (Supplementary Fig. 1B). It reflected that the sequencing depth approached stable, and elevated sequencing data amount will generate only a few new OTUs.

The richness of observed OTUs was comparable in participants before and after vaccines treatment (Fig. 1A). In addition, parameters related to alpha-diversity including chao1 richness, shannon diversity, pielou evenness, ACE index, simpson index and invsimpson index were accessed (Fig. 1B–G). We found that both chao1 and ACE were similar between the two groups prior to and post receiving vaccines. Interestingly, the results showed that intestinal microbial alpha diversity in shannon diversity, pielou evenness, simpson index as well as invsimpson index was remarkably reduced in recipients after vaccination than before vaccination, indicating that bacterial homeostasis was affected by vaccines. The beta diversity results based on NMDS, PCA and PCoA plots showed that there was no significant distinction in fecal microbial community distribution between groups (Supplementary Fig. 2). Consequently, it was suggested that BBIBP-CorV vaccines impacted gut microbiome by eliciting depressed microbial alpha-diversity.

Fig. 1.

Alterations of the bacterial diversity and richness indexes of the fecal microbiota in exactly individuals before and after vaccination. (A) The number of observed OTUs in gut microbiome of healthy volunteers was examined before and after vaccination. (B–G) Box plots depicted changes of alpha diversity by BBIBP-CorV vaccines treatment. B, Chao1; C, shannon; D, pielou; E, ACE; F, simpson; G, Invsimpson. The horizontal lines within the box plots represent median values; the upper and lower ranges of the box is 75% and 25% quartiles, and each dot represents a sample. p values were obtained from Wilcoxon matched-pairs signed rank test.

Having established that microbial alpha-diversity differed significantly by vaccines taxonomically, we next examined phylogenetic profiles in the gut microbiome and investigated the taxonomic composition alterations associated with vaccines treatment. Relative abundances of bacterial phylum and primary genus in the fecal samples from subjects who underwent vaccination were revealed. Firmicutes, Bacteroidetes, Verrucomicrobia and Proteobacteria were the most dominant phylum identified. Before vaccination, the proportions of Bacteroidetes and Proteobacteria were higher in their relative bacterial abundances. By contrast, the proportions of Firmicutes were more enriched after vaccination (Fig. 2A,B). In addition, at the genus level, Faecalibacterium, Bacteroides, Blautia, Prevotella, Bifidobacterium, Lachnospira and Roseburia were the most abundant entities in the gut microbiota. After vaccination, the subjects exhibited enriched Faecalibacterium, Bifidobacterium, Lachnospira and Roseburia in terms of relative bacterial abundance, but less abundance of Bacteroides and Blautia (Fig. 2C,D).

Fig. 2.

Characterization of the phylogenetic profiles in the gut microbiome of recipients before and after vaccination. (A) Average composition of the bacterial community at the phylum level was compared before and after vaccines treatment. (B) Relative abundance of the top 10 bacterial phylum in the fecal samples from each participant who underwent vaccination. (C) Average abundance of the top 10 genera was compared before and after receiving vaccines. (D) Relative composition for the top 10 bacterial genus in the fecal samples of each participant before and after vaccination.

Further, we performed differential expression analysis of bacteria abundance between Before and After vaccines inoculation groups through Wilcoxon matched-pairs signed rank test. At the phylum level, Verrucomicrobia was observed to be decreased with marginal significant after vaccination (Fig. 3A). We found that the classes of Coriobacteriia, Verrucomicrobiae and order of Coriobacteriales, Verrucomicrobiales were reduced marginal significantly, whereas Betaproteobacteria class and Burkholderiales order elevated after receiving the vaccines (Fig. 3B,C). In addition, decreased abundance of the family Verrucomicrobiaceae, Coriobacteriaceae, Peptostreptococcaceae, Mogibacteriaceae, but enhanced Alcaligenaceae was detected following vaccines treatment (Fig. 3D). The Actinomycetales order and two families including Oxalobacteraceae and Actinomycetaceae were prominently decreased after vaccination (Fig. 3C,D). Among them, 12 taxa were depleted and 3 were enhanced by vaccines (Fig. 3E). At the genus level, Collinsella and Akkermansia genera were reduced, while Roseburia and Sutterella genera were overgrowth after vaccines treatment (marginal significant) (Fig. 4A). It was intriguingly that Ruminococcus and Actinomyces were strikingly deficient after vaccination, with fold change at 0.52 and 0.31, respectively (Fig. 4A,B). And Faecalibacterium was dramatically augmented with 1.58-fold changes in recipients of vaccines (Fig. 4A,B).

Fig. 3.

Comparison of the relative abundance in fecal taxa between samples obtained from recipients before and after vaccination. (A) Gut phylum Verrucomicrobia displayed marginal statistical difference after vaccination. (B) Classes with marginal statistical difference following vaccination. (C,D) Orders and families shifted in gut microbiome after receiving vaccines treatment. The horizontal lines within the box plots were medians; the upper and lower ranges of the box represent 75% and 25% quartiles, and dot represents each sample. p values were obtained from Wilcoxon matched-pairs signed rank test. (E) The enrichment and depletion of gut taxa, and fold changes for their alterations following vaccination were illustrated in heat map.

Fig. 4.

Comparison for the abundance of genera in fecal taxa between samples obtained from recipients before and after vaccination. (A) Intestinal bacteria at genus level were identified to be altered following vaccines treatment. The horizontal lines within the box plots were medians; the upper and lower ranges of the box represent 75% and 25% quartiles, and dot represents each sample. p values were obtained from Wilcoxon matched-pairs signed rank test. (B) The enrichment and depletion of gut genera, and corresponding fold changes after vaccination.

To predict the microbial community function profiles, all the observed OTUs were aligned into the PICRUSt built-in reference database. Therefore, a total of 5661 KOs were parsed and mapped into 276 KEGG pathways. Functionally, we again failed to identify obviously distinct clusters between individuals before and after vaccination based on NMDS, PCA and PCoA plots of Bray-Curtis pairwise dissimilarities (Supplementary Fig. 3). Nevertheless, the top 30 most significantly affected KOs and a total of 17 prominently altered KEGG pathways by vaccination were presented in the heat map with abundance (Fig. 5A,B). There were 6 KOs suppressed and 24 KOs activated by vaccines, including branched-chain amino acid transport system and glutamate synthase (NADPH) small chain (Fig. 5A). Among the 17 KEGG pathways identified to be with significant differences between Before and After group, 10 functions, such as amino acid metabolism, lipid biosynthesis proteins and Steroid biosynthesis, were remarkably increased, while 7 functions, such as renin-angiotensin system and biotin metabolism, were remarkably decreased following vaccines treatment.

Fig. 5.

Functional capability of gut microbiome varied dramatically in individuals after receiving vaccines. (A) Heat map showing the top 30 most significantly fluctuated KEGG orthologs (KOs) by vaccination, with 6 KOs reduced and 24 KOs enhanced after vaccination. (B) There were seventeen KEGG pathways at level 3 identified to be prominently affected by BBIBP-CorV vaccines. The abundance profiles are transformed into Z scores by subtracting the average and dividing the standard deviation of all the samples. Z score is negative in blue when the row abundance is lower than the mean, and positive in orange when the row abundance is higher than the mean.