From January 2020 to March 2021, a total of 36 volunteers were recruited from the first People’s Hospital of Lanzhou, Gansu Province, including six who did not develop an immune response (NR) and 30 who displayed a high immune response (HR) to the recombinant hepatitis B vaccine. All the volunteers met all the following criteria: (1) at least 18 years of age, standard body mass index and transaminase, and no current or past history of the disease; (2) negative for serum markers of hepatitis B and had a history of hepatitis B vaccination (all serum markers are described in Supplementary Methods ); (3) willing to be inoculated with 20 μg/ml recombinant hepatitis B vaccine (Shenzhen Kangtai Biological Products Co., Ltd, Shenzhen, China) following a 0-, 1-, and 6-month schedule. One month after the third dose of the hepatitis B vaccine, participants whose serum concentrations of the hepatitis B surface antibody (HBsAb) were lower than 10 mIU/ml were considered non-responders and those whose HBsAb titer was greater than or equal to 1,000 mIU/ml were considered high responders (HR). The instruments and reagents used for the detection of serum markers are provided in Supplementary Methods .

One month after the administration of the third dose of the vaccine, 8 ml of peripheral venous blood was collected from each volunteer. Whole blood was diluted with an equal volume of phosphate-buffered saline (PBS) (Tianjin Hao Yang Biological Manufacture Co., Ltd, Tianjin, China), transferred to a centrifuge tube containing Human Lymphocyte Separation Medium (Dakewe Biotech Co., Ltd, Shenzhen, China), and centrifuged at 700×g for 20 min. PBMCs, located between the plasma layer and the cell separation medium, were carefully transferred into another centrifuge tube and washed twice with PBS.

PBMCs obtained from three high responders and three non-responders matched for age, sex, and ethnicity were used for single-cell RNA-seq. The PBMCs were resuspended in 1 ml of Iscove’s modified Dulbecco’s medium supplemented with 10% fetal bovine serum and a single-cell suspension was prepared. Cell counts and cell survival (greater than 85%) were determined using a Luna-FL automatic cell counter. Dead cells were removed using a dead cell removal kit. Finally, the cell concentration was adjusted to 700–1,200 cells/μl.

Single-cell RNA-seq libraries were constructed using a Single Cell 3′ Library and Gel Bead Kit (v3.1; 10X Genomics, 1000121) and a Chromium Single Cell A Chip Kit (10X Genomics, 1000120) according to the manufacturer’s instructions. Single cells were suspended in PBS containing 0.04% BSA (700–1,200 living cells per ml as determined using the LUNA-FL cell counter) and loaded onto a Chromium single-cell controller (10X Genomics) to generate single-cell gel bead-in-emulsion (GEMs) according to the manufacturer’s protocol. Approximately 10,000 cells were added to each channel and approximately 10,000 target cells were captured from each sample. The captured cells were lysed, and the released RNA was barcoded through reverse transcription in each GEM. Reverse transcription was performed using a C1000 Touch Thermal Cycler (Bio-Rad) at 53°C for 45 min, followed by 85°C for 5 min, and held at 4°C. The resulting cDNA was amplified and then assessed for quality using an Agilent 4200 TapeStation. The libraries were paired-end (150 bp) sequenced using the Illumina NovaSeq 6000 system. The above experiments were carried out by BioMiao Biological Technology Co., Ltd, Beijing, China.

Raw gene expression matrices were generated for each sample using the Cell Ranger (v.4.0.0) pipeline coupled with the human transcriptome reference version GRCh38-2020-A. The generated filtered gene expression matrices were analyzed by R software (v.4.1.1) with the Seurat 4 package (v.4.0.0). In brief, genes that were expressed in more than 1% of cells were discarded, while cells were selected for further analyses if they met the following criteria: (1) 500 < unique molecular identifiers (UMIs) < mean + 2 standard deviations; (2) 400<genes<mean + 2 standard deviations; (3) <10% mitochondrial genome. The gene expression matrices were then normalized using LogNormalize in Seurat and 2,000 features with high cell-to-cell variation were calculated using the FindVariableFeatures function. To reduce the dimensionality of the datasets, the RunPCA function was run with default parameters on linear-transformation scaled data generated by the ScaleData function. Next, the true dimensionality of each dataset was identified using the functions ElbowPlot, DimHeatmap, and JackStrawPlot according to the recommendations of the Seurat developers. Finally, the top 50 principal components were clustered using the FindNeighbors and FindClusters functions, and non-linear dimensional reduction was performed using the RunUMAP function with default settings. All the details of the Seurat analyses that are relevant to this study can be found in the website tutorial (https://satijalab.org/seurat/articles/pbmc3k_tutorial.html).

To compare cell types and the proportions of each cell type between the NR group and the HR group, six datasets were integrated using canonical correlation analysis plus mutual nearest neighbor (CCA+MNN) (6). Finally, six distinct single-cell RNA-seq datasets were assembled into an integrated and unbatched dataset.

The FindAllMarkers function was used to find markers for each cluster. Clusters were then classified and annotated based on the expression of canonical markers. Clusters that expressed two or more canonical cell-type markers were classified as doublet cells and excluded from further analysis.

The FindMarkers function and the Wilcox rank-sum test in Seurat were used to analyze gene expression levels in each cell cluster. Finally, DEGs were screened according to the following conditions: (1) Avg_log₂ (fold change) > 0.25; (2) min. PCT > 0.1; and (3) P. value _ adjust < 0.05. Functional enrichment analysis of the DEGs was conducted using the clusterProfiler package (7), Metascape (8), and the STRING database (9).

The UCELL package (10) was used to calculate the scores of the expression of HLA-related, cytotoxic activity-related, MAPK pathway-related, and cell activity-related gene sets. According to GSEA (Molecular Signature Database) and related literature (11–15), we selected 10 cytotoxicity-related genes, NKG7, CCL4, CST7, PRF1, GZMA, GNLY, GZMH, KLRB1, KLRD1, and CTSW, and 17 MAPK pathway-related genes, DDIT4, DUSP1, FOS, FOSB, JUN, JUNB, JUND, DDX5, PDCD4, DNAJA1, BLK, BLNK, UBA52, CXCR4, HSPA8, NFKBIA, and NFKBID. The cell activity-related genes were derived from the Gene Ontology (GO) Biological Process category and comprised ACTB, ACTG1, ANXA1, RHOA, B2M, ZFP36L2, CD8B, CD44, KLF6, FCGR3A, GNAS, HLA-A, HLA-E, HMGB1, IL7R, ITGB1, ITGB2, JUN, MSN, MYH9, PRKCB, PTGDS, RPS6, CCL3, VAMP2, TLN1, YY1, ITM2A, KLRK1, KMT2E, DOCK8, JAML, and NCR3.

Total RNA was extracted from PBMCs using a GeneJET RNA Purification Kit (Thermo Fisher Scientific, USA) according to the manufacturer’s instructions. RNA concentration (ng/μl) and purity (260/280 ratio) were determined using a NanoDrop One UV-Vis Spectrophotometer (Thermo Fisher Scientific, USA). The β-actin gene served as the internal control. All primers were designed by Primer-Blast and spanned exon–exon junctions. Extracted RNA was reverse transcribed to cDNA using a Goldenstar RT6 cDNA Synthesis Kit Ver.2 (Tsingke Biotechnology Co., Ltd, Beijing, China). qPCR was performed using TSINGKE Master SYBR Green I qPCR Mix (Tsingke Biotechnology Co., Ltd, Beijing, China) with a QuantStudio™ 5 Real-Time PCR System according to the manufacturer’s instructions. Relative gene expression levels were calculated using the 2^−ΔΔCt method. The qPCR results were plotted and analyzed in GraphPad Prism v7.0 and are expressed as mean ± standard deviation. The Mann–Whitney U Test was used to compare differences between independent groups and a P-value <0.05 was considered significant. The sequences of the primers used for qPCR were (5′–3′) β-actin, forward 5′-CCTTCCTGGGCATGGAGTC-3′ and reverse 5′-TGATCTTCATTGTGCTGGGTG-3′; HLA-B, forward 5′-AGCCATCTTCCCAGTCCA-3′ and reverse 5′-AGCTCCGATGACCACAA-3′; TXNIP, forward 5′-CCAGCACTTGGTCAGTCA-3′ and reverse 5′-CCATCTTCAGCCCACACT-3′; and JUN, forward 5′-GAGCATGACCCTGAACCTG-3′ and reverse 5′-CCGTTGCTGGACTGGATT-3′. Details of the RT-qPCR assay are described in Supplementary Methods .

To characterize the immunological characteristics of non-responders and high responders to the hepatitis B vaccine, we performed single-cell RNA-seq on PBMCs from the three volunteers in each group ( Figure 1A ). The characteristics of the six samples are detailed in Supplementary Table 1 . After quality control of the sequencing data, 52,299 high-quality cells were obtained, including 26,434 from the three non-responders and 25,865 from the three high responders. All these cells were grouped in 28 clusters according to the transcription characteristics of the top 50 principal components ( Supplementary Figures 1A–F ). Based on the expression of canonical marker genes, we identified 17 major cell types among the 28 clusters ( Figures 1B–F ; Supplementary Table 2 ). No difference in the percentage of each cell type was found between the two groups (P>0.05) ( Supplementary Figure 1G ).

As previously shown in Figure 1 , six types of APCs were identified, including naive B cells, memory B cells, CD14⁺ monocytes, CD16⁺ monocytes, classical dendritic cells (cDCs), and plasmacytoid dendritic cells (pDCs). We first examined the expression of HLA-I and HLA-II molecules in the APCs and found significant heterogeneity in the expression of these molecules between the APCs of the two groups ( Figure 2A ). Differential analysis indicated that, in the NR group, the expression level of HLA-DRB5 was markedly downregulated in several APCs, while HLA-B was significantly underexpressed in CD14⁺ monocytes, CD16⁺ monocytes, and pDCs relative to the HR group ( Supplementary Figure 2 ; Supplementary Table 3 ).

We further found that the expression score for HLA-II molecules in the various APCs was significantly lower in the NR group than in the HR group ( Figure 2B ). Except for pDCs, the expression score of HLA-I molecules in the APCs was also significantly lower in non-responders than in high responders ( Figure 2C ). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis showed that the downregulated genes in the various APCs of the NR group were significantly enriched in the antigen processing and presentation pathway ( Figure 2D ).

To screen for the key cells and molecules that regulate the humoral immune response, we analyzed the DEGs in the naive B cell and memory B cell subsets between the NR and HR groups. A total of 155 DEGs were identified in the naive B cell subset, accounting for 3.95% (155/3,925) of the total number of captured genes; this number was higher than that detected in the memory B cell subset, in which 72 DEGs were identified, accounting for 2.17% (72/3479) of the total number of captured genes ( Figure 3A ; Supplementary Table 4 ). Additionally, in the protein–protein interaction (PPI) network constructed with the DEGs, the average node degree and the number of nodes were higher, and the edges were more numerous, in the naive B cell subset than in the memory B cell subset ( Supplementary Table 5 ). Meanwhile, 15 MAPK pathway-related DEGs were identified in the naive B cell subset, 14 of which were downregulated and one upregulated in the NR group relative to that in the HR group; eight DEGs were identified in the memory B cell subset, including seven that were downregulated and one that was upregulated in the NR group ( Figure 3B ; Supplementary Figures 3A, B ; Supplementary Table 6 ). Furthermore, the MAPK pathway scores of the naive B cell and memory B cell subsets were lower in the NR group than in the HR group, with the lowest score being observed in the NR-derived naive B cell subset ( Figure 3C ; Supplementary Figure 3C ). Additionally, we found that several molecules (BLNK, BLK, IL4R, SCIMP, NFKBIA, JUN, and JUND, among others) that act downstream of the B cell and IL4R receptor signaling pathways and participate in Ca²⁺, MAPK, and NF-κB signaling were downregulated in the naive B cell subset in the NR group ( Figure 3D ; Supplementary Figure 3D ).

As previously shown in Figure 1 , seven T cell subsets were identified in this study. To investigate the features of adaptive immunity in non-responders, we compared the gene expression profile of each T cell subset between the two groups ( Supplementary Figures 4A–G ). Intersection analysis of the DEGs indicated that the three cell clusters with the highest numbers of DEGs and cell type-specific DEGs were CD8⁺ effector memory T (Tem) cells, CD4⁺ effector T (Teff) cells, and CD8⁺ Teff cells ( Figure 4A ). Meanwhile, 12 DEGs were common to all the T cell subsets, including nine that were downregulated and three that were upregulated ( Figures 4A, B ). A PPI network of the 12 DEGs demonstrated that TXNIP, JUN, and FOSB are involved in calcium ion reactions and are associated with oxidative phosphorylation ( Figure 4C ).

Next, we performed a GO term enrichment analysis on the three above-mentioned T cell subsets ( Figure 4D ). The results showed that all the downregulated genes in these T cell subsets from the NR group were involved in biological processes associated with immune cell activity, including cell activation, calcium ion reaction, interferon-γ, and oxidative phosphorylation. Two genes downregulated in CD4⁺ Teff and CD8⁺ Teff cells—CEBPB and NDFIP1—were involved in the production of interleukin-4 (IL-4). Other downregulated genes identified ( Supplementary Table 7 ) were involved in pathways such as the T cell receptor signaling pathway, the antigen receptor-mediated signaling pathway, and T cell proliferation. Finally, all the cell activity scores of the three T cell subsets were significantly lower in the NR group than in the HR group ( Figure 4E ).

According to their transcriptional characteristics, two subsets of NK cells were characterized as CD56^lowCD16^high and CD56^highCD16^low. Three subsets of NKT cells identified were categorized as CD3^lowCD16^high, CD3^highCD16^lowKLRG1^high, and CD3^lowCD16^lowMKI67^high ( Figures 5A, B ). No difference in the proportion of each subset was observed between the NR and HR groups ( Supplementary Figure 5A, B ).

To understand the role of NK and NKT cells in the immune response to the hepatitis B vaccine, the transcriptional characteristics of each cell subset were compared between the two groups ( Supplementary Figure 5C ). CD56^lowCD16^high NK cells and CD3^highCD16^lowKLRG1^high NKT cells were the subtypes with the greatest number of cell-type specific DEGs ( Supplementary Figure 5D ). Next, we undertook a comprehensive analysis of the DEGs from the two subtypes. Notably, we found that the relative expression levels of cytotoxicity-associated markers (such as GNLY, NKG7, GZMB, GZMM, KLRC1, KLRD1, PRF1, CST7, and CTSW) in the CD3^highCD16^lowKLRG1^high NKT cell subset were higher in the NR group than in the HR group ( Figure 5C ). Additionally, the cytotoxicity scores of the cluster from the NR group were significantly higher than those in the HR group ( Figure 5D ). Furthermore, we observed that the upregulated DEGs in the NR group were involved in biological processes such as cell killing and NK cell-mediated cytotoxicity ( Figure 5E ).

Given the difference in sensitivity between RT-qPCR and single-cell RNA-seq, we selected three genes—JUN, HLA-B, and TXNIP—that were noticeably downregulated in most cell subtypes in our single-cell RNA-seq analysis for verification by RT-qPCR ( Supplementary Figure 6 ). Consistent with the sequencing results, the RT-qPCR analysis showed that these three genes were significantly downregulated in the NR group compared with the HR group ( Figure 6 ).