Transcriptome Analysis of Peripheral Blood Mononuclear Cells Reveals Distinct Immune Response in Asymptomatic and Re-Detectable Positive COVID-19 Patients

The existence of asymptomatic and re-detectable positive coronavirus disease 2019 (COVID-19) patients presents the disease control challenges of COVID-19. Most studies on immune responses in COVID-19 have focused on moderately or severely symptomatic patients; however, little is known about the immune response in asymptomatic and re-detectable positive (RP) patients. Here we performed a comprehensive analysis of the transcriptomic profiles of peripheral blood mononuclear cells (PBMCs) from 48 COVID-19 patients which included 8 asymptomatic, 13 symptomatic, 15 recovered and 12 RP patients. The weighted gene co-expression network analysis (WGCNA) identified six co-expression modules, of which the turquoise module was positively correlated with the asymptomatic, symptomatic, and recovered COVID-19 patients. The red module positively correlated with symptomatic patients only and the blue and brown modules positively correlated with the RP patients. The analysis by single sample gene set enrichment analysis (ssGSEA) revealed a lower level of IFN response and complement activation in the asymptomatic patients compared with the symptomatic, indicating a weaker immune response of the PBMCs in the asymptomatic patients. In addition, gene set enrichment analysis (GSEA) analysis showed the enrichment of TNFα/NF-κB and influenza infection in the RP patients compared with the recovered patients, indicating a hyper-inflammatory immune response in the PBMC of RP patients. Thus our findings could extend our understanding of host immune response during the progression of COVID-19 disease and assist clinical management and the immunotherapy development for COVID-19.


INTRODUCTION
Coronavirus disease 2019 (COVID- 19) is an infectious disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), which was declared a pandemic by the WHO on 11 March 2020 (1). As of May 1, 2021, 153 million cases have been reported across 188 countries and territories with more than 3.20 million deaths. Disease in most COVID-19 patients is classified as mild; however, about 20% become seriously ill and require hospitalization due to pneumonia (2). The primary symptoms of COVID-19 are listed as fever, dry cough, and shortness of breath but also include other symptoms such as diarrhea, loss of taste and smell, and rashes (1). However, increasing evidence has shown that individuals infected with SARS-CoV-2 can be asymptomatic and simultaneously silent spreaders of SARS-CoV-2 (3,4). Identification and isolation of asymptomatic patients as early as possible is critical in controlling COVID-19 outbreaks. According to a recent study, four medical workers aged from 30 to 36 years were re-detectable positive (RP) for SARS-CoV-2 within 5-13 days after recovery, suggesting that these four recovered patients may still have been virus carriers (5). All four patients continued to be asymptomatic upon clinical examination, and chest CT findings showed no change from previous images. RP cases have also been reported by several other studies. For example, one report showed that 10.99% of patients (20/182) were SARS-CoV-2 RNA re-positive, and none of them showed any clinical symptomatic recurrence (6). Both asymptomatic and RP patients are receiving more attention because of their potential infectivity.
As there are no effective drugs available at this moment against SARS-CoV-2 and most patients receive symptomatic treatment, it is essential to understand the host-pathogen biology of COVID-19, which will provide the novel insights into more targeted therapeutic approaches for this disease. Over the past few months, studies have reported the characteristics of innate and adaptive immune responses in SARS-CoV-2 infection, which have helped us to understand the potential pathogenesis of COVID-19 (7,8). Most of these studies focused on the peripheral immune response of symptomatic patients, particularly those with severe COVID-19 (9,10). A dysfunctional immune response was observed in patients with severe COVID-19, which triggered a cytokine storm causing severe lung and even systemic pathology (11). However, the immunological features and the molecular mechanisms involved in asymptomatic and RP patients still remains elusive. Here we report the transcriptome profiles of peripheral blood mononuclear cells (PBMCs) from COVID-19 patients including asymptomatic, symptomatic, recovered, and RP groups. We performed a comprehensive bioinformatics analysis of these transcriptome profiles using multiple methods including weighted gene coexpression network analysis (WGCNA) (12), single sample gene set enrichment analysis (ssGSEA) (13), and gene set enrichment analysis (GSEA) (14). Our study revealed a weaker peripheral immune response in the asymptomatic group and a hyperinflammatory immune response in the RP group.

Patients
This study was reviewed and approved by the ethics committees of the four hospitals including Guangzhou Eighth People's Hospital, Shunde Hospital of Guangzhou University of Chinese Medicine, Fourth People's Hospital of Foshan, and First People's Hospital of Foshan. Written informed consent was obtained from all participants enrolled in this study. The study enrolled 48 patients with confirmed SARS-CoV-2 admitted to the above four hospitals from March through May in 2020 and classified into four groups. The asymptomatic group (8 cases) included patients who did not have self-perceived or clinically recognizable symptoms but were diagnosed as having COVID-19 according to the Protocol for Prevention and Control of COVID-19 (Edition 6) of the National Health Commission of China. According to the Guidelines for the Diagnosis and Treatment of COVID-19 (7th edition), the symptomatic group (13 cases) was defined as patients with evident clinical symptoms, and the recovered group (15 cases) was defined as patients meeting the discharge standards. The RP group (12 cases) was composed of those with re-positive results of SARS-CoV-2 nucleic acid during the follow-up period after discharge from hospital. The time between the re-positive results of the PCR test and discharge from hospital ranged from 2 weeks to 1 month. In addition, we enrolled 22 healthy donors from Shunde Hospital of Guangzhou University of Chinese Medicine. All patients had routine laboratory investigations, including complete blood count, liver function tests, blood gas analysis, and coagulation tests. The PBMC samples were collected within 4 days of admission from asymptomatic, symptomatic and RP groups to maintain consistent timing for comparison between groups.

RNA Extraction and Library Construction
Total RNA was isolated and purified using TRIzol (Life Services Inc., Glendale, CA, USA) following the manufacturer's procedure. After the quality inspection of Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and NanoPhotometer ® (IMPLEN GmbH, Munich, Germany), poly (A) RNA was purified from 1mg total RNA using VAHTS ® mRNA Capture Beads with Oligo (dT) (Vazyme, Nanjing, China) using two rounds of purification. Then the poly(A) RNA was fragmented into small pieces using VAHTS ® Universal V6 RNA-seq Library Prep Kit (Vazyme, Nanjing, China) at 94°C for 8 min. Then the cleaved RNA fragments were reverse-transcribed to create the cDNA with reverse transcription reagent, which was next used to synthesize Ulabeled second-stranded DNAs. An A-base was then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contained a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. After the heat-labile UDG enzyme treatment of the U-labeled secondstranded DNAs, size selection was performed with VAHTS ® DNA Clean Beads (Vazyme, Nanjing, China). Then the ligated products were e amplified with PCR. Finally, we performed the 2×150bp paired-end sequencing (PE150) on an Illumina Novaseq ™ 6000 (Illumina, Inc., San Diego, CA, USA) following the manufacturer's recommended protocol.

Data Processing and Analysis
Cutadapt software (version 1.16) was used to remove the reads containing adaptor contamination. Then we prepared CleanData by removing the low quality bases and undetermined bases. The clean reads were mapped to the hg19 UCSC transcript set using Bowtie2 version 2.1.0 and the gene expression level was estimated using RSEM v1.2.15. 24 immunologically relevant gene sets were collected from the ImmPort database and molecular signature database. The detailed lists of the gene sets are shown in the Supplementary Table 1. Gene sets were scored using single-sample gene set enrichment (ssGSEA) analysis, as implemented in the GSVA R package. The statistical analyses were performed using the unpaired, two-sided Mann-Whitney U test. All statistical analyses were performed using R software (version 4.0.2).

Construction of Co-Expression Modules
First, flashClust function was used to carry out cluster analysis of the samples with the appropriate threshold value to detect and remove the outliers. Then, the soft thresholding power b−value was screened during module construction by the pickSoftThreshold function in the WGCNA package. Once the appropriate power value had been determined when the degree of independence was 0.8, the WGCNA algorithm was performed to construct co −expression networks (modules); the minimum module size was set to 30. Co−expression networks or modules were defined as branches of a hierarchical clustering tree, and each module was assigned a unique color label. Subsequently, the information pertaining to the corresponding genes in each module was extracted. The ClusterProfiler R package was used to perform the GO enrichment analysis of the genes in the modules.

Construction of Module-Trait Relationships
The WGCNA algorithm utilizes module eigengenes (MEs) to assess the potential correlation of gene modules with clinical traits. The module eigengene was defined as the first principal component of the expression matrix for a given module. The module eigengene can be considered an average gene expression level for all genes in each module. The disease conditions of COVID-19 patients were converted into numerical values and a regression analysis was performed between the module eigengene values and the disease conditions of COVID-19 patients.

Differences in Immune Cell Abundance Across Disease Conditions
To investigate the immunological features of asymptomatic and RP patients with COVID-19, we performed mRNA sequencing to study the transcriptomic profiles of PBMCs from 48 COVID-19 patients which include 8 asymptomatic, 13 symptomatic, 15 recovered and 12 RP patients. There were 22 healthy donors included as the control group. The demographic, clinical, and laboratory characteristics of enrolled patients are summarized in Table 1. The laboratory results revealed a significant lower level of lymphocytes/leukocytes in symptomatic patients compared with the healthy donors (Figure 1), which is consistent with a previous report of lymphocytopenia/leukopenia in COVID-19 patients (2). In addition, eosinophil and basophil counts were also significantly lower than in symptomatic patients compared with the healthy donors, although no significant differences were identified for monocyte and neutrophil counts ( Figure 1). As shown in Figure 1, the lymphocyte counts in asymptomatic patients were higher than in symptomatic patients, although the difference was not significant because of the small sample size. The basophils counts in asymptomatic patients were significantly higher compared with the symptomatic patients. Interestingly, a significantly lower level of monocytes in RP patients was observed when compared with the recovered patients and it was comparable for other blood cell counts between RP and recovered patients.

Identification of Gene Co-Expression Modules
To clarify the critical modules and hub genes that were associated with different status of COVID-19 patients, we performed a weighted gene co-expression network construction analysis (WGCNA) in the 70 samples. Because non-varying genes usually represent noise, we selected the top 3000 varied genes for constructing the weighted gene co-expression network. The cluster analysis by flashClust revealed no outlier samples, and all 70 samples were included in further analysis (Supplementary Figure 1A). The power of b was set at 7 to ensure a scale-free network (Supplementary Figure 1B). The topological overlap measure (TOM) for each gene pair was then calculated. Hierarchical clustering analysis based on the TOM dissimilarity measure (1-TOM) revealed six modules (Supplementary Figure 1C). Each module contained a group of genes with high TOM that were coordinately expressed and potentially involved in similar biological processes. Each module was assigned a unique color to distinguish the individual modules.

Association of Modules With Disease Conditions
To investigate the association of co−expression modules with the disease conditions of in COVID19 patients, we calculated the correlations between the module eigengenes and the disease conditions in COVID19 patients. As shown in Figure 2A, the turquoise module displayed strong positive correlation with the asymptomatic, symptomatic and recovered group, whereas the red module showed a significant positive correlation with the symptomatic group only. Furthermore, both blue and brown modules had a significant positive correlation with the RP group. Consistent with the correlation data, a heatmap revealed higher expression of the turquoise module genes in the asymptomatic, symptomatic and recovered groups and higher expression of the red module genes in the symptomatic group ( Figure 2B).

Functional Enrichment Analysis of Genes in Modules of Interest
To determine the biological function of the four modules of interest, we employed GO enrichment analysis of the genes in each module. As shown in Figure 3A, genes in the turquoise module were enriched in response to the virus, response to type I interferon and histone binding. For the red module, the genes were mainly enriched in antigen binding, humoral immune response, and complement activation ( Figure 3B). Finally, genes in blue module were enriched in response to bacterium and cell chemotaxis and genes in the brown module were enriched in response to lipopolysaccharide, hemoglobin complex and cytokine activity (Figures 3C, D).
To investigate in more detail the interplay between genes in the above four modules, we used cytoscape software to visualize the network of targeted modules and the intra-modular connectivity. The top 40 genes with the highest intra-modular connectivity were used for the network visualization. Most of these genes in the turquoise module are related to response to virus and response to type I interferon, which include TLR7, TLR8, CXCL10, IL27, IFI27, IFI44, and TRIM6, et al. (Figure 3A). These top genes in the red module are related to humoral immune response, including IGHG1, IGHA1, IGLC2, IGLV1-47, and IGKV3-20, et al. (Figure 3B). For the blue module, the top 40 genes include the CXCL1, CSF3R, EGR3, IL1R2, and TREM1, et al, which are related to cell chemotaxis and cytokine receptor activity ( Figure 3C). For the brown module, the top 40 genes include JUNB, TNF, CXCL2, CXCL3, CCL3, CXCL16, and FOSL1, et al, which are related to response to lipopolysaccharide and response to TNFa ( Figure 3D).

Differences in Immune Responses Across Disease Conditions
To identify the difference of immune-related function among these five groups, we collected 24 immunologically relevant gene sets including 17 from the ImmPort database (15) and 7 from the   Table 1). We used the ssGSEA score (13,17) to quantify the activity or enrichment levels of immune-related functions or pathways in these PBMCs samples. Multiple immune-related functions were dysregulated in PBMCs from COVID-19 patients. The ssGSEA analysis revealed a dynamic change in the interferon response among these five groups. As expected, both the interferon-alpha response and interferongamma response activity are significantly higher in the symptomatic group compared with the healthy donors ( Figures 4A, B). Interestingly, we observed that the levels of interferon response in the asymptomatic group were comparable with those of the healthy group and significantly lower than in the symptomatic group, indicating a lower level of IFN in the serum of asymptomatic patients. Natural killer (NK) cells are an important arm of the innate lymphocytic antiviral response. The level of NK cell cytotoxicity was significantly lower in both asymptomatic and symptomatic patients when compared with the healthy group ( Figure 4C). The RP patients had a significantly higher level of NK cell cytotoxicity than the recovered patients, indicating a restoration of the dysregulated immune system. In addition, both the symptomatic and asymptomatic patients had a significantly lower level of antigen processing and presentation than the healthy donors, suggesting a defective antigen presentation in the COVID-19 patients ( Figure 4D). Interestingly, we also observed the RP patients had antigen presentation restored to the normal level. Furthermore, the level of cytokines, chemokines, and interleukins in the PBMCs showed a similar pattern as the antigen processing and presentation ( Figures 4E, F and Supplementary Figure 2). The lower expression of cytokines in symptomatic patients is consistent with previous reports of undetectable expression of most cytokines in the PBMCs of COVID-19 patients (7,18). The lower level of cytokines in PBMC indicates the elevated serum cytokines largely arise from the local infection site instead of the peripheral blood. On the contrary, the levels of cytokines, chemokines, and interleukins in the RP group are significantly higher than in other groups of COVID-19 patients, indicating a hyper-inflammatory immune response in the PBMC of RP patients.
To further elucidate the up-and down-regulated biological pathways in asymptomatic verse symptomatic, RP verse recovered and RP verse healthy groups, we performed the GSEA analysis (14) to analyze the pathways from KEGG, Reactome and Biocarta as well as the hallmark gene sets from MSigDB. We observed the downregulation of IFN response and complement activation ("Creation of C4 and C2 activators") in the asymptomatic patients ( Figures 5A,  D), indicating a weaker immune response of the PBMCs in asymptomatic patients. Interestingly, some cell cycle related pathways were down-regulated and NGF related pathways upregulated in the asymptomatic patients. The potential role of these pathways needs further investigation. In addition, both TNFa/NF-kB and influenza infection activity were enriched in the RP patients compared with the recovered patients ( Figures 5B, E). Further, as shown in Figures 5C, F, the TNFa/NF-kB and inflammatory response activity were significantly enriched in the RP patients compared with the healthy donors. NF-kB is one of the hallmark signaling factors activated by influenza infection (19). Moreover, a comparative study of COVID-19 and influenza showed the activation of STAT3/NF-kB in PBMCs of patients infected with influenza A virus (IAV)vs. the activation of STAT1/IRF3 in COVID-19 patients (18). Taken together, these findings suggest PBMCs in RP patients have a hyper-inflammatory, flu-like immune response.

DISCUSSION
COVID-19 patients are characterized by a broad spectrum of disease severity ranging from asymptomatic to critically severe (20). The heterogeneous manifestation of the disease may be due to the patients' different immune response to SARS-CoV-2. Dysregulated immune responses have been described in symptomatic COVID-19 patients, particularly those with severe disease (9,21). However, little is known about the immune response in the asymptomatic and RP patients.
Here, we performed the transcriptome analysis on PBMCs from 48 COVID-19 patients in different phases including asymptomatic and RP patients, as well as 22 healthy donors. We found the level of IFN response and complement activation in asymptomatic patients was lower than in the symptomatic patients. This observation is agreement with the data showing lower levels of 18 cytokines including IFN in the serum of the asymptomatic compared with the symptomatic patients (22). It was reported that plasma IFN levels were associated with the COVID-19 disease severity (23) and a recent longitudinal analysis showed that IFNa in peripheral blood was sustained at high levels in patients with severe COVID-19 (24). A robust type I interferon response could exacerbate hyperinflammation in the progression to severe COVID-19 through diverse mechanisms (25).
WGCNA analysis revealed the genes in the red co-expression module were enriched in the humoral immune response. The red module was found to be positively correlated with the symptomatic patients but not asymptomatic patients, suggesting a high level of humoral immune response in the symptomatic patients. Consistently, two recent studies of serological responses in  COVID-19 patients revealed that asymptomatic patients did not respond or have lower antibody levels upon SARS-CoV-2 infection (26,27). Moreover, the asymptomatic patients were reported to exhibit reduced proportions of SARS-CoV-2 specific T cells, suggesting a weaker cell-mediated immune response (27). Collectively, the well-controlled immune response in the asymptomatic patients may protect the patients from progressing to the inflammatory secondary phase of the disease. Several recent studies have reported the existence of RP patients and the underlying mechanism of RP occurrence remains unknown (5,28). The potential reasons might be related to some factors including virology, immunology and sampling methodology. Factors such as false negatives in detection (29), viral residuals (30), intermittent viral release (31) and viral distribution (32) are usually considered major reasons. However, to the best of our knowledge, there is no report on the peripheral immune response in the RP patients. In this study, we identified a hyper-inflammatory immune response in the peripheral blood of RP patients. Several limitations should be taken into consideration for interpreting this study. First, this study focused on the transcriptome analysis of PBMCs in blood. Adding data on immune cells from lesion sites such as the lung and bronchoalveolar lavage fluid will make the analysis more comprehensive and conclusive. Second, the PBMC samples were collected within 4 days of admission from asymptomatic and symptomatic groups to maintain uniformity of timing for comparison between groups. Future studies with longitudinal samples from COVID-19 patients may help to determine the cause-and-effect relationships between immune response and disease outcome.
In summary, this study will help us extend our understanding of host immune response during the progression of COVID-19 disease, and may help elucidate the COVID-19 infection and provide a basis for rationally designed immuno therapies.

DATA AVAILABILITY STATEMENT
The datasets generated and analyzed during the current study are available from the corresponding author on request. RNA-seq data will be stored at NCBI/GEO, Accession No. GSE179627.

ETHICS STATEMENT
Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.