Stratification of COVID-19 patients based on quantitative immune-related gene expression in whole blood

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) causes mild symptoms in the majority of infected individuals, yet in some cases it leads to a life-threatening condition. Determination of early predictive biomarkers enabling risk stratification for coronavirus disease 2019 (COVID-19) patients can inform treatment and intervention strategies. Herein, we analyzed whole blood samples obtained from individuals infected with SARS-CoV-2, varying from mild to critical symptoms, approximately one week after symptom onset. In order to identify blood-specific markers of disease severity status, a targeted expression analysis of 143 immune-related genes was carried out by dual-color reverse transcriptase multiplex ligation-dependent probe amplification (dcRT-MLPA). The clinically well-defined subgroups of COVID-19 patients were compared with healthy controls. The transcriptional profile of the critically ill patients clearly separated from that of healthy individuals. Moreover, the number of differentially expressed genes increased by severity of COVID-19. It was also found that critically ill patients can be distinguished by reduced peripheral blood expression of several genes, which most likely reflects the lower lymphocyte counts. There was a notable predominance of IFN-associated gene expression in all subgroups of COVID-19, which was most profound in critically ill patients. Interestingly, the gene encoding one of the main TNF-receptors, TNFRS1A, had selectively lower expression in mild COVID-19 cases. This report provides added value in understanding COVID-19 disease, and shows potential of determining early immune transcript signatures in the blood of patients with different disease severity. These results can guide further explorations to uncover mechanisms underlying immunity and immunopathology in COVID-19.


Introduction
The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) causes mild symptoms in the majority of infected individuals, albeit in some cases it requires hospitalization and may lead to a life-threatening condition (Gorbalenya et al., 2020;Wang et al., 2020a). This single-stranded RNA virus' S-protein binds to angiotensin-converting enzyme 2 (ACE2), which is highly expressed by epithelial cells in the lungs but also present in other organs (Hoffmann et al., 2020). The coronavirus disease 2019 (COVID-19) pandemic has so far led to 4.4 million confirmed deaths (reported to the World Health Organization by August 16, 2021), and has put a strain on all layers of societies (UN, 2020;WHO, 2021). Critical illness and deaths are mainly caused by severe pneumonia with respiratory distress (Gorbalenya et al., 2020;. A dysregulated immune response with excessive inflammation is believed to play an important role in pathogenesis of severe COVID-19 (Laing et al., 2020;Qin et al., 2020;Zhou et al., 2020).
Older age quickly emerged as the most important risk factor for poor outcome, but also male sex, obesity, diabetes mellitus and several other co-morbidities are associated with increased risk for severe disease and mortality (Rawshani et al., 2021;Williamson et al., 2020). Moreover, objective clinical and laboratory parameters, like hypoxemia, low lymphocyte counts and high levels of C-reactive protein (CRP), aid in the diagnostic and prognostic of the disease . Determination of early predictive biomarkers enabling risk stratification for COVID-19 patients can inform treatment and intervention strategies.
Herein, we analyzed whole blood samples obtained from individuals infected with SARS-CoV-2, with varying degree of severity, around one week after symptom debut for the expression of 143 immune-related genes by dual-color reverse transcriptase multiplex ligation-dependent probe amplification (dcRT-MLPA) (Haks et al., 2015). The underlying premise was to determine early blood transcriptional signatures of disease severity in COVID-19 patients.

Patients and study design
A total of 50 whole blood samples were selected from a cohort of adult individuals at the Department of Infectious Diseases, Sahlgrenska University Hospital, Gothenburg, Sweden (Patel et al., 2021). All samples were collected during the period March to May, 2020. Out of these, 35 displayed symptoms compatible with COVID-19, and were subsequently PCR-verified as SARS-CoV-2 infected. The peak severity COVID-19 symptoms varied from mild (score 2-3, neither treatment nor in-patient hospital care) to moderate/severe (score 4-6, requiring low-flow nasal oxygen to high-flow nasal oxygen) to critical (score ≥ 7, requiring invasive mechanical ventilation or deceased) according to the WHO Clinical Progression Scale (Marshall et al., 2020). All blood samples for this study were retrieved before the start of treatment, e.g. glucocorticoid administration. The remaining 15 individuals were healthy and were tested negative for SARS-CoV-2 by PCR. All clinical data were retrieved from the patient records.

Sample collection
Whole blood samples were collected in PAXgene Blood RNA tubes (PreAnalytiX, Qiagen/BD) according to the manufacturer's recommendations. Each tube was filled with a volume of 2.5 ml whole blood, and stored at room temperature for a minimum of 2 h before transferred to − 70 • C.

RNA isolation
Total RNAs were isolated from whole blood using the PAXgene Blood RNA kit (PreAnalytiX, Qiagen/BD) according to the manufacturer's recommendations. Briefly, the samples were thawed, the cells were pelleted and then lysed. The cell content was treated with proteinase K. After filtering out cell debris, the flow-through was mixed with ethanol to precipitate nucleic acids and added to a spin column, designed to bind nucleic acids. The membrane in the spin column was treated with DNase I, and the remaining RNAs were eluted in a total volume of 80 μl. The whole procedure included several washing steps. Quantification of RNA (mean 5.8 μg/sample) and purity (A 260 /A 280 ratio 1.8-2.2) were carried out using a NanoDrop-1000 instrument.

dcRT-MLPA
For each target sequence, a specific RT primer was designed located immediately downstream of the left-and right-hand half-probe target sequence. RNA was reverse transcribed to cDNA by incubating at 37 • C for 15 min, using RT-primer mix and the Moloney Murine Leukemia Virus reverse transcriptase kit (Promega). The enzyme was inactivated by heating at 98 o C for 2 min. The probes were hybridized to the cDNA at 60 • C overnight and annealed half-probes were ligated at 54 • C for 15 min using ligase-65. The ligase was subsequently inactivated by heating at 98 • C for 5 min Ligated probes were amplified by PCR (33 cycles at 95 • C for 30 s, 58 • C for 30 s and 72 • C for 60 s, followed by one cycle at 72 • C for 20 min). To monitor assay performance, substituting RNA with nuclease-free water (Thermo Fisher) was used as negative control while Human Universal References RNA (Clontech, Palo Alto, CA, USA) and synthetic template oligonucleotides as hybridization templates were used as positive controls. Primers and probes were from Sigma-Aldrich -7 (5-9) 8 (5-9) 7 (4-9) Low-flow oxygen, no.
-0 8 0 High-flow oxygen, no.  Whole blood samples were collected from healthy control subjects (n = 13) and COVID-19 patients (n = 29) with varying degree of symptoms. The COVID-19 symptoms were classified, based on clinical presentation, as either mild (n = 7), moderate/severe (n = 10) or critical (n = 12). Total RNA samples were isolated from the blood, converted into cDNA, and profiled by a dual-color Reverse Transcriptase Multiplex Ligation-dependent Probe Amplification (dcRT-MLPA) assay that targeted 143 immune genes. A supervised global analysis of the mRNA expression data was performed with the multivariate analysis tool Orthogonal Partial Least Square Discriminant Analysis (OPLS-DA). The OPLS-DA score plot explained 68% (R 2 X) of the variability. The model displays inter-group variation for mild, moderate and severe COVID-19 in relation to healthy controls (xaxis) as well as the variability within each group (y-axis). Each data point represents one sample and the circles define respective group cluster. and the SALSA MLPA reagent kit from MRC-Holland (Amsterdam, The Netherlands). RT primers and half-probes were designed by Leiden University Medical Centre (LUMC, Leiden, The Netherlands) (Joosten et al., 2012), and comprised sequences for four housekeeping genes and 143 selected genes to profile the innate, adaptive and inflammatory immune responses. PCR products were diluted 1:10 in highly deionized formamide (Thermo Fisher) containing 400HD Rhodamine X fluorophore size standard (Thermo Fisher). PCR products were denatured at 95 • C for 5 min, stored immediately at 4 • C and analyzed on an Applied Biosystems 3730 capillary sequencer in GeneScan mode (BaseClear, Leiden, The Netherlands). Trace data were analyzed using GeneMapper software 5 (Applied Biosystems). The areas of each assigned peak (arbitrary units) were exported for analysis in R (version 3.5.1). The housekeeping genes included in the analyzes were beta-2-microglobulin (B2M), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), glucuronidase beta (GUSB), and RhoGEF and GTPase activating protein (ABR). Data were normalized to GAPDH as it was shown to be the most reliable housekeeping gene for whole blood derived RNAs in the assay (Haks et al., 2015;Joosten et al., 2012). GeneMapper (Log2 transformed peak area 7.64) were used to assign the threshold value for noise cutoff. Finally, the normalized data were Log2-transformed for statistical analysis. Eight samples were removed from data analysis due to quality issues (two from control, four from mild and two from critical). Sample number, per group, according to Table 1.

Orthogonal Partial Least Square Discriminant Analysis (OPLS-DA) and Principal Component Analysis (PCA)
OPLS-DA is a discriminant version of orthogonal projections to latent structures (OPLS) and a modification of traditional projection to latent structures (PLS-DA) (Bylesjö et al., 2006;Trygg and Wold, 2002). OPLS-DA was used to model and visualize the separation between the COVID-19 subgroups and the controls. The models were evaluated using the OPLS-DA performance measures R 2 and Q 2 (the closer to 1, the better model performance). Cross-validation ANOVA (CV-ANOVA) was used for significance testing for the OPLS models, and as a diagnostic tool for assessing the reliability of the OPLS models. ANOVA is a method to test if two models are significantly different when fitted to the data used for cross-validation. The default SIMCA cross-validation procedure is a 7-fold cross-validation. CV-ANOVA for significance testing of PLS and OPLS1 models (Bylesjö et al., 2006;Eriksson et al., 2008;Trygg and Wold, 2002). The Simca v17 (Sartorius/Umetrics, Goettingen, Germany) was used for the analyzes. Further, gene expression variances were displayed as PCA plots to highlight similarities and differences depending on COVID-19 disease status, age classes (according to the United Nation definition: young 18-39 years; middle age 40-60 years; and older than 60 years) as well as gender (female/male), using the R-packages prcomp, ggfortify and ggplot2 on Log2 transformed data.

Differential expression analysis
To illustrate the distribution of gene expression, a heatmap was

Fig. 2. Pattern of expression for individual immune-related genes in COVID-19 patients and in healthy controls.
A heatmap displaying color-coded expression levels of genes in healthy individuals and COVID-19 patients with mild, moderate/severe or critical symptoms. The experimental setup is described in Fig. 1. The heatmap was created with normalized Log2 transformed dcRT-MLPA data. Only genes with a pre-specified call rate > 80% were included. The heatmap was generated by hierarchical clustering on both rows (genes) and columns (groups). The color intensity corresponds to the Z-score that indicates the number of standard deviations by which the normalized raw value was below (green) or above (red) the mean value row by row. Black color denotes that the normalized raw value coincided with the mean value of that specific gene, while white gaps represent no detectable signal in the dcRT-MLPA assay.
created with normalized Log2 transformed dcRT-MLPA data. Of the 143 genes, 101 showed a pre-specified call rate > 80% and were included. The heatmap was generated by hierarchical clustering on both genes and groups, using the R-package gplots. Differentially expressed genes in the three groups (mild, moderate/severe and critical COVID-19) were identified using an unpaired t-test analysis on Log2 normalized data compared with controls, using a false discovery rate (FDR) correction. All 143 genes were considered for the FDR correction and the significance cut-off was set to 0.05. The differentially expressed genes were visualized using a volcano plot, that displays Log2 fold change against -Log10 p-value from the t-tests, using the R-package ggplot2. A Venn diagram of the overlapping genes among the groups (only genes with FDR-adjusted p-value below p < 0.05 were included) were created using the R-package VennDiagram. Differentially expressed genes were then visualized in Circos plots generated using CIRCOS v 0.69-6, to demonstrate the similarities and differences between the COVID-19 subgroups compared to controls. The R package ggpubr was used to add significance levels to box plots comparing the gene expression between COVID-19 groups. FDR-adjusted P-values < 0.05 were considered statistically significant. The R software version 4.0.2 (https://www.r-project.org/, The R project, Vienna, Austria) was used for the analyzes.

Results and discussion
To determine blood-specific markers of disease severity status, a targeted expression analysis of 143 immune-related genes was carried out on total RNA extracted from whole blood of patients with acute COVID-19 and healthy control samples. The disease status varied in the COVID-19 cohort, classified as either mild, moderate/severe or critical. The blood samples were retrieved approximately one week from the onset of illness, which coincides with the onset of clinical aggravation of the disease . The demographics and clinical characteristics of the study cohort are summarized in Table 1. The patients with a critical condition were approximately ten years older (mean age 57) than the other two COVID-19 subgroups. A majority of the individuals suffering from either moderate/severe or critical illness were male (80% and 76%, respectively). The observed increased vulnerability linked to older age and male gender as described by others (Docherty et al., 2020;Peckham et al., 2020). In line with previous reports, the subgroups moderate/severe and critical COVID-19 had relatively high CRP (109 respectively 257 mg/L) and low blood lymphocyte counts (1.1 respectively 1.0 × 10 9 /L) (Guan et al., 2020;Herold et al., 2020;Yang et al., 2020). Of note, the CRP value and the lymphocyte count denote the highest and the lowest recorded value during the hospital care, respectively (Table 1).
The dcRT-MLPA data were first modeled using OPLS-DA, a supervised statistical multivariate data analysis technique that uses labeled groups. The model explained 68% of the variability (R 2 X = 68%, R 2 Y = 32%, Q 2 = 25%) and was significant according to cross-validation ANOVA (p-value = 0.045). The model was significant but of modest Q2, most likely due to heterogeneity in study groups. The inter-group variation, visualized in the horizontal direction of the score scatter plot, revealed that the transcriptional profile of critical COVID-19 patients can be separated from that of healthy controls, while those of mild and moderate/severe COVID-19 patients showed close to overlap in the model (Fig. 1). The intra-group variation, displayed in the vertical dimension, was large for the three COVID-19 patient subgroups. When projecting the data in an unsupervised manner (Fig. S1), PCA axis 1 and 2 explain 68% of the variance. The PCA plot did not show a clear separation between the subgroups, even though a tendency of separation could be observed, especially for the critical subgroup. The intra-group variation is highlighted in Fig. 2, a heatmap displaying normalized expression values for 101 genes that had available data from at least 80% of the samples. The clustering, based on both genes and subgroups, identified varying transcriptional response types within each subgroup. This heterogeneity in the transcriptional changes of immune-related  Fig. 1). The vertical axis (y-axis) corresponds to the Log 10 FDR-adjusted pvalues, and the horizontal axis (x-axis) displays the Log 2 fold change value. The vertical dashed line indicates 0 Log 2 fold change, and the horizontal dashed line represent 1.3-Log 10 FDR-adjusted p-value (corresponding to the significant threshold of 0.05 FDR-adjusted p-value). Positive x-values represent higher expression and negative x-values represent lower expression. The color of the dots indicates whether the gene had higher expression (red), lower expression (blue) or was unregulated (gray). The abbreviation next to each dot corresponds to the gene with significantly changed expression. genes in each clinically well-defined subgroup can, at least in part, be explained by the varying impact of SARS-CoV-2 infection in different individuals (Mathew et al., 2020). No confounding co-variation was observed with regards to age or gender of the infected subjects ( Fig. S2A and B).
Despite the variance observed within each subgroup, numerous differentially expressed genes showed a statistically significant difference (FDR-adjusted p < 0.05) in COVID-19 patients relative to control samples. The number of differentially expressed genes increased by severity of COVID-19 (Fig. 3). More specifically, out of the 143 analyzed immune-related genes, we found ten genes (three down-and seven upregulated) for mild, eleven genes (nine down and two up) for moderate/severe and 27 genes (13 down, 14 up) for critical, to be significantly altered in comparison with the control group.
As shown in Fig. 4, only three differentially expressed genes were shared between all of the COVID-19 subgroups, namely IFI6 (interferon alpha inducible protein 6), IFITM3 (interferon-induced transmembrane protein 3) and NLRP1 (nod-like receptor family pyrin domain containing 1). The two interferon (IFN) signaling-related genes, IFI6 and IFITM3, were found to be expressed higher in all COVID-19 subgroups. These genes have previously been suggested as markers for COVID-19 infection (Hachim et al., 2020;Shaath et al., 2020b). NLRP1, which is involved in pattern recognition of innate immunity as part of an inflammasome, had lower expression. It has recently been reported that SARS-CoV-2 can activate double-stranded (ds) RNA-mediated innate immune responses (Bauernfried et al., 2021;Li et al., 2021b). It is therefore reasonable to speculate that reduced expression of NLRP1, as one of the sensors of dsRNA, may provide an innate immune escape mechanism for SARS-CoV-2.
The subgroup displaying mild symptoms had one downregulated gene, HCK (hematopoietic cell kinase), in common with the moderate/ severe subgroup (Fig. 4). This kinase is specifically expressed by myeloid and B-lymphocyte cell lineages, activated by certain inflammatory stimuli and can initiate several downstream signaling pathways that are involved in immunity to viruses (reviewed in Poh et al., 2015). In addition, mild cases had five genes with higher expression, IFI35 (interferon induced protein 35), IFI44, IFIT2 (IFI with tetratricopeptide repeats 2), OAS2 (2'-5'-oligoadenylate synthetase 2) and OAS3, in common with the critical subgroup (Fig. 4). The moderate/severe subgroup showed a trend towards upregulation of these five genes, albeit not statistically significant (FDR-adjusted p-values 0.066-0.405). This can be explained by the observed intra-group variation and the small sample size. In order to identify expression patterns within the three COVID-19 subgroups more clearly, genes with FDR-adjusted p-values < 0.1 were also marked in the Circos plots (Fig. 5) that have been based on gene functionality. The interferon-stimulated genes (IFI35, IFI44 and IFIT2) and antiviral restriction enzyme activators (OAS2 and OAS3) have recently been reported in association with SARS-CoV-2 infection (Li et al., 2021a;Pairo-Castineira et al., 2021;Shaath et al., 2020a;Ziegler et al., 2020). Hence, significantly higher expression of interferon-related genes dominates transcriptional changes that occurred within a week after COVID-19 symptoms onset (Fig. 5).
The moderate/severe and the critical subgroup had seven commonly down-regulated genes: CCL5 (C-C motif chemokine ligand 5); CD3E (cluster of differentiation 3E); CD4; CTLA4 (cytotoxic T-lymphocyteassociated protein 4); GATA3 (GATA binding protein 3); PTPRCv1 (protein tyrosine phosphatase receptor type C transcript variant 1, encoding CD45RA); and TAGAP (T cell activation RhoGTPase activating protein) (Fig. 4). This gene set is preferentially expressed by T cells (Fig. 5) (Connelly et al., 2014;Mao et al., 2004;Takahashi et al., 2000;Tian et al., 2017;Van de Walle et al., 2016). It is hence plausible that the observed down-regulation of these genes stemmed from the decline in peripheral blood lymphocytes, indicated by the reduced CD3E expression which correlated with the observed lymphopenia. There are several reports on lymphocytopenia in serious cases of COVID-19 (Guan et al., 2020;Yang et al., 2020), with profound T cell losses (Laing et al., 2020;Liu et al., 2020b). The reduction of blood lymphocytes following SARS-CoV-2 infection may for instance be due to Fig. 4. Similarities and differences in gene expression patterns of COVID-19 patients with different disease severity. A Venn diagram summarizing genes with significantly changed mRNA expression (FDRadjusted p-value < 0.05) in COVID-19 patients with either mild, moderate/severe or critical symptoms. The data presented are based on the results displayed in Fig. 3. Each circle represents one group and any unique regulated genes can be found in its outer portion, while shared regulated genes appear in the intersections between the circles. Red color indicates genes with higher expression, while blue color denotes genes with lower expression.
infiltration and sequestration of these cells in target organs, death due to being target cells of the virus and/or because of the burst of pro-inflammatory cytokines (Huang and Pranata, 2020;Tan et al., 2020). It is noteworthy that overexpression of inhibitory checkpoint molecules, such as CTLA-4 has been reported during SARS-CoV-2 infection (Kong et al., 2020;Zheng et al., 2020a). However, the expression can apparently be transient in CD4+ T cells and unchanged in CD8+ T cells, during the first week . It has also been shown that regulatory T cell (FOXP3+) mechanisms are impaired in the lung, which may lead to T-cell hyperactivation (Kalfaoglu et al., 2020).
Likewise, the lower expression of BCL2 (B-cell lymphoma 2), CCL4, CD8A, GNLY (granulysin) and IL7R (interleukin-7 receptor), uniquely observed for critical COVID-19, most likely relates to the changes of circulating lymphocytes (Figs. 4 and 5). Low expression of the antiapoptotic BCL2 can have a direct effect on cell numbers (Youle and Strasser, 2008). It has been described that the SARS-CoV-2 ORF3a protein can induce apoptosis, in cell lines, without any effect on BCL2 expression levels Ren et al., 2020). As lymphocytes produce the chemokine CCL4 (Eberlein et al., 2020) as well as granulysin, a pore-forming protein involved in chemotaxis (Sparrow and Bodman-Smith, 2020), the reduced expression of these genes may be a

Fig. 5. Functional profile of immune-related genes in the blood of COVID-19 patients.
Circos plots displaying expression of all analyzed genes in blood samples from COVID-19 samples compared with healthy individuals, as described in Fig. 1. Each quarter of the plot represents one group, controls as well as mild, moderate/ severe and critical COVID-19, as indicated. The immune-focused genes were categorized into four main groups: (A) Innate pattern recognition and response; (B) Cell functions; (C) Cytokine signaling and immune cell markers; and (D) T cell cytokines and markers. Each of the main gene categories were further subdivided according to the outer circle but single genes do not have a subheading. Genes with lower expression are represented by green (p < 0.1) and blue (p < 0.05) lines, while genes with higher expression are represented by orange (p < 0.1) and red (p < 0.05) lines. All unregulated genes (p > 0.1) have black lines. direct effect of the lymphopenia. Similarly, the IL7R expression could be a direct consequence of the lymphocyte reduction in the blood. This may also be due, however, to the state of lymphocyte differentiation during COVID-19 as the IL-7R expression varies quite a bit depending on whether the T cells are naïve, activated or of memory type (reviewed in Mazzucchelli and Durum, 2007). IL-7 treatment has been used to restore normal lymphocyte levels in both septic shock (Francois et al., 2018) and COVID-19 (Laterre et al., 2020).
The critical COVID-19 subgroup also uniquely displayed higher expression of a few genes: FCGR1A (Fc-gamma receptor 1A); GBP1 (guanylate-binding protein 1); GBP2, STAT1 (signal transducer and activator of transcription 1); and TAP1 (antigen peptide transporter 1) (Fig. 4). The expression of the above-mentioned genes is increased in response to IFN stimuli (Fig. 5). It is believed that IFNs play a crucial role in COVID-19 (Lee et al., 2020;Wilk et al., 2020), as part of the antiviral defense, and are, at least in part, responsible for the aggravation during the second week of illness (Tincati et al., 2020). The transcription factor STAT1 is located downstream in the IFN signaling pathway, and thereby involved in the induction of IFN-responsive genes (Darnell et al., 1994). The high affinity FcγRI (CD64) is continuously expressed on the surface of most myeloid cells and some granulocytes. Upregulation of FcγRI on monocytes, macrophages and NK cells enables efficient phagocytosis, binding to immune-complexes, antibody-dependent cellular cytotoxicity and cytokine secretion during infection (Kårehed et al., 2007;Pearse et al., 1991). High levels and/or dysregulation of peripheral pro-inflammatory monocytes have been described in COVID-19 patients suffering from severe disease (Schulte-Schrepping et al., 2020;Tincati et al., 2020). There are also reports of neutrophilia and/or highly active NK cells in severely ill COVID-19 patients (Maucourant et al., 2020;Qin et al., 2020;Tincati et al., 2020;Wang et al., 2020b). TAP1 is engaged in peptide transport for antigen presentation, forming a complex with MHC class I molecules (Eggensperger and Tampé, 2015), and hence facilitating activation of cytotoxic T cells. However, several studies described functional exhaustion of cytotoxic lymphocytes during COVID-19 (Kusnadi et al., 2021;Liu et al., 2020a;Mathew et al., 2020;Zheng et al., 2020b).
The gene TNFRSF1A (tumor necrosis factor receptor superfamily member 1A) had selectively lower expression in the subgroup displaying mild symptoms (Fig. 4). This finding is of particular interest as TNFRSF1A is one of the major receptors for TNF-α, which is believed to Fig. 6. Comparative analysis of differentially expressed genes in patients with varying degree of illness. Box plots comparing the expression levels of a set of immune-related genes between clinically defined COVID-19 subgroups, according to Fig. 1. The selection was based on the genes found to differ significantly between critically ill cases (right) relative to individuals with mild symptoms (left) and/or the subgroup with moderate/severe symptoms (middle). The x-axis displays Log 2 transformed GAPDH-normalized relative gene expression. Each box contains the interquartile range (IQR), quartile 1-quartile 3, and the center line indicates the median value. The whiskers mark 1.5 × IQR and data points outside this range represent outliers. The comparisons were performed pairwise with asterisks denoting significance level. Differences were considered statistically significant at p values of p < 0.05 (*), p < 0.01(**) and p < 0.001(***). play a pivotal role in the hyperinflammatory response associated with more-severe COVID-19 cases . TNF-α has been identified to independently predict patient outcomes (Del Valle et al., 2020), and TNF-α inhibitors were reported to protect against severe forms of COVID-19 (Brito et al., 2021;Robinson et al., 2020). It is to be expected that the level of TNF-α signaling would be influenced by receptor availability. Several viruses actually have the capacity to downregulate TNFRAF1A (Rahman and McFadden, 2006). There is no previous report, to our knowledge, that describes stratification of COVID-19 patients based on TNF-α receptor expression.
In order to pinpoint what distinguishes the individuals in need of hospitalization from outpatients, gene expression profile was compared between the COVID-19 patient subgroups. A set of 15 immune-related genes was found to differ significantly between the two extremes, namely critically ill individuals and those with mild symptoms (Fig. 6). All but three genes (FCGR1A, TAP1 and TLR2) displayed lower expression in the critical subgroup. For most of the genes, the expression pattern in the moderate/severe subgroup resembled that of critical cases, with three exceptions (CCL5, TAP1 and TLR2). The genes expressed to a lesser extent in more-severe disease can presumably be explained by the decline in circulating lymphocytes. The high expression of FCGR1A and TLR2 may probably be linked to the level/activation states of monocytes and NK cells found in peripheral blood of patients with critical SARS-CoV-2 infection as a feature of the documented imbalanced immune response to SARS-CoV-2 (Maucourant et al., 2020;Qin et al., 2020;Schulte-Schrepping et al., 2020;Tincati et al., 2020;Wang et al., 2020b;Zhou et al., 2020).
In summary, we found numerous differentially expressed immunerelated genes in acute COVID-19 patients compared with healthy controls, although heterogeneity in the transcriptional profile was observed within each of the clinically well-defined subgroup. Our results propose that critical COVID-19 patients can be distinguished by reduced peripheral blood expression of several T-cell-associated genes, which most likely relates to the lower lymphocyte counts, observed in this study and reported in other studies (Guan et al., 2020;Herold et al., 2020;Yang et al., 2020). There was also a notable dominance of IFN-associated genes in critically ill COVID-19 patients, a majority displayed higher expression, but several of these gene expression patterns could also be observed in less severely ill patients. Interestingly, the expression of the gene encoding one of the main TNF-receptors, TNFRS1A, was selectively lower in mild COVID-19 patients. Although it needs validation in larger COVID-19 patient cohorts, TNFRS1A could potentially help stratify COVID-19 patients. We are however aware that the limited sample size in the patients' subgroups could represent a limitation of this study.
Whole genome transcriptomics analyses of blood samples from COVID-19 patients have recently been reported (Aschenbrenner et al., 2021;Monaghan et al., 2021;Ng et al., 2021;Pairo-Castineira et al., 2021). However, our highly quantitative assessment of 143 immune-related genes by dcRT-MLPA method provides added value and can help inform determination of early immune transcript signatures in the blood of COVID-19 patients with different disease severity. These results warrant further exploration to uncover mechanisms underlying immunity and immunopathology in COVID-19.

Declaration of Competing Interest
The authors report no conflict of interest.

Ethics declarations
All study participants were included after written informed consent and agreed on the use of anonymized information as per the Swedish law on the General Data Protection Regulation (GDPR). The study was approved by the Swedish Ethical Review Authority (2020-01771) and participants were included after written informed consent.

Appendix A. Supporting information
Supplementary data associated with this article can be found in the online version at doi:10.1016/j.molimm.2022.03.004.