Screening the hub genes and analyzing the mechanisms in discharged COVID-19 patients retesting positive through bioinformatics analysis.

BACKGROUND: After encountering COVID-19 patients who test positive again after discharge, our study analyzed the pathogenesis to further assess the risk and possibility of virus reactivation.
METHODS: A separate microarray was acquired from the Gene Expression Omnibus (GEO), and its samples were divided into two groups: a "convalescent-RTP" group consisting of convalescent and "retesting positive" (RTP) patients (group CR) and a "healthy-RTP" group consisting of healthy control and RTP patients (group HR). The enrichment analysis was performed with R software, obtaining the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). Subsequently, the protein-protein interaction (PPI) networks of each group were established, and the hub genes were discovered using the cytoHubba plugin.
RESULTS: In this study, 6622 differentially expressed genes were identified in the group CR, among which RAB11B-AS1, DISP1, MICAL3, PSMG1, and DOCK4 were up-regulated genes, and ANAPC1, IGLV1-40, SORT1, PLPPR2, and ATP1A1-AS1 were down-regulated. 7335 genes were screened in the group HR, including the top 5 up-regulated genes ALKBH6, AMBRA1, MIR1249, TRAV18, and LRRC69, and the top 5 down-regulated genes FAM241B, AC018529.3, AL031963.3, AC006946.1, and FAM149B1. The GO and KEGG analysis of the two groups revealed a significant enrichment in immune response and apoptosis. In the PPI network constructed, group CR and group HR identified 10 genes, respectively, and TP53BP1, SNRPD1, and SNRPD2 were selected as hub genes.
CONCLUSIONS: Using the messenger ribonucleic acid (mRNA) expression data from GSE166253, we found TP53BP1, SNRPD1, and SNRPD2 as hub genes in RTP patients, which is vital to the management and prognostic prediction of RTP patients.

INTRODUCTION

The coronavirus disease 2019 (COVID‐19), induced by the novel coronavirus, severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), has emerged as a global public health crisis. Typical symptoms of the disease are fever, sore throat, fatigue, cough, or dyspnea ; droplet transmission is the main route of transmission, with close contact providing an opportunity for such virus transfer. , As of October 1, 2021, there are 233,136,147 confirmed cases throughout the world, including 4,771,408 deaths reported by the World Health Organization (WHO).

So far, a large number of patients have been cured and discharged from hospitals. However, some, described as “retesting positive” (RTP) patients, were retested positive for SARS‐CoV‐2 after discharge. According to a systematic review, the RTP cases accounted for about 12% of discharged patients with the median duration of viral shedding in RTP patients after admission being 14.0. Most of these patients' clinical results revealed moderate severity conditions along with clinical symptoms including cough, sputum production, headache, and dizziness. As for the cause of RTP, the related researches focus on the false‐negative RT‐PCR testing at discharge or the persistence of the virus and periodic shedding in the gastrointestinal tract. However, there still remain concerns that discharged patients may be at risk of viral reactivation and should be considered potential sources of SARS‐CoV‐2 infection, which may become a public health crisis. Therefore, further exploration of the impaired immune function and pathogenesis is required.

Based on the previous study in immunology, a higher level of white blood count (WBC) was independently relevant to RTP by RT‐PCR. Besides, patients infected with COVID‐19 also exhibited higher leukocyte numbers and elevated plasma pro‐inflammatory cytokine levels. Additionally, compared with convalescent patients, they presented lower levels in lactate dehydrogenase (LDH) at the first week after admission but higher levels in eosinophil count. According to Yao MQ et al., patients presenting with decreased leukocyte, platelet, and CD4+ T counts were at a higher risk of RTP. Xiangying Ren et al. revealed that ribonucleic acid (RNA) of viral could still be detected in patients with SARS‐CoV‐2‐specific immunoglobulin M (IgM) and immunoglobulin G (IgG) present. Besides, a systematic review presented that it's predictive to recognize RTP through a combination of lower high‐sensitivity C‐reactive protein (hs‐CRP) and higher WBC. Overall, little is known about RTP in immunology and pathogenesis, and many results are in conflict with the previous study.

Therefore, we believe that it is necessary to analyze the biological information of the samples. The research results may further expand the insight into RTP events in COVID‐19 by analyzing the RTP immunological characteristics and molecular pathogenesis genes of peripheral blood mononuclear cells containing RTP patients, convalescent patients, and healthy control subjects.

MATERIALS AND METHODS

Data collection

Aiming at analyzing the biological information of RTP patients, we performed a search of the Gene Expression Omnibus (GEO) data using SARS‐CoV‐2 as the keyword (ncbi.nlm.nih.gov/geo/). Dataset GSE166253, contributed by Wang D et al., was found and used in our study (Table 1). The sequencing platform employed was the GLP20795 HiSeq X Ten (Homo sapiens). The dataset consisted of 10 retesting positive (RTP) samples, 6 convalescent, and 10 healthy ones, with the sample source containing the peripheral blood mononuclear cells (PBMCs) of human beings. To analyze the pathogenesis of RTP occurrence in COVID‐19, we divided the samples into two groups: a “convalescent‐RTP” group comprising of convalescent patients and RTP patients (group CR) and a “healthy‐RTP” group of healthy controls and RTP patients (group HR).

TABLE 1

Details of the data sources for this study

Gene expression profile	Sample collection	Sample genetic data included	Platform	Disease description
GSE166253	Peripheral blood mononuclear cells	GSM5066812 GSM5066813 GSM5066814 GSM5066815 GSM5066816 GSM5066817 GSM5066818 GSM5066819 GSM5066820 GSM5066821 GSM5066822 GSM5066823 GSM5066824 GSM5066825 GSM5066826 GSM5066827 GSM5066828 GSM5066829 GSM5066830 GSM5066831 GSM5066832 GSM5066833 GSM5066834 GSM5066835 GSM5066836 GSM5066837”	GPL20795 HiSeq X Ten (Homo sapiens)	Retesting positive(RTP): patients who have met the discharge criteria from COVID‐19 subsequently tested positive again for SARS‐CoV‐2

Screening of DEGs

After acquiring the gene expression data, the R pack (R Foundation for statistic computing) was utilized, and the fold change (FC) and p‐value were calculated to identify the differentially expressed genes (DEGs) from the two groups (CR and HR). Following this, in the R language program, p < 0.05 and |log2FC| ≥ 1 were designed as the threshold, and the genes that met the criteria were sifted out as DEGs. Finally, we constructed volcano plots and heatmaps to visualize all the DEGs of GSE166253.

Functional enrichment analysis

The GO annotation can be divided into three categories: biological process (BP), cellular component (CC), and molecular function (MF). In our study, we execute the GO function and KEGG pathway enrichment analysis for the DEGs, with the p‐value <0.05 deemed statistically significant.

Establishment of PPI and screening of Hub genes

Exploring functional interactions between proteins is crucial to understanding the molecular mechanisms of RTP for the COVID‐19 infection. The Search Tool for the Retrieval of Interacting Genes (STRING), an online tool to assess and integrate PPI information, contains physical and functional associations. To find the potential correlations of the DEGs, we employed STRING to map their PPI network. If the interaction score is >0.4, it will permit protein interaction. Then, the Cytoscape software was used to show the interaction of DEGs. To examine the protein interaction further, we employed this based on the PPI associations to screen out hub genes. Finally, the immune infiltration analysis was performed for the two groups to further confirm the results of our study. The results of the immunoblotting analysis are attached to the Appendix S1.

RESULTS

Downloading the GSE166253 dataset from the GEO database, we conducted differential expression analysis using the R language. In group CR, 6622 DEGs were identified, containing 5003 up‐regulated and 1619 down‐regulated genes (Figure 1A,B). RAB11B‐AS1, DISP1, MICAL3, PSMG1, and DOCK4 were the top 5 up‐regulated genes, and ANAPC1, IGLV1‐40, SORT1, PLPPR2, and ATP1A1‐AS1 were the top 5 down‐regulated genes. Simultaneously, 7335 DEGs were found in group HR, consisting of 4323 up‐regulated and 3012 down‐regulated genes (Figure 2A,B). The top 5 up‐regulated genes were ALKBH6, AMBRA1, MIR1249, TRAV18, and LRRC69, and the top 5 down‐regulated genes were FAM241B, AC018529.3, AL031963.3, AC006946.1, and FAM149B1. Subsequently, volcano plots and heatmaps were performed to present the DEGs in both groups (Figures 1 and 2). The volcano plots depicted up‐regulated (blue) and down‐regulated (red) genes, while the heatmaps presented the expression levels of the DEGs.

FIGURE 1

(A) Volcano plot of DEGs in group CR; (B) The heatmap of DEGs in group CR. The colors in (A) indicate different gene expressions in the corresponding group (red for up‐regulated DEGs and green for down‐regulated DEGs, gray indicates no difference); for (B), the abscissa axis represents the sample types and the ordinate axis represents the gene names

FIGURE 2

(A) Volcano plot of DEGs in group HR; (B) The heatmap of DEGs in group HR. The colors in (A) indicate different gene expressions in the corresponding group (red for up‐regulated DEGs and green for down‐regulated DEGs, gray indicates no difference); for (B), the abscissa axis represents the sample types and the ordinate axis represents the gene names

Functional enrichment

With regard to the CR group, the GO analysis in the biological processes demonstrated that the DEGs were significantly enriched in the processes that neutrophils are involved in, such as the activation of neutrophil, neutrophil‐mediated immune, and neutrophil activation participated in immune response. As for the cellular components, the mitochondrial inner membrane, and mitochondrial protein complex were predominantly recognized. Besides, the terms of the most enriched in molecular functions were associated with transcription coactivator activity as well as magnesium ion binding (Figure 3A). Moreover, according to the KEGG enrichment analysis results of group CR, three significant pathways of DEGs were identified, containing endocytosis, nucleotide‐binding oligomerization domain‐like receptors (NOD‐like receptor) signal path, and tumor necrosis factor (TNF) signal path (Figure 3B).

FIGURE 3

Results of functional analysis in group CR and group HR. (A) The GO analysis of group CR; (B) the KEGG analysis of group CR; (C) the GO analysis of group HR; (D) the KEGG analysis of group HR. In (A) and (C), the abscissa axis represents the gene ratio while the ordinate axis represents term names. The size of a single bubble represents the degree of enrichment; the color variety represents the different q‐values (those in red are considered to be of significance). For (B) and (D), the abscissa axis represents counts; the ordinate axis represents KEGG pathways; color represents the same meaning as the bubble diagram

Similarly, the results in group HR were shown as follows: the GO enrichment in biological processes revealed that DEGs were primarily located on neutrophil activation, neutrophil‐mediated immunity, and neutrophil activation participating in immune response. In terms of the category of cellular components, the mitochondrial inner membrane, mitochondrial matrix, and nuclear speck were chiefly enriched. For molecular functions, the mitochondrial inner membrane, mitochondrial matrix, and nuclear speck were significantly identified (Figure 3C). In the KEGG pathway enrichment analysis, a total of 10 key pathways were found, especially the pathways associated with prion disease, Huntington's disease, and Parkinson's disease, which take part in the corresponding neurodegenerative disorder (Figure 3D).

We constructed the PPI networks for both the CR and HR groups, after which the top 10 genes with the highest degree of connectivity were determined using the Cytoscape and cytoHubba plugin (Figures 4 and 5). In group CR, 1246 nodes and 5665 edges were included in the PPI network for up‐regulated DEGs (Figure 4A), while 1060 nodes and 613 edges were found in the PPI network of down‐regulated DEGs (Figure 4B). Afterward, TP53BP1, SNRPD1, SNRPD2, SF3B1, SNRNP200, MRPS16, MRPS9, CALM1, PPP2R1A, and YWHAZ were identified as hub genes (Figure 4C). As for group HR, 1223 nodes and 1257 edges were found in the PPI network of up‐regulated DEGs (Figure 5A), whereas for down‐regulated DEGS, 732 nodes and 228 edges were found in the PPI network (Figure 5B). With the same methods, 10 hub genes, including TP53BP1, RPS15, EFTUD2, MRPL16, MRPL17, MRPS14, RPL35A, MRPL32, MRPS6, and POLR2G, were screened out and their interactions were displayed in Figure 5C.

FIGURE 4

(A) PPI network of up‐regulated genes in group CR; (B) PPI network of down‐regulated genes in group CR; (C) Hub genes in group CR. Nodes represent genes; lines represent the interactions between gene‐encoded proteins; the redder the color, the more significant it is

FIGURE 5

(A) PPI network of up‐regulated genes in group HR; (B) PPI network of down‐regulated genes in group HR; (C) Hub genes in group HR. Nodes represent genes; lines represent the interactions between gene‐encoded proteins; the redder the color, the more significant it is

DISCUSSION

With the spread of the COVID‐19 epidemic, the situation of RTP for SARS‐CoV‐2 after discharge renders the management of discharged patients a difficult problem, potentially leading to a public health crisis. At present, there is still no definitive conclusion about the cause of RTP, and the possibility of their re‐infection and disease transmission cannot be ruled out. In this study, the peripheral blood samples of RTP patients, convalescent patients, and healthy people were divided into two groups—CR and HR groups, and the immune mechanism of COVID‐19 RTP was explored based on the results of the enrichment analysis.

Through the GO analysis for both the groups, the immune processes' neutrophils involved and the mitochondrial function and contents were significantly enriched. As Cynthia Magro et al. studied, patients with COVID‐19 had markedly increased neutrophils. The extracellular bactericidal networks (NETs) of SARS‐CoV‐2‐activated neutrophils are considered important in the progress of the COVID‐19 infection and releasing NETs in lung tissue could result in the death of lung epithelial cells, indicating that the immune processes that neutrophils are involved in RTP patients may be more activated in lysing infected cells and release the virus. The mitochondria play a variety of interrelated functions in cellular metabolism, producing adenosine triphosphate (ATP) and many biosynthesis intermediates, as well as participating in stress responses such as autophagy and apoptosis. Studies have pointed out that the monocytes of patients with COVID‐19 accumulate dysfunctional mitochondria, resulting in their metabolic defects, decreased OXPHOS and glycolysis, impaired oxidative burst, altered immune response, and possibly increased susceptibility to secondary infections. RTP patients may experience more intense mitochondrial dysfunction. Additionally, mitochondrial dysfunction was a key factor in neurodegenerative diseases, which also corresponded to the KEGG‐rich pathway in the HR group (prion disease, Huntington's disease)—this may be related to the prognosis of RP patients.

Apart from neurodegenerative diseases, the KEGG analysis of both the CR and HR groups also revealed that endocytosis, the NOD‐like receptors and TNF pathways were the main enrichment pathways. Endocytosis plays a vital role in the infection of disease and restoration of health. Recent studies have suggested that the key event supporting the SARS‐CoV‐2 infection may be the endocytosis of viral particles, which may be more intense in RTP patients. NOD‐like receptors (NLRs), a class of pattern recognition receptors (PRRs), can recognize ligands from microbial pathogens, host cells, and environmental sources, then activate an inflammatory response. Similar to mitochondria, NLRs play a significant part in regulating apoptosis, acting as a key regulator of apoptosis, and aiding in early development and directing autophagy by recruiting ATG16L1 to the plasma membrane at the site of the bacterial entry sites. The TNF superfamily is a multifunctional pro‐inflammatory cytokine that participates in the processes of cell survival, apoptosis induction, inflammatory response, and cell differentiation. In the SARS‐CoV‐2 infection, innate immune cells will produce a variety of inflammatory cytokines, but only the combination of TNF‐α and IFN‐γ can damage vital organs by inducing the death of inflammatory cells. In summation, the NOD‐like receptors and TNF pathways may contribute extensively to RTP patients through apoptosis and inflammatory response.

Since the results of the HR group cannot exclude the interference with the COVID‐19 infection, we selected the hub genes of the CR group for discussion next, with the HR group taken as verification. TP53BP1 was enriched in both groups, and its expression in the CR group was down‐regulated. TP53BP1, one of the regulator genes of TP53, is known to be highly associated with the DNA damage signal pathway, checkpoint signal, and DNA repair capacity. Moreover, it has been associated with cancer and neurological diseases. , In the immune system, DNA repair factors TP53 can support immunoglobulin class switch recombination. B cells are critical for the production of antibodies and protective immunity to viruses, and Tan et al. discovered a stronger antibody response is related to delayed viral clearance and disease severity. According to their research, the down‐regulation of TP53BP1 resulted in a lower antibody response, indicating a higher level of viral clearance in RTP patients. In a study on the mechanisms of the maintenance of foot‐and‐mouth disease (FMD) virus persistence, TP53BP1 in the persistently positive group was also significantly down‐regulated. It was reported that the down‐regulation of TP53BP1 could inhibit the Th17 response, natural killer cell (NK cell) cytotoxicity, and apoptosis. The main effector of Th17 cells is IL‐17, while the lymphocytes in COVID‐19 patients can produce more IL‐17, and blocking IL‐17 could be a novel therapeutic strategy for COVID‐19. The lower expression of IL‐17 in RTP patients may lead to milder symptoms. Preliminary studies of patients with severe COVID‐19 indicate that the decrease in the number or function of NK cells will reduce the clearance of infected and activated cells. Moreover, the limitation of apoptosis and reduction in clearing infected and activated cells may lead to the slower death of infected cells, resulting in the longer discharge of RNA of SARS‐CoV‐2 in RTP patients.

SNRPD1 and SNRPD2 are both up‐regulated in group CR, with the main regulatory pathway of the former being the mRNA splicing. The study confirmed that spliceosome machinery is a replicative machinery of SARS‐CoV‐2 for evading host challenges, which may prove that the replication activity of SARS‐CoV‐2 is more active in RTP patients. According to Xiaofeng Dai et al., SNRPD1 over‐expression is observably connected to genes involved in the cell cycle, cell mitosis and chromatin replication, and the expression of SNRPD1 in malignant or the number of hyperproliferative cells is higher than in normal cells. The high correlation with the cell cycle in SNPRD1 expression emphasizes that the cell proliferation of RTP patients may be more active. Similar to SNRPD1, SNRPD2 regulates mRNA splicing and is highly associated with Alzheimer's disease, which is similar to the KEGG results in group HR.

In our study, many enriched pathways and the effects of hub genes are related to apoptosis, which could be evidence of the hypothesis that the cause of reactivation leads to delayed virus clearance. The immunity in RTP patients, the immune processes' neutrophils involved, endocytosis, NOD‐like receptors and TNF pathways, as well as the immune cell function regulated by the hub genes, still require further research to explore their characters in RTP. Besides, the pathway of neurodegenerative diseases was significantly enriched, which makes the management and close attention of the RTP patients crucial.

However, due to the discrepancy in the theory of RTP and the small number of samples, our study has limitations. To make up these limitations, we expect further studies with larger samples to explore the association between each hub gene and the role of the pathways in RTP.

CONCLUSION

In general, using the biological analyses, we have explored the pathogenesis of RTP occurrence in COVID‐19 from the immune mechanism and molecular level. With the analysis of DEGs, we found TP53BP1, SNRPD1, and SNRPD2 as hub genes that are considered vital to apoptosis, changes in immune cell function, cancer and neurodegenerative disease. We regard apoptosis as a key pathway in the persistence of the virus and its periodic release. Besides, the prognosis of RTP patients should be of concern because the pathway is regulated by hub genes in cancer and neurodegenerative disease. With these findings, the mechanism in RTP and related management measures for RTP patients would induce further clarity. However, aiming at deepening the understanding of the immune mechanism in RTP patients, we expect future studies to further explore these at the genetic and immunity levels.

ETHICAL APPROVAL

The GSE166253 dataset used during the current study is available from GEO database (https://www.ncbi.nlm.nih.gov/gds/). Funding statement, conflict of interest, ethics approval statement, patient consent statement, permission to reproduce material from other sources, and clinical trial registration do not applicable to this article.

Appendix S1

Click here for additional data file.