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

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 b-actin, using the 2 -DDCt 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://stringdb.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 b 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.

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.

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

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

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

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.

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 -DDCt method, the difference in expression levels of the 10 genes relative to b-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 b = 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 coexpression module is shown in Figure 6E

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,

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 upregulated 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-kB 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 a subunits and the 14 b subunits, in which b1 (coded by the PSMB6 gene), b2 (coded by PSMB7 gene), and b5 (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).

ACKNOWLEDGMENTS
This article was submitted as a preprint to bioRxiv as .