Maturation signatures of conventional dendritic cell subtypes in COVID‐19 suggest direct viral sensing

Abstract Growing evidence suggests that conventional dendritic cells (cDCs) undergo aberrant maturation in COVID‐19, which negatively affects T‐cell activation. The presence of effector T cells in patients with mild disease and dysfunctional T cells in severely ill patients suggests that adequate T‐cell responses limit disease severity. Understanding how cDCs cope with SARS‐CoV‐2 can help elucidate how protective immune responses are generated. Here, we report that cDC2 subtypes exhibit similar infection‐induced gene signatures, with the upregulation of interferon‐stimulated genes and interleukin (IL)‐6 signaling pathways. Furthermore, comparison of cDCs between patients with severe and mild disease showed severely ill patients to exhibit profound downregulation of genes encoding molecules involved in antigen presentation, such as MHCII, TAP, and costimulatory proteins, whereas we observed the opposite for proinflammatory molecules, such as complement and coagulation factors. Thus, as disease severity increases, cDC2s exhibit enhanced inflammatory properties and lose antigen presentation capacity. Moreover, DC3s showed upregulation of anti‐apoptotic genes and accumulated during infection. Direct exposure of cDC2s to the virus in vitro recapitulated the activation profile observed in vivo. Our findings suggest that SARS‐CoV‐2 interacts directly with cDC2s and implements an efficient immune escape mechanism that correlates with disease severity by downregulating crucial molecules required for T‐cell activation. This article is protected by copyright. All rights reserved


Introduction
Clinical outcomes of COVID-19 are highly variable. Patients may show either no or mild symptoms (such as mild fever and cough) or severe respiratory involvement requiring hospitalization. In the most severe cases, Acute Respiratory Distress Syndrome (ARDS) can develop, with high levels of inflammatory molecules in the blood [1], [2] and diffuse intravascular coagulation (DIC) [3], [4]. COVID-19 is lethal in a non-negligible number of cases [5]. Patients presenting severe symptoms show immune dysregulation, characterized by excessive release of type 1 and type 2 cytokines [2] and alterations of lymphoid and myeloid populations in the peripheral blood [6]. Severely ill patients, in contrast to patients with mild disease, also show alterations in both Th17 and Th1 cell activation, with defects in the acquisition of effector functions [7].
Cells of myeloid origin play a pivotal role during infections by sensing pathogens, producing inflammatory mediators, and contributing to the activation of adaptive immunity. In this context, dendritic cells (DCs) are particularly relevant, as they are specialized in antigen presentation and T-cell priming [8]. The differences observed in the activated T-cell compartments of patients with severe versus mild disease suggest abnormal activation in the conventional DC (cDC) compartment of patients presenting with more severe disease.
Myeloid cDC2s express various Pattern Recognition Receptors (PRRs) and can promote a wide range of immune responses, especially CD4 + T-cell responses [15]. Recently, cDC2s have been divided into two subsets, DC2s and DC3s [16]- [18]. DC3s have been described as a heterogeneous population that expands in inflammatory conditions [16].
Functional impairment of cDCs has been described in COVID-19 patients, with decreased numbers in the blood [19], [20] and reduced functionality, measured in terms of cytokine production and T-cell priming capacity [21] upon in vitro restimulation. Nevertheless, a defect in maturation upon in vitro restimulation does not necessarily indicate functional impairment, as activated DCs may not further respond to PRR agonists. No specific information is available concerning the impact of SARS-CoV-2 infection on the maturation of DC subtypes. A deeper understanding is crucial, given the important role of cDC subtypes in the activation and skewing of adaptive immune responses that ultimately contribute to COVID-19 pathogenesis [22], [23].
Here, we characterized the transcriptional signatures reflecting the functional state of cDC subpopulations by high-throughput single-cell RNA sequencing (scRNA-seq).

Results and Discussion
We analyzed peripheral blood cDCs from COVID-19 patients with severe and mild disease, according to the World Health Organization (WHO) classification, to better understand the impact of SARS-CoV-2 infection on cDC subtypes. Patients were enrolled from the STORM cohort (see Table S1 for clinical data of the patients) of San Gerardo Hospital in Monza, Italy.
Amongst CD11c + MHCII + peripheral blood mononuclear cells (PBMCs), cDC1s were identified as CLEC9A + and cDC2s were identified among the CD1c + FcεRIα + cells, excluding cells expressing markers for T or B lymphocytes (CD3 and CD19, respectively) or monocytes (CD88 and CD89) [16]. CD14 was included in the analysis to identify DC3s [17] (Supplementary Fig. 1 for gating strategies). Consistent with previous studies, we found a decreasing trend in the frequency of cDC1s and DC2s [19], [20] and an increasing trend in DC3s in the blood of COVID-19 patients relative to that of healthy donors (HDs) (Fig. 1A). The decreased number of cDCs in circulation may be due to both their recruitment to the respiratory tract and to apoptotic death caused by exposure to an inflammatory environment or direct viral encounter.
We systematically characterized the transcriptional response of cDCs to SARS-CoV-2 infection by analyzing three different single-cell transcriptomic datasets, two that were publicly available and one that was newly generated. Analysis of three independent datasets allowed us to identify consistently altered signaling pathways, minimizing the effects of possible biases in the single datasets.
The new dataset (dataset 1) [24] was generated using a droplet-based single-cell platform (10X Chromium) and contains scRNA-seq data of CD11c + MHCII + cells isolated from PBMCs of three COVID-19 patients (two with mild and one with severe disease) and two HDs ( Table S1). The second dataset [22] (dataset 2) contains cellular indexing of transcriptomes and epitopes by sequencing (CITE-seq) data of PBMCs and enriched DCs obtained from seven COVID-19 patients (three with mild and four with severe disease) and five HDs, and the third dataset [25] (dataset 3) contains scRNA-seq data of PBMCs obtained from 18 COVID-19 patients (8 with mild and 10 with severe disease) and 21 HDs.
Single-cell data from datasets 1 and 2 were first visualized using non-linear dimensionality reduction through uniform manifold approximation and projection (UMAP) and graph-based clustering algorithms (Supplementary Fig. 2A, 3A). Clusters containing myeloid DCs were identified based on the expression of markers that discriminate cDC2s and cDC1s from all other cell populations. Specifically, CD1C, FCER1A, and CLEC10A were used to identify cDC2s, whereas CLEC9A was used to identify cDC1s (Supplementary Fig. 2B, 3B). Myeloid DCs already This article is protected by copyright. All rights reserved. annotated by the authors were considered for dataset 3 [25]. Clusters corresponding to myeloid DCs in the three datasets were re-clustered in further iterations to separate cDC1s from cDC2s, to discriminate cDC2 subpopulations and to exclude possible contaminants or doublets.
Specifically, DC3s were distinguished from DC2s based on the expression of CD14, CD163, and S100A8. This approach allowed us to clearly identify cDC subsets (Fig. 1B, Supplementary Fig.   3C,D,E and Supplementary Fig. 4A).
Next, we aggregated cell-level counts into sample-level pseudo-bulk counts, mitigating singlecell mRNA measurement noise, to unravel the transcriptional response of each cDC subset during SARS-CoV-2 infection and identified differentially expressed genes (DEGs) between COVID-19 patients and HDs ( Table S2). The low numbers of cDC1s allowed their analysis solely in dataset 2.
Comparison of the expression profiles of COVID-19 patients with those of HDs showed most of the genes upregulated in COVID-19 to be interferon (IFN) stimulated genes (ISGs) in all cDC subsets ( Fig. 1C and Supplementary Fig. 5A,B). On the other hand, genes encoding MHCII molecules were among the most significantly downregulated genes in cDCs from COVID-19 patients (Fig. 1C), indicating that these cells have an impaired antigen presentation capacity.
We performed gene set enrichment analysis (GSEA) to better understand the biological signaling pathways that are differentially regulated in cDCs from COVID-19 patients relative to HDs using two gene sets: the Hallmark collection from the Molecular Signatures Database  Table S3). As expected, maturation was dominated by ISGs in all cDC subtypes from the three datasets, consistent with the identified DEGs, whereas we were unable to detect the upregulation of signatures containing classical activation markers and cytokines for T-cell priming ( Fig. 1D and Supplementary Fig. 6A). Along with the IFN-induced pathways, IL-6 pathways (IL-6-JAK-STAT3 and PI3K-AKT-mTOR [27]) were consistently upregulated in cDC2s in all datasets (Fig.   1D). This is in accordance with the relevance of IL-6 in COVID-19 pathogenesis and the expansion of activated Th17 cells in COVID-19 patients [28].
The absence of a conventional maturation signature (absence of upregulation of genes encoding MHCI-II, costimulatory molecules, and cytokines) in circulating DCs prompted us to determine whether it was, in fact, possible to identify activated DCs in the blood.
We, therefore, investigated the transcriptional responses of circulating DC2 and DC3 subsets at single-cell resolution under different clinical conditions. Two distinct publicly available datasets were analyzed: the dataset from Reyes et al. [29], containing scRNA-seq data of PBMCs and  [30], containing CITE-seq data of PBMCs obtained from healthy volunteers who received an adenovirus-based vaccine. As already described, we performed dimensionality reduction and unsupervised clustering to identify cDC subpopulations. Our approach clearly identified cDC subsets in both datasets ( Fig.   2A,B and Supplementary Fig. 7). We then determined DEGs in infected or vaccinated donors with respect to the corresponding HDs (Table S4) and performed GSEA.
The results are in stark contrast to those obtained from COVID-19 patients. Indeed, circulating   Fig. 10). Interestingly, we observed a stronger activation response for circulating cDC2s from patients with localized bacterial infections (Leuk-UTI group) and transient organ dysfunction (Int-URO group) than those with bacterial sepsis and persistent organ dysfunction (URO group) (Fig. 2E, left panel). This was expected, as sepsis induces functional impairment of myeloid cells.
These results show that conventionally activated DCs can be detected in the blood in response to both vaccination and infection (localized and systemic). Therefore, the absence of a classic activation signature in cDC subpopulations from the blood of COVID-19 patients is not a general phenomenon.
Recent studies have shown potential functional differences between DC2s and DC3s [36] in inflammatory diseases, such as Systemic Lupus Erythematosus (SLE), in which type I IFNs play a major role [16]. We investigated the potential specific role of DC3s with respect to DC2s by determining the genes that are differentially induced/downmodulated by these two subpopulations in response to SARS-CoV-2 stimulation and compared them to bacterial infection. The resolution was increased by pooling cDCs from the three COVID-19 datasets and performing Harmony integration [37], followed by graph-based clustering. After integration, we obtained 2,663 cDCs ( Fig. 3A) and clearly identified cDC1s, DC2s, and DC3s (Fig. 3B).
Only 30 genes (p-value < 0.05 and absolute log2 fold change > 1) were differentially expressed in DC3s relative to DC2s in response to COVID-19 infection, of which 19 were upregulated and 11 downregulated (Fig. 3C, left panel and Table S5). On the other hand, 146 genes (p-value < 0.05 and absolute log2 fold change > 1) were identified as being differentially regulated in DC3s relative to DC2s in response to intermediate urosepsis (Int-URO condition), of which 60 were upregulated and 86 downregulated (Fig. 3C, right panel and Table S5).
The diversity in the responses of DC3s and DC2s during bacterial infections can be explained, at least partially, by the differential expression of certain receptors, such as CD14, exclusively expressed by DC3s. CD14 is a component of the receptor complex of lipopolysaccharide (LPS), a major factor of the outer membrane of Gram-negative bacteria, and contributes to LPS recognition and internalization of the receptor complex [38]. Thus, DC3s can respond more efficiently to Gram-negative bacteria than DC2s. CD14 also has an important role as a chaperon for ligands of endosomal and cytosolic PRRs [38]. Therefore, the differences between DC2 and DC3 responses observed in Gram-negative bacterial infections may also occur after Grampositive bacterial recognition.
In conclusion, these findings suggest that DC2s and DC3s respond similarly to SARS-CoV-2 infection, whereas they show more diversified responses to bacterial infections.
Among the small number of genes differentially expressed between DC3s and DC2s in response to COVID-19 were those encoding complement factors (C1QA) and, most importantly, antiapoptotic genes, such as AXL and CLU, which were among the most significantly upregulated (Fig. 3C, left panel). This suggests that DC3s are less susceptible to apoptosis than DC2s and may explain why they tend to increase during SARS-CoV-2 infection (Fig. 1A). Comparison of these results with those obtained from bacterial infections showed genes associated with cell-cycle progression and cell proliferation (RGCC, SENP5, SMC6, SERTAD3, MAD2L1BP) to be specifically upregulated in DC3s (Fig. 3C, right panel). Therefore, DC3s may proliferate during inflammatory responses or circulating DC3s may contain proliferating progenitors that expand the DC3 population during bacterial infections.
Coherent with these observations and the frequency of cDCs found in the blood of COVID-19 patients (Fig. 1A), we also observed an alteration in the relative abundance of cDCs by scRNAseq. Specifically, DC3s showed higher frequencies in COVID-19 patients with respect to HDs (Fig. 3D).
Overall, these observations may explain why DC3s increase in number in acute and chronic inflammatory conditions [16], [39].
DC3s have been described as a heterogeneous population [16], [18], which was confirmed by our FACS analysis ( Fig. 3E and Supplementary Fig. 1). CD1c + CD5 -cells included CD14 low/+ CD163 low/+ cells, for which the abundance differed between HDs and COVID-19 patients. CD14 + CD163 + cells were present in COVID-19 patients but almost absent in HDs, who instead showed only CD14 low/+ CD163 low cells (Fig. 3E). Consistent with these findings, there was a significant increase in the frequency of CD14 + CD163 + DC3s, but not of CD14 low/+ CD163 low DC3s, in the blood of COVID-19 patients (Fig. 3E). We then investigated whether this heterogeneity was also measurable at the transcriptional level. Hence, we retained clusters annotated as DC3s ( Fig. 3B) and performed a re-clustering. We obtained a continuum of cells in which two subpopulations could be identified based on the expression of CD14 and CD163 (CD163 + CD14 + and CD163 low CD14 low cells) (Fig. 3F). Mirroring previous FACS analysis (Fig. 3E), we found an increase in the relative abundance of CD163 + CD14 + DC3s in COVID-19 patients relative to HDs (Fig. 3G).
Finally, we investigated cDC2 gene expression profiles in COVID-19 patients with severe versus mild disease to seek specific alterations in the innate immune signature between these two groups of patients and to link variations in the immune response to disease severity. As already described, we aggregated cell-level counts into sample-level pseudo-bulk counts and identified DEGs between COVID-19 patients with severe and mild disease (Table S6).
We identified a large number of DEGs in both DC2s and DC3s (101 for DC2s and 203 for DC3s, pvalue < 0.05 and absolute log2 fold change > 1) between patients with severe and mild disease, indicating relevant differences in the transcriptional response of these two groups (Fig. 4A).
Inflammatory genes not directly related to the activation of adaptive immunity, such as complement factors (C1QB) and complement receptors (C5AR1), those involved in the production of leukotrienes known to exacerbate respiratory syndromes (ALOX5AP) and the coagulation cascade (THBS1, THBD), those coding for factors involved in vasodilation (ADM), and other inflammatory genes, such as CD14, S100A8/A9/A12 and ADAM9, were significantly upregulated in patients with severe versus mild disease in DC2s and/or DC3s (Fig. 4A). In addition, genes that negatively interfere with the maturation of DCs for T-cell activation, such as TMEM176B [40] and MT1E [41], were upregulated in patients with severe versus those with mild disease (Fig. 4A). Moreover, CLU was among the significantly upregulated genes in DC3s of patients with severe disease, further supporting a potential role of anti-apoptotic genes in DC3s during infections (Fig.4A, right panel).
Strikingly, genes encoding MHCII molecules, the costimulatory molecule CD86, and cytokines, such as IL-1β, CCL3, and CCL4, showed progressive downregulation in DC2s and/or DC3s from HDs to patients with mild disease and, finally, those with severe disease (Fig. 4B).
These observations were consolidated by pathway analysis, which showed clear upregulation of pathways involved in metabolism, coagulation, angiogenesis, and reactive oxygen species in patients with severe versus mild disease (Fig. 4C). Among the leading edge genes of the allograft rejection pathway, which was found to be downregulated in patients with severe versus mild disease, were many critical for DC-mediated T-cell activation, such as those coding for proteins involved in antigen presentation in both the MHCI and MHCII pathways (TAP1, TAP2, HLA-DMB, HLA-DRA) and genes encoding molecules relevant for T-cell recruitment and activation (IL16, IL1B, CCL4, CCL5, CCR1, CCR2) (Fig. 4D). The specific downregulation of these genes in severely ill patients emphasizes the alteration of cDC functions in these individuals, which may be associated with worse disease progression.
Overall, these observations indicate that cDC2s progressively skew toward inflammatory activities and lose their antigen presenting function as disease severity increases. This could explain the alteration of the activated T-cell compartment observed in severely ill patients.
The results shown until now indicate that DC3s and DC2s respond similarly to the virus during COVID-19 infections, with three main features: i) the upregulation of ISG and IL-6 pathways, ii) progressive downregulation from patients with mild to severe disease of genes encoding signal 1 and signal 2 molecules associated with antigen presentation, and iii) the upregulation of an inflammatory signature, mainly represented by complement and coagulation factors, in severely ill patients.
It is possible that these features are due to exposure to mediators released during SARS-CoV-2 infection or to the direct interaction of cDCs with the virus. The first cDC characteristic observed in our study is compatible with both direct interaction with the virus and exposure to paracrine cytokines, such as IFNs and IL-6 produced by bystander cells. The lack of expression of IFN and IL6 genes in circulating DCs does not necessarily mean that cDCs cannot be a source of these cytokines, as their expression is acutely regulated and may be repressed in DCs early after activation. By contrast, the systematic downregulation of genes encoding MHCII molecules is more likely explained by a direct interaction of cDCs with the virus. This prediction is also supported by evidence that the virus can directly activate monocyte-derived DCs following abortive infection [42].
Therefore, we investigated whether the direct interaction of cDC2s with the virus could induce a similar response to that observed at single-cell resolution. We measured the response to the virus of cDC2s (CD1c + CD19 -cells) freshly isolated from HDs using IL-6 and MHCII as readouts.
As predicted, we found that SARS-CoV-2 directly induced significant downregulation of MHCII surface expression and the upregulation of IL-6 in both DC2s and DC3s (Fig. 5A,B). Conversely, the exposure of cDC2s to sera of patients with mild and severe disease, which contain inflammatory cytokines and other mediators, did not induce any modification in MHCII or IL-6 expression (Fig. 5C). This suggests that at least part of the peculiar response of cDCs induced in vivo by SARS-CoV-2 infection may be directly mediated by the virus.
cDCs are likely to detect the virus through C-type lectins, receptors that recognize Nglycosylation on the Spike protein and identified as being responsible for myeloid-cell activation [43]. C-type lectins can bind to many different viruses. However, in addition to playing a positive role in antiviral responses by inducing the activation of DCs, they also contribute to viral spread through the migratory capacity of DCs [44], [45]. It would be informative to determine whether and where the activated DCs present in the circulation encountered the virus. In the case of cytomegalovirus infections in the mouse, it has been shown that infected DCs reach the lymph nodes and then exit them through high endothelial venules to reach the circulation [46]. Therefore, it is possible that DCs encounter the virus in the lungs and then reach the circulation, with migration favored by the strong inflammatory condition. Alternative possibilities are also plausible, such as the encounter of the virus by immature DCs directly in the blood or even in the bone marrow in severely ill patients.
In conclusion, SARS-CoV-2, like other viruses, can be directly detected by DCs and induces the downregulation of signals necessary for T-cell activation, a phenomenon that is accentuated with disease severity. This allows the virus to evade control of the adaptive immune system, while the host attempts to counteract the viral infection through innate immunity.

Understanding how DCs manage SARS-CoV-2 infection will help in identifying ad hoc
This article is protected by copyright. All rights reserved.
interventions to achieve optimal adaptive responses, a prerequisite for a good prognosis [23], [47]. Committee approval was received for the study and the informed consent of all participating subjects was obtained.

Flow cytometric analysis
PBMCs from COVID-19 patients were extracted from peripheral blood by density gradient

Single-cell RNA sequencing datasets analyzed in the study
Three different single-cell datasets from COVID-19 patients and healthy controls were analyzed in this study.
Dataset 1 was newly generated. Myeloid cells were sorted (CD11c + MHCII + ) from three COVID-19 patients (two with mild and one with severe disease, enrolled from the STORM cohort) and 2 HDs using a MACSQuant Tyto instrument (Miltenyi) (Supplementary Fig. 1B). After sorting, Dataset 3 [25] was generated from a scRNA-seq experiment with PBMCs from 18 COVID-19 patients (8 with mild and 10 with severe disease) that was integrated with publicly available 10x scRNA-seq data of healthy controls from [29]. Seurat objects were downloaded from

Single-cell data processing and analysis
Data processing and analysis for all single-cell datasets was performed using the Seurat package First, filters were applied to remove low-quality cells. These were based on the number of genes and UMIs detected in each cell and the percentage of reads mapping to mitochondrial genes (cells with < 500 genes and > 10% of reads mapping to mitochondrial RNA were removed). For COVID-19 dataset 1, 16,325 genes and 15,400 cells were available before quality control (QC).
This article is protected by copyright. All rights reserved.
After QC filtering, we retained 16,325 genes and 13,364 cells. For COVID-19 dataset 2, 33,538 genes and 45,547 cells were available before QC. After QC filtering, we retained 33,538 genes and 33,430 cells. COVID-19 dataset 3 was pre-processed by the authors and 46,584 genes and 99,049 cells were available.
Counts were then normalized and log-transformed using sctransform [49], while regressing out UMI counts and the percentage of mitochondrial counts.
PCA was performed to reduce dimensionality. Principal components (PCs) were fed to Harmony [37] for batch correction and/or the integration of datasets from both disease and healthy conditions. UMAP was used for 2D visualization. Clusters were identified using the shared nearest neighbor (SNN) modularity optimization-based clustering algorithm, followed by Louvain community detection. Cell type assignment was manually performed using marker genes, as detailed in the figures. cDCs were retained and again re-clustered to identify subsets.

Doublet identification
We observed a number of clusters likely comprised of primarily doublets, based on mixed lineage markers, when re-clustering the cDCs in dataset 2. This was confirmed by two distinct doublet-calling algorithms: the computeDoubletDensity method from scDblFinder package [50] and the cxds methods from the scds package [51].

Pseudo-bulk differential gene expression analysis
After the identification of cDC subsets, we aggregated cell-level counts into sample-level pseudo-bulk counts. For each DC subset, only donors with at least 10 cells were retained. For the dataset from Reyes et al., only samples from the primary cohort were considered for differential analysis. Differential expression analysis was performed using the quasi-likelihood framework of the edgeR package [52], using each donor as the unit of independent replication. This article is protected by copyright. All rights reserved.

Integration between COVID-19 datasets
cDCs identified in the three COVID-19 datasets were pooled, integrated using the Harmony algorithm [37], and further subclustered using the shared nearest neighbor (SNN) modularity optimization-based clustering algorithm, followed by Louvain community detection to identify cDC1s, DC2s, and DC3s. For dataset 3, only COVID-19 samples were retained for the integration.

Data and Code Availability
For 54. Sergushichev AA. An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation. bioRxiv. 2016.  Table S2 for the full list of DEGs This article is protected by copyright. All rights reserved. This article is protected by copyright. All rights reserved. This article is protected by copyright. All rights reserved. (B, lower panel) Quantitative analysis of the percentage of IL-6-producing cells. Statistical significance was determined using the unpaired student t-test. **p < 0.01; n = 4 NT donors and n