APC subset distribution in patients with COVID-19
To characterize the molecular profile of circulating APCs, we performed scRNAseq on freshly sampled APC-enriched PBMCs from five patients with moderate COVID-19 (non-mechanically ventilated, oxygen supply <10 l min−1) and ten patients with severe COVID-19 (mechanically ventilated or oxygen supply ≥10 l min−1), at day 1 and day 4 following hospital and/or intensive care unit (ICU) admission, as well as four elderly healthy controls (HC) (Fig. 1a and Supplementary Tables 1 and 2). To obtain single-cell suspensions and minimize DC–DC and DC–T cell clusters and clumps, EDTA-containing medium was used for the enrichment steps in the first set of samples, which we further define as the ‘discovery set’. This set is composed of a total of 12 samples from two HCs, three patients with moderate COVID-19 and four patients with severe COVID-19 from both day 1 and day 4 time points (results are presented in the figures and Extended Data Figs. 1 and 2). However, EDTA is known to decrease reverse transcription (RT) efficiency through RT deactivation and ion chelation, resulting in reduced amounts of complementary DNA (cDNA) during amplification. We therefore validated the main results derived from the discovery set by using RPMI (EDTA-free) medium for the enrichment steps in a second set of samples, including a total of 15 samples (defined as the ‘validation set’) from two HCs, two patients with moderate COVID-19 and six with severe COVID-19 (Extended Data Figs. 3–7).
For each fresh sample, 25,000 cells (~20,000 monocytes and 5,000 total DCs) were loaded onto the 10X lane (10X Genomics technology) (Fig. 1a). As expected, more cells per sample were effectively sequenced in the EDTA (discovery set) than the RPMI (validation set) dataset (mean: 3,360 versus 2,528 cells, respectively) (Extended Data Fig. 3a), confirming that EDTA optimizes single-cell suspension efficiency for rare DC types. This retrospectively justified the importance of using two complementary experimental protocols, split into two independent datasets, to avoid biasing the results. All main findings were validated in both datasets, indicating the reproducibility and robustness to experimental procedures. Altogether, we analysed a total of 81,643 APCs, split into 42,784 cells in the discovery set and 38,859 cells in the validation set. The two sets were analysed separately after sample integration using Harmony38. Graph-based clustering (SNN-based), community detection and nonlinear dimension reduction, using uniform manifold approximation and projection (UMAP), were independently applied to both sets for cell cluster visualization. Manual annotation of the cell clusters using canonical gene signature markers for each APC subset established a comprehensive map of APCs in HCs and patients with COVID-19 in both sets (Fig. 1b,c and Extended Data Fig. 3b,c). In our discovery set, among the 42,784 APCs, we recovered six subsets: 22,690 CD14+ monocytes; 866 CD16+ monocytes; 13,252 CD1c+ DCs; 1,754 CLEC9a+ DCs; 3,538 pDCs; 684 Axl+Siglec6+ AS-DCs (Extended Data Fig. 3a). The validation set included 29,409 CD14+ and 1,021 CD16+ monocytes, 5,754 CD1c+, 197 CLEC9a+ DCs, 1,602 pDCs and 876 AS-DCs (Extended Data Fig. 3a). In both sets, APC populations were captured across all the collected samples (Supplementary Table 3).
The accurate identification of all six APC populations was confirmed by the expression of canonical markers defining each subset (Extended Data Fig. 3d). All DC populations expressed higher levels of human leukocyte antigen HLA-DR and CD86 compared to monocytes (Extended Data Fig. 3d). None of the cells expressed CD19 (B-cell marker), GNLY (natural killed (NK) marker) or CD3E (T-cell marker), validating the pure APC populations. CD14+ monocytes expressed lineage-defining CD14, whereas CD16+ monocytes expressed FCGR3A. AXL expression distinguished AS-DCs from pDCs, whereas CD1c and CLEC9a characterized the respective cDC subsets39,40. In both sets (discovery and validation), UMAP embeddings coloured by severity revealed the heterogeneity of APC distribution between the three groups (Fig. 1c). This was confirmed by splitting the UMAP embeddings per severity (Fig. 1d). Overall, our enrichment strategy allowed the efficient identification of all APC populations including the rare pDCs, AS-DCs and CLEC9A+ DCs, enabling further molecular and phenotypic characterization.
Inflammation-related pathways are hallmarks of COVID-19 APCs
We performed differential expression and pathway enrichment analyses among APC severity groups, revealing 368 differentially expressed genes (DEGs) among the three groups (absolute fold change > 1.4). Among them, 101 genes were upregulated in HCs (as compared to patients with moderate and severe COVID-19), 109 in patients with moderate COVID-19 and 134 in patients with severe COVID-19 as compared to the two other groups, respectively (Fig. 2a). The top 50 DEGs upregulated in severe APCs as compared to HCs and patients with moderate COVID-19 included pro-inflammatory molecules (IL1B, CXCR4), surface markers (CD36, CD83, AREG, ITGAM), enzymes (CTSD, CTSB) and secreted molecules (RETN, EREG, ANXA2) (Fig. 2b). Next, we sought to identify enriched pathways discriminating each severity group from HCs. We found enriched IFN-γ and IFN-α response pathways in APCs from patients with moderate COVID-19, whereas hypoxia and TNF-α signalling were enriched in patients with severe COVID-19 (Fig. 2c).
We next compared the enriched pathways upregulated in severe versus moderate COVID-19 and in moderate versus severe, respectively. We found that IFN-γ and IFN-α pathways could be used to discriminate moderate from severe APCs at the global level (Fig. 2d).
To allow for an accurate comparison between the two transcriptional signatures, we ranked the DEGs of the pairwise comparison according to decreasing fold change. Severe APCs significantly upregulated AREG (amphiregulin), IL1R2 (IL-1 receptor), NRGN (calmodulin binding protein) and pro-inflammatory molecules (S100A12) (Fig. 2e). However, moderate APCs overexpressed interferon-stimulated genes (ISGs; IFITM2, ISG15 and IFI27) and HLAII molecules (HLA-DRB5 and HLA-DQA2), suggesting decreased Ag presentation and antiviral programs in severe as compared to moderate APCs (Fig. 2e). Similar observations were recovered from our validation set (Extended Data Fig. 4a–d). Additional upregulated genes in severe as compared to moderate APCs were found in the validation set, including CXCL8, NAMPT and G0S2 (Extended Data Fig. 4e).
Defective IFN responses in COVID-19 APCs
Increases in inflammatory cytokines have been reported in COVID-19. We addressed the global contribution of APCs to the expression of inflammatory cytokines and their receptors. As compared to APCs derived from HCs, IL1B, CXCL2, CXCL8 and CCL3 were significantly increased, whereas IL18 was decreased in both severity groups (Fig. 3a and Extended Data Fig. 1a). TGFB1 and IL10RA expression decreased in severe, but not in moderate subsets, as compared to HCs (Fig. 3a and Extended Data Fig. 1a), whereas IL6 was not detected in our discovery set (Extended Data Fig. 1a). Despite the low expression levels of most cytokines, we explored downstream biological pathways associated with inflammatory cytokine signalling (mainly IL1B, IL6 and TNF-α). In comparison to APCs from HCs, both moderate and severe APCs showed higher score levels for hallmark inflammatory pathways, including ‘IL6_JAK_STAT3’, ‘TGF-β’, ‘P53’, ‘TNFa_SIGNALLING_VIA_NFKB’ and ‘KRAS_SIGNALLING’ (Fig. 3b). The IFN family of cytokines is one of the most important for innate and adaptive antiviral responses. We showed the expression levels of IFNL1, IFNL1R, IFNAR1, IFNAR2, IFNA1, IFNGR1 and IFNGR2 and further explored their distribution in the three severity groups via a scaled heatmap (Fig. 3c). Both IFN receptor types (IFNAR1, IFNAR2, IFNGR1 and IFNGR2) were broadly expressed in the APC subsets, whereas detection of IFNL1 and IFNLR1 was patchy in our discovery dataset (Extended Data Fig. 3b). The heatmap representation indicated that severe APCs expressed lower levels of IFN molecules, suggesting a potential defect in IFN signalling (Fig. 3c). To further validate this hypothesis, we investigated the expression levels of ISGs. We observed higher expression levels of ISGs (MX2, ISG15, IRF7, BST2, IFITM2 and ADAR) in moderate APCs, but lower levels in severe APCs, supporting the hypothesis of defective antiviral programs contributing to the severity of COVID-19 (Fig. 3c). We further stratified a more exhaustive ISG signature according to their respective functions related to ‘antiviral’ and ‘regulators of IFN signalling’. Moderate APCs displayed higher levels of these two ISG families compared to both severe and HC groups (Fig. 3d,e). These results suggest a global perturbation of IFN downstream functions in severe COVID-19 APCs.
Multi-process effector defects in severe COVID-19 pDCs
After having analysed COVID-19 APCs at the global level, we sought to decipher alterations occurring in specific APC subsets. To depict the alterations occurring in pDC subsets, we isolated and sub-clustered pDCs (Fig. 4a) and performed pairwise differential expression among the three severity groups. Pathway enrichment analysis using MsigDB hallmark signatures was conducted on the upregulated genes in each subset. Compared to both moderate and HC pDCs, severe pDCs were enriched for the ‘TNFa_SIGNALING’, ‘IL2_STAT5’ and ‘HYPOXIA’ signalling pathways. In parallel, compared to pDCs from patients with moderate COVID-19, pDCs from patients with severe COVID-19 were enriched in the ‘IL6_JAK_STAT3’, ‘P53’ and ‘MTORC’ signalling pathways (Fig. 4b). When comparing pDCs between patients with moderate and severe COVID-19, the most notably enriched pathways were related to IFN signalling (IFNG and IFNA response), along with MYC targets signalling pathways (Fig. 4b). We asked whether apoptosis and pro-inflammatory signalling signatures would be associated with changes in pDC innate sensing receptors, including TLR9, DHX36, IFNAR1 and IFNAR2. We imputed the expression values to recover the signal from dropped-out features using Nebulosa (https://github.com/powellgenomicslab/Nebulosa), and plotted the density estimation values on UMAP embeddings (Fig. 4c). We observed zero-value density levels for TLR9, along with decreased density levels for DHX36, IFNAR1 and IFNAR2, in pDCs from patients with severe COVID-19 (Fig. 4c). To explore whether these modulations may impact pDC functions, we defined four original functional modules using a literature-driven manual curation: ‘immune cell attraction’ (hereafter ‘attraction’) (18 genes), ‘innate sensing’ (12 genes), ‘antiviral effector molecules’ (23 genes) and ‘cytotoxicity’ (12 genes) (Fig. 4d). Each of these modules was crossed with the pDC expression matrix, and detected genes were depicted for each patient group (Fig. 4d). No major differences between groups were detected within the ‘attraction’ module. On the contrary, many genes in the ‘innate sensing’, ‘antiviral effector molecules’ and ‘cytotoxicity’ modules were detected in the three groups, and followed the same pattern: baseline in HCs, increased in patients with moderate COVID-19 and decreased in patients with severe COVID-19 (Fig. 4d,e and Extended Data Fig. 5). This was particularly striking for the viral sensors TLR7, DHX9 and DHX36, the cytotoxic molecule TNFSF10 and the antiviral effector IRF7. These results were supported by the downregulation of antiviral ISGs and innate sensors in pDCs from patients with severe COVID-19, including BST2 and PYCARD (Fig. 4e), in both experimental datasets (Extended Data Fig. 5).
Coordinated transcriptional adaptation in monocyte subsets
Monocytes have been implicated in the physiopathology of severe sepsis and COVID-19. We performed dimensionality reduction through independent component analysis (ICA) and highlighted cells according to their severity group. We observed that IC1 clearly separated moderate from severe and HC CD14+ monocytes, whereas IC2 distinguished HC from COVID-19 CD14+ monocytes (Fig. 5a). The top 50 genes contributing to either IC1 or IC2 revealed distinct transcriptional signatures for the CD14+ monocyte subsets identified in each severity group: the severe subset expressed higher levels of complement (C1GC and C1GB), B7 family (VSIG4) and CD163, which may function as an innate immune sensor and inducer of local inflammation. The moderate monocyte subset expressed increased levels of antiviral ISGs (IFITM1, IFITM3, IFI27, MZB1 and IFI6) and the HLA-II gene (HLA-DRB5), suggesting an efficient antiviral program (Fig. 5b). Compared to HCs, several transcription factors (TFs) were downregulated in both moderate and severe groups, including the AP-1 superfamily (FOS, JUNB and ZFP36) and DUSP1, involved in MAPK dephosphorylation (Fig. 5b). Pathway enrichment analysis on the top 50 genes contributing to IC1 and IC2 identified key pathways that segregated COVID-19 CD14+ monocytes from HCs (Fig. 5c). The ‘complement’, ‘TNF-α’, ‘KRAS’ and ‘hypoxia’ signalling pathways were upregulated in COVID-19 monocytes, whereas ‘IFN-α’ and ‘IFN-γ’ response signalling were decreased in the severe subset, as compared to the HC and moderate subsets (Fig. 5c). To estimate antiviral effector functions, we used our manually curated gene functional module across patient groups (Fig. 5d). We observed a decrease of almost all antiviral effector molecules in patients with severe COVID-19, as compared to either HCs or patients with moderate COVID-19, in both experimental datasets (Fig. 5d and Extended Data Fig. 6). In parallel, we subclustered CD16+ monocytes and reduced the data dimension using UMAP projection to depict the corresponding clusters for each severity group (Fig. 5e). Differential expression between the three severity groups of this subset indicated similar trends as described in CD14+ monocytes (Fig. 5b,f). This included overexpression of ‘complement’-related genes (C1QA, C1QB and C1GC) by the severe subset, upregulation of antiviral ISGs (ISG15, IFI6 and IFI44L) in the moderate subset, as compared to the HC subset (Fig. 5f). Overall, these disease-associated changes in CD16+ paralleled those observed in CD14+ monocytes, suggesting common adaptation mechanisms.
CLEC9A+ DC- and AS-DC-specific transcriptional alterations
Thanks to our APC enrichment protocol, we could recover rare CLEC9a+ DC and AS-DC subsets. Differential expression of AS-DC severity groups revealed significant upregulated genes in severe AS-DCs (SEPT7 and AREG), compared to the moderate and HC subsets. We could also observe a significant downregulation of the HLA-DQA2 gene and antiviral IFI27 gene in severe, compared to moderate AS-DCs (Fig. 6a). In the search for upstream regulatory mechanisms, we inferred TF activity using the Dorothea algorithm41 and scored the activity of each regulon using the Viper inference tool42. This identified a large number of highly variant TF activity scores (Fig. 6b). In moderate AS-DCs, we observed a higher activity scored for IRF1, IRF9 and STAT2, reported to be involved in the ISG transcription cycle (Fig. 6b). In AS-DCs from patients with severe COVID-19, we found increased TF activities for RELA, NFKB1, STAT5 and STAT3, indicative of a higher activation of NFKB/STAT signalling, potentially induced by the pro-inflammatory cytokines described in the ‘APC subset distribution in patients with COVID-19’ section, along with hypoxia activation, indicated by a higher activity of HIF1A (Fig. 6b).
DEGs among the CLEC9a+ DC subclusters included specific transcriptional signatures segregating patients with moderate and severe COVID-19 from HCs (Fig. 6c). We remarkably observed a downregulation of HLA-II genes, including HLA-DQB1 and HLA-DPB1, in severe as compared to HCs, along with a significant upregulation of a larger subset of ISGs, including IRF1, IFI44L, IFI6, IFI27, IFITM2, IFITM3, IFI44L, ISG15 and ISG20, in moderate as compared to both HC and severe subsets (Fig. 6c). Expression values representation indicated a significant increase of AREG and SEPT7 genes, which were also upregulated by severe AS-DCs (Fig. 6a,d). Most importantly, we noted a significant decrease of the IFNGR1 CLEC9a+ DC subset in patients with moderate and severe COVID-19 as compared to HCs (Fig. 6d), supporting a defective antiviral program.
Downregulation of MHC-II and CIITA activity in CD1C+ DCs
We next focused on disease-induced alterations in CD1c+ DCs. We first explored the gene expression levels of MHC-II-related genes (Fig. 7a). We noted a global decrease of HLA-II genes (mainly HLA-DQA2 and HLA-DRB5) in patients with severe COVID-19 (Fig. 7a). Similar findings were reported for the validation set (Extended Data Fig. 7). We then grouped the expression values of these MHC-II genes (HLA-DRB1, HLA-DMA, HLA-DQA2, HLA-DRB5, HLA-DPB1, HLA-DQB1 and HLA-DMB), constructed a signature that we named the ‘HLAII’ module, and scored CD1c+ DCs using the ‘AddModuleScore’ Seurat function. The ‘HLAII’ signature was significantly reduced in severe as compared to HC CD1c+ DCs (Fig. 7b). This was associated with decreased expression of upstream MHC-II regulators (including RFX5, RFXANK and CIITA) in the severe group (Fig. 7b). Comparison of the scaled values revealed a reduction of IRF1 and RFX5 TF activities, mainly described to be involved in MHC-II gene synthesis, whereas C/EBP family member (CEBPB and CEBPD) TF activities, known to be involved in myeloid fate differentiation, were increased in severe subsets (Fig. 7c). We also noted higher TF activities for RELA (NFKB superfamily) and the AP-1 family in patients with severe COVID-19, including FOSL1, FOSL2 and JUN, which regulate a large range of cellular processes, including cell survival, death and proliferation (Fig. 7c).
To further decipher the transcriptional changes occurring in DCs when transitioning from healthy to moderate and severe conditions, we conducted pseudo-temporal inference using Monocle3, using the UMAP embedding two-dimensional space of the DC subsets (Fig. 7d). The pseudotime tree revealed a continuous trajectory from healthy to moderate, and a marked transition to the severe subsets. This trajectory was correlated to pseudotime values (Fig. 7d, right). To recover the genes contributing to this transition tree, we conducted a graph-based test to assess the most significant genes. The top genes were associated with Ag presentation, including B2M and HLA-DPA1, along with genes related to ISG expression (GABARAP and IFITM3; Fig. 7e).
Given that MHC-II genes are involved in the DC–T cell interaction, we hypothesized that a more global dysfunction of DC–T cell communication may occur in COVID-19 APCs. To test this hypothesis, we applied our cell communication inference computational framework ICELLNET43. Using our ‘reference partner cell’ methodology, we inferred potential communication between each of the APC subsets and CD4+ T cells in each of the disease groups (Fig. 8a). Cell connectivity networks revealed a global decrease in APC–T cell communication in patients with severe, as compared to moderate COVID-19 and HCs, predominantly in CD1c+ DCs, CLEC9a+ DCs and CD14+ monocytes (Fig. 8a,b). We then explored the various molecular families that may explain this decrease. This revealed a dominant contribution of immune checkpoint molecules and cytokines for CD1c+ DC–T cell communication (Fig. 8b), in particular decreased JAG-NOTCH, CD80-CD28 and CD48-CD2 interactions (Fig. 8c). Cytokines were mostly underlying the decrease in CD14+ monocytes–T cell communication (Fig. 8b). As expected, signalling through HLA-II-related genes (HLA-II/LAG3 pairs) was significantly decreased in moderate and severe subsets as compared to HCs. Among the cytokines, IL10-, CCL5– and TGFb-mediated interactions were predominantly damped in patients with severe COVID-19, which may contribute to immunopathology through excessive Th1 responses (Fig. 8c).
Persistent defects in severe COVID-19 APCs across samples
We also asked whether functional pathway alterations observed across APC subtypes were sustained over time. In parallel, we wanted to ensure that our main findings were not driven by a single patient and/or time point. We compared scRNAseq datasets generated at day 1 versus day 4 post hospital admission for each patient, in all severity groups. We used a focused approach, by selecting genes involved in previously identified altered functions, in APCs from patients with severe COVID-19, and systematically compared day 1 and day 4 expression levels. Most of the day 1 defects were sustained at day 4, in particular the low score of the HLA-II module in CD1c+ DCs (Extended Data Fig. 2a) together with the decreased expression of the antiviral effector molecules in CD14+ monocytes (Extended Data Fig. 2b) and increased ‘apoptosis p53 pathway’ in pDCs (Extended Data Fig. 2c). At a patient level, we observed that HC samples displayed similar score (or expression) levels for these three biological processes, whereas both moderate and severe samples displayed slight differences due to inter-individual heterogeneity. Overall, we confirmed that our findings were not associated with either a specific time point or a dominant single patient effect. However, this does not exclude changes in APC molecular profiles at later times in the course of moderate and severe COVID-19.