Single-cell RNA sequencing reveals the transcriptomic characteristics of peripheral blood mononuclear cells in hepatitis B vaccine non-responders

The emergence of a vaccine against hepatitis B has proven to be an important milestone in the prevention of this disease; however, 5%–10% of vaccinated individuals do not generate an immune response to the vaccine, and its molecular mechanism has not been clarified. In this study, single-cell RNA sequencing was performed on peripheral blood mononuclear cells (PBMCs) from three volunteers with a high immune response (HR) and three with no immune response (NR) to the hepatitis B vaccine. We found that the antigen-presenting activity scores of various antigen-presenting cells, the mitogen-activated protein kinase (MAPK) pathway activity scores of naive B cells, and the cell activity scores of three types of effector T cells were significantly decreased, whereas the cytotoxicity scores of CD3^highCD16^lowKLRG1^high natural killer T (NKT) cells were significantly increased in the NR group compared with those in the HR group. Additionally, the expression levels of some classical molecules associated with distinct signaling pathways—including HLA-B, HLA-DRB5, BLNK, BLK, IL4R, SCIMP, JUN, CEBPB, NDFIP1, and TXNIP—were significantly reduced in corresponding subsets of PBMCs from the NR group relative to those of the HR group. Furthermore, the expression of several cytotoxicity-related effector molecules, such as GNLY, NKG7, GZMB, GZMM, KLRC1, KLRD1, PRF1, CST7, and CTSW, was significantly higher in CD3^highCD16^lowKLRG1^high NKT cells derived from non-responders. Our study provides a molecular basis for the lack of response to the hepatitis B vaccine, including defective antigen presentation, decreased T cell activity, and reduced IL-4 secretion, as well as novel insight into the role of NKT cells in the immune response to the hepatitis B vaccine.

1 Introduction

The emergence of the hepatitis B vaccine is an important milestone in the prevention of hepatitis B; however, approximately 5%–10% (1) of individuals do not produce an immune response to the existing recombinant hepatitis B vaccine (2, 3). The reason for the lack of an immune response and the characteristics of peripheral blood mononuclear cells (PBMCs) in non-responders are still not fully understood.

To date, studies investigating the absence of an immune response to the hepatitis B vaccine have mainly been carried out at the population level. However, the heterogeneity among cells that play an important role in regulating cellular function is often masked by the population signal (4). Single-cell RNA sequencing (Single-cell RNA-seq) technology has been applied in the field of immunology and revealed the heterogeneity among immune cells (5). In this study, using this technology, we generated a comprehensive transcriptomic landscape of PBMCs from non-responders who received enhanced immunization with recombinant hepatitis B vaccine and identified several key cell type-specific molecules that may serve as targets for enhancing the immune response to the hepatitis B vaccine.

2 Subjects and methods

2.1 Subjects

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.

2.2 Sample collection and PBMC isolation

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.

2.3 Preparation of single-cell suspensions

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.

2.4 Single-cell RNA-seq

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.

2.5 Single-cell RNA-seq data processing

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).

2.6 Integration of different datasets

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.

2.7 Identification of cell type

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.

2.8 Identification and functional enrichment analysis of differentially expressed genes

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).

2.9 Expression enrichment scores of cell-state genes

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.

2.10 RT-qPCR

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.

3 Results

3.1 Single-cell atlas of PBMCs from non-responders and high responders

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).

FIGURE 1

Figure 1 The transcriptional map of peripheral blood mononuclear cells (PBMCs) from non-responders (NR) and high responders (HR) to the hepatitis B vaccine. (A) Flowchart showing the overall study design; 3′ single-cell RNA sequencing was undertaken for PBMCs from three NR and three HR volunteers. (B) A UMAP showing the 18 cell clusters identified from 52,299 high-quality cells derived from both groups; each dot represents one cell and dots are colored according to the cell type. (C) A balloon plot showing the number of cells that were captured from the various cell types in each sample. (D) A UMAP plot displaying the classical marker genes used to identify each cell type; genes are colored according to their expression level in the cells. (E) A dot plot showing the expression level of the selected classical marker genes in the 18 cell clusters; the shade of the dot represents the expression level, and the size of the dot represents the proportion of cells expressing the gene. (F) A heat map showing the three most highly expressed genes in each of the 18 cell clusters.

3.2 Transcriptional characteristics of antigen-presenting cells from non-responders

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).

FIGURE 2

Figure 2 Comparison of the characteristics of antigen-presenting cells (APCs) between the non-responder (NR) and high-responder (HR) groups. (A) A dot plot showing the proportions of cells expressing HLA-related genes and the average expression levels of HLA-related genes in APCs from the NR and HR groups. (B, C) A violin diagram showing the expression scores of HLA-I- and HLA-II-related genes in APCs from the NR and HR groups; the horizontal line dividing the box represents the median value and the two ends of the box represent the upper quartile (Q3) and lower quartile (Q1); the difference between Q1 and Q3 is the interquartile range (IQR). A P-value <0.05 was considered significant. *P<0.05; **P<0.01; ****P<0.0001; ns, not significant (Student’s t-test). B indicates the expression score of HLA-II-related genes and C indicates the expression score of HLA-I-related genes. (D) Enrichment analysis of genes with low expression in APCs; the Kyoto Encyclopedia of Genes and Genomes (KEGG) terms of interest are shown in red.

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).

3.3 Features of humoral immunity in non-responders

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).

FIGURE 3

Figure 3 The characteristics of humoral immunity in the non-responder (NR) and high-responder (HR) groups. (A) A volcano plot showing the differentially expressed genes (DEGs) as determined by the Wilcox rank-sum test (NR group vs. the HR group); adjusted P-values <0.05 and |avg_log2FC| > 0.25 indicate that the difference is significant. The blue dots and red dots, respectively, represent downregulated and upregulated genes in naive B cell or memory B cell subsets from the NR group; the genes of interest are marked with the respective symbol. (B) A Venn diagram showing the overlap between MAPK-related genes and the DEGs of two B cell subsets from the NR group. (C) A violin diagram showing the expression scores of MAPK-related DEGs in memory B cells and naive B cells from the NR group. A P-value <0.05 was considered significant. ***P<0.001 (Student’s t-test). (D) A Kyoto Encyclopedia of Genes and Genomes (KEGG) network showing the regulatory information of some downregulated genes in naive B cells from the NR group; the blue squares represent downregulated genes.

3.4 Characteristics of adaptive immunity in non-responders

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).

FIGURE 4

Figure 4 The features of T cell subsets from the non-responder (NR) and high-responder (HR) groups. (A) An upset diagram depicting the overlapping and non-overlapping differentially expressed genes (DEGs) in seven T cell subtypes. The horizontal bar chart on the left represents the number of elements in each set, the colored points in the middle and the lines between the points represent the intersection of different T cell subsets, and the vertical bar chart at the top represents the number of corresponding intersection elements. (B) A boxplot showing the log2 fold changes of the unique molecular identifier (UMI) counts of 12 DEGs common to the seven T cell subtypes. (C) Protein–protein interaction (PPI) network showing the GO terms associated with the 12 common DEGs. (D) Enrichment analysis of the downregulated genes in three T cell subtypes with marked differential characteristics; the Gene Ontology (GO) terms of interest are shown in red. (E) Cell activity scores of three T cell subtypes. A P-value <0.05 was considered significant. ****P<0.0001 (Student’s t-test).

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).

3.5 Transcriptional characteristics of natural killer and NKT cells of non-responders

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).

FIGURE 5

Figure 5 Features of natural killer (NK) and natural killer T (NKT) cell subsets from the non-responder (NR) and high-responder (HR) groups. (A) A UMAP displaying five cell subsets identified from 5,615 high-quality NK cells and 6,989 high-quality NKT cells in the NR and HR groups; each dot represents one cell and dots are colored according to the cell subset. (B) A dot plot showing the expression level of the selected classical marker genes in five cell subsets; the color of the dots represents the expression level, the size of each dot represents the proportion of cells expressing the corresponding gene, and the squares along the x-axis represent the cell subsets corresponding to the subsets in A. (C) A violin plot showing the expression levels of nine cytotoxicity-related marker genes with statistical significance in CD3^highKLRG1^high NKT cells; the Wilcox rank-sum test was used for analysis. The horizontal line dividing the box represents the median value and the two ends of the box represent the upper quartile (Q3) and lower quartile (Q1); the difference between Q1 and Q3 represents the interquartile range (IQR). (D) The cytotoxicity scores of CD3^highKLRG1^high NKT cells are shown; a P-value <0.05 was considered significant. ****P<0.0001 (Student’s t-test). (E) Enrichment analysis of upregulated genes in CD3^highKLRG1^high NKT cells from the NR group. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms of interest are shown in red.

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).

3.6 RT-qPCR verification of the results

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).

FIGURE 6

Figure 6 The relative expression levels of JUN, HLA-B, and TXNIP in peripheral blood mononuclear cells (PBMCs) of the non-responder (NR) (n=6) and high-responder (HR) (n=30) groups; the Mann-Whitney U Test was used to compare differences between independent groups and a P-value <0.05 was considered significant. *P<0.05, ***P<0.001.

4 Discussion

The mechanism underlying the absence of a response to the hepatitis B vaccine is still not fully understood. In this study, through a comparative analysis of single-cell transcriptional characteristics of non-responders and high responders to the vaccine, we have provided potential targets that may contribute to improving the intensity of the immune response to hepatitis B vaccination. Seventeen immune cell types were identified from the peripheral blood of volunteers. Although no differences in cell type composition were found between the HR and NR groups in our study, several studies have reported that the proportion of CD4⁺ Tregs is higher in non-responders than in high responders (16, 17), which may be associated with differences in the protocols used between these studies and ours.

In this study, we demonstrated the heterogeneity of expression of HLA molecules in different APCs at the single-cell level. In addition, except for pDCs, the expression enrichment score of HLA-I molecules in APC sub-clusters was significantly lower in non-responders than in high responders, while the expression levels of HLA-II molecules in the various APC sub-clusters were significantly lower in the NR group than in the HR group. This implied that antigen presentation activity was weaker in the former group than in the latter one. Additionally, the downregulated DEGs in the various APCs derived from the NR group were significantly enriched in the antigen presentation and processing signaling pathways, which again suggested that the antigen processing and presentation functions of APCs are impaired in non-responders. Our results also supported the hypothesis that the peptides of the hepatitis B surface antigen may not be degraded into smaller antigenic peptides because of deficient antigen presentation and processing or, alternatively, that a lack of HLA-I and HLA-II molecules connected to these peptides may lead to functional deficits in APCs (18, 19). Furthermore, our data suggested that HLA-DRB5 and HLA-B may be the key molecules responsible for hepatitis B virus (HBV) antigen presentation, and the low expression of HLA-DRB5 and HLA-B may be one of the main factors influencing the lack of response to the vaccine. Similar to this result, it was reported that HLA-DRB5 is significantly downregulated in the liver tissue of patients with hepatitis B who showed a poor response to antiviral therapy (20). A different study demonstrated that the concentrations of HLA-DR and HLA-ABC antigens, markers for activated mature monocyte-derived dendritic cells (moDCs), were significantly lower in non-responders than in high responders (21).

The number of B cells producing HBsAb was reported to be significantly reduced in non-responders; however, the molecular mechanism underlying this effect remains unclear (22). In our study, the scores for MAPK pathway activity in naive B cells from non-responders were especially low. The MAPK pathway is the hub of the eukaryotic signal transmission network and may play an important role in defending against virus infection (12, 23). Another study showed that the inhibition of the NF-κB and p38 MAPK pathways suppressed the HBsAg-induced secretion of cytokines in moDCs (24). We further found that molecules such as B cell linker protein (BLNK), B lymphoid tyrosine kinase (BLK), interleukin-4 receptor (IL-4R), SLP adaptor and CSK interacting membrane protein (SCIMP), NFKBIA, and JUN were underexpressed in naive B cells from non-responders. BLNK acts as a central linker protein downstream of the B-cell receptor (BCR) and connects spleen-associated tyrosine kinase (SYK) to several signaling pathways (25). Additionally, BLNK plays a crucial role in ERK/EPHB2, p38MAPK, JNK, NF-κB, NFAT, and AP-1 activation, the regulation of B cell function and development (26), and BCR-mediated activation of PLCG1, PLCG2, and Ca2⁺ mobilization (27). BLK, a member of the Src tyrosine kinase family, can exert both positive and negative regulatory effects on B cells (28, 29). IL-4R regulates the expression of numerous genes and promotes cell proliferation and differentiation through the common cytokine receptor γ chain-mediated activation of STAT proteins (30). The low expression of IL-4R in naive B cells in the NR group observed in the present study is in agreement with the results of a study demonstrating that IL-4 secretion is reduced in non-responders to the hepatitis B vaccine (31). SCIMP is expressed in B cells and other professional APCs and undergoes tyrosine phosphorylation under the stimulation of the major histocompatibility complex (MHC) type II signal. Phosphorylated SCIMP then recruits SLP65/SLP76 and GRB2 complexes, thereby activating the PLCγ1/2 and RAS pathways, among others (32–34). NFKBIA can promote its own synthesis by binding to the NF-κB transcription factor, which can block the binding of NF-κB to DNA, thus inhibiting the NF-κB signaling pathway (35). One study reported that the expression of IκBα (encoded by NFKBIA) is increased in gastric cancer after EBV infection (36). Combined, these observations suggest that the reduced expression of these molecules in naive B cells in non-responders reduces the activity of the MAPK and NF-κB pathways, consequently weakening the expression of AP-1 family members. These effects may underlie the lack of response to the hepatitis B vaccine.

The recognition of foreign antigens by both CD4⁺ and CD8⁺ T cells is MHC-restricted (37, 38). Like that seen with HLA-I and HLA-II molecules, we found that the expression levels of genes associated with Ca²⁺ signaling and oxidative phosphorylation were reduced in seven subsets of T cells, with JUN being the representative gene. c-Jun N-terminal phosphorylation was reported to play an important role in the progression of HBV infection (39). Furthermore, two genes downregulated in CD4⁺ Teff cells in the NR group—CEBPB and NDFIP1—were associated with IL-4 production, which was in line with a previous report that showed the secretion of IL-4 was reduced in non-responders to the hepatitis B vaccine (31). CCAAT enhancer binding protein beta (CEBPB) can influence the activation of IL-4/IL-13 and NF-κB signaling through transcriptional repression during ER stress (40, 41). Nedd4 family interacting protein 1 (NDFIP1) plays a crucial role in CD4⁺ T cell proliferation and can maintain immune homeostasis by inhibiting the secretion of IL-4 (42). Our results further demonstrated that the expression scores of cell activity-related genes in CD4⁺ Teff cells were significantly reduced in non-responders. Additionally, several studies have suggested that changes in cell activity may affect the immune response to HBV (43, 44).

In this study, we found that CD3^highCD16^lowKLRG1^high NKT cells from the NR group had higher cytotoxicity-related gene set expression scores than those of the HR group, and the upregulated genes were enriched in cytotoxicity-related pathways. The above characteristics are similar to those reported for CD56^dim NK cells and KLRG1⁺ NKT cells from patients with chronic hepatitis B, which had specific memory for and stronger cytotoxic effects against HBsAg-pulsed moDCs, as well as proliferative responses to HBV antigen relative to unvaccinated people (45). In addition, it has been shown that the expression scores for these representative cytotoxicity-related gene sets in a subset of NKT cells were significantly higher in patients infected with SARS-CoV-2 than in healthy people (11). There is evidence to support that, on the one hand, vaccines can induce specific immune memory in NK cells, while on the other, NK cells can also indirectly regulate vaccine-induced humoral and adaptive immunity (46, 47). These observations largely support our findings that the increased cytotoxicity of CD3^high CD56^low KLRG1^high NKT cells is a key factor associated with non-response to the hepatitis B vaccine.

In addition to JUN and HLA-B, we found that TXNIP, one of the identified DEGs between the NR and HR groups, was downregulated in most cells. One study showed that reducing the expression of TXNIP can inhibit hepatitis C virus replication; however, the mechanism responsible for this phenomenon is not clear (48). Another study showed that TXNIP can maintain cells in a higher state of oxidative stress and promote immunoglobulin production (49).

In summary, some key cells and molecules that may regulate the intensity of the immune response to the hepatitis B vaccine were revealed via a comparative single-cell analysis of the transcriptional characteristics of PBMCs from non-responders and high responders. A key reason for the small sample in this study is that it is extremely difficult to identify hepatitis B vaccine non-responders. Accordingly, the mechanism relating to the production of HBsAb by these key cells, the roles played by the genes identified in this study in the absence of response to the hepatitis B vaccine, and the ability of CD3^highCD56^highKLRG1^high NKT cells to fight HBV infection should be further explored through cell- and animal-based experiments.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE229431.

Ethics statement

The studies involving human participants were reviewed and approved by Ethics Committee of Lanzhou First People’s Hospital. The patients/participants provided their written informed consent to participate in this study.

Author contributions

All the authors contributed to sample collection and the isolation of peripheral blood mononuclear cells. MZ analyzed the single-cell data and performed RT-qPCR. MZ, CW, PL, TS, and XDX participated in the design of the study. MZ and XDX wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded by the Gansu Province Science and Technology Innovation Service Platform Construction Project (18JR2TA024) and the Youth Science and Technology Foundation of Gansu Province (20JR5RA141).

Acknowledgments

We thank all the volunteers who actively participated in this project for their contribution to blood sample collection as well as the relevant specialized personnel from Bomiao Biotechnology Co., Ltd for their technical support. We are extremely grateful to Dr. Jianming Zeng (University of Macau) and all the members of his bioinformatics team for generously sharing their experience and codes.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1091237/full#supplementary-material

References

1. Zuckerman JN. Protective efficacy, immunotherapeutic potential, and safety of hepatitis B vaccines. J Med Virol (2006) 78:169–77. doi: 10.1002/jmv.20524

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Nishida N, Sugiyama M, Sawai H, Nishina S, Sakai A, Ohashi J, et al. Key HLA-DRB1-DQB1 haplotypes and role of the BTNL2 gene for response to a hepatitis B vaccine. Hepatology (2018) 68:848–58. doi: 10.1002/hep.29876

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yue X, Ge C, Zhuge S, He H, Yang H, Xu H, et al. Changes and analysis of anti-HBs titres after primary immunization in 1- to 16-year-old Chinese children: A hospital-based study. J Viral Hepat (2018) 25:373–80. doi: 10.1111/jvh.12818

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Hedlund E, Deng Q. Single-cell RNA sequencing: Technical advancements and biological applications. Mol Aspects Med (2018) 59:36–46. doi: 10.1016/j.mam.2017.07.003

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Waickman AT, Victor K, Li T, Hatch K, Rutvisuttinunt W, Medin C, et al. Dissecting the heterogeneity of DENV vaccine-elicited cellular immunity using single-cell RNA sequencing and metabolic profiling. Nat Commun (2019) 10:3666. doi: 10.1038/s41467-019-11634-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Stuart T, Butler A, Hoffman P, Hafemesiter C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell (2019) 177:1888–902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhsi AH, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun (2019) 10:1523. doi: 10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res (2021) 49:D605–d12. doi: 10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Andreatta M, Carmona SJ. UCell: Robust and scalable single-cell gene signature scoring. Comput Struct Biotechnol J (2021) 19:3796–98. doi: 10.1016/j.csbj.2021.06.043

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Zhang JY, Wang XM, Xing X, Xu Z, Zhang C, Song JW, et al. Single-cell landscape of immunological responses in patients with COVID-19. Nat Immunol (2020) 21:1107–18. doi: 10.1038/s41590-020-0762-x

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Huang L, Shi Y, Gong B, Jiang L, Zhang Z, Liu X, et al. Dynamic blood single-cell immune responses in patients with COVID-19. Signal Transduct Target Ther (2021) 6:110. doi: 10.1038/s41392-021-00526-2

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pleiman CM, Clark MR, Gauen LK, Winitz S, Coggeshall KM, Johnson GL, et al. Mapping of sites on the Src family protein tyrosine kinases p55blk, p59fyn, and p56lyn which interact with the effector molecules phospholipase C-gamma 2, microtubule-associated protein kinase, GTPase-activating protein, and phosphatidylinositol 3-kinase. Mol Cell Biol (1993) 13:5877–87. doi: 10.1128/mcb.13.9.5877-5887.1993

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Lv H, Qi Z, Wang S, Feng H, Deng X, Ci X. Asiatic Acid Exhibits Anti-inflammatory and Antioxidant Activities against Lipopolysaccharide and d-Galactosamine-Induced Fulminant Hepatic Failure. Front Immunol (2017) 8:785. doi: 10.3389/fimmu.2017.00785

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yang Z, Wang L, Niu W, Wu Y, Zhang J, Meng G. Increased CD4+CD25+FoxP3+ regulatory T cells in the blood of nonresponders after standard hepatitis B surface antigen vaccine immunization. Clin Immunol (2008) 127:265–6. doi: 10.1016/j.clim.2007.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Li J, Tan D, Liu H, Li K. CD4(+) CD25(+) FoxP3(+) T regulatory cells in subjects responsive or unresponsive to hepatitis B vaccination. Zhong Nan Da Xue Xue Bao Yi Xue Ban (2011) 36:1046–51. doi: 10.3969/j.issn.1672-7347.2011.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Milich DR, Leroux-Roels GG. Immunogenetics of the response to HBsAg vaccination. Autoimmun Rev (2003) 2:248–57. doi: 10.1016/s1568-9972(03)00031-4

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Mineta M, Tanimura M, Tana T, Yssel H, Kashiwagi S, Sasazuki T. Contribution of HLA class I and class II alleles to the regulation of antibody production to hepatitis B surface antigen in humans. Int Immunol (1996) 8:525–31. doi: 10.1093/intimm/8.4.525

PubMed Abstract | CrossRef Full Text | Google Scholar

20. He Y, Zhou Y, Wang H, Yin J, Chang Y, Hu P, et al. Identifying potential biomarkers in hepatitis B virus infection and its response to the antiviral therapy by integrated bioinformatic analysis. J Cell Mol Med (2021) 25:6558–72. doi: 10.1111/jcmm.16655

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Verkade MA, van Druningen CJ, Op de Hoek CT, Weimar W, Betjes MGH. Decreased antigen-specific T-cell proliferation by moDC among hepatitis B vaccine non-responders on haemodialysis. Clin Exp Med (2007) 7:65–71. doi: 10.1007/s10238-007-0127-x

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Valats JC, Tuaillon E, Funakoshi N, Hoa D, Brabet MC, Bolloré K, et al. Investigation of memory B cell responses to hepatitis B surface antigen in health care workers considered as non-responders to vaccination. Vaccine (2010) 28:6411–6. doi: 10.1016/j.vaccine.2010.07.058

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Manley GCA, Parker LC, Zhang Y. Emerging regulatory roles of dual-specificity phosphatases in inflammatory airway disease. Int J Mol Sci (2019) 20:678. doi: 10.3390/ijms20030678

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Jan RH, Lin YL, Chen CJ, Lin TY, Hsu YC, Chen LK, et al. Hepatitis B virus surface antigen can activate human monocyte-derived dendritic cells by nuclear factor kappa B and p38 mitogen-activated protein kinase mediated signaling. Microbiol Immunol (2012) 56:719–27. doi: 10.1111/j.1348-0421.2012.00496.x

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Fu C, Turck CW, Kurosaki T, Chan AC. BLNK: a central linker protein in B cell activation. Immunity (1998) 9:93–103. doi: 10.1016/s1074-7613(00)80591-9

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Minegishi Y, Rohrer J, Coustan-Smith E, Lederman HM, Pappu R, Campana D, et al. An essential role for BLNK in human B cell development. Science (1999) 286:1954–7. doi: 10.1126/science.286.5446.1954

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Taguchi T, Kiyokawa N, Takenouch H, Matsui J, Tang WR, Nakajima H, et al. Deficiency of BLNK hampers PLC-gamma2 phosphorylation and Ca2+ influx induced by the pre-B-cell receptor in human pre-B cells. Immunology (2004) 112:575–82. doi: 10.1111/j.1365-2567.2004.01918.x

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Tretter T, Ross AE, Dordai DI, Desiderio S. Mimicry of pre-B cell receptor signaling by activation of the tyrosine kinase Blk. J Exp Med (2003) 198:1863–73. doi: 10.1084/jem.20030729

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Bewarder N, Weinrich V, Budde P, Hartmann D, Flaswinkel H, Reth M, et al. In vivo and in vitro specificity of protein tyrosine kinases for immunoglobulin G receptor (FcgammaRII) phosphorylation. Mol Cell Biol (1996) 16:4735–43. doi: 10.1128/MCB.16.9.4735

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Lin JX, Leonard WJ. The common cytokine receptor γ Chain family of cytokines. Cold Spring Harb Perspect Biol (2018) 10:a028449. doi: 10.1101/cshperspect.a028449

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Qiu S, He P, Fang X, Tong H, Lv J, Liu J, et al. Significant transcriptome and cytokine changes in hepatitis B vaccine non-responders revealed by genome-wide comparative analysis. Hum Vaccin Immunother (2018) 14:1763–72. doi: 10.1080/21645515.2018.1450122

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Draber P, Vonkova I, Stepanek O, Hrdinka M, Kucova M, Skopcova T, et al. SCIMP, a transmembrane adaptor protein involved in major histocompatibility complex class II signaling. Mol Cell Biol (2011) 31:4550–62. doi: 10.1128/MCB.05817-11

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kralova J, Fabisik M, Pokorna J, Skopcova T, Malissen B, Brdicka T. The transmembrane adaptor protein SCIMP facilitates sustained dectin-1 signaling in dendritic cells. J Biol Chem (2016) 291:16530–40. doi: 10.1074/jbc.M116.717157

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Vanshylla K, Bartsch C, Hitzing C, Krümpelmann L, Wienands J, Engels N. Grb2 and GRAP connect the B cell antigen receptor to Erk MAP kinase activation in human B cells. Sci Rep (2018) 8:4244. doi: 10.1038/s41598-018-22544-x

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Hoffmann A, Levchenko A, Scott ML, Baltimore D. The IkappaB-NF-kappaB signaling module: temporal control and selective gene activation. Science (2002) 298:1241–5. doi: 10.1126/science.1071914

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang Y, Liu W, Zhang W, Wang W, Song Y, Xiao H, et al. Constitutive activation of the canonical NF-κB signaling pathway in EBV-associated gastric carcinoma. Virology (2019) 532:1–10. doi: 10.1016/j.virol.2019.03.019

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Andersen MH, Schrama D, Thor Straten P, Becker JC. Cytotoxic T cells. J Invest Dermatol (2006) 126:32–41. doi: 10.1038/sj.jid.5700001

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Caminero F, Iqbal Z, Tadi P. Histology, cytotoxic T cells. In: StatPearls. Treasure Island (FL) StatPearls Publishing LLC (2022). Available at: https://www.ncbi.nlm.nih.gov/books/NBK559279/.

Google Scholar

39. Benn J, Su F, Doria M, Schneider RJ. Hepatitis B virus HBx protein induces transcription factor AP-1 by activation of extracellular signal-regulated and c-Jun N-terminal mitogen-activated protein kinases. J Virol (1996) 70:4978–85. doi: 10.1128/JVI.70.8.4978-4985.1996

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Arensdorf AM, Rutkowski DT. Endoplasmic reticulum stress impairs IL-4/IL-13 signaling through C/EBPβ-mediated transcriptional suppression. J Cell Sci (2013) 126:4026–36. doi: 10.1242/jcs.130757

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Hayakawa K, Nakajima S, Hiramatsu N, Okamura M, Huang T, Yukinori S, et al. ER stress depresses NF-kappaB activation in mesangial cells through preferential induction of C/EBP beta. J Am Soc Nephrol (2010) 21:73–81. doi: 10.1681/ASN.2009040432

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Layman AAK, Deng G, O’Leary CE, Tadros S, Thomas RM, Dybad JM, et al. Ndfip1 restricts mTORC1 signalling and glycolysis in regulatory T cells to prevent autoinflammatory disease. Nat Commun (2017) 8:15677. doi: 10.1038/ncomms15677

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Goncalves L, Albarran B, Salmen S, Borges L, Fields H, Montes H, et al. The nonresponse to hepatitis B vaccination is associated with impaired lymphocyte activation. Virology (2004) 326:20–8. doi: 10.1016/j.virol.2004.04.042

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Goncalves L, Barboza L, Albarrán B, Salem S, Montes H, Hernández M, et al. [Pattern of T cell activation in absence of protective immunity against hepatitis B virus. Review]. Invest Clin (2006) 47:83–96.

PubMed Abstract | Google Scholar

45. Wijaya RS, Read SA, Truong NR, Han S, Chen D, Shahidipour H, et al. HBV vaccination and HBV infection induces HBV-specific natural killer cell memory. Gut (2021) 70:357–69. doi: 10.1136/gutjnl-2019-319252

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Hatfield SD, Daniels KA, O'Donnell CL, Waggoner SN, Welsh RM. Weak vaccinia virus-induced NK cell regulation of CD4 T cells is associated with reduced NK cell differentiation and cytolytic activity. Virology (2018) 519:131–44. doi: 10.1016/j.virol.2018.04.012

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Rydyznski C, Daniels KA, Karmele EP, Brooks TR, Mahl SE, Moran MT, et al. Generation of cellular immune memory and B-cell immunity is impaired by natural killer cells. Nat Commun (2015) 6:6375. doi: 10.1038/ncomms7375

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Blackham S, Baillie A, Al-Hababi F, Remlinger K, You S, Hamatake R, et al. Gene expression profiling indicates the roles of host oxidative stress, apoptosis, lipid metabolism, and intracellular transport genes in the replication of hepatitis C virus. J Virol (2010) 84:5404–14. doi: 10.1128/JVI.02529-09

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Huang Y, Mao Z, Zhang X, Yang X, Sawada N, Takeda M, et al. Connexin43 is required for the effective activation of spleen cells and immunoglobulin production. Int J Mol Sci (2019) 20:5789. doi: 10.3390/ijms20225789

PubMed Abstract | CrossRef Full Text | Google Scholar