The plasma profiling of 62 cytokines, chemokines and other soluble factors involved in several immune responses was first performed in 23 patients with severe COVID-19, responding or not to tocilizumab. The drug was given twice, 12 hours apart, by intravenous or subcutaneous route (see Methods for details). Plasma molecules were quantified at the baseline (T0, before the beginning of treatment), two days (T2) and seven days (T7) after the first dose of the drug. Most patients were sampled during the acute phase of the infection: 13 patients responded to therapy, whereas 10 patients did not. Demographic information, clinical data, blood parameters and arterial haemogas analysis, collected at T0, are reported in Supplementary Table 1 .

Figure 1 shows plasma levels of several cytokines/chemokines and other soluble factors in responder and non-responder patients. Concerning pro-inflammatory cytokines, we found that after 2 days of therapy IL-6 and IL-27 plasma concentrations strongly increased in non-responders. IL-11 levels did not differ between responders and non-responders, but increased in responders after therapy. Non-responders had lower basal levels of interleukin IL-1α, IL-1β, while in responders at T2 IL-18 decreased and IL-17E increased. Level of TGF-α increased over time in responders. Leptin level increased from T2 to T7 in non-responders, while PD-L1 decreased. IL-2 level was lower in non-responders at T0 and significant differences were observed at T7 between responders and non-responders for IL-12p70.

The concentrations of IL-6RA reached the minimum value at T2 in responders. At T2, Leptin receptor (Leptin R) was different between responders and non-responders, with higher level in the latter. IL-10 decreased during therapy in non-responders. IL-2 and IL-13 plasma levels were lower in non-responders at T0, while IL-4 level was lower at T2. Significant differences were observed at T7 between responders and non-responders for IL-12p70. We also noted that IL-23 increased after one week of treatment in both groups of patients. Levels of optineurin (OPN) decreased during therapy in patients responding to tocilizumab. The level of B-cell activating factor (BAFF) decreased in non-responders during therapy.

Apoptosis-related molecules have been investigated. FAS levels increased during therapy in responder and non-responder patients, but higher levels have been detected in non-responders at T2. FAS ligand and TRAIL levels increased in responders, at T7 and T2, respectively.

Plasma concentrations of C-C Motif Chemokine Ligand (CCL)-11 were lower in non-responders at T7, while C-X-C motif chemokine ligand (CXCL)-2 was higher at T2 in responders. CX3CL1 was higher at T2 in non-responders. CXCL10 plasma levels decreased over time in both groups of patients, but non-responders displayed higher levels at T0. This difference was maintained at T2 and at T7.

Regarding growth factors, in responders epidermal growth factor (EGF) was higher at T2 whereas GM-CSF decreased with time. Vascular endothelial growth factor (VEGF) decreased in non-responders over time. Levels of platelet-derived growth factor (PDGF)-AA were lower at the baseline in non-responders. (PDGF)-AB was lower in non-responders at T2. Fms related receptor tyrosine kinase 3 (FLT-3 ligand, also known as CD135) decreased after seven days in non-responders, but increased in responders. All the other soluble molecules remained unchanged in responders and non-responders ( Supplementary Figure 1 ).

To investigate how immune responses vary in the presence of tocilizumab, we profiled B cells and T cells by using polychromatic flow cytometry. After unsupervised analysis, nine distinct cell clusters have been identified (UMAP and heatmap are reported in Figures 2A, B ) on the basis of the surface protein expression ( Supplementary Figure 2 ). Besides a population of cells negative for all the activation/differentiation markers analyzed, naïve B cells are defined as IgM⁺IgD⁺CD24⁺CD21⁺CD38⁻CD27⁻ while transitional B cells are IgM⁺IgD⁺CD24⁺CD21⁺CD38⁺CD27⁻ memory unswitched B cells are defined as IgM⁺IgD⁺CD24⁺CD21⁺CD38⁻CD27⁺ while memory switched B cells are IgM⁻IgD⁻CD24⁺CD21⁺CD38⁻CD27⁺. Memory IgM only B cells are IgM⁺IgD^-CD24⁺CD21⁺CD38⁻CD27⁺. Plasmablasts are IgM^−/+IgD⁻CD24⁻CD21⁻CD38⁺CD27⁺. Finally, exhausted B cells are those CD21⁻ and CD24⁻ (2).

Before therapy (T0), responders showed a higher percentage and absolute number of total B cells and a higher absolute number of memory IgM-only subpopulation compared to non-responders. After two days of therapy, plasmablasts not expressing IgM increases in non-responders compared to responder patients. Responders show an increase in both absolute number and percentage of memory switch B cell population compared to non-responders. After seven days of therapy, we observed an increase in absolute number of total B cells and naïve B cell subpopulation in both responder and non-responder patients, an increase in percentage of plasmablast IgM⁺ in responder patients, a decrease in percentage of plasmablast IgM⁻ and an increase in absolute number of exhausted B cell and surface- subpopulation in both groups. Finally, non-responder patients display a decrease of absolute number of transitional B cells ( Figure 2C ).

T cell compartment of patient with COVID-19 pneumonia displays marked T cell activation, senescence, exhaustion and skewing towards Th17 (3). Given that also T cells expressed CD126 and could be tuned by Tocilizumab administration, we investigated the T cells landscape of patients undergoing anti-IL6R-therapy.

Eighteen clusters of CD4⁺ T cells were recognized after unsupervised analysis on the basis of different expression of differentiation (CD45RA, CCR7, CD25, CD95, CD27, CD28, CD57) and activation markers (HLA-DR, CD28, PD1) ( Figures 3A, B ; Supplementary Figure 3, 4 ). Naïve T cells (N) were defined as expressing CD45RA, CCR7, CD27, CD28, CD127, CD25 and lack of expression of CD95, CD38 and HLA-DR. Moreover, a population of recently activated naïve T cells has been recognized as those expressing CD38. CD38 is a multifunctional molecule, belonging to the family of ectoenzymes and it induces intracellular calcium release. The expression of CD38 identifies a hypo-proliferative CD4⁺ T-cell subset that, following TCR stimulation, retains expression of naïve cell surface markers including CD45RA and CCR7 (20). We were able also to identify T cells stem memory (T_SCM) as naïve cells expressing CD95 and CD38. Central memory T cells (CM) were those expressing CCR7, CD95 and not CD45RA, then four subsets of central memory were found: one expressing CD38, one expressing PD-1, one co-expressing CD38 and PD-1 and one activated (CD38⁺HLA-DR⁺) expressing also PD-1. Three populations of transitional memory (TM) were found as they do not express CD45RA and CCR7, but they express CD28. Moreover, one population of TM is phenotypically resting, one population is expressing CD38 and PD1 and another subset is highly activated as it expresses high level of HLA-DR, CD38 and PD-1. In particular, PD-1 is expressed during the early phase of T cell activation when naïve T cells differentiate into effector cells: it is rapidly expressed after antigen stimulation of naïve T cells with a kinetics of expression similar to the expression of early activation markers (CD69 and CD25) (21, 22). Effector memory T cells (EM) are characterized by the lack of the expression of CD28, CD45RA and CCR7. Two different populations of EM have been identified: one has been identified as a population of phenotypically resting cells and the other as exhausted cells as it expresses CD57 and PD-1. A population of terminally differentiated T cells (TE) were found as effector memory re-expressing CD45RA. Finally, putative T regulatory cells (T_reg) were defined as those expressing CD25 and not expressing CD127. Four populations were found: naïve T_reg expressing CD45RA and CCR7, central memory T_reg as those not expressing CD45RA and expressing CCR7, central memory T_reg expressing also CD38 and effector memory T_reg as those characterized by the lack of expression of both previous markers.

The percentage of CD4⁺ T cells was influenced by therapy, as in responders it increased at T2 and decreased at T7. The absolute number increased at T2 and T7 in responders, and only from T0 and T2 in non-responders. At T7, the absolute number of CD4⁺ T cells was lower in non-responders if compared to responders. The absolute number of naïve T cells increased in responders (until T7) and non-responders (until T2). Naïve CD38⁺ T cells increased in responders from T0 to T2, while the absolute number of naïve Treg at T7 was higher in responders vs non-responders. Absolute number of CM T cells increased after treatment in responders, but not in non-responders. The same trend was observed for the absolute number of CM CD38⁺ T cells and CM CD38⁺ T cells expressing PD-1 for responders and non-responders. The absolute number of activated CM expressing PD-1 increased in non-responders, but did not change in responders. Absolute numbers of CM T cells expressing PD-1 were higher in responders at each time points. Absolute number of TM cells was higher at T0 and T7 in responders and increased during therapy for responder patients. The absolute number of TM expressing CD38 and PD-1 increased in non-responders after therapy, but did not change in responders. The percentage of EM T cells decreased in non-responders from T0 to T7 and in responders from T0 to T2. The absolute number of EM Treg increased in both responders and non-responders after therapy (from T0 to T7).

To sum up, even if the percentage of total CD4⁺ T is modulated by the therapy, the proportion of different CD4⁺ T cell subpopulations remained unchanged. However, the increased of the absolute number of CD4⁺ T cells in responders was likely responsible for all the differences found in the main subpopulations ( Figure 3C , exact p-values are reported in the figures).

Regarding CD8⁺ T cell subsets, twenty-two clusters of CD8 T cells have been identified by unsupervised analysis on the basis of different expression of differentiation (CD45RA, CCR7, CD25, CD95, CD27, CD28, CD57) and activation markers (HLA-DR, CD28, PD-1) ( Figures 4A, B ; Supplementary Figures 5, 6 ). These populations span the entire spectrum of differentiation and activation status. As for CD8⁺ T cell subsets, naïve cells, naïve cells expressing CD38, T_SCM, T_SCM expressing CD38, CM, CM co-expressing CD38 and PD-1, TM, EM and TE were identified. Of the TM population, we found four different subsets: a subset of phenotypically resting cells, a subset of TM expressing HLA-DR, a subset of TM expressing CD38 and PD-1, a subset of activated TM also expressing PD1 and a subset of exhausted TM (CD57⁺PD1⁺). Moreover, we found four populations of EM: EM phenotypically resting, EM expressing only CD57, EM expressing CD38 and CD57, EM co-expressing CD57 and HLA-DR. Finally, seven populations of TE cells have been identified: TE phenotypically resting, TE expressing HLA-DR, TE expressing CD38⁺, TE expressing CD57, TE activated (CD38⁺HLA-DR⁺), TE co-expressing CD57 and HLA-DR, and TE co-expressing CD38 and CD57.

The percentage of CD8⁺ T cells decreased from T0 at T2 in both groups of patients, even if in responders the percentage increased at T7, reaching a similar level of T0. As far as the absolute number is concerned, responders were characterized by an increase of CD8⁺ T cells at T7. Similar absolute numbers of CD8⁺ T cells were observed in non-responders. At T7, responders were characterized by increased absolute number of naïve, naïve CD38⁺, CM, TM and EM T cells. The absolute number of T pedix SCM CD38⁺ was different at T0 and T2 between responders and non-responders, with higher level in responders. Non-responders had higher absolute number of CM CD38⁺PD-1⁺ and TM CD38⁺PD-1⁺, at T7. Moreover, at T7, responders showed higher levels of TM and EM if compared to non-responders. The percentages of TE decreased at T7 in non-responders and at T2 in responders. The percentage of TE CD38⁺CD57⁺ decreased at T2 in non-responders ( Figure 4C , exact p values are reported in the figures).

We wondered whether the two groups of patients could be identified based on immunological and clinical parameters before therapy initiation (T0), predicting who will benefit from therapy and who will not. By using the principal component analysis (PCA), we saw that non-responders clusterize on the left side of the first dimension (PC1, representing 17% of variance), while responder patients clusterize on the right side ( Figure 5A ), meaning that before starting therapy they were extremely different from both immunological and clinical point of view. Besides the level of pO2/FiO2, the immunological parameters that were the main responsible for the clusterization in responders were plasma levels of IL-2, IL-13, IL-1α, the absolute number of B cells as well as the memory unswitched B subset and among CD4⁺ T cells those TM expressing CD38 and PD1. Plasma levels of CXCL10, CXC3CL1, BAFF and the absolute number of different subpopulations of CD8⁺ T cells such as those EM expressing CD38, CD57 and HLA-DR were the major contributors for the clusterization of non-responders along PC1 ( Figure 5B ). Moreover, as previously reported, age, levels of D-dimer, SOFA score and levels of creatine are the major demographic/clinical features that are mostly associated to worse disease outcome and lack of response to therapy (2–6).

This is a case-control, longitudinal, single-center study, approved by the local Ethical Committee (Comitato Etico dell’Area Vasta Emilia Nord, protocol number 177/2020, March 11th, 2020) and by the University Hospital Committee (Direzione Sanitaria dell’Azienda Ospedaliero-Universitaria di Modena, protocol number 7531, March 11th, 2020). Each participant provided informed consent according to Helsinki Declaration, and all uses of human material have been approved by the same Committees. A total of 23 patients with severe COVID-19 was included in the study; they were all treated with tocilizumab. Thirteen patients were responders to therapy, whereas ten patients were non-responders. All patients had severe pneumonia, defined as at least one of the following: (i) presence of a respiratory rate of 30 or more breaths per minute, (ii) peripheral blood oxygen saturation (SaO₂) of less than 93% in room air, (iii) a ratio of arterial oxygen partial pressure (PaO2) to fractional inspired oxygen (FiO₂) of less than 300 mm Hg in room air, and (iv) lung infiltrates of more than 50% within 24–48 h, according to Chinese management guidelines for COVID-19 (version 6.0) (39). All patients received standard of care treatment at the time of hospital admission according to the regional COVID-19 guidelines of Emilia Romagna. In addition to the standard of care treatment, patients also received tocilizumab treatment. Depending on the availability of specific formulation at time of treatment, twelve patients received tocilizumab by intravenous administration at 8 mg/kg bodyweight (up to a maximum of 800 mg) given twice, 12 hours apart, while eleven patients received it by subcutaneous route at a dose of 162 mg administered in two simultaneous doses, one in each thigh (i.e., 324 mg in total). This particular subcutaneous approach was used to mimic, as much as possible, the pharmacokinetic activity of the intravenous formulation in order to achieve similar levels of drug exposure. After the administration of tocilizumab, patients were divided into 2 groups based on its clinical effects: 13 did not required invasive mechanical ventilation and were classified as responder, meanwhile 10 required invasive mechanical ventilation and were classified as non-responder. Three patients of the non-responder group died between 13 and 29 days after admission. We recorded demographic data, medical history, and main laboratory findings from each patient. For details, see Supplementary Table 1 . These data are referred to hospital admission, before starting therapy (T0).

Up to 20 mL of blood were collected from each patient in vacuettes containing ethylenediamine-tetraacetic acid (EDTA). Blood was immediately processed. Isolation of peripheral blood mononuclear cells (PBMC) was performed using ficoll-hypaque according to standard procedures. PBMC were stored in liquid nitrogen in fetal bovine serum (FBS) supplemented with 10% dimethyl sulfoxide (DMSO). Plasma was stored at -80°C until use. Measurements were taken from individual patients; in the case of plasma, each measurement was performed in duplicate and only the mean was considered and shown.

The plasma levels of 62 molecular species were quantified using a Luminex platform (Human Cytokine Discovery, R&D System, Minneapolis, MN) for the simultaneous detection of the following molecules: G-CSF, PDGF-AA, EGF, PDGF-AB/BB, VEGF, GM-CSF, FGF, GRZB, IL-1A, IL-1RA, IL-2, IL-27, IL-4, IL-6, IL-10, IL-13, TNF, IL-17C, IL-11, IL-18, IL-23, IL-6RA, IL-19, IFN-α, IFN-β, IL-3, IL-5, IL-7, IL-12p70, IL-15, IL-33, TGF- β, IFN- γ, IL-1B, IL-17, IL-17E, CCL3, CCL11, CCL20, CXCCL1, CXCL2, CCL5, CCL2, CCL4, CCL19, CXCL1, CXCL10, PD-L1, FLT-3, TACI, FAS, LEPTIN R, APRIL, OPN, BAFF, LEPTIN, BMP4, CD40 LIGAND, FAS LIGAND, BMP7, BMP2, TRAIL, according to the manufacturer’s instruction. Data represent the mean of two technical replicates.

All data obtained by flow cytometric assays have been produced according to the state-of-the-art technologies, as described (40).

Thawed PBMC were washed twice with RPMI 1640 supplemented with 10% fetal bovine serum and 1% each of l-glutamine, sodium pyruvate, nonessential amino acids, antibiotics, 0.1 M HEPES, 55 μM β-mercaptoethanol and 0.02 mg/ml DNAse. Thawed PBMC were stained with viability marker Promokine IR-840 (PromoCell GmbH, Heidelberg, Germany) for 20 min at room temperature in PBS. One million PBMC were washed with FACS buffer and stained with DuraClone IM B cells containing the following lyophilized directly conjugated mAbs: IgD-FITC, CD21-PE, CD19-ECD, CD27-PC7, CD24-APC, CD38-AF750, IgM-PB, CD45-KrO. Cells were washed with FACS buffer and acquired at Cytoflex LX flow cytometer (Beckman Coulter, Hialeah, FL). A minimum of 500,000 cells was acquired on a CytoFLEX LX flow cytometer (Beckman Coulter).

Thawed PBMC were washed twice with RPMI 1640 supplemented with 10% fetal bovine serum and 1% each of l-glutamine, sodium pyruvate, nonessential amino acids, anti- biotics, 0.1 M HEPES, 55 μM β-mercaptoethanol and 0.02 mg/ml DNAse. Up to 1 million PBMC were stained with the Duraclone IM T cell panel (Beckman Coulter) added with another five fluorescent mAbs and a marker of cell viability. Along with side and forward scatter signals, signals were obtained from different fluorochrome-labeled mAbs, i.e., CD45 conjugated with Krome Orange, CD3 APC-A750, CD4-APC, CD8-AF700, CD27-PC7, CD57-Pacific Blue, CD279(PD1)-PC5.5, CD28-ECD, CCR7-PE, CD45RA-FITC, HLA-DR-BUV661, CD127-BV650, CD25-BV785, CD95-BUV395, CD38-BUV496, and PromoFluor-840 (Promokine, PromoCell, Heidelberg, Germany). A minimum of 500,000 cells per sample were acquired on a CytoFLEX LX flow cytometer (Beckman Coulter).

Flow Cytometry Standard (FCS) 3.0 files were imported into FlowJo software version X (Becton Dickinson, San Joseè, CA), and analyzed by standard gating to eliminate aggregates and dead cells, and to identify CD3⁺CD4⁺ T cells, CD3⁺CD8⁺ T cells ( Supplementary Figures 7 ) and CD19⁺ B cells ( Supplementary Figures 8 ). Data were exported for further analysis in R, by following a script that makes use of Bioconductor libraries and R statistical packages (CATALYST 1.10.1). The script is available at: https://github.com/HelenaLC/CATALYST). The selection of cofactor for data transformation was checked on Cytobank premium version (see: cytobank.org). Metaclustering (K=25) was performed by using FlowSOM algorithm. Dimensional reduction was performed using UMAP algorithm.

B cell UMAP graphs stratified by patient at time T0, T2 and T7 are reported in Supplementary Figures 9A–C , respectively. B cell projection of UMAP graphs stratified by patient showing the FlowSOM clusters at time T0, T2, T7 are reported in Supplementary Figures 10A–C , respectively. CD4⁺ T cell UMAP graphs stratified by patient at time T0, T2 and T7 are reported in Supplementary Figures 11A–C , respectively. CD4⁺ T cell projection of UMAP graphs stratified by patient showing the FlowSOM clusters at time T0, T2, T7 are reported in Supplementary Figures 12A–C , respectively. CD8⁺ T cell UMAP graphs stratified by patient at time T0, T2 and T7 are reported in Supplementary Figures 13A–C , respectively. CD8⁺ T cell projection of UMAP graphs stratified by patient showing the FlowSOM clusters at time T0, T2, T7 are reported in Supplementary Figures 14A–C , respectively.

We have then investigated the role of B and T cells along with several clinical and biochemical parameters. For this purpose, we used the principal component analysis (PCA), a dimension reduction method that retains the characteristics of a data set that contribute most to its variance, by keeping lower order principal components (PCs) and ignoring the others (41). PCA uses an orthogonal transformation to collapse the dataset containing correlated parameters to a smaller set of linearly uncorrelated variables known as PCs, such that each PC is a weighted combination of all the markers. We performed this analysis to test whether subject classification was possible based on the B and T cell profile, cytokines and/or clinical data. Thus, PCA was carried out on the clinical events (responder or non-responder), and a dataset comprising 137 parameters that included age, blood pH, pCO₂ (mmHg), pO₂ (mmHg), sO₂ (%), pO₂/FiO₂ ratio, ALT (U/L), total bilirubin (mg/dL), CK (U/L), creatinine (mg/dL), D-Dimer (ng/mL), Hb (g/dL), white blood cells, WBC (N/μL), red blood cells, RBC (10⁶/μL), INR (ratio), LDH (U/L), CRP (mg/dL), platelets (10⁹/L), respiratory rate (breaths/minute), systolic blood pressure (mmHg), heart rate (beats/minute), SOFA score, different absolute numbers of B-cell populations (naïve, memory switched, memory IgM-only, transitional, memory unswitched, plasmablasts, exhausted, surface-) and T-cell subpopulations.

High-dimensional cytometric analysis was performed by using differential discovery in high-dimensional cytometry via high-resolution clustering. Quantitative variables were compared using Two-way ANOVA with Bonferroni correction. Data are represented as individual values, means, and standard errors of the mean. Statistical analyses were performed using Prism 8.4.3 (GraphPad Software Inc., La Jolla, USA).