Transcriptome Analysis of Peripheral Blood Mononuclear Cells in SARS-CoV-2 Naïve and Recovered Individuals Vaccinated With Inactivated Vaccine

The urgent approval of the use of the inactivated COVID-19 vaccine is essential to reduce the threat and burden of the epidemic on global public health, however, our current understanding of the host immune response to inactivated vaccine remains limited. Herein, we performed serum IgG antibody detection and transcriptomics analysis on 20 SARS-CoV-2 naïve individuals who received multiple doses of inactivated vaccine and 5 SARS-CoV-2 recovered individuals who received single dose of inactivated vaccine. Our research revealed the important role of many innate immune pathways after vaccination, identified a significant correlation with the third dose of booster vaccine and proteasome-related genes, and found that SARS-CoV-2 recovered individuals can produces a strong immune response to a single dose of inactivated vaccine. These results help us understand the reaction mechanism of the host’s molecular immune system to the inactivated vaccine, and provide a basis for the choice of vaccination strategy.

Introduction

Since December 2019, a new severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has swept the world, causing a variety of clinical syndromes termed coronavirus disease 2019 (COVID-19) (Wu et al., 2020). The clinical manifestations of COVID-19 include fever, dry cough, fatigue, sore throat, pneumonia, diarrhea and other symptoms, and may even develop into severe pneumonia, acute respiratory distress syndrome (ARDS) or multiple organ failure (Huang et al., 2020). The World Health Organization declared a pandemic in March 2020. COVID-19 has caused considerable impacts on the global economy and public health.

Although for a long time, people have relied on social distancing, hygiene measures, and repurposed drugs to coping with it, now many researchers are committed to developing safe and effective vaccines to establish herd immunity to prevent SARS-CoV-2. In view of the turmoil caused by the COVID-19 pandemic and the urgent need for effective vaccine, vaccine development can be accelerated by combining originally requested phases. The vaccine does not go through a complete approval process, but may be approved for emergency use (Sharma et al., 2020). Early clinical trial results of the inactivated vaccines produced by Sinopharm and Sinovac showed a low incidence of adverse reactions and good immunogenicity (Palacios et al., 2020; Xia et al., 2020; Xia et al., 2021a; Xia et al., 2021b; Zhang et al., 2021). On July 22, 2020, the above two candidate inactivated vaccines were approved for use (Zhao et al., 2020a). These two vaccines are widely promoted and vaccinated, but the reaction mechanism of the host’s molecular immune response to the inactivated vaccine is not yet fully understood, and the implementation of the third booster dose is also being actively discussed recently (Yigit et al., 2021) (https://www.cdc.gov/coronavirus/2019-ncov/vaccines/booster-shot.html). In addition, the impact of prior SARS-CoV-2 infection status on vaccination response is also worthy of further analysis. These insights may provide a theoretical basis for the determination of vaccination strategies and the allocation of vaccine resources.

The application of high-throughput technology to systematically scan the transcriptome response and evaluate changes in gene expression levels is very suitable for identifying immune response dynamics and gene regulatory networks. Previously, transcriptome analysis of Hantavax vaccine (Khan et al., 2019), influenza vaccine (Alcorn et al., 2020), VSV-EBOV vaccine (Menicucci et al., 2019) and BNT162b mRNA vaccine (Lee et al., 2021) have fully revealed the dynamics of the host immune response after vaccination. In this study, we characterized the PBMC transcriptome changes of SARS-CoV-2 recovered individuals receiving one dose of vaccine and healthy SARS-CoV-2 naïve individuals receiving one to three doses of vaccine respectively. This real-world study showed the changes of various cytokines and the regulation of immune pathways induced after vaccination, reveal the indispensable role of innate immune pathways, and reflect the key modules information of vaccine response in different individuals and different doses of vaccine.

Methods

Study Population and Recruitment

From January 2021, through August 2021, we recruited five SARS-CoV-2 recovered and twenty healthy SARS-CoV-2-naïve individuals to participate in this study in Linyi City, Shandong Province. Under different doses of the COVID-19 inactivated vaccine, anticoagulant and procoagulant blood samples were collected from participants. This study was approved by the Ethical Approval Committee of Shandong Center for Disease Control and Prevention.

Serum and PMBCs Isolation

The venous blood was collected from each participant to separate serum or isolate PBMCs. We separated sera by centrifugation at 2500 rpm/min for ten minutes and preserved at −80°C until testing. PBMCs were isolated by density-gradient sedimentation using Lymphoprep™ density gradients (Axis-Shield, Norway), frozen in cell saving media and stored in liquid nitrogen.

Qualitative SARS-CoV-2 IgG Detection

The IgG antibodies were detected using an indirect ELISA kit (Beijing Wantai Biological Pharmacy Enterprise Co, China) (Zhao et al., 2020b; Jiang et al., 2021) based on a recombinant nucleoprotein of SARS-CoV-2. The specific operation is carried out in strict accordance with the instructions. All serum samples were diluted to a concentration gradient of 1:10, 1:20, 1:40, 1:80, 1:160, 1:320, 1:640, 1:1280. The cut-off value for IgG is the mean OD value of three negative controls (if the mean absorbance value for three negative calibrators is < 0.03, take it as 0.03) + 0.16. A serum sample with an OD value ≥cut-off OD value was considered to be an anti-N IgG antibody positive. The geometric mean titer (GMT) of the antibody was calculated after three repetitions.

RNA-Seq

Total RNA was extracted from PBMCs by using the RNeasy Mini Kit (Qiagen, Germany) according to the manufacturer’s instructions. The concentration and integrity of total RNA were checked using the Qubit RNA Assay Kit in Qubit 4.0 Fluorometer (Life Technologies, USA) and the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 System (Agilent Technologies, USA) respectively. A total amount of 100 ng total RNA per sample was used to prepare for the rRNA-depleted cDNA library by Stranded Total RNA Prep Ligation with Ribo-Zero Plus Kit (Illumina, USA). The final library size and quality was evaluated using an Agilent High Sensitivity DNA Kit (Agilent Technologies, CA), and the fragments were found to be between 250 and 350 bp in size. The library was sequenced using an Illumina NextSeq 2000 platform to generate 100 bp paired-end reads.

Differentially Expressed Genes (DEGs) Identification

The quality control of RNA-seq raw fasta data is completed in the “prepare sequencing data” of CLC Genomics Workbench, and trimming of the base imbalance area is performed according to the quality control results. Processed data were mapped to the human reference genome hg38 using CLC Genomics Workbench. Gene expression level was measured based on the transcripts per million (TPM). We calculated normalization factors using iterative edgeR (Robinson et al., 2010) and limma (Ritchie et al., 2015) package, and the standardization and filtering of gene expression is accomplished by the voom, lmFit and eBayes functions in the limma package (Supplementary Figure 1). DEGs were filtered based on a p-value < 0.05 and a 2^logFC_cutoff. DEGs are visualized in the form of heat maps through the pheatmap (https://CRAN.R-project.org/package=pheatmap) package in R software.

Reverse Transcription Quantitative Polymerase Chain Reaction (RT−qPCR)

To validate the accuracy of DEGs analysis results, 2 µg RNA was transcribed into cDNA by Evo M-MLV RT Mix Kit with gDNA Clean for qPCR (Accurate Biotechnology). Relative expression levels of mRNA were quantitatively normalized and analyzed by RT-qPCR using SYBR^® Green Premix Pro Taq HS qPCR Kit (Accurate Biotechnology) against the expression of β-actin, using the 2^-ΔΔCt method. The primer sequences were derived from PrimerBank (https://pga.mgh.harvard.edu/primerbank/) (Table SI). The thermocycling conditions for RT-qPCR are presented in Table SII.

PPI Network Construction

The initial PPI network for the protein products of identified up and down regulated DEGs was constructed using the STRING Database (Szklarczyk et al., 2015) (STRING v11.5; https://string-db.org/), and then the network was visualized and analyzed with Cytoscape software (Shannon et al., 2003). The CentiScape and MCODE plugins in the Cytoscape software were used to extract the characteristic genes from the DEGs.

Pathway Enrichment Analysis

To be aware of the prospective functions of characteristic genes identified by the PPI network analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified using CluGO from Cytoscape software. P<0.05, subjected to Bonferroni adjustment, were defined as the cut−off criterion.

Weighted Gene Co-Expression Network Analysis (WGCNA)

We constructed the co-expression network through WGCNA package (Langfelder and Horvath, 2008) in R software. The best β value was confirmed with a scale-free fit index bigger than 0.85 as well as the highest mean connectivity by performing a gradient test from 1 to 30. Subsequently, we transformed the adjacency matrix into a topological overlap matrix (TOM). The TOM obtained was then clustered by dissimilarity between genes, and we performed hierarchical clustering to identify modules, each containing at least 30 genes (minModuleSize=30). Some modules were combined according to the correlation coefficient. We calculated the correlation between the modules and the sample characteristics to identify significant modules. The genes in significant modules were analyzed using cytoHubba and GeneMANIA in Cytoscape software.

Results

Characteristic of Samples

The median age of twenty-five participants was 40 years (interquartile range [IQR], 24–61), of which twelve (48%) were female and thirteen (52%) were male. They are divided into five groups according to their previous SARS-CoV-2 infection status and the difference in the number of doses of the COVID-19 inactivated vaccine. The Healthy-unvaccinated (H-un) group contains SARS-CoV-2 naïve individuals who have not been vaccinated with an inactivated vaccine; the Patient-1st (P-1st) group consists of COVID-19 patients who have received a single dose of vaccine and have recovered for 18 months; the Healthy-1st (H-1st) group includes SARS-CoV-2 naïve individuals who have received one dose of inactivated vaccine; the Healthy-2nd (H-2nd) group consists of SARS-CoV-2 naïve individuals who have received two doses of inactivated vaccine and the Healthy-3rd (H-3rd) group is SARS-CoV-2 naïve individuals who received the third booster vaccine half a year after receiving two doses of inactivated vaccine. We collected anticoagulant and procoagulant blood samples of participants fourteen days after vaccination. The specific details are shown in Figure 1.

FIGURE 1

Figure 1 The details of this study design. Participants’ previous SARS-CoV-2 infection status, vaccine dose information and blood sample time points are shown.

Detection of SARS-CoV-2 IgG in Serum Samples

In order to study the antibody response of serum samples to SARS-CoV-2 after vaccination, we used a fully validated commercial diagnostic ELISA kit to qualitatively measure IgG antibodies against N in serum samples diluted 1:10. The critical value for judging that the antibody level is positive is OD>0.19. The IgG antibodies of unvaccinated healthy participants were all negative. The GMT in healthy participant who received one dose of the vaccine was 1.58, whereas the GMT of recovered patients who received one dose of vaccine and healthy participants who received two or three doses of vaccine were 121.26, 80.00 and 226.27 respectively (Figure 2).

FIGURE 2

Figure 2 The antibody titer of IgG in the five groups of serum samples. Dots are antibody titres in serum samples. Numbers above the bars are geometric mean titre(GMT).

Identification of DEGs

The transcriptome data of five groups of PBMCs under different vaccination conditions were analyzed by limma and edgeR, and a total of 613 DEGs were found, of which 304 DEGs were up-regulated and 309 DEGs were down-regulated in the vaccination group. These results are shown in heatmap which is a simple yet effective way to compare the content of multiple major gene lists (Figure 3).

FIGURE 3

Figure 3 Hierarchical cluster analysis of DEGs in five groups. (A) The expression profile of 304 up-regulate DEGs. (B) The expression profile of 309 down-regulate DEGs. Each row represents mRNA and each column represents a sample. Orange indicates higher expression and blue indicates low expression in vaccination groups as compared with that in H-un group.

Prioritization of DEGs by PPI Network Analysis

The PPI network of these up-regulated and down-regulated DEGs were established using STRING database, and the node and edge relationships were imported into Cytoscape software for visualization, and then the DC evaluation was completed. The degree results of 304 up-regulated DEGs and 309 down-regulated DEGs were demonstrated by the size of the point and the thickness of the edge (Figures 4A and 5A). The first 61 up-regulated genes and 76 down-regulated genes were selected as characteristic genes (Figures 4B and 5B) as a focus for subsequent examination. Using MCODE’s k-core decomposition function, a variety of sub-networks were obtained, and two new networks were created, in which 23 up-regulated and 35 down-regulated differentially expressed genes were retained, respectively (Figures 4C and 5C). The results of DC and MCODE of up-regulated and down-regulated characteristic DEGs analysis are summarized in Tables 1 and 2. The differential expression of these characteristic DEGs in the four vaccination groups is shown in the heatmap (Figures 4D and 5D).

FIGURE 4

Figure 4 PPI networks of the up-regulated DEGs in vaccination groups. (A) The PPI network of all up-regulated DEGs, consisting of 304 nodes. (B) PPI network extracted from A by Degree Centrality, consisting of 61 nodes. (C) The PPI network extracted from B by MCODE, consisting of 32 nodes. PPI, protein-protein interaction; DEGs, differentially expressed genes. (D) The differential expression of 23 characteristic up-regulated DEGs in the four vaccination groups. Each row represents mRNA and each column represents a group. (E) A total of 35 KEGG pathways of the 23 characteristic up-regulated DEGs. (F–J) Expression levels of IL1B, CXCL8, IL10, JUN and VEGFA.

FIGURE 5

Figure 5 PPI networks of the down-regulated DEGs in vaccination groups. (A) The PPI network of all down-regulated DEGs, consisting of 309 nodes. (B) PPI network extracted from A by Degree Centrality, consisting of 76 nodes. (C) The PPI network extracted from B by MCODE, consisting of 35 nodes. PPI, protein-protein interaction; DEGs, differentially expressed genes. (D) The differential expression of 35 characteristic down-regulated DEGs in the four vaccination groups. Each row represents mRNA and each column represents a group. (E) A total of four KEGG pathways of the 35 characteristic down-regulated DEGs. (F–J) Expression levels of CXCL10, SELL, TLR9, CCR5, and HIST1H4F.

TABLE 1

Table 1 Results of Degree Centrality and MCODE analysis of up-regulated characteristic genes.

TABLE 2

Table 2 Results of Degree Centrality and MCODE analysis of down-regulated characteristic genes.

KEGG Pathway Enrichment Analyses of Characteristic Genes

The KEGG term enrichment analyses were performed using CluGO. There were 19 KEGG pathways with P<0.05 (including Bonferroni adjustment) identified within the analysis of the 23 characteristic up-regulated genes, including TNF signaling pathway, IL-17 signaling pathway, Viral protein interaction with cytokine and cytokine receptor, Cytokine-cytokine receptor interaction, Pertussis, NF-kappa B signaling pathway and other pathways (Figure 4E). In these pathways, genes that play important roles include chemokine, interleukin, macrophage colony-stimulating factor, interferon regulatory factor 1 and tumor necrosis factor alpha-induced protein 3. There were four KEGG pathways with P<0.05 (including Bonferroni adjustment) identified within the analysis of the 35 characteristic down-regulated genes, including Systemic lupus erythematosus, Alcoholism, Neutrophil extracellular trap formation, Viral protein interaction with cytokine and cytokine receptor (Figure 5E). Among these pathways, it was hypothesized that genes which serve a significant role include histone cluster, chemokine and interleukin-6 receptor subunit alpha (IL6R). Detailed enrichment information for the KEGG pathways is summarized in Tables 3, 4.

TABLE 3

Table 3 KEGG pathway of the 23 characteristic up-regulated genes.

TABLE 4

Table 4 KEGG pathways of the 35 characteristic down-regulated genes.

Verification of the Relative Expression Levels of the Top Five Up-Regulated and Down-Regulated Genes in Degree Centrality Value

A total of 10 genes, IL1B, CXCL8, IL10, JUN, VEGFA, CXCL10, SELL, TLR9, CCR5, and HIST1H4F, were selected from the important genes identified, and their differential expression levels between the vaccination and unvaccination groups was verified by RT-qPCR. Using the 2^-ΔΔCt method, the difference in expression levels of the 10 genes relative to β-actin was calculated, and graphpad prism 8 software was used to visually demonstrate the comparison between RT-qPCR results and transcriptome analysis results. It was identified that the experimental verification results were consistent with bioinformatics analysis data (Figures 4F–J and 5F–J). These genes were differentially expressed in vaccination group and are associated with the response mechanism of the host to inactivated vaccine, which requires additional study.

WGCNA: Identify Important Feature Modules Related to Sample Features

To explore the impact of different vaccine doses and SARS-CoV-2 infection status on the transcriptome, we made an analysis of the transcriptome using WGCNA. In this study, the power of β = 5 (scale-free R 2 = 0.85) was selected as the soft-thresholding parameter to ensure a scale-free network (Figures 6A, B). A total of 37 modules were identified through hierarchical clustering, and a dendrogram of all DEGs was clustered based on a dissimilarity measure (1-TOM) (Figures 6C, D). The correlation between the sample characteristics and the co-expression module is shown in Figure 6E, where the Brown (eigengene value = 0.69) module is significantly positively correlated with the vaccine dose. The darkorange (eigengene value = 0.91), yellow (eigengene value = 0.81) and steelblue (eigengene value = 0.71) modules are significantly correlated with the previous SARS-CoV-2 infection status (Figures 6F–I).

FIGURE 6

Figure 6 WGCNA of the PBMCs transcriptome. (A, B) Analysis of network topology for various soft-thresholding powers. The left panel shows the scale-free fit index (y-axis) as a function of the soft-thresholding power (x-axis). The right panel displays the mean connectivity (degree, y-axis) as a function of the soft-thresholding power (x-axis). (C) Clustering dendrogram of genes, with dissimilarity based on topological overlap, together with assigned module colors. (D) Heatmap depicts the Topological Overlap Matrix (TOM) of genes selected for weighted co-expression network analysis. Light color represents lower overlap and red represents higher overlap. (E) 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. (F–I) Scatter diagram for MM vs GS in the brown, darkorange, yellow and steelblue module.

Analysis of Important Modules and Identification of Genes

In order to analyze the network in the corresponding modules and determine the importance and function of genes, we combined CytoHubba and GeneMANIA in Cytoscape to analyze and identify the hub genes of the above four modules. The functions of the hub gene of the brown module mainly focus on activation of innate immune response, response to interleukin-1, and cellular response to tumor necrosis factor, etc. The hub genes of darkorange module act on regulation of tyrosine phosphorylation of STAT protein, tyrosine phosphorylation of STAT protein, pattern specification process, postsynaptic neurotransmitter receptor activity, neurotransmitter receptor activity pathway. The hub gene of the yellow module performs adaptive immune response, regulation of T cell activation, and B cell mediated immunity, etc. The hub gene of steelblue module can participate in humoral immune response, detection of biotic stimulus, antimicrobial humoral response, glycosaminoglycan binding, peptidyl-arginine modification and other pathways (Figures 7A–D).

FIGURE 7

Figure 7 PPI network construction. The PPI network of brown (A), darkorange (B), yellow (C) and steelblue (D) module. Different colors in each node represent different functions as indicated. Function of genes in dark node was not revealed through the analysis. Different colors in each edge represent different connection.

Discussion

The SARS-CoV-2 pandemic poses an imminent threat to humans, and the development and application of vaccines have brought hope to the effective prevention and control of the epidemic. In the field of rapid development and promotion of vaccines, the real data of vaccinated persons is crucial to the management and control of COVID-19 by government agencies and public health departments. At present, there is no research on the transcriptome response of COVID-19 recovered individuals and SARS-CoV-2- naïve individuals after vaccination with COVID-19 inactivated vaccine. In this study, we systematically revealed the molecular immune response and response mechanism of the vaccinated individuals to the inactivated vaccine.

Compared with the unvaccinated population, the up-regulated characteristic genes with substantial transcriptional changes in the vaccinated groups included IL1B, CXCL8, IL10, JUN, VEGFA, PTGS2, IL1A, CXCL1, CXCL2, CD80, etc. They participate in the activation of many innate immune pathways, including TNF signaling pathway, IL-17 signaling pathway, Viral protein interaction with cytokine and cytokine receptor, Cytokine-cytokine receptor interaction, Pertussis, NF-kappa B signaling pathway, etc. These transcriptional responses are highly consistent with the transcriptional regulation induced by inactivated vaccines of other pathogens. In research on inactivated influenza vaccines, the expression of IL-17 and NF-κB pathway genes and related innate immune activation have been observed (Alcorn et al., 2020); Rubella vaccine significantly affects the host’s antigen presentation and innate/inflammatory genome (Haralambieva et al., 2013); Hantavax vaccine can mediate Th17 cell differentiation, antigen processing and presentation, NF-kappa B signaling pathway, phenylalanine metabolism, phagosome, Fc gamma R-mediated phagocytosis (Khan et al., 2019). In the study of BNT162b2 mRNA vaccine, it was found that vaccination not only stimulated antiviral and interferon response modules, but also induced a broader innate immune response such as Toll-like receptor signaling, monocyte and neutrophil modules (Arunachalam et al., 2021). The above research results emphasize the key role of innate immune response in vaccination.

The characteristic genes that were significantly down-regulated in the vaccinated group included Histone cluster, chemokine and IL6R that participated in Systemic lupus erythematosus, Alcoholism, Neutrophil extracellular trap formation, Viral protein interaction with cytokine and cytokine receptor pathways. In a systems biology analysis of PBMC transcriptome data from COVID-19 patients, it was also found that key gene modules are not only involved in infection-related pathways, but also significantly enriched in Systemic lupus erythematosus and Alcoholism (Auwul et al., 2021). This coincidental consistency suggests that the role of these two pathways in COVID-19 needs to be further explored. In a study on Live-Attenuated Francisella tularensis vaccine, it was found that while a variety of immune pathways were activated early, some innate immune signaling pathways were inhibited, especially cytokine-cytokine receptor interactions. This result is related to the host escape mechanism of Francisella tularensis (Goll et al., 2020). Recent investigations suggested that in addition to the virus itself, dysregulated host immune response to SARS-CoV-2 may also contribute to the pathogenesis of COVID-19 (Blanco-Melo et al., 2020; Mehta et al., 2020; Qin et al., 2020). The inhibition of Neutrophil extracellular trap formation and Viral protein interaction with cytokine and cytokine receptor pathway caused by COVID-19 inactivated vaccine may also be related to the immune escape of SARS-CoV-2.

In the early days, a two-dose vaccination program was implemented in the SARS-CoV-2 naïve individuals, with an interval of 14 days between the two doses. Now we are trying to carry out a third booster vaccine half a year after the second vaccine. Our research data found that although the first dose of vaccine played a partial regulatory role at the transcriptome level, the GMT of IgG antibodies in serum samples was very low. After the second and third doses of vaccine, the GMT of IgG antibodies increased to 80.00 and 226.27 respectively. WGCNA analysis shows that with the increase of vaccine doses, it is significantly related to the brown modules involved in activation of innate immune response, response to interleukin-1, cellular response to tumor necrosis factor, Fc receptor signaling pathway, antigen processing and presentation of exogenous peptide antigen via MHC class I and innate immune response activating cell surface receptor signaling pathway. These functions are mainly performed by hub genes related to the 20S proteasome subunit, including PSMA1, PSMA4, PSMA5, PSMB1, PSMB10, PSMB2, PSMB3, PSMB4, PSMB5, ​​PSMB7, PSMB8, PSMB9, PSMD11, PSMD4, PSMD6, PSMD7, PSME4. The proteasome is a multicatalytic proteinase complex with a highly ordered ring-shaped 20S core structure. The 20S proteasome is formed by the 14 α subunits and the 14 β subunits, in which β1 (coded by the PSMB6 gene), β2 (coded by PSMB7 gene), and β5 (coded by PSMB5 gene) subunits possess protease activities (Arimochi et al., 2016). It is the catalytic core of the 26S proteasome and is distributed throughout eukaryotic cells at a high concentration (Voges et al., 1999). An essential function of a modified proteasome is the processing of class I MHC peptides. It can continue to decompose intracellular proteins that may appear in the process of viral infection, and produce antigen peptides, which can be combined with MHC1 molecules in the endoplasmic reticulum and transport the complex to the cell surface to play a role in antigen presentation (Lecker et al., 2006). Obviously, the function of the proteasome in the immune system is significant, and the third booster vaccine is beneficial and harmless.

A pivotal question in the process of vaccination promotion is whether individuals who have previously been infected with SARS-CoV-2 need to be vaccinated with multiple doses, because these people have developed a primary immune response to the virus during natural infection (Juno et al., 2020; Ahluwalia et al., 2021; Dan et al., 2021). This issue is particularly important in environments where vaccine supply is limited and vaccine deployment is challenging. The GMT levels of serum IgG antibodies in recovered patients who received one dose of the vaccine were between those of healthy participants who received two and three doses of the vaccine. Co-expression analysis of transcriptome data showed that darkorange, yellow, and steelblue modules were significantly related to SARS-CoV-2 infection status in people who received one dose of vaccine. The hub genes of these modules perform many adaptive immune regulation functions, including cellular immune responses such as the activation and regulation of T cells, and humoral immune responses mediated by B cells and antibodies. It suggests that a single dose of the inactivated vaccine elicits a strong immune response in population vaccinated 18 months after recovering from COVID-19. In some study of mRNA vaccines on SARS-CoV-2 naïve and recovered individuals, it was found that after a single dose of vaccine, recovered individuals can produce a strong immune response, including serum antibody response, B cell response, and transcriptome changes (Ebinger et al., 2021; Goel et al., 2021; Lee et al., 2021). In view of the relatively short duration of this study, future studies will be necessary to evaluate the impact of multiple doses of vaccine on the long-term immune response of recovered individuals.

Our data provides quantitative and global insights into the dynamic changes of immune-related genes, reveals the differential expression of cytokines and inflammatory factors, as well as the regulation of antigen presentation and innate immune pathways, and helps us explain the effects of molecular immune mechanism induced by inactivated vaccines. The transcriptome characteristics of individuals with different SARS-CoV-2 infection and different doses may provide a new reference for the determination of vaccination strategies in the future. However, the results of this study rely on bioinformatics analysis, the sample size is small, and the long-term evaluation of the inoculation effect is lacking. We intend to further expand the sample size and conduct long-term follow-up to explore comprehensive insights into the inactivated COVID-19 vaccine.

Data Availability Statement

The original contributions presented in the study are publicly available in the GEO database (NCBI) under accession number GSE189263.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethical Approval Committee of Shandong Center for Disease Control and Prevention. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

XJ and YZ conceived the study. YZ devised the statistical methodology. CL and LL recruited participants and collected samples. YZ, XG, LL, MY, and BP did the formal analysis. XJ, YZ, and ZK were project administrators. XZ and QD curated and validated the data. XG, XT, and YX did the literature review. XJ acquired the funding. YZ and XG wrote the original draft of the manuscript. XJ reviewed and edited the manuscript. XJ and YZ have shared responsibility for the decision to submit for publication. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Major Scientific and Technological Innovation Project in Shandong Province (grant numbers: 2020SFXGFY02-1).

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.

Acknowledgments

This article was submitted as a preprint to bioRxiv as Zhang et al. (2021).

Supplementary Material

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

Supplementary Figure 1 | The standardization and filtering processes of gene expression. (A) The fitting relationship between the average expression of unprocessed raw data (red line) and Sqrt. (B) The fitting relationship between the average expression of normalized and filtered data (red line) and Sqrt. (C) The linear fitting model constructed by Bayesian method. The horizontal blue line represents the average log2 residual standard deviation.

References

Ahluwalia, P., Vaibhav, K., Ahluwalia, M., Mondal, A. K., Sahajpal, N., Rojiani, A. M., et al. (2021). Infection and Immune Memory: Variables in Robust Protection by Vaccines Against SARS-CoV-2. Front. Immunol. 12, 660019. doi: 10.3389/fimmu.2021.660019

PubMed Abstract | CrossRef Full Text | Google Scholar

Alcorn, J. F., Avula, R., Chakka, A. B., Schwarzmann, W. E., Nowalk, M. P., Lin, C. J., et al. (2020). Differential Gene Expression in Peripheral Blood Mononuclear Cells From Children Immunized With Inactivated Influenza Vaccine. Hum. Vaccines Immunother. 16 (8), 1782–1790. doi: 10.1080/21645515.2020.1711677

CrossRef Full Text | Google Scholar

Arimochi, H., Sasaki, Y., Kitamura, A., Yasutomo, K. (2016). Dysfunctional Immunoproteasomes in Autoinflammatory Diseases. Inflamm. Regen. 36, 13. doi: 10.1186/s41232-016-0011-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Arunachalam, P. S., Scott, M. K. D., Hagan, T., Li, C., Feng, Y., Wimmers, F., et al. (2021). Systems Vaccinology of the BNT162b2 mRNA Vaccine in Humans. Nature 596 (7872), 410–416. doi: 10.1038/s41586-021-03791-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Auwul, M. R., Rahman, M. R., Gov, E., Shahjaman, M., Moni, M. A. (2021). Bioinformatics and Machine Learning Approach Identifies Potential Drug Targets and Pathways in COVID-19. Briefings Bioinf. 22 (5), bbab120. doi: 10.1093/bib/bbab120

CrossRef Full Text | Google Scholar

Blanco-Melo, D., Nilsson-Payant, B. E., Liu, W. C., Uhl, S., Hoagland, D., Moller, R., et al. (2020). Imbalanced Host Response to SARS-CoV-2 Drives Development of COVID-19. Cell 181 (5), 1036–1045.e9. doi: 10.1016/j.cell.2020.04.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Dan, J. M., Mateus, J., Kato, Y., Hastie, K. M., Yu, E. D., Faliti, C. E., et al. (2021). Immunological Memory to SARS-CoV-2 Assessed for Up to 8 Months After Infection. Science 371 (6529), eabf4063. doi: 10.1126/science.abf4063

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebinger, J. E., Fert-Bober, J., Printsev, I., Wu, M., Sun, N., Prostko, J. C., et al. (2021). Antibody Responses to the BNT162b2 mRNA Vaccine in Individuals Previously Infected With SARS-CoV-2. Nat. Med. 27 (6), 981–984. doi: 10.1038/s41591-021-01325-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Goel, R. R., Apostolidis, S. A., Painter, M. M., Mathew, D., Pattekar, A., Kuthuru, O., et al. (2021). Distinct Antibody and Memory B Cell Responses in SARS-CoV-2 Naive and Recovered Individuals Following mRNA Vaccination. Sci. Immunol. 6 (58), eabi6950. doi: 10.1126/sciimmunol.abi6950

PubMed Abstract | CrossRef Full Text | Google Scholar

Goll, J. B., Li, S., Edwards, J. L., Bosinger, S. E., Jensen, T. L., Wang, Y., et al. (2020). Transcriptomic and Metabolic Responses to a Live-Attenuated Francisella Tularensis Vaccine. Vaccines 8 (3), 412. doi: 10.3390/vaccines8030412

CrossRef Full Text | Google Scholar

Haralambieva, I. H., Oberg, A. L., Ovsyannikova, I. G., Kennedy, R. B., Grill, D. E., Middha, S., et al. (2013). Genome-Wide Characterization of Transcriptional Patterns in High and Low Antibody Responders to Rubella Vaccination. PloS One 8 (5), e62149. doi: 10.1371/journal.pone.0062149

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C., Wang, Y., Li, X., Ren, L., Zhao, J., Hu, Y., et al. (2020). Clinical Features of Patients Infected With 2019 Novel Coronavirus in Wuhan, China. Lancet 395 (10223), 497–506. doi: 10.1016/S0140-6736(20)30183-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, X. L., Wang, G. L., Zhao, X. N., Yan, F. H., Yao, L., Kou, Z. Q., et al. (2021). Lasting Antibody and T Cell Responses to SARS-CoV-2 in COVID-19 Patients Three Months After Infection. Nat. Commun. 12 (1), 897. doi: 10.1038/s41467-021-21155-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Juno, J. A., Tan, H. X., Lee, W. S., Reynaldi, A., Kelly, H. G., Wragg, K., et al. (2020). Humoral and Circulating Follicular Helper T Cell Responses in Recovered Patients With COVID-19. Nat. Med. 26 (9), 1428–1434. doi: 10.1038/s41591-020-0995-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Shin, O. S., Na, J., Kim, J. K., Seong, R. K., Park, M. S., et al. (2019). A Systems Vaccinology Approach Reveals the Mechanisms of Immunogenic Responses to Hantavax Vaccination in Humans. Sci. Rep. 9 (1), 4760. doi: 10.1038/s41598-019-41205-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., Horvath, S. (2008). WGCNA: An R Package for Weighted Correlation Network Analysis. BMC Bioinf. 9, 559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

Lecker, S. H., Goldberg, A. L., Mitch, W. E. (2006). Protein Degradation by the Ubiquitin-Proteasome Pathway in Normal and Disease States. J. Am. Soc. Nephrol.: JASN 17 (7), 1807–1819. doi: 10.1681/ASN.2006010083

CrossRef Full Text | Google Scholar

Lee, H. K., Knabl, L., Knabl, L., Kapferer, S., Pateter, B., Walter, M., et al. (2021). Robust Immune Response to the BNT162b mRNA Vaccine in an Elderly Population Vaccinated 15 Months After Recovery From COVID-19. MedRxiv: Preprint Server Health Sci. 09.08.21263284. doi: 10.1101/2021.09.08.21263284

CrossRef Full Text | Google Scholar

Mehta, P., McAuley, D. F., Brown, M., Sanchez, E., Tattersall, R. S., Manson, J. J., et al. (2020). COVID-19: Consider Cytokine Storm Syndromes and Immunosuppression. Lancet 395 (10229), 1033–1034. doi: 10.1016/S0140-6736(20)30628-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Menicucci, A. R., Jankeel, A., Feldmann, H., Marzi, A., Messaoudi, I. (2019). Antiviral Innate Responses Induced by VSV-EBOV Vaccination Contribute to Rapid Protection. mBio 10 (3), e00597–619. doi: 10.1128/mBio.00597-19

PubMed Abstract | CrossRef Full Text | Google Scholar

Palacios, R., Patino, E. G., de Oliveira Piorelli, R., Conde, M., Batista, A. P., Zeng, G., et al. (2020). Double-Blind, Randomized, Placebo-Controlled Phase III Clinical Trial to Evaluate the Efficacy and Safety of Treating Healthcare Professionals With the Adsorbed COVID-19 (Inactivated) Vaccine Manufactured by Sinovac - PROFISCOV: A Structured Summary of a Study Protocol for a Randomised Controlled Trial. Trials 21 (1), 853. doi: 10.1186/s13063-020-04775-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, C., Zhou, L., Hu, Z., Zhang, S., Yang, S., Tao, Y., et al. (2020). Dysregulation of Immune Response in Patients With Coronavirus 2019 (COVID-19) in Wuhan, China. Clin. Infect. Dis.: Off. Publ. Infect. Dis. Soc. Am. 71 (15), 762–768. doi: 10.1093/cid/ciaa248

CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). Edger: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data. Bioinformatics 26 (1), 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 13 (11), 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Sharma, O., Sultan, A. A., Ding, H., Triggle, C. R. (2020). A Review of the Progress and Challenges of Developing a Vaccine for COVID-19. Front. Immunol. 11, 585354. doi: 10.3389/fimmu.2020.585354

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING V10: Protein-Protein Interaction Networks, Integrated Over the Tree of Life. Nucleic Acids Res. 43 (Database issue), D447–D452. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Voges, D., Zwickl, P., Baumeister, W. (1999). The 26S Proteasome: A Molecular Machine Designed for Controlled Proteolysis. Annu. Rev. Biochem. 68, 1015–1068. doi: 10.1146/annurev.biochem.68.1.1015

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Zhao, S., Yu, B., Chen, Y. M., Wang, W., Song, Z. G., et al. (2020). A New Coronavirus Associated With Human Respiratory Disease in China. Nature 579 (7798), 265–269. doi: 10.1038/s41586-020-2008-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, S., Duan, K., Zhang, Y., Zhao, D., Zhang, H., Xie, Z., et al. (2020). Effect of an Inactivated Vaccine Against SARS-CoV-2 on Safety and Immunogenicity Outcomes: Interim Analysis of 2 Randomized Clinical Trials. Jama 324 (10), 951–960. doi: 10.1001/jama.2020.15543

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, S., Zhang, Y., Wang, Y., Wang, H., Yang, Y., Gao, G. F., et al. (2021a). Safety and Immunogenicity of an Inactivated SARS-CoV-2 Vaccine, BBIBP-CorV: A Randomised, Double-Blind, Placebo-Controlled, Phase 1/2 Trial. Lancet Infect. Dis. 21 (1), 39–51. doi: 10.1016/S1473-3099(20)30831-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Xia, S., Zhang, Y., Wang, Y., Wang, H., Yang, Y., Gao, G. F., et al (2021b). Safety and Immunogenicity of an Inactivated COVID-19 Vaccine, BBIBP-CorV, in People Younger Than 18 Years: A Randomised, Double-Blind, Controlled, Phase 1/2 Trial. Lancet Infect. Dis. S1473-3099(21), 00462–X. doi: 10.1016/S1473-3099(21)00462-X

CrossRef Full Text | Google Scholar

Yigit, M., Ozkaya-Parlakay, A., Cosgun, Y., Ince, Y. E., Bulut, Y. E., Senel, E. (2021). Should a Third Booster Dose be Scheduled After Two Doses of CoronaVac? A Single-Center Experience. J. Med. Virol. 94 (1), 287–290. doi: 10.1002/jmv.27318

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zeng, G., Pan, H., Li, C., Hu, Y., Chu, K., et al. (2021). Safety, Tolerability, and Immunogenicity of an Inactivated SARS-CoV-2 Vaccine in Healthy Adults Aged 18-59 Years: A Randomised, Double-Blind, Placebo-Controlled, Phase 1/2 Clinical Trial. Lancet Infect. Dis. 21 (2), 181–192. doi: 10.1016/S1473-3099(20)30843-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. W., Guo, X. Y., Li, C. B., Kou, Z. Q., Lin, L. F., Yao, M. X., et al. (2021). Transcriptome Analysis of SARS-CoV-2 Naïve and Recovered Individuals Vaccinated With Inactivated Vaccine. bioRxiv 2021.11.05.467537. doi: 10.1101/2021.11.05.467537

CrossRef Full Text | Google Scholar

Zhao, J., Yuan, Q., Wang, H., Liu, W., Liao, X., Su, Y., et al. (2020b). Antibody Responses to SARS-CoV-2 in Patients With Novel Coronavirus Disease 2019. Clin. Infect. Dis.: Off. Publ. Infect. Dis. Soc. Am. 71 (16), 2027–2034. doi: 10.1093/cid/ciaa344

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, J., Zhao, S., Ou, J., Zhang, J., Lan, W., Guan, W., et al. (2020a). COVID-19: Coronavirus Vaccine Development Updates. Front. Immunol. 11, 602256. doi: 10.3389/fimmu.2020.602256

CrossRef Full Text | Google Scholar