key: cord-0760815-hlubxphj
authors: Julià, Antonio; Bonafonte-Pardàs, Irene; Gómez, Antonio; López-Lasanta, María; López-Corbeto, Mireia; Martínez-Mateu, Sergio H.; Lladós, Jordi; Rodríguez-Nunez, Iván; Myers, Richard M.; Marsal, Sara
title: Targeting of the CD80/86 proinflammatory axis as a therapeutic strategy to prevent severe COVID-19
date: 2021-06-01
journal: Sci Rep
DOI: 10.1038/s41598-021-90797-0
sha: b71f48a22e0d8a729781de21aca451973a7e8bfd
doc_id: 760815
cord_uid: hlubxphj

An excessive immune response known as cytokine storm is the hallmark of severe COVID-19. The cause of this cytokine rampage is yet not known. Based on recent epidemiological evidence, we hypothesized that CD80/86 signaling is essential for this hyperinflammation, and that blocking this proinflammatory axis could be an effective therapeutic approach to protect against severe COVID-19. Here we provide exploratory evidence that abatacept, a drug that blocks CD80/86 co-stimulation, produces changes at the systemic level that are highly antagonistic of the proinflammatory processes elicited by COVID-19. Using RNA-seq from blood samples from a longitudinal cohort of n = 38 rheumatic patients treated with abatacept, we determined the immunological processes that are significantly regulated by this treatment. We then analyzed available blood RNA-seq from two COVID19 patient cohorts, a very early cohort from the epicenter of the pandemic in China (n = 3 COVID-19 cases and n = 3 controls), and a recent and larger cohort from the USA (n = 49 severe and n = 51 mild COVD-19 patients). We found a highly significant antagonism between SARS-CoV-2 infection and COVID-19 severity with the systemic response to abatacept. Analysis of previous single-cell RNA-seq data from bronchoalveolar lavage fluid from mild and severe COVID-19 patients and controls, reinforce the implication of the CD80/86 proinflammatory axis. Our functional results further support abatacept as a candidate therapeutic approach to prevent severe COVID-19.

| (2021) 11:11462 | https://doi.org/10.1038/s41598-021-90797-0 www.nature.com/scientificreports/ than CD28 on the surface of T cells, thereby preventing the necessary costimulatory signal required by T cells to progress to activation. It was originally developed to reduce inflammation in rheumatoid arthritis (RA), where T cell activation is central to the disease pathology, and has been recently approved also for the treatment of psoriatic arthritis and juvenile idiopathic arthritis 9 .

In severe COVID-19, macrophages in the lungs activate and produce large amounts of proinflammatory cytokines, mainly IL6. Recent evidence suggests that during SARS-CoV-2 infection, resident alveolar macrophages are gradually substituted by monocyte-derived proinflammatory macrophages 10, 11 . However, how these macrophages are brought to drastically increase the proinflammatory cytokine production in severe COVID-19 is yet not known 5 . An increase in proinflammatory signals by hyperactivated T cells is a possible explanation 12 . Co-stimulation through the CD80/86 axis is a key mechanism for T-cell activation. We hypothesize that blocking this co-stimulation with abatacept could be an effective way to dampen the excessive T cell hyperactivation and subsequent macrophage cytokine rampage that leads to severe COVID-19. This mechanism could explain the lower prevalence of COVID-19 symptomatology observed in abatacept-treated patients. Here we provide additional evidence supporting the utility of blocking the CD80/86 axis to prevent COVID-19 associated hyperinflammation. Using RNA-seq, we have characterized the longitudinal response to abatacept at the systemic level. Analyzing blood RNA-seq data from two different COVID-19 cohorts, we demonstrate that this therapeutic response is highly antagonistic to that elicited by COVID-19. Our results provide additional evidence at the functional level that blocking the CD80/86 axis could be a useful therapeutic approach to prevent severe COVID-19.

In our first analysis, we evaluated how the biological features associated with COVID-19 pathology at the systemic level are affected by treatment with abatacept. A total of 22 processes representing previously described anti-viral and hyperinflammatory responses during COVID-19 were evaluated (Supplementary Table S1 ). Using a longitudinal cohort of n = 38 rheumatoid arthritis (RA) patients treated with abatacept, we tested whether these COVID-19-related processes are modified during CD80/86 axis blocking. Using this approach, we found that 16 COVID-19 associated processes were significantly modified by treatment with abatacept at the systemic level ( Fig. 1A ,B, Supplementary  Fig. S1A ,B, Supplementary Table S2 ). These pathways include processes crucial for COVID-19 severity including T cell mediated immunity, macrophage activation and cytokine production (IL1, IL6, IL8 and TNF, among others). Importantly, all the functional changes occurred in the opposite direction to that previously described in COVID-19 pathology.

To provide additional evidence supporting the antagonism of abatacept to COVID-19, we directly analyzed the perturbation of the systemic gene expression caused by . For this objective, we determined the biological pathways differentially activated by COVID-19 using two different COVID-19 datasets. First, we used RNA-seq data from a very early cohort of PBMC data from 3 COVID-19 patients and 3 healthy controls obtained at the epicenter of the pandemic in Wuhan, China 13 . After adjusting for granulocyte percentage, we found that, from a total of 6,187 biological pathways, 260 were differentially activated due to COVID-19 (FDR < 0.05, Supplementary Table S3 ). The same analysis approach in our longitudinal cohort identified 109 biological pathways significantly associated to the treatment with abatacept (FDR < 0.05, Supplementary Table S4 ). We next determined the level of overlap between the two exposures and found a total of 49 pathways significantly modified by both COVID-19 and abatacept ( Fig. 2A ). From these, 47 (96%) were modified by abatacept in an opposite direction to that elicited by COVID-19. The probability that this level of antagonism occurred by chance, is very low (P value < 2.5e−18, binomial test). Figure 2B shows the biological pathways that are more antagonistically modified by abatacept in relation to COVID-19. Among the pathways most significantly antagonized, B cell activation and Fc-mediated immune response were found to be highly activated during COVID-19 and significantly repressed by treatment with abatacept. The complete list of overlapping pathways is available in Supplementary Table S5. In October 2020, a larger dataset on COVID-19 patients became available (late cohort), consisting on a transcriptomic, metabolomic, lipidomic and proteomic data from moderate and severe COVID-19-positive patients as well as COVID-negative patients admitted to the hospital 14 . The whole blood RNA-seq data from n = 51 moderate COVID-19 patients and n = 49 severe patients from this study allowed us not only to corroborate our findings in the early cohort, but also to specifically test for association with clinical severity. From 6187 tested biological pathways, a total of 99 pathways were significantly enriched when comparing severe to moderate COVID-19 patients (FDR < 0.05, Supplementary Table S6 ). Using the same approach, 259 biological pathways were associated with the response to abatacept in our cohort (FDR < 0.05, Supplementary Table S7 ). Comparing the two exposures, we found that 67 pathways were significantly associated with both COVID-19 severity and the response to therapy (Fig. 3A) . Importantly, all of these common pathways were modified by abatacept in the opposite direction to that of severe COVID-19. This level of antagonism was highly unlikely to occur by chance (P value < 3.07e−35). Figure 3B shows the biological pathways that are more antagonistically modified by abatacept in relation to severe COVID-19. Key biological processes associated with severity, like platelet activation and interleukin-6 production, were among the most significantly antagonized biological pathways. The complete list of overlapping pathways is available in Supplementary Table S8 . At the single-gene level, from the 44 genes differentially expressed in response to Abatacept therapy (FDRadjusted P-value < 0.05), 15 were found to be expressed also in the COVID-19 late cohort. From these, five genes were differentially expressed in severe COVID-19 patients compared to mild patients (P-value < 0.05, Supplementary Table S9 ). Importantly, all five genes-HJURP, GTSE1, TYMS, CDC20 and BIRC5-are down-regulated by abatacept and up-regulated in severe COVID-19 patients. www.nature.com/scientificreports/ Survival analysis: antagonism to abatacept signature as a predictor of severity. To further corroborate the antagonism of the abatacept response signature to COVID-19, we explored its utility as a biomarker to predict severity. For this objective, we built a score (abatacept score, ABS) using the genes differentially expressed by treatment with abatacept (n = 15, Supplementary Information). This set was enriched in immunoglobulin encoding genes, suggesting a predominant gene expression change in the B cell compartment. Using this gene expression signature, we identified the COVID-19 patients in the late cohort with more similar (ABS-High) and dissimilar (ABS-Low) expression. A total of 11 COVID-19 patients were classified as ABS-High and 20 were classified as ABS-Low. Using this patient categorization, survival analyses were performed using two critical and complementary measures of severity: number of days without mechanical ventilation and number of days in the hospital. In both cases, we found that ABS-High patients (i.e. patients with a gene expression profile more similar to the abatacept-response signature) were significantly associated with a better prognosis (Fig. 4A ,B, Supplementary Fig. S2A ,B). ABS-Low patients, instead, had a higher probability to require mechanical ventilation (P = 0.047, log-rank test) and to spend more days in the hospital (P = 0.032, log-rank test).

After identifying a significant antagonism between the response to abatacept and COVID-19 pathology at the systemic level, we next sought to corroborate the relevance of this immune-signaling axis in the target tissue of COVID-19 severity. While COVID-19 can affect multiple organs, the lungs are the predominant target of the severity of the immune response and where the cytokine storm originates 3 . To evaluate the relevance of CD80/86 co-stimulation in the lungs during COVID-19 we took advantage of an available single-cell RNA-seq dataset generated from the bronchoalveolar lavage fluid obtained from mild (n = 3) and severe (n = 6) COVID-19 patients and from healthy controls (n = 3) 10, 11 . Figure 5A ,B show the cell clustering analysis and identification of the marker genes within each cell type, respectively. Based on this clustering, Fig. 6A -C display the expression of gene markers of the key elements of the C80/86 activation axis, respectively. Circle size is also proportional to the NES. P-value significance levels are indicated as *< 0.05, **< 0.005, ***< 0.0005, ****< 0.00005 and *****< 0.000005. The two middle columns show the change and statistical significance of these processes in the early COVID-19 cohort (comparing COVID-19 patients to healthy controls), and the late COVID-19 cohort (comparing severe to mild cases). www.nature.com/scientificreports/ According to our hypothesis, an excess of CD80/86 co-stimulation of antigen-presenting cells (APCs) to CD4+ T cells in SARS-CoV-2 infected individuals leads to a hyperactivation of macrophages that subsequently ramp up their production proinflammatory cytokines, mainly IL6. Treatment with abatacept would dampen this signaling axis, thereby reducing the likelihood of developing a cytokine storm. In this disease severity model, we would expect to find an association between the abundance of the three key elements of the CD80/86 axis: CD80/86 + APCs, activated CD4+ T cells and IL-6 producing macrophages. Using the single-cell BALF dataset, we found a significant and positive association between several of the pairwise combinations of the three elements of the CD80/86 axis ( 

Since the beginning of the pandemic, SARS-CoV-2 has infected more than 67 million people and caused more than 1,544,000 deaths world-wide as of December 8 2020 (https:// coron avirus. jhu. edu/ map. html). There is a pressing need to identify therapies that can prevent COVID-19 patients from reaching the hyperinflammatory life-threatening stage known as cytokine storm. Based on previous epidemiological evidence from rheumatic patients treated with several targeted immunomodulators, we hypothesized that blocking of the CD80/86 proinflammatory axis could be a useful therapeutic approach to prevent severe COVID-19. Here we provide additional functional evidence to support this hypothesis. To do so, we have characterized the systemic response to abatacept in a longitudinal cohort of rheumatic patients and subsequently tested for its antagonism to COVID-19 pathological features. Using two patient cohorts, we have found a high level of antagonism of abatacept-mediated response to COVID-19 pathology. The results of this study support abatacept as a candidate therapy to prevent severe COVID-19.

Our results suggest that T cell stimulation by CD80/86-expressing macrophages is a driver of COVID-19 severity. In non-pathological conditions, the lungs are the tissue expressing the highest levels of CD80 mRNA and the second tissue expressing more CD86 mRNA ( Supplementary Fig. S3A,B) . This high basal expression of the key elements of the CD80/86 axis might reflect the need for a rapid activation of the immune system in an (96%) were found to be modified antagonistically by abatacept while only 2 (4%) processes showed changes in the same direction. (B) 20 top biological pathways that show a higher level of antagonism between COVID-19 exposure and treatment with abatacept. FDR-adjusted P-values from the two analyses were combined using Fisher's method, and the resulting combined P-value are used here as a ranking measure of antagonism. Blue indicates upregulation and red downregulation of the pathway as measured by the GSEA Normalized Enrichment Score, which is also reflected in the circle size. FDR-corrected P-value levels are indicated as *< 0.05, **< 0.005, ***< 0.0005, ****< 0.00005 and *****< 0.000005. www.nature.com/scientificreports/ organ that is systematically exposed to multiple environmental cues 15 . However, this enhanced potential for T cell activation could become a highly detrimental factor in the presence of a rapidly propagating immune insult like SARS-CoV-2. A large population of CD80/86 APCs activated by a fast-spreading infection will lead to a higher activation of T cells which, in turn, will increase the activation proinflammatory cell programs, including IL6 production by lung macrophages 16 (Fig. 8 ). This high signal amplification potential could therefore explain the sudden transition into the extremely high cytokine production stage observed during the cytokine storm of most COVID-19 fatalities 17 .

We have found that the treatment with abatacept significantly down-regulates myeloid leukocyte differentiation and macrophage activation at the systemic level. The hyperinflammatory stage in COVID-19 is characterized by the activation of myeloid cells, mostly monocyte-derived macrophages that massively infiltrate the lung 5 . Therefore, the specific anti-inflammatory effect of this therapy could be a useful property to avoid entering into the life-threatening massive cytokine release stage. Our results also show that treatment with abatacept downregulates the pathways associated with both the production of and the response to several key cytokines associated to the severe stage in COVID-19, including TNF 18 , IL1 19 , IL8 20 , IL7 18 and IL6 18 . In particular, the downregulation by abatacept of systemic IL6 -the proinflammatory cytokine that has been shown to be predominant in patients progressing to severity 17 -was found to be highly significant. Our study suggests that IL6 production can also be effectively reduced by the upstream regulation of the immune system mediated by abatacept.

Our longitudinal analysis of the systemic response to abatacept showed that blocking of the CD80/86 axis leads to a strong downregulation of complement activation and B cell mediated immunity. Excessive complement activation has been associated to damage in the lungs of COVID-19 patients, and antibody therapies against elements of the complement system are currently being evaluated to treat severe cases 21 . An over-activation of the B cell mediated response could also be a harmful mechanism in severe COVID-19. During the SARS-CoV1 epidemic, deceased patients were found to have higher levels of neutralizing antibodies 22 , which raised the possibility of disease exacerbation through antibody-dependent enhancement (ADE) 23, 24 . Experimental evidence with rhesus macaques infected with SARS-CoV-1 showed that IgG antibodies against the spike protein contributed to the massive influx of monocytes to the lung and subsequent severe injury 25 . The observed downregulation of the B cell immunity by abatacept might reduce the likelihood of low-affinity antibodies and the probability of ADE in patients 26 . To our knowledge, however, solid evidence that ADE contributes to COVID-19 severity has not been yet provided 27 . www.nature.com/scientificreports/ A limitation of our study is that the longitudinal response to abatacept was measured in RA patients and not directly in COVID-19 patients. Therefore, our analysis of the systemic changes due to abatacept could have been confounded by the specific immune activation of RA. Like severe COVID-19 28 , however, RA is a disease where increased production of proinflammatory cytokines IL-6, IL-1β, and TNF is central to the pathological hyperinflammatory process 29 . Our longitudinal study design was key to capturing the effect and direction of abatacept-mediated regulation of these key pathways associated to severity. Nonetheless, a direct analysis of the transcriptional impact abatacept on SARS-CoV-2 infected patients or animal models was not evaluated. This is a limitation, because other aspects more specific to viral immune response will not have been represented. A direct analysis of the clinical and transcriptional impact abatacept on SARS-CoV-2 infected patients would be the ideal design to confirm our findings. In October 20th 2020, the US National Institutes of Health launched the ACTIV-1 IM protocol (#NCT04593940), a randomized clinical trial to evaluate immune modulators for the treatment of moderately or severely ill hospitalized COVID-19 patients. From a pool of more than 130 immune modulators, three agents, including abatacept, were selected to be evaluated in a randomized clinical trial with > 2000 COVID-19 patients admitted to the hospital. Our results show that there is potential for stratification of patients according to the gene signature at baseline. If the clinical trial shows a protective role for abatacept in severe COVID-19, our abatacept score could be tested for validation and be the basis of a biomarker profile to prioritize patients for preventative treatment.

There is increasing evidence that abatacept is an efficacious therapy to treat interstitial lung disease (ILD) in RA 30 . ILD is a comorbidity of RA that is characterized by the inflammation of the lung tissue and subsequent fibrotic development, which are also two canonical features of severe COVID-19. Treatment of this comorbidity with abatacept was motivated by the finding that blocking T cell co-stimulation is highly effective in a mouse model of hypersensitivity pneumonitis 31 . Much like in severe COVID-19, hypersensitivity pneumonitis is characterized by a massive influx of activated T cells in the lungs and a high expression of B7 (i.e. mouse CD80/ CD86 molecules) by alveolar macrophages. Evidence from other diseases also support the role of the CD80/86 activation mechanism in excessive inflammatory responses in the lung 32, 33 . Although the triggering immunogenic insults will be different in each disease, the importance of the CD80/86 axis in the immune response of the lungs might lead to a convergent hyperinflammatory response. www.nature.com/scientificreports/ In support of the safety and utility of CD80/86 axis blocking to prevent severe COVID-19, there is substantial epidemiological evidence showing that abatacept does not increase the risk of serious infections [34] [35] [36] . While other immunosuppressive agents have been shown to increase this susceptibility to infections in rheumatic patients, abatacept seems to render the immune system fully capable of fighting off viral infections. In accordance with this previous evidence, in our abatacept-treated cohort of patients, we have found that basic antiviral response mechanisms like the type I interferon response were not significantly altered by the therapeutic blocking of the CD80/86 axis. Together, these observations support the safety of using this therapy to prevent the development of severe COVID-19.

An extraordinary international effort and unprecedented collaborations have recently led to the identification of vaccines for COVID-19 that appear to be safe and efficacious. Despite this formidable advance, there is a clear consensus that it is still critical to identify therapies that can reduce the high mortality rates of COVID-19 patients. The results of our study provide additional functional evidence that blocking of the CD80/86 axis could www.nature.com/scientificreports/ board and informed consent was obtained in all cases. All research was performed in accordance with the Declaration of Helsinki, and informed consent was obtained from all participants. Supplementary Table S10 shows the clinical features of the RA patient cohort. Whole blood samples were obtained from all patients at the beginning of the therapy with abatacept and at week 12. Blood was collected using RNA-stabilizing PaxGene tubes (PreAnalytiX, Switzerland), to preserve total RNA integrity from the time of venipuncture.

Total RNA was extracted using the PaxGene blood isolation kit (Qiagen). All samples had a RIN index > 7 and were included for RNA-seq analysis. RNA-seq libraries were performed using with the KAPA RNA HyperPrep Kit, with RiboErase (HMR) Globin (globin and rRNA depletion protocol) for Illumina sequencing platforms. In order to improve the accuracy of the gene expression quantification, we added Unique Molecule Identifiers (UMIs). This was performed during the library construction, in which xGen UDI-UMI Adapters (IDT) were ligated after the second strand synthesis. These full-length adapters include P5 and P7 sequences that are compatible with Illumina sequencing platforms. Sequencing was performed using the NovaSeq platform (Illumina) at average of 250 paired-end reads. A total of 3,376,481,954 sequencing reads were obtained in the two conditions (weeks 0 and 12). FASTQ files were subsequently aligned to the GRCh37 human reference genome assembly using STAR 2.7.6a 38 . The global alignment ratio (mean of all samples) was 95.2% ( Supplementary Fig. S4 shows the read mapping distribution in all the samples and Supplementary Table S11 the per-sample QC). After alignment, we performed deduplication of the PCR duplicates using UMItools 39 to mark duplicates and samtools 40 to remove them. Gene-level read counts of the deduplicated bam files were obtained with featureCounts 41 for gene annotations downloaded from Ensembl Genes GRCh37v75. This resulted in a total of 3,212,779,029 (95.15%) of overall reads mapping as pairs to annotated genes. The RA longitudinal dataset is available through the NCBI Gene Expression Omnibus database with id GSE151161.

To characterize the gene expression profile associated with COVID-19, we used RNA-seq data from two COVID-19 cohorts. The first cohort 13 (early cohort) consists on the transcriptional profile on peripheral blood mononuclear cells (PBMCs) from three COVID-19-positive patients and three healthy controls from Wuhan, China and made available in March 2020. The raw sequencing data was processed as described by the authors. Abatacept is a recombinant CTLA4Ig fusion protein that selectively blocks CD80 or CD86 co-stimulation of T cells by antigen presenting cells (APCs). Reducing the level of CD4+ T cell activation with abatacept could be a useful therapeutic strategy to prevent reaching the hyperinflammatory state that defines severe COVID-19. Once an individual becomes infected with SARS-CoV-2, an immune response is mounted (first stage of COVID-19). In some individuals, inadequate clearance of the coronavirus by the anti-viral immune response leads to a progressive hyperinflammation (second stage of . In this stage, CD80/86-expressing APCs in the lungs-alveolar and monocyte-derived macrophages-activate CD4+ T cells which, in turn, stimulate macrophages to produce large quantities of proinflammatory cytokines, predominantly IL6. The large pool of CD80/86 + APCs in the lungs together with the rapidly spreading capacity of SARS-CoV-2 produce a dangerous combination that can lead to an abrupt and life-threatening overexpression of proinflammatory cytokines known as the cytokine storm. Blocking of the CD80/86 + axis with abatacept would reduce the probability of reaching this hyperinflammatory stage and therefore prevent severe COVID-19. www.nature.com/scientificreports/ The second cohort 14 (late cohort) consists on the whole blood RNA-seq profile from n = 100 COVID-19-positive patients admitted to the Albany Medical Center (NY, USA). Severe patients were defined as COVID-19 patients having less than 25 hospital free days (HFD45 < 25) in a 45 day follow-up.

To determine the gene expression changes induced at a systemic level by the treatment with abatacept, we performed a longitudinal analysis using the linear regression framework implemented in limma 42 . First, we selected genes with a minimal evidence of gene expression in blood, defined as > 0.6 counts per million in at least 20 samples. TMM-normalized values were then log-transformed using voom 43 and a multivariate linear model was fitted accounting for repeated measures per individual (i.e. baseline and week 12 measures), and also adjusting for sex, age and batch effects (i.e. sequencing library plate number). In the early COVID-19 cohort we identified an unequal percentage of granulocytes between cases and controls ( Supplementary Fig. S5 ). For this reason, for analysis involving this cohort, adjustment was also performed for granulocyte proportions (estimated using DeconCell 44 ). A total of n = 44 genes were found to be differentially expressed (FDR < 0.05) without granulocyte-adjustment, and n = 101 when adjusting. The genes differentially expressed between COVID-19 patients and healthy controls were determined using edgeR 45 with adjustment for granulocyte percentage (sex and age information were not available in this dataset). Genes with low expression values (defined as < 2 counts per million in at least 2 samples) were previously filtered out. In the late COVID-19 cohort, we used the RSEM counts provided by the authors 14 and differential expression analysis was performed like in the abatacept dataset. A total of n = 1065 genes were found to be differentially expressed (FDR < 0.05) between severe and mild COVID-19 patients. Testing for up-or down-regulation of biological processes due to the treatment with abatacept and by COVID-19 was performed using the Gene Set Enrichment Analysis (GSEA) 46 method. The GSEA is a commonly used method to identify coordinated differences in expression in functionally-related sets of genes. In this method, a normalized enrichment score (NES) is calculated for each pathway. The NES measures if the genes in a particular pathway are over-represented at the top or at the bottom of a ranked list of genes. In our analysis, we ranked genes based on the corresponding differential expression results. In order to be able to test for antagonism, the directionality of the differential gene expression was preserved by multiplying the -log (P-value) of the differential expression by the sign of the fold change, as described previously 47 . This resulted in a ranked list of n = 18,370 genes, n = 14,829 genes, n = 12,389 genes, in the Abatacept, early COVID-19 and late COVID-19 datasets respectively, that are used as input in the GSEA analysis. For each biological pathway, P-values for association were calculated using the adaptive Monte-Carlo method implemented in the fgsea package with n = 1e6 permutations. GO terms of size < 10 or > 300 where excluded, and the GSEA parameter value was kept to 1 (default). Multiple test p-value adjustment was performed using Benjamini and Hochberg false discovery rate (FDR) method 48 .

To determine the level of antagonism induced by abatacept to COVID-19 pathology in the two COVID-19 cohorts we used two alternative approaches. In the first approach we analyzed exclusively biological processes that have been consistently associated with COVID-19 disease pathology and, principally, with disease severity. After manually identifying the most cited studies and more replicated evidence, we identified a total of 22 relevant processes and its representing GO term. Supplementary Table S1 describes the most relevant bibliographical evidences for each process. In order to avoid a potential bias in considering pathways with a very similar gene composition, we measured the overlap between each pair of pathways using the Jaccard Index ( Supplementary Fig. S6 ). We found a low degree of overlap in most of the pathway annotations, ruling out the bias of redundant process annotation. The list of pathology-associated processes can be divided according to the two immune response stages of COVID-19: (i) anti-viral response mechanisms (i.e. first stage of COVID-19, n = 6) 49 , and (ii) hyperinflammatory process and disease severity (i.e. second stage of COVID-19, n = 16) 3, 19, 21, 50 . The change in activity of these processes induced by abatacept was subsequently tested using the GSEA approach. Antagonistic processes were determined as those were the pathway was significantly modified by abatacept (P-value < 0.05) and in the opposite direction to that described for COVID-19 patients.

In the second approach, we determined the level of antagonism of abatacept to COVID-19 using an agnostic analysis, and therefore not limited to the current knowledge of the disease. In this approach, we determined the biological pathways associated to COVID-19 as those differentially activated by the disease at the transcriptional level using the case-control RNA-seq dataset. We next determined all the biological pathways differentially activated by treatment with abatacept. In both cases, pathway annotation was defined using the GO biological process annotation. Pathways with less than 10 genes or more than 300 genes were excluded. The FDR method was applied to correct for multiple testing; an FDR < 0.05 was considered significant. Since GO terms that are functionally close can have very similar gene composition, we used the Jaccard index to identify redundant annotations. For this objective, a distance matrix between all pairs of BPs was built using 1 − Jaccard Index. Using this Jaccard distance matrix, clusters of redundant BPs were then defined using hierarchical clustering and cutting at height h = 0.6. For each cluster, the GO term with the most significant association was selected for downstream analysis, leading to the selection of n = 6187 GO terms. Similar to the previous approach, antagonistic processes were determined as those that are significantly modified by abatacept (FDR < 0.05) and in the opposite direction to that found in the COVID-19 case-control group (FDR < 0.05). In order to test for the significance of the observed antagonistic processes a binomial test was used. Details of this analysis are included in the supplementary methods.

Survival analysis. In order to investigate the potential benefits of abatacept to prevent severity, we developed a transcriptional-based score that captures the similarity to abatacept-induced signature. www.nature.com/scientificreports/ to build this score (abatacept score, ABS) is described in the supplementary material. Using the score distribution in the complete patient cohort (COVID-19-positive and negative patients, Supplementary Fig. S7) , we selected the patients showing the more extreme ABS values for downstream analysis.

To evaluate the association between the ABS patient categorization with severity, we performed two survival analyses. In both analyses, Cox proportional hazards models by Kaplan-Meier estimation were carried out using outcomes of interest. In the first analysis, Cox proportional hazards were estimated using days without mechanical ventilation, considering intubation (yes or no) as censoring status. In the second analysis, a Cox proportional hazards model was conducted to evaluate the differences between the two ABS groups with regards to the number of days spent in the hospital. As censorship status, we used the categorization of patients into severe and mild COVID-19 groups based on the HFD-45 described in the original manuscript 14 (i.e. < 25 days as severe and ≥ 25 as mild patients).

Single cell RNA-seq analysis of COVID-19 BALF cells. In order to corroborate the relevance of the CD80/86 signaling axis in COVID-19, we analyzed single cell RNA-seq data generate from the lungs, the target tissue of severity. For this objective we analyzed a single cell RNA-seq recently generated from bronchoalveolar lavage fluid (BALF) cells from COVID-19 patients and healthy controls (GEO database accession GSE145926) 10, 11 . Single-cell RNA-seq raw data from BALF samples was generated using the 10× Genomics platform and included 9 COVID-19 patients (n = 6 severe and n = 3 mild). The raw data was downloaded and processed using Seurat (v3) 51,52 and scTransform 53 . scTransform performs normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression, which has shown to effectively reduce the effect of technical variation, without affecting biological heterogeneity. This dataset initially consisted of a total of 62,460 cells. Cells with < 200 or > 6000 unique feature counts, > 10% mitochondrial counts or < 1000 UMIs were subsequently filtered out. After this quality control, a total 54,420 cells (87.1%) were available for downstream analyses. Samples were log-normalized and scaled for the number of genes, number of UMIs and percentage of mitochondrial reads. Cell type clustering was performed using the FindClusters function from Seurat using the default resolution (0.8). This clustering method uses a shared nearest neighbour (SNN) modularity optimization-based clustering algorithm to identify cell clusters of cells based on their PCs. The number of PCs used for each clustering round (k = 50) was estimated by the elbow criterion in the PCA scree plot (Supplementary Fig. S8 ). Using this approach, we identified a total of 34 cell clusters (Fig. 5A,B) . The cell clusters aggregated into four major regions representing (i) T and NK cells (n = 8014 cells, 14.7%, clusters 9, 11, 5, 25 and 18), (ii) two regions of FABP4 + alveolar macrophage cells (n = 17,399 cells, 32%, clusters 1, 2, 3, 30, 7, 32 and 12), and (iii) a large cluster aggregate of active, IL6-expressing macrophages both from FABP4 + alveolar and monocyte-derived FCN1 + macrophages (n = 24,742 cells, 45.5%, clusters 6, 20, 13, 22, 0, 8, 4, 15, 10, 16, 14, 26) . As described by the original study, FABP4 + alveolar macrophages show a marked transition to FCN1 + monocytederived macrophages from healthy to infected and as COVID-19 severity increases (data not shown). Cluster 5 aggregated markers of different types of CD4+ T cells including CCR7 (naïve TCD4), IL2RA (Treg), FOXP3 (Treg), IL7R (naive), LTB (naive), CXCL13 (T peripheral helper) (Fig. 6C ). This cluster also was characterized by the upregulation of CTLA4, the protein that is recombined in abatacept, and also a marker that is highly upregulated in active T cells 54 . Consequently, this cell cluster was classified as activated CD4+ T cells. COVID-19 patients had a much higher number of CD80 + and CD86 + macrophages compared to controls (19.6% vs 1.7% and 45.4% vs 26.5%, respectively, P < 1e−16 in both comparisons) (Fig. 6A ). IL6-expressing macrophages, were present almost exclusively in COVID-19 patients (Fig. 6B) . Association between the key cell types (CD80 + , CD86 + macrophages, active T cells and IL-6 expressing macrophages) was performed using linear regression. Log-transformed cell counts were normally distributed (i.e., Shapiro-Wilk test P > 0.05) and were used in this analysis. To account for potential confounding due to differences in total cell numbers, these were included as a covariate in the regression model.

Characteristics of and important lessons from the coronavirus disease 2019 (COVID-19) outbreak in China: Summary of a report of 72 314 cases from the Chinese Center for Disease Control and Prevention

Estimating clinical severity of COVID-19 from the transmission dynamics in Wuhan, China

Cytokine release syndrome in severe COVID-19

Respiratory support for patients with COVID-19 infection

Pathological inflammation in patients with COVID-19: A key role for monocytes and macrophages

Incidence of COVID-19 in a cohort of adult and paediatric patients with rheumatic diseases treated with targeted biologic and synthetic disease-modifying anti-rheumatic drugs

Clinical predictors of mortality due to COVID-19 based on an analysis of data of 150 patients from Wuhan

Tocilizumab treatment in severe COVID-19 patients attenuates the inflammatory storm incited by monocyte centric immune interactions revealed by single-cell analysis

Abatacept for rheumatoid arthritis refractory to tumor necrosis factor α inhibition

The landscape of lung bronchoalveolar immune cells in COVID-19 revealed by single-cell RNA sequencing

Single-cell landscape of bronchoalveolar immune cells in patients with COVID-19

Deep immune profiling of COVID-19 patients reveals distinct immunotypes with therapeutic implications

Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood mononuclear cells in COVID-19 patients

Large-scale multi-omic analysis of COVID-19 severity

The development and function of lung-resident macrophages and dendritic cells

Targets of T cell responses to SARS-CoV-2 coronavirus in humans with COVID-19 disease and unexposed individuals

COVID-19: Consider cytokine storm syndromes and immunosuppression

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Immunology of COVID-19: Current state of the science

Correlation analysis between disease severity and inflammation-related parameters in patients with COVID-19 pneumonia

Highly pathogenic coronavirus N protein aggravates lung injury by MASP-2-mediated complement over-activation

Antibody responses against SARS coronavirus are correlated with disease outcome of infected individuals

Antibody-dependent SARS coronavirus infection is mediated by antibodies against spike proteins

Anti-severe acute respiratory syndrome coronavirus spike antibodies trigger infection of human immune cells via a pH-and cysteine protease-independent FcγR pathway

Anti-spike IgG causes severe acute lung injury by skewing macrophage responses during acute SARS-CoV infection

The potential danger of suboptimal antibody responses in COVID-19

Antibody-dependent enhancement and SARS-CoV-2 vaccines and therapies

Longitudinal analyses reveal immunological misfiring in severe COVID-19

Effect of interleukin-6 receptor inhibition with tocilizumab in patients with rheumatoid arthritis (OPTION study): A double-blind, placebo-controlled, randomised trial

Abatacept in patients with rheumatoid arthritis and interstitial lung disease: A national multicenter study of 63 patients

Blockade of T cell costimulation by CTLA4-Ig inhibits lung inflammation in murine hypersensitivity pneumonitis

Upregulation of CD40, CD80, CD83 or CD86 on alveolar macrophages after lung transplantation

CD80/CD86 signaling contributes to the proinflammatory response of Staphylococcus aureus in the airway

Abatacept initiation in rheumatoid arthritis and the risk of serious infection: A population-based cohort study

Risk of hospitalized infection and initiation of abatacept versus tumor necrosis factor inhibitors among patients with rheumatoid arthritis: A propensity score-Matched Cohort Study

Long-term safety of abatacept in patients with rheumatoid arthritis

ASCORE, a 2-year, observational, prospective multicentre study of subcutaneous abatacept for the treatment of rheumatoid arthritis in routine clinical practice: 1-year interim analysis

STAR: Ultrafast universal RNA-seq aligner

Elimination of PCR duplicates in RNA-seq and small RNA-seq using unique molecular identifiers

The sequence alignment/map format and SAMtools

featureCounts: An efficient general purpose program for assigning sequence reads to genomic features

limma powers differential expression analyses for RNA-sequencing and microarray studies

Precision weights unlock linear model analysis tools for RNA-seq read counts

Deconvolution of bulk blood eQTL effects into immune cell subpopulations

edgeR: A Bioconductor package for differential expression analysis of digital gene expression data

Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

Etiology matters-Genomic DNA methylation patterns in three rat models of acquired epilepsy

Controlling the false discovery rate: A practical and powerful approach to multiple testing

COVID-19 infection: The perspectives on immune responses

COVID-19 coagulopathy in Caucasian patients

Integrating single-cell transcriptomic data across different conditions, technologies, and species

Comprehensive integration of single-cell data

Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression

Co-inhibitory molecules of the B7-CD28 family in the control of T-cell immunity

The code used to generate the results and figures of this study is openly available at https:// github. com/ Rheum atolo gy-Resea rch-Group/ COVID-19_ Abata cept.

The PACTABA project was funded Bristol-Myers Squibb. We thank all participants from the PACTABA study for their collaboration. AJ and SM are supported by the DoCTIS project funded by the European Union's H2020 programme (Grant #848028). This work was supported by funds from the Vall d'Hebron Hospital Research Institute and from IMIDomics S.L. We thank Dr Ariel Jaitovich (Albany Medical Centre, USA) for providing additional clinical data on the late COVID-19 cohort.

SM and RMM are co-founders of IMIDomics, a biotech company in precision medicine for immune-mediated inflammatory diseases. AJ is also affiliated with IMIDomics. The other authors declare no competing interests.

Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-90797-0.Correspondence and requests for materials should be addressed to A.J. or S.M.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.