Immunomodulatory mechanisms of abatacept: A therapeutic strategy for COVID-19

Coronavirus disease 2019 (COVID-19) caused by coronavirus-2 (SARS-CoV-2) infection has rapidly spread throughout the world and become a major threat to human beings. Cytokine storm is a major cause of death in severe patients. Abatacept can suppress cytokines used as antirheumatic drugs in clinical applications. This study analyzed the molecular mechanisms of abatacept treatment for COVID-19. Differentially expressed genes (DEGs) were identified by analyzing expression profiling of abatacept treatment for rheumatoid arthritis (RA) patients and SARS-CoV-2 infection patients. We found that 59 DEGs were upregulated in COVID-19 patients and downregulated following abatacept treatment. Gene set enrichment analysis (GSEA) and Gene Ontology (GO) analysis showed that immune and inflammatory responses were potential regulatory mechanisms. Moreover, we verified 8 targeting genes and identified 15 potential drug candidates for the treatment of COVID-19. Our study illustrated that abatacept could be a promising property for preventing severe COVID-19, and we predicted alternative potential drugs for the treatment of SARS-CoV-2 infection.


Introduction
Coronavirus disease 2019 (COVID- 19) is an infectious and global pandemic disease caused by coronavirus-2 (SARS-CoV-2) (1). The clinical symptoms were characteristically severe acute respiratory syndrome and hyperinflammatory immune response (2). T-cell immune responses are at the forefront of eliciting potent antiviral responses (3). Uncontrolled viral infection in advanced diseases resulting from insufficient T-cell responses may lead to systemic inflammation and severe lung damage (4,5). Clinical manifestation is represented by an immune defense-based protective .
/fmed. . phase first and is characterized by broad inflammation subsequently, which may lead to multiple organ failure (MOF) in severe patients (6). The drastic cytokine storm in severe patients is one of the causes leading to death (7). The widespread use of biological therapies in rheumatoid arthritis has shown a rapid resolution. There was evidence that the combination of anti-cytokine agents and systemic corticosteroid therapy had lower mortality rates (8). While there are no effective drugs for COVID-19 presently, the current treatment option for severe patients is symptomatic supportive therapy. Many rheumatoid drugs have shown potential for the treatment of COVID-19 based on their pharmacological properties. Abatacept, a T-cell selective co-stimulation modulator specifically binding to CD80 and CD86, was approved for use in rheumatoid arthritis (RA) with good efficacy in clinical application. Abatacept modulates Tcell activation and prevents the production of cytokines and downstream immune responses in RA (9)(10)(11). Therefore, neutralizing inflammatory factors in cytokine release syndrome (CRS) based on the pathophysiology of rheumatic diseases will be of great value in preventing disease progression in severe COVID-19 patients (12). Steroid hormone is widely used in COVID-19 treatment, but its side effects are also obvious, such as femoral head necrosis. Drug repurposing is a strategy for identifying new indications for approval. Evidence has been provided for abatacept as a candidate therapeutic approach to prevent severe COVID-19 (13). However, the mechanism of abatacept acting on COVID-19 remains unclear. Meanwhile, foreknowledge of the similar drug-target network providing alternative treatment strategies for COVID-19 patients is still challenging.
In this study, we conducted a bioinformatics analysis to explore the role of abatacept in the treatment of COVID-19 and provide more alternative anti-COVID-19 therapeutic drugs of similar pharmacological effects in the absence of sophisticated drugs for COVID-19. We also performed gene set enrichment analysis (GSEA), Gene Ontology (GO) analysis, and miRNA-mRNA network construction to reveal potential molecular mechanisms of abatacept preventing excessive inflammation and MOF. We also identified 8 abatacept targeting genes in the COVID-19 treatment. Finally, we identified the top 15 drug candidates combined with 8 targeting genes, offering a therapeutic strategy for severe COVID-19.

Data acquisition
The RNA-seq data (GSE151161, GSE152418, and GSE157103) were downloaded from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database. The GSE151161 dataset consists of 76 whole blood samples from 38 RA patients before and post abatacept treatment. The GSE152418 dataset included 34 peripheral blood mononuclear cell (PBMC) samples from 17 COVID-19 patients and 17 normal controls, which were used for analysis. The GSE157103 dataset involved 126 whole blood leukocytes from 100 COVID-19 patients and 26 normal controls, which were used for validation ( Table 1).

Enrichment analysis of DEGs
Gene Ontology analysis of DEGs (biological processes, cellular components, and molecular functions) was performed using the R software (version 3.6.3) and the R package clusterProfiler (version 3.14.3). R package org.Hs.eg.db (version 3.10.0) was used for ID conversion. GO for the immune system process was analyzed using the CluGO (version 2.5.7) (16) apps of Cytoscape Software (version 3.8) (17). Padj < 0.05 was considered as the threshold. The top 15 GO terms with the smallest Padj value were presented.

Di erential protein and mRNA expression of hub genes in multiple tissues
We further explore the RNA and protein expression of 8 hub genes in multiple tissues in the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) (20). The RNA expression data were in transcripts per kilobase of exon model per million mapped reads (TPM) format from the Consensus dataset.

Single-cell mRNA expression of hub genes in lung
In the HPA database, we explored the single-cell RNA expression of 8 hub genes in lung tissues. The cell types included macrophages, alveolar cells type 2, T cells, granulocytes, fibroblasts, club cells, respiratory ciliated cells, endothelial cells, and alveolar cells type 1.
Di erential mRNA expression of hub genes in multiple immune cells

Prediction of the miRNA-mRNA interaction
The NetworkAnalyst database (https://www. networkanalyst.ca/) was used to predict target miRNAs of hub genes (21). In addition, the R ggalluvial package (version 0.12.3) was used to construct the mRNA-miRNA co-expressed interaction networks based on the interaction information.

Relationship between hub gene and immunocyte-related gene expression
The Spearman's rank correlation coefficient was used to analyze the relationship between the expression of hub genes and immunocyte-related genes, including CD3, CD4, CD8, CD19, CD20, CD27, CD28, CD56, CD11b, CD66b, CCL4, and CCL5. R (version 3.6.3) was used for analysis and the R ggplot2 package (version 3.3.3) was used for visualization. Threshold values were considered as follows: * p < 0.05 indicates a mild correlation and * * p < 0.01 indicates a moderate correlation.

Identifying drug candidates
We used the Drug Signatures database DSigDB tool of Enrichr (https://maayanlab.cloud/Enrichr/) to identify drug candidates targeting hub genes for COVID-19 treatment (22). The top 15 drug candidates were selected according to combined scores, from highest to lowest.

Statistical analysis
All statistical analysis was carried out using R (version 3.6.2). The Spearman's rank correlation coefficient was used for correlation analysis. A Wilcoxon rank-sum test was used to compare the expression of hub genes in two groups. The hypothetical test was two-sided in all tests, and a p-value of <0.05 was considered statistically significant.

Identification of DEGs after SARS-CoVinfection and abatacept treatment
The volcano map showed 36 upregulated genes and 93 downregulated genes in GSE151161 after the abatacept treatment for 12 weeks in patients with RA ( Figure 1A). COVID-19 caused significant gene expression changes in the samples of

PBMC.
A total of 4,182 DEGs were identified from GSE151161, including 3,881 upregulated genes and 301 downregulated genes ( Figure 1B). The Venn plot showed that 59 DEGs were both upregulated in COVID-19 patients and downregulated post abatacept treatment ( Figure 1C). Conversely, there were no DEGs downregulated in COVID-19 patients and upregulated post abatacept treatment ( Figure 1D).

GSEA showed the immune-inflammatory response activity
To obtain insight into the effect of abatacept treatment and SARS-CoV-2 infection, GSEA analysis indicated several biological processes. The top 20 significant gene sets are shown in Tables 2, 3. We showed representative 5 gene sets mainly participated in immunity and inflammation pathways enriched in whole blood with abatacept treatment, such as interleukin (IL)-6 JAK STAT3 signaling, interferon-alpha response, and complement (Figures 2A-E), and representative 10 gene sets in PBMC with SARS-CoV-2 infection such as anti-inflammatory response favoring leishmania parasite infection, antigen-activated B cell receptor BCR leading to the generation of the second messenger, and complement cascade ( Figures 2F-O).

GO analysis of DEGs
We used GO analysis to identify 58 GO terms of 59 DEGs upregulated in COVID-19 patients and downregulated post abatacept treatment, which included biological process (BP), cellular component (CC), and molecular function (MF) (Figure 3). Many immune-mediated pathways were enriched, such as the immune response-activating cell surface receptor signaling pathway, B cell-mediated immunity, lymphocytemediated immunity, and humoral immune response. The circular diagrams displayed the corresponding relationship between DEGs and GO terms ( Figure 4). The predominantly related pathways of immune response mainly include the classical pathway of complement activation, immunoglobulinmediated immune response, Fc-gamma receptor signaling pathway, immune response-activating signal transduction, and cell surface receptor signaling pathway. The Cluego tool of Cytoscape software was used to enrich immune-related pathways of 59 DEGs, such as the classical pathway of   Table 4).

PPI network construction and hub gene identification
We obtained the top 14 hub genes with the highest interaction degrees and constructed a PPI network using the STRING database ( Figure 5A), which included DAAM2, KIF20A, PCSK9, MIXL1, CDC20, FOXC1, SDC1, CAV1, GPRC5D, IGF1, TSHR, RIMS2, NTRK3, and ADAMTS2. PPI network and function analyses were further analyzed for 14 hub genes. The results illustrate that the complex PPI network was with physical interactions of 70.90%, co-expression of 16.01%, co-localization of 3.22%, genetic interactions of 2.63%, predicted of 4.96%, shared protein domains of 0.55%, and a pathway of 1.74%. Receptor-mediated endocytosis, insulin-like growth factor binding, insulin-like growth factor receptor signaling .
/fmed. . pathway, growth factor binding, regulation of cellular protein catabolic process, regulation of nuclear division, and regulation of MAP kinase activity were identified as the main functions of those hub genes ( Figure 5B).

Validation of hub gene expression
We then verified the expression of 14 hub genes in the GSE157103 dataset. The expression of eight genes, including CAV1, CDC20, GPRC5D, IGF1, KIF20A, MIXL1, SDC1, and TSHR, was significantly upregulated in whole blood leukocyte samples of COVID-19 patients than normal controls. All of them were consistent with our previous analysis in the GSE152418 dataset ( Figures 6A-H). The characteristics of the eight genes are described in Table 5. The expression of ADAMTS2, DAAM2, FOXC1, NTRK3, PCSK9, and RIMS2 genes has no significant difference between COVID-19 patients and normal controls (Figures 6I-N).

Construction of mRNA-miRNA network and immune correlation analysis
To further explore the potential mechanisms and regulating axis of hub genes in the process of abatacept-treated COVID-19, we established mRNA-miRNA co-expressed interaction networks using the NetworkAnalyst database and analyzed the correlation between hub genes and immunocyte-related molecular markers. In mRNA-miRNA co-expressed networks, 230 target miRNAs of seven hub genes without GPRC5D  were included ( Figure 7A). The expression of 8 hub genes was related to inflammatory factors and lymphocyte-related gene markers, including CD3, CD4, CD8, CD19, CD20, CD27, CD28, CD56, CD11b, CD66b, IL-6, CCL4, and CCL5. The correlations indicated the importance of immune cells participating in the immune-inflammation response ( Figure 7B). The noticeable correlation between 8 hub genes and immune cell-related genes, including IL-6, CD3, CD8, CD19, CD20, CD27, CD28, CD56, CD66b, and CCL5, is also displayed in the annular chart ( Figure 7C).

Identification of drug candidates
To identify potential drug candidates for the treatment of COVID-19, we predicted the top 15 drug candidates targeting 8 hub genes such as deferoxamine, monobenzone, bicalutamide, trifluridine, and raloxifene (Supplementary Table 1).

Tissue-specific and cell-specific expression of eight hub genes
The expression of 8 hub genes in multiple tissues was analyzed at RNA and protein levels. Remarkably, CDC20 RNA expression was highest in bone marrow, KIF20A was highest in the thymus, and MIXL1 was highest in the tonsil (Supplementary Figure 1). At the protein level, CAV1 was highly expressed in the lung, KIF20A was highly expressed in lymph nodes and bone marrow, GPRC5D was highly expressed in the spleen, and SDC1 was highly expressed in the tonsil (Supplementary Figure 2). CAV1 was highly expressed in lung tissue both at RNA and protein levels. The single-cell expression of 8 hub genes in lung tissues showed that CAV1, KIF20A, and SDC1 genes were highly expressed in alveolar cells. The CDC20, IGF1, and TSHR genes were highly expressed in endothelial cells (Supplementary Figure 2). The expression of 8 hub genes in immune cell subtypes was mainly in basophil, regulatory T cells, B cells, and NK cells (Supplementary Figure 3).

Discussion
The SARS-CoV-2 is responsible for the COVID-19 pandemic, which has rapidly spread throughout the world and become a major threat to human beings (23). COVID-19 continues to expand in the pandemic form and is responsible for significant morbidity and mortality. The initial symptoms of COVID-19 are mainly fever, cough, myalgia, fatigue, and dyspnea. In addition, acute respiratory distress syndrome or multiple organ failure may gradually occur in terminal COVID-19 patients (24). COVID-19 patients requiring hospitalization were generally aged with multiple comorbidities. There is a relatively high rate of severe mortality in the early days of the epidemic, at 49 and 33%, respectively (25). Approximately, 14% of patients were severe cases that required ventilation in an intensive care unit (ICU), 5% were critical, and around 2.3% died in a report of 72,314 cases in China (26). While the treatment for severe patients is mainly symptomatic supportive therapy. Drug repositions are necessary to be accelerated for severe patients.
Angiotensin-converting enzyme-2 (ACE2) is a receptor for the binding and entry of the SARS-CoV-2 virus into cells (27). On binding to epithelial cells in the respiratory tract, .
/fmed. . SARS-CoV-2 starts replicating and migrating down to the airways and finally enters alveolar epithelial cells in the lung tissues. The rapid replication of SARS-CoV-2 in the lungs may trigger strong immune and inflammatory responses and produce large amounts of inflammatory cytokines, which leads to the cytokine storm, a major cause of patient death (28,29). These uncontrolled inflammatory responses may lead to local and systemic tissue damage (30). Elevated pro-inflammatory cytokine is a typical profile in patients with COVID-19, such as IL-6, IL-1, and tumor necrosis factor (TNF)-α (31). Tocilizumab, as an anti-IL-6 receptor antibody, has performed its superiority in preventing severe outcomes such as mechanical ventilation and reducing death risk (32). Anakinra, Canakinumab, and Rilonacept are IL-1 blockades, which are currently used not only for therapy of RA and many other autoimmune rheumatic diseases but also for the treatment of the cytokine storm (33). TNF-α is an important mediator of other cytokines and chemokine production, so anti-TNFα therapy could be useful in . Similarly, these antirheumatic drugs might play a protective role in the development of the exaggerated  immune-mediated inflammatory response associated with SARS-CoV-2 infection. The T-cell immune response is essential for protecting from COVID-19 and participates in abating innate immune responses involved in cytokine syndrome (31). In terms of abatacept, it can compete with CD28 for CD80/CD86 receptors, inhibiting its downstream inflammation reaction and suppressing the expression of other costimulatory molecules on antigen-presenting cells, resulting in a decrease in immune responses (35). Hence, viral inhibition is expected to be the most effective in the early disease course, while immunosuppressive treatment may be useful in the later stages to prevent severe disease. Multiple biological agents have shown their huge potential in the COVID-19 treatment (12).
We used bioinformatics to deeply analyze the RNA-Seq datasets GSE151161 and GSE152418, which included the samples of abatacept-treated RA and SARS-CoV-2-infected patients. First, we screened 129 DEGs in the GSE151161 dataset and 4,182 DEGs in the GSE152418 dataset. Venn diagram results identified 59 DEGs, upregulating in COVID-19 patients and downregulating in RA patients post abatacept treatment for 12 weeks. We then performed a GSEA analysis of all detected genes.  Previous studies have shown that COVID-19 induced a poor immune response, leading to virus-induced pathology, or a hyperactive immune response that leads to cytokine storms associated with uncontrolled inflammation, severe pulmonary tissue damage, and even death in severe COVID-19 patients (36). Similar to previous studies, our GSEA analysis showed immune-related signal pathways, which were consistent with the pathogenesis of COVID-19 and RA.
The top genes with the highest degree of interaction in the PPI network were considered hub genes, which may be critical for the treatment of COVID-19 patients. In our study, functional analysis of 14 hub genes showed their associations with receptor-mediated endocytosis, insulin-like growth factor binding, insulin-like growth factor receptor signaling pathway, growth factor binding, regulation of cellular protein catabolic processes, regulation of nuclear division, and regulation of MAP kinase activity, which were related to virus entry, cell growth, regulation of protein catabolism, and mitosis. We confirmed the expression of 8 hub genes in the GSE152418 dataset, including CAV1, CDC20, GPRC5D, IGF1, KIF20A, MIXL1, SDC1, and TSHR. A previous study found that a high level of SDC1 and a low level of IGF1 increased mortality, which are potential biomarkers of severe COVID-19 (37,38). In addition, CAV1, CDC20, and KIF20A were also identified as hub genes in COVID-19 by other researchers (39)(40)(41). Our study confirmed and complemented previous studies to some extent.
We identified 15 drug molecules targeted at 8 hub genes. Consistent with our results, deferoxamine, bicalutamide, trifluridine, raloxifene, etoposide, methotrexate, and progesterone were considered potential drugs for COVID-19 in previous studies among the top 15 drug candidates (42-48). However, the therapeutic effect of these candidate drugs on COVID-19 needs further study.
It is generally accepted that increased pro-inflammatory cytokines are associated with disease severity. Blocking immunity-induced inflammation was essential for reversing immunopathology.
Our mRNA-miRNA co-expressed interaction networks revealed the potential regulatory mechanism of targeting genes, which guides the follow-up study of COVID-19 treatment mechanisms. A cytokine storm is potentially fatal and is characterized by the high-level activation of immune cells and the excessive production of massive inflammatory cytokines. A comparison between ICU and non-ICU patients showed that plasma concentrations of IL2, IL7, IL10, and TNFα were higher in ICU patients than in non-ICU patients (49). A retrospective study suggested elevated IL-6 was possibly fatality contributors and that mortality might be due to viral-driven hyperinflammation (50). Cytotoxic T cells, B cells, NK cells, and neutrophils can trigger SARS-CoV-2 infection-mediated CRS (51). Increased neutrophil number and neutrophil/lymphocyte ratio are usually accompanied by advanced severity and poor clinical outcomes (52). In some COVID-19 patients, the feature of severe SARS-CoV-2 infection is lymphopenia with severely exhausted CD4+ T cells, CD8+ T cells, B cells, and NK cell counts (51). Some studies also analyzed histological changes induced by SARS-CoV-2 infection. Tissue samples obtained from COVID-19 patients showed marked infiltration of various T-cell subclasses in the lungs (53). Thus, a dysregulated and excessive immune response not only initiates immune infiltration but also results in extensive inflammation and immunopathology through the induction of proinflammatory cytokines (51). Our study found a high correlation between 8 target genes and immune cell markers (CD8, CD19, CD20, CD27, CD56, and CD66b), which further demonstrated the crucial role of 8 target genes in COVID-19 treatment.
Overall, we provided insights into the SARS-CoV-2 infection and the immunomodulatory mechanisms of abatacept treatment. Our study suggested that abatacept was a potential strategy for severe COVID-19 based on the intersection of DEGs, the resemblance of antirheumatoid mechanisms and immune-inflammatory responses in COVID-19. Abatacept treatment could be a critical way to avoid an inflammatory storm in COVID-19 (13). Compared with the study of Julia et al., we conducted an analysis using larger sample size, identified potential biomarkers, and explored deeper mechanisms through which abatacept treats CODVID-19. More clinical trials of drug candidates should be tried to exert efficacy and tolerable safety in clinical therapy. The current study has several limitations. Only bioinformatics analysis was conducted in this study. We have further planned randomized and controlled trials (RCT) to support our conclusions.

Data availability statement
Publicly available datasets were analyzed in this study. This data can be found at: GEO database, under accessions GSE151161, GSE152418, GSE157103.

Author contributions
DY, HL, and QJ: conception and design. QJ: administrative support. DY, YC, and QJ: provision of study materials. HL, YC, and WR: collection and assembly of data. DY, WR, MD, and CL: data analysis and interpretation. DY and HL: manuscript writing. All authors approved the final manuscript.

Funding
This study was supported by Natural Science Foundation of Shanxi Province: 20210302123287.