Characteristics of participants and samples

A total of 71 participants who had completed two doses of primary immunization with the inactivated vaccine for more than 6 months were recruited to receive a third dose of either the BBIBP-CorV or ZF2001 vaccine. All participants were divided into four groups according to age and booster dose, including group A (boosted by BBIBP-CorV, < 70 years old), group B (boosted by BBIBP-CorV, ≥ 70 years old), group C (boosted by ZF2001, < 70 years old), and group D (boosted by ZF2001, ≥ 70 years old) (Fig. 1). At days 7 and 28 after vaccination, we collected a total of 131 serum and 57 PBMC samples. The detailed characteristics of the participants are shown in Table S1 and Table S2.

Fig. 1

The details of this study design, including age characteristics of participants, booster vaccine type, and samples collection time in the 4 booster groups

Antibody response

First, we estimated the titers of binding antibodies from sera collected at 7 and 28 days after vaccination of the participants. After 7 days of homologous boost with BBIBP-CorV, the mean value of IgG antibody in the younger group was 221.09 BAU/ml, and that in the elderly group was 104.35 BAU/ml. The elderly group was significantly lower than the younger group (P = 0.0282). A similar phenomenon was observed in different age groups heterologously boosted by ZF2001, with 366.63 BAU/ml in the younger group and 89.24 BAU/ml in the elder group, and the difference was statistically significant (P = 0.0014) (Fig. 2a). After 28 days of the third booster dose, the mean IgG antibody levels in the four groups increased to 610.6, 733.9, 1631, and 1617 BAU/ml, respectively. Remarkably, there was no age-associated difference in antibody responses(P = 0.8237, P = 0.0644) (Fig. 2b).

Fig. 2

The level of IgG antibodies against SARS-CoV-2 of serum samples 7 (a) and 28 (b) days after booster vaccination. P values were computed by unpaired t-tests with log-transformation using GraphPad Prism 8.0

Differential expression analysis of transcriptome profiles of participants with the third dose of COVID-19 booster vaccine

To determine the effect of aging on the transcriptional program of the third dose of COVID-19 vaccine, we performed a transcriptome systematic scan of 57 PBMC samples from different groups. After 7 days of homologous boost with BBIBP-CorV, the number of up-regulated and down-regulated DEGs in the elder group compared with the younger group was 213 and 65, respectively(Fig. 3a). After 7 days of ZF2001 heterologous boost, 136 and 83 genes were up-regulated and down-regulated in the elder group compared with the younger group, respectively(Fig. 3b). The up-regulated DEGs of the elder group after 28 days with the BBIBP-CorV booster were 100, and the down-regulated DEGs were 134(Fig. 3c). There were 131 up-regulated DEGs and 84 down-regulated DEGs in the elder group after 28 days of the ZF2001 booster(Fig. 3d). We used the heatmap to show the expression of the top 20 up-regulated and down-regulated differentially expressed genes in each sample(Fig. 4a-d).

Fig. 3

Volcano plot of global gene expression changes in the elderly group induced by the third dose of BBIBP-CorV (a) and ZF2001 (b) at 7 days. Volcano plot of global gene expression changes in the elderly group induced by the third dose of BBIBP-CorV (c) and ZF2001 (d) at 28 days. Red: up-regulated DEGs in the elderly group compared with the younger group; Blue: down-regulated DEGs in the elderly group compared with the younger group; Grey: non-DEGs. Criteria: |Fold Change|= 2^logFC_cutoff is marked out by two black dotted lines. P-value = 0.05 is marked out by a horizontal black line y =  − Log10.(0.05)

Fig. 4

a The expression profile of the top 20 up-regulated and down-regulated DEGs in the elderly group compared with the younger group boosted by BBIBP-CorV after 7 days. b The expression profile of the top 20 up-regulated and down-regulated DEGs in the elderly group compared with the younger group boosted by ZF2001 after 7 days. c The expression profile of the top 20 up-regulated and down-regulated DEGs in the elderly group compared with the younger group boosted by BBIBP-CorV after 28 days. d The expression profile of the top 20 up-regulated and down-regulated DEGs in the elderly group compared with the younger group boosted by ZF2001 after 28 days. Each row represents mRNA and each column represents a sample. Red indicates higher expression and blue indicates low expression

Impairment of the immune transcriptional response to the vaccine in aging on Day 7

In order to analyze the interaction relationships between DEGs, the STRING database was used to create a PPI network. Using the MCODE plug-in of Cytoscape software, we extracted 4 subnetworks from the PPI network at the nodes and edges of all DEGs, including the characteristic induced up-regulated and depleted down-regulated DEG sets identified from the third dose of BBIBP-CorV(Figure S1a, b) and ZF2001(Figure S1c, d) in elderly individuals. These characteristic gene subnetworks contained 21, 25, 34, and 42 DEGs, respectively.

For the purpose of clarifying the functional pathways of the signature DEGs mediated by the third booster vaccine in the elderly group, we conducted GO-term analysis, focusing on the biological process entries of these DEGs. First, we showed the top 20 GO-terms enriched by BBIBP-CorV and ZF2001 in the elderly group in the GO dot-plot (Fig. 5a, b). Coincidentally, DEGs were found to be involved in multiple innate immune pathways, such as lymphocyte proliferation, mononuclear cell proliferation, and leukocyte proliferation. The connection of GO functions with specific DEGs was shown in Figures S1e, f. The chordal diagram was used to further demonstrate the degree of differential expression changes of genes involved in each functional project (Figs. 5c, d). We found that most of the DEGs involved in immune-related functional items were down-regulated, implying that these immune-related pathways were impaired to some extent in the elderly group. This may be a potential reason of early humoral immune deficiency in elderly people.

Fig. 5

a The GO analysis on DEGs induced by the third dose of BBIBP-CorV in elderly groups and the top 20 enriched terms in the biological process were shown. b The GO analysis on DEGs induced by the third dose of ZF2001 in elderly groups and the top 20 enriched terms in the biological process were shown. c The GO chordal diagram showed the degree of differential expression changes of genes involved in each functional project. Criteria: p-value < 0.05. d The GO chordal diagram showed the degree of differential expression changes of genes involved in each functional project. Criteria: p-value < 0.05

The booster vaccine activated PBMC gene expression patterns related to immune responses on day 28

We similarly used the STRING database and MCODE plug-in of Cytoscape software to prioritize DEGs 28 days after the booster dose. A total of 4 characteristic up-regulated and down-regulated DEG sets were identified induced by the third dose of BBIBP-CorV and ZF2001 in the elderly on Day 28 (Figs. 6a, b and 7a, b).

Fig. 6

The PPI network extracted from initial PPI networks for the protein products of up and down-regulated DEGs induced by the third dose of BBIBP-CorV in elderly groups after 28 days, consisting of 36 nodes (a) and 30 nodes (b). c The GO analysis on DEGs induced by the third dose of BBIBP-CorV in elderly groups and only the top 20 enriched terms in the biological process were shown. d The connection of GO functions with specific DEGs. e The GO chordal diagram showed the degree of differential expression changes of genes involved in each functional project. Criteria: p-value < 0.05

Fig. 7

The PPI network extracted from initial PPI networks for the protein products of up and down-regulated DEGs induced by the third dose of ZF2001 in elderly groups after 28 days, consisting of 37 nodes (a) and 32 nodes (b). c The GO analysis on DEGs induced by the third dose of ZF2001 in elderly groups and only the top 20 enriched terms in the biological process were shown. d The connection of GO functions with specific DEGs. e The GO chordal diagram showed the degree of differential expression changes of genes involved in each functional project. Criteria: p-value < 0.05

These characteristic DEGs were also enriched in immune-related biological process GO-terms. Specifically, 28 days after the BBIBP-CorV and ZF2001 booster vaccine, DEGs in the senile group were enriched in several pathways such as the chemokine-mediated signaling pathway and humoral immune response (Figs. 6c and 7c). In the BBIBP-CorV booster vaccine group, these functional entries were mainly executed by CD83, CR2, CXCL1, CXCL2, CXCL3, CXCL8, IL1B, JCHAIN, PF4, PF4V1, PPBP, and TNF, etc. (Fig. 6d). In the ZF2001 booster vaccine group, C2, CXCR1, CXCR2, CX3CR1, CXCL5, CXCR1, CXCR2, LTF, PF4, PF4V1, and PPBP participated in immune regulation (Fig. 7d). Moreover, the chordal diagram further demonstrated that most of the DEGs involved in immune-related functional items were up-regulated, meaning that these immune-related pathways were widely activated 28 days after the booster vaccine in the elderly group (Figs. 6e and 7e).

The various immune landscapes of age-related gene modules between 7 and 28 days after booster immunization

To explore the impact of age characteristics on the transcriptome, we performed an analysis of the transcriptome using weighted gene co-expression network analysis(WGCNA). In order to ensure the construction of a scale-free network, the most suitable β was determined as a soft threshold parameter. A total of four co-expression modules were identified through hierarchical clustering, and a dendrogram of all DEGs was clustered based on a dissimilarity measure (1-TOM) at day 7. These four co expression modules based on differences in topological overlap were displayed in the 3D clustering map and assigned as brown, blue, turquoise, and grey module colors(Fig. 8a, b). Heatmap depicted the Topological Overlap Matrix (TOM) of genes selected for weighted co-expression network analysis(Fig. 8c). For each module, the gene co-expression was summarized by the eigengene and we calculated the correlations of each eigengene with clinical traits, such as vaccine type, gender and age. The correlation between the vaccination signatures and the co-expression module is shown in Fig. 8d, of which the brown module (eigengene value = -0.75, p = 0.03 × 10–6) was significantly negatively correlated with the age of the participants, but not significantly correlated with vaccine type and gender. The brown module includes 56 genes involved in many immune function-related genes, such as AK5, BTLA, CCR5, CD24C, CR2, DPP4, GCSAM, IL7R, IL23R, LEF1, LTB, RORC, and STAP1, etc. The GO enrichment analysis showed that these genes were mainly involved in "lymphocyte costimulation", "T cell differentiation involved in immune response", "CD4-positive, alpha–beta T cell differentiation involved in immune response", and "T-helper cell differentiation" (Fig. 10a). The above co-expression modules analysis, which is significantly negatively correlated with age, once again indicates that the immune response of the elderly to booster vaccines at 7 days was impaired and delayed.

Fig. 8

WGCNA of the PBMCs transcriptome after 7 days of the third booster dose. a, b The 3D cluster map of genes, based on differences in topological overlap matrix (TOM), and assigned module colors. c Heatmap depicts the TOM of genes selected for weighted co-expression network analysis. Light color represents lower overlap and red represents higher overlap. d Module-trait associations: Each row corresponds to a module eigengene and each column to a trait. Each cell contains the corresponding correlation and p-value

A similar WGCNA was performed on transcriptome data at 28 days. Five co-expression modules were identified on day 28, of which the brown module (eigengene value = 0.47, p = 0.01) and yellow module (eigengene value = 0.46, p = 0.01) were significantly positively correlated with the age of the participants (Fig. 9a-d). The co-expressed genes of these two modules are mainly involved in "fever generation", "chemokine activity", "chemokine receptor binding", "response to chemokine", and "chemokine-mediated signaling pathway" (Fig. 10b). This means that after 28 days of the booster vaccine, the immune transcriptional landscape associated with the elderly significantly improved.

Fig. 9

WGCNA of the transcriptome after 28 days of the third booster dose. a, b The 3D cluster map of genes, based on differences in topological overlap matrix (TOM), and assigned module colors. c Heatmap depicts the TOM of genes selected for weighted co-expression network analysis. Light color represents lower overlap and red represents higher overlap. d Module-trait associations: Each row corresponds to a module eigengene and each column to a trait. Each cell contains the corresponding correlation and p-value

Fig. 10

a The GO functional enrichment analysis of brown gene module after 7 days of the third booster dose. b The GO functional enrichment analysis of brown and yellow gene module after 28 days of the third booster dose