The study protocol was approved by the Human Biomedical Research Ethic Committee of Shandong Provincial Hospital. Written informed consent was obtained from the patients. Inclusion criteria for septic patients were based on the Third International Consensus Definitions of Sepsis and Septic Shock (Sepsis 3.0) and required standardized anti-infective therapy (21). Patients under 18 years old, pregnant or breastfeeding, those with chronic liver or kidney diseases, or those with incomplete medical records were excluded from the study.

For bulk RNA-Seq, peripheral blood samples were collected from 5 healthy controls and 20 septic patients who were admitted to the ICU of Shandong Provincial Hospital in China. The baseline characteristics of patients in the study are shown in Table S1 . A volume of 2.5 ml of blood was collected from each participant into PAXgene Blood RNA tubes (BD Biosciences) to ensure the stability of intracellular RNA. Total RNA was isolated from PAXgene Blood RNA tubes using the PaxGene Blood miRNA kit (Qiagen). RNA sequencing libraries were prepared using the NEBNext UltraTM RNA library prep kit. The libraries were then analyzed on an Illumina NovaSeq 6000 platform.

For flow cytometry analysis, mononuclear cells were isolated from peripheral blood using density centrifugation. Each replicate included a total of 10 healthy controls and 20 septic patients. The experiment was performed in triplicate, resulting in a total of three replicates.

TrimGalore was used to trim raw reads, and the trimmed reads were mapped to the hg19 genome using HISAT2, generating sam files that were then converted to bam files by SAMtools. HTSeq was used to calculate the read count of each gene, and R package edgeR was employed to identify differentially expressed genes (DEGs) using a cutoff of adjust.p < 0.05 and |log2FC| > 1 (22).

The scRNA-seq data were downloaded from the Broad Institute Single Cell Portal and analyzed using the R package Seurat (23). Cell quality control was applied based on four metrics: total unique molecular identifiers (UMI) counts, the number of genes detected, the expression ratio of hemoglobin genes and the expression ratio of mitochondrial genes. Cells were filtered out if they met any of the following criteria: (1) more than 25000 UMI counts; (2) more than 4000 detected genes or less than 500 detected genes; (3) more than 20% of mitochondrial genes; (4) more than 1% of hemoglobin genes. R package DoubletFinder was used to remove doublets on a per-sample basis, and then all ribosomal and mitochondrial genes were removed to avoid unexpected noise (24).

To correct for batch effects in the dataset, the Harmony algorithm was used, and 3,000 variable genes were identified using the ‘SCTransform’ function (25). Principal component analysis (PCA) was performed using the ‘RunPCA’ function, and the PCA matrix was fed into the ‘RunHarmony’ function. Clustering and dimensionality reduction were then performed using the ‘FindClusters’ and ‘RunUMAP’ functions, respectively, based on the batch-corrected matrix.

Following the first round of clustering with a resolution of 0.8, 13 major cell types were identified, including B cells, plasma cells, NK cells, CD4+ T cells, CD8+ T cells, CD14+ monocytes, CD16+ monocytes, megakaryocytes, AXL+ SIGLEC6+ dendritic cells (AS-DCs), pDCs, cDC1s, cDC2s and cycling cells. Clusters expressing marker genes of more than one cell type were excluded and a total of 82273 cells expressing 22661 genes were retained for the downstream analysis. In the analysis for DCs and monocytes, each DC type and monocytes was extracted and underwent a second-round clustering with the same procedure as the first round to improve the annotation results.

Differential gene expression analysis was conducted using the ‘FindAllMarkers’ or ‘FindMarkers’ functions in R based on the ‘RNA’ assay, with a Wilcoxon signed-rank test used as the statistical method. A cutoff of adjusted P < 0.05 and |log2FC| > 0.25 was applied identify genes that were significantly differentially expressed.

For module analysis of scRNA-seq data, we utilized the R package high dimensional WGCNA (hdWGCNA), which is an extension of WGCNA (26). Specifically, we focused on cDC2 and monocytes and applied the analysis pipeline using the R packages Harmony and Seurat. The analysis pipeline included pooling cells within the same group and cell type to create metacells using the ‘MetacellsByGroups’ function, identifying an appropriate soft threshold using the ‘TestSoftPowers’ function, and selecting the module with the highest correlation to sepsis for further evaluation of its function using GO enrichment analysis. We also calculated the preservation of cDC2-related modules in pDC and cDC1 using the ‘ModulePreservation’ function and projected the cDC2-related modules and monocytes-related modules into other APC types using the ‘ProjectModules’ function.

GO enrichment analysis of DEGs was performed using R package clusterProfiler (27).

To analyze transcription factor (TF) activity, we utilized the python single-cell regulatory network inference and clustering (PYSCENIC) tool on all single cells (28). Differentially expressed TFs were calculated using R package Limma (29). Only significant differentially expressed TFs with   adjusted p-value < 0.05 were involved in downstream analysis.

For prediction of differentially expressed TFs in bulk RNA-seq data between healthy controls and septic patients, we used the R package DoRothEA (30).

For bulk RNA-seq data, gene set score was calculated using the R package gene set variation analysis (GSVA), and in the case of scRNA-seq data, we calculated the gene set score using the ‘AddModuleScore’ function (31). Score changes between two groups were evaluated using the R package ggpubr based on the Wilcoxon signed-rank test or Student’s t test.

To compare the three newly identified anergic DC subtypes with the previously reported DC3, we conducted pairwise comparisons using gene set enrichment analysis (GSEA). The GSEA signature list used in the analysis was obtained from the Supplementary Tables of Villani et al. (18). ‘DC3 > DC2’ gene set contained genes highly expressed in DC3 (characterized by low MHC II expression), while ‘DC2 > DC3’ gene set contained genes highly expressed in DC2 (similar to traditional cDC2).

To analyze heterogeneity in septic patients, we merged cells by calculating the sum of monocyte gene expression for each sample. Then, we corrected the gene expression of each sample using the ‘SCTransform’ function and selected the top 3000 genes for downstream analysis. For dimensionality reduction, we performed PCA and utilized the Louvain clustering method to cluster samples. We also applied reversed graph embedding to uncover the potential pseudostages that could link the clusters of patients (32).

Pseudotime analysis was performed on cDC1, cDC2 and pDC subtypes in sepsis using the R package Monocle3 (32–34). Top highly differential genes were plotted along the inferred developmental trajectories.

Intercellular communication analysis was conducted using the R package CellChat (35). The ‘CellChatDB.human’ database was used for analysis. The control and sepsis groups were analyzed separately and then merged for downstream analysis. We did not account for the effect of cell proportion, as the scRNA-seq data were enriched for DCs.

Bulk-data deconvolution was performed using CIBERSORTX (36). The reference signature matrix was constructed based on the count matrix in ‘RNA’ assay of scRNA-seq data.

A weighted co-expression network was constructed using the R package Weighted Correlation Network Analysis (WGCNA), followed by a selection of genes with the top 75% variance (37). The appropriate soft power β was selected using the ‘pickSoftThreshold’ function. The expression matrix is then converted to adjacency matrix to identify modules based on the topological overlap. A hierarchical clustering dendrogram was further built and the genes were clustered into different modules. The modules with the highest correlation to sepsis or control were selected for functional evaluation through gene ontology (GO) enrichment analysis.

All antibodies used for flow cytometry were mouse anti-human mAbs. Anti-human CD45 (Cat No: 340910), CD56 (Cat No: 345811), CD123 (Cat No: 564195), HLA-DR (Cat No:756414) were purchased from BD Biosciences. Anti-human CD3 (Cat No: 317321), CD19 (Cat No: 302239), CD88 (Cat No: 344304), CD89 (Cat No: 354120), CD14 (Cat No: 325618), CD45RA (Cat No: 304142), CD1c (Cat No: 331520), CD141 (Cat No: 344112) were purchased from BioLegend.

Samples were run on LSRFortessa (BD Biosciences) flow cytometer. Initially, all DCs, including pDCs, cDC1s and cDC2s, were gated based on CD45⁺CD3-CD56^-CD19^-CD88^-CD89^-. Subsequently, pDCs were specifically identified by gating on CD123⁺CD45RA⁺. Following the exclusion of pDCs, cDC1s were identified based on the expression of CD141, while cDC2s were identified based on the expression of CD1c. Data were analyzed with FlowJo v.10.8.1.

All analyses were performed using R software v4.2.1. The Student's t-test is used to assess the statistical differences between continuous variables that follow a normal distribution, while the Wilcoxon signed-rank test is used to assess the statistical differences between continuous variables that do not follow a normal distribution. The p values for differential gene expression analysis, GSEA, and GO analyses were adjusted using the Benjamini & Hochberg method. In Table S1 , continuous variables were presented as mean ± standard deviation or median with interquartile range depending on their normality. Categorical variables were presented as counts (percentages). The level of significance was set at 0.05.

We performed extensive literature mining at scientific search engines such as PubMed, Google Scholar, BASE, CORE to search related human sepsis scRNA-seq studies. Ultimately, we selected a study published in Nature Medicine with the largest sample size and DC enrichment for further analysis (9). The dataset was divided into 5 groups, including healthy control (n=15), leukocytosis (n=10, patients with blood WBC ≥ 12000 per mm³, but no organ dysfunction), Int-Sepsis (n=7, septic patients with mild or transient organ dysfunction), NoICU-Sepsis (n=14, septic patients in hospital wards or in emergency department), ICU-Sepsis (n=8, septic patients admitted to ICU). It is worth noting that all the groups except for leukocytosis and healthy control fulfilled the diagnostic criteria for sepsis ( Figure 1A ). After quality control and removal of doublets, we identified 13 cell types based on canonical annotation marker genes for immune cells used by many studies ( Figures 1B, C ). Our results showed that the proportion of monocytes increased gradually from the healthy control to the Int-Sepsis group, and then decreased as the disease progressed, while lymphocytes showed an opposite trend. Moreover, the proportion of DCs consistently decreased with disease progression ( Figure 1D ).

We next identified three cell subtypes, cDC2 A, cDC2 B and cDC2 C changed dynamically in sepsis disease progression in cDC2 population ( Figure 2A ). The populations of cDC2 A and cDC2 C gradually decreased, whereas the population of cDC2 B increased, particularly in the ICU-Sepsis group ( Figure 2B ). Additional analysis showed that cDC2 A and cDC2 C exhibited high gene expression levels of CD74 and HLA-DRA, which are associated with antigen presentation, the primary function of DCs. Conversely, cDC2 B was characterized by the expression of monocyte-associated genes such as S100A8, CD14, VCAN, and FCGR3A ( Figure 2C ). The marker genes of cDC2 B exhibited similarities to the previously reported DC3, whereas cDC2 A and cDC2 C were found to be more closely associated with DC2 (18). GSEA provided additional support for the idea that cDC2 B aligns with DC3 ( Figure 2D ) (18). Given that cDC2 B was previously identified as an inflammatory DC subtype (7), we calculated the inflammatory score of the three cDC2 subtypes and performed pathway enrichment analysis on the marker genes of cDC2 B. Interestingly, cDC2 B exhibited higher inflammatory score and enriched pathways related to the inflammatory response, such as ‘positive regulation of cytokine production’ and ‘positive regulation of inflammatory response’ ( Figures 2E and S1A ). Despite cDC2 B had a higher inflammatory score, we found that its highly expressed genes were monocyte-associated genes and it lost its primary function as a dedicated APC, specifically activating T cells via the MHC II pathway. Consistent with our findings, previous studies have demonstrated that cDC2 B has a reduced capacity to activate T cells in comparison to cDC2 A (7, 38, 39). Hence, we designated cDC2 B as anergic DCs, as previously described (15).

Next, we sub-clustered cDC1 population into two subtypes, cDC1 A and cDC1 B ( Figure 2F ). Surprisingly, cDC1 B, which shares similar marker genes with cDC2 B, was predominantly observed in the ICU-Sepsis group, suggesting the possible existence of an anergic subtype within the cDC1 population ( Figures 2G , S1B ). GSEA confirmed the similarity between cDC1 B and cDC2 B ( Figure 2H ). Moreover, analysis of inflammatory scores and pathway enrichment revealed that cDC1 B exhibited heightened inflammatory characteristics and enrichment in inflammatory pathways ( Figures S1C, D ).

Subsequently, we focused on pDC population to investigate whether the anergic subtype exists in all three traditional DC populations. We identified two clusters, named pDC A and pDC B ( Figure 2I ). Notably, pDC B, characterized by high expression of S100A8, S100A9, FCN1 and CD14, was the dominant subtype in the ICU-Sepsis group ( Figure 2J ). Further analysis using GSEA revealed a significant enrichment of the ‘DC3>DC2’ gene set in pDC B compared to pDC A ( Figure 2K ). Moreover, pDC B showed a higher inflammatory score and enrichment in inflammatory pathways ( Figures S1F, G ). These findings suggest that pDC B and cDC1 B represent two novel DC subtypes with gene expression profiling similar to the previously reported DC3, which we have designated as cDC2 B.

To further demonstrate the transcriptomic similarity of the three anergic DC subtypes, we performed hdWGCNA to identify co-expression gene modules in cDC2 and analyzed the preservation of these modules across other two DC populations (26). We identified a total of 15 modules in cDC2, of which five modules were found to be specific to cDC2 B, including the salmon, turquoise, magenta, red and greenyellow modules ( Figures S2A–C ). Of these five modules, only the red module was highly preserved in both pDC and cDC1 ( Figures S2D, E ). GO enrichment analysis revealed that this module was associated with inflammatory response and myeloid cell differentiation and activation ( Figure S2F ). We next projected the cDC2 gene modules onto pDC and cDC1 and found that the red module was also highly expressed in cDC1 B and pDC B, as well as in sepsis ( Figures S2G–J ). Furthermore, we performed trajectory and pseudotime analysis on all DC populations and found significant differences in gene expression between normal and anergic DC subtypes. The normal and anergic DCs were located at opposite ends of the pseudotime trajectory, with the expression of monocyte-related genes such as CD14, S100A8, S100A9, VCAN, and FCN1 increasing significantly along the pseudotime trajectory, while the expression of antigen-presentation-related genes such as CD74 and MHC II-related genes decreased significantly ( Figures S3A–F ). These findings further illustrate the consistency of the transcriptome among the three anergic DC subtypes and their close association with sepsis.

After demonstrating the transcriptome-level consistency of the three anergic DC subtypes, we then analyzed their TFs and pathway activity changes. Our analysis using PYSCENIC revealed that these anergic DCs had a high degree of similarity in the top upregulated TFs ( Figures 3A, B ) (28). Given the important role of the C/EBP family in inflammation-induced myelopoiesis and the high expression of monocyte-associated genes by anergic DCs, we evaluated the activity of C/EBP TFs (40). Notably, we observed a significant increase in the activity of almost all C/EBP family members in anergic DC subtypes. Furthermore, we discovered that cDC1 B and pDC B also exhibited low IRF8 activity, similar to the IRF8^low developmental pathway observed in cDC2 B (DC3), whereas cDC2 A/C (DC2) differentiated follows an IRF8^hi trajectory (41) ( Figure 3C ). These findings suggest that the three DC subtypes may share a similar developmental trajectory. We also examined antigen-presenting capacity changes between all DC subtypes and found that all three anergic subtypes had significantly lower ability to present external antigens via the MHC II pathway, indicating their anergic nature ( Figure 3D ). We then calculated the pathway activities based on 50 hallmark gene sets from MSigDB database to identify the pathway activity changes of all DC subtypes. Comparing the three anergic DCs with their corresponding normal subtypes, we observed similarities in pathway activity, including high activities in hallmark inflammatory pathways such as ‘INFLAMMATORY_RESPONSE’, ‘TNFA_SIGNALING_VIA_NFKB’, ‘TGF_BETA_SIGNALING’, ‘IL6_JAK_STAT3_SIGNALING’, ‘IL2_STAT5_SIGNALING’, and ‘COMPLEMENT’. Moreover, these anergic DC subtypes exhibited similarities in metabolic activity, including low oxidative phosphorylation, lipid metabolic activity and glycolytic activity ( Figure 3E ).

To investigate the distinct characteristics of anergic DCs in sepsis, we performed additional analysis to compare the phenotypic differences between septic patients and controls. The enrichment analysis showed that pathways related to antigen presentation were mainly enriched in the control and leukocytosis groups. In contrast, the sepsis group showed significant enrichment in inflammation-related signaling pathways, such as “positive regulation of cell activation,” “positive regulation of cytokine production,” and “NF-κB signaling” ( Figures S4A–F ). Similarly, a subsequent assessment of the antigen-presenting capacity of anergic DCs revealed that their ability to present antigens was enhanced in the leukocytosis group, but significantly reduced during sepsis ( Figure S4G ). Overall, our findings suggest that among the traditional DC types (cDC1, cDC2, pDC), there exist anergic subtypes with similar transcriptional regulation and enriched pathways. These subtypes may potentially contribute to the development of the observed immunosuppression in the context of sepsis.

Monocytes have been extensively implicated in the development of sepsis, and the heterogeneity between patients is considered a major contributing factor to the ineffective immune modulatory interventions in sepsis (42). To investigate whether monocyte subtypes are associated with heterogeneity in septic patients, we examined the expression profiling of monocytes from patients. We isolated all monocytes and identified six subtypes, including Mono1-CD14-RETN, Mono2-CD14-HLA, Mono3-CD14-CTNNB1, Mono4-CD14-CCL3, Mono5-CD16, and Mono6-CD16-C1QA ( Figures 4A , S5A ). Notably, the proportion of Mono1, which shares similar marker genes with previously reported immunosuppressive MS1( Table S2 , Figure S5B ) (9), was significantly elevated in septic patients ( Figure S5C ). To investigate the covariation in gene expression patterns between patients, we merged the single-cell expression matrix into a sample-level matrix and used unsupervised analysis techniques, including dimensionality reduction, clustering, and pseudostage analysis, to identify three distinct clusters (cluster 1, cluster 2, and cluster 3) of samples that varied along a pseudostage. A three-stage model associated with different immune states was then developed ( Figure 4B ). Samples from cluster 2 exhibited high expression of pro-inflammatory cytokines, including CCL2, CCL3, IL1B, IL8 and TNF, while samples from cluster 3 showed high expression of immunosuppressive cytokines such as IL10 and IL1RN ( Figure 4C ). These findings suggest that peripheral inflammation in septic patients is highly heterogeneous. We then examined the composition of monocyte subtypes in the three clusters and observed a significant enrichment of different cell subtypes. Specifically, Mono1 and Mono6 were enriched in cluster 3, while the proportion of Mono2 and Mono4 decreased progressively in clusters 2 and 3 ( Figures 4D, E ). These results indicate that the differential abundance of certain monocyte subtypes may contribute to the observed differences in inflammation levels among patients.

Given the significant expansion of Mono1 in septic patients and its enrichment in immunosuppressive cluster3, Mono1 should be the one of the main contributors of sepsis. To further investigate the co-varying genes among Mono1, we employed hdWGCNA and identified nine modules ( Figure S5D ). Notably, we observed that the blue module exhibited high expression specifically in Mono1 and was exclusive to sepsis ( Figure 4F ). Surprisingly, GO enrichment analysis showed that the blue module was associated with myeloid cell activation and inflammatory response, which is highly similar to the cDC2 B-specific red module ( Figure S5E ). We therefore projected the Mono1-related blue module into all DCs and found that it was specifically expressed in anergic DC subtypes ( Figure 4G ). Furthermore, the cDC2 B-specific red module was highly expressed in Mono1 ( Figures S5F–G ). We then performed GSEA on monocyte subtypes and found that the ‘DC3 > DC2’ gene set was also enriched in Mono1. This suggests that Mono1 share a similar gene expression pattern with anergic DCs ( Figure 4H ). PYSCENIC analysis showed that the top 10 highly expressed TFs were also similar to those in anergic DCs ( Figure 4I ). In addition, Mono1 and anergic DCs shared similar transcriptional characteristics, such as low expression of MHC II-related genes. The onset of sepsis leads to the expansion of anergic subtypes within all APCs, which may be one of the main reasons for the development of immunosuppression in sepsis.

Based on our observation of reduced expression of the MHC II genes in APCs and their specific role in interactions with T cells, we hypothesized that sepsis may lead to a generalized dysfunction in APC-T cell communication. We then utilized CellChat to analyze communication patterns among different cell subtypes separately for each condition (35). We observed that CD4 T cells and CD8 T cells exhibited a significant increase in the strength of interactions in septic patients. Specifically, CD4 T cells showed a greater increase in outgoing signals than CD8 T cells, while CD8 T cells had a higher increase in incoming signals. Interestingly, normal DC subtypes (cDC1 A, cDC2 A, cDC2 C and pDC A) showed no significant enhancement in outgoing signals or even a decrease, while anergic DC subtypes (cDC1 B, cDC2 B and pDC B) exhibited significantly enhanced interactions with T cells ( Figure 5A ). Moreover, DCs were found to be the primary receivers and senders of signals regardless of the condition ( Figure 5B ). We then focused on the interactions between all APC subtypes and T cells. The strength of the interaction between these anergic DCs and T cells was weak in healthy controls, but significantly enhanced in septic patients, even exceeding that of the other DCs ( Figure 5C ). Additionally, the interaction between Mono1 and T cells was also significantly enhanced in septic patients ( Figure 5D ). Although the anergic DCs are less likely to act via ligand-receptor pairs in the MHC I and MHC II pathway compared to other DC subtypes, they showed enhanced communication probability to use these pathways in sepsis ( Figure 5E ).

To confirm the findings from the single-cell dataset, bulk RNA-seq analysis was conducted on 20 septic patients and 5 healthy controls. After quality control based on PCA and heatmap, we performed differential expression and GO enrichment analysis. The upregulated genes were chiefly enriched in cytokine release and myeloid activation, while the downregulated genes were mostly enriched in T cell activation and MHC II pathway, which is in alignment with the changes in cell proportions demonstrated in the scRNA-seq dataset ( Figure 6A , Figures S6A–C ). Furthermore, correlation analysis revealed a significant positive correlation between the pathway enrichment scores of the scRNA-seq dataset and the bulk RNA-seq dataset ( Figure 6B ). To assess the potential diagnostic value of these anergic DC subtypes in sepsis, deconvolution analysis was carried out. The analysis indicated higher proportions of the three anergic DC subtypes in septic patients compared to healthy controls. Additionally, the proportions of these anergic DC subtypes inferred by CIBERSORTX for each subject individually can be used as a classifier of sepsis (area under the curve (AUC), cDC1 B = 0.90, cDC2 B = 0.85, pDC B = 0.65) ( Figure 6C ). The expression level of the anergic DC-related red co-expression module and Mono1-related blue module in the bulk RNA-seq data was also evaluated. Consistent with cell proportions, the expression of both modules was significantly increased in septic patients, and further augmented when the patients presented with shock ( Figures 6D, E ). We next conducted WGCNA on bulk RNA-seq data to identify modules specific to healthy controls and septic patients (37). The analysis identified 13 distinct modules, of which three (purple, magenta, and tan) were specific to sepsis, and two (blue and turquoise) were highly expressed in healthy individuals. The purple module was strongly associated with exocytosis, myeloid cell differentiation, and endoplasmic reticulum stress, and exhibited significantly higher expression levels in sepsis patients. Additionally, its expression was observed to further increase in the presence of shock. The magenta and tan modules were linked to cell division and humoral immunity, respectively. Their expression levels were enhanced in septic patients but did not appear to be further increased at the onset of shock, and in fact, their expression levels even seemed to be reduced. The expression of the blue module showed a clear decrease with increasing disease, although this module was not enriched in any specific pathway. Finally, the turquoise module, related to antigen presentation and lymphocyte differentiation, displayed a consistent decrease in expression with increasing disease, in line with the results obtained from the differential analysis ( Figures S6F, G ). Given that the purple module demonstrated the strongest correlation with sepsis and its expression was further enhanced with disease progression, the expression of hub genes within this module was evaluated in different APC subtypes. Notably, these hub genes were significantly overexpressed in all three anergic DCs, as well as in Mono1, further supporting the relevance of these cells to sepsis ( Figure S6H ).

Following the insights gained from the single-cell analysis, we proceeded to perform flow cytometry analysis on cohort comprising both healthy individuals and septic patients. The detailed gating strategy steps employed in the analysis are provided in the Methods section. Our findings unveiled a notable elevation in the relative proportion of CD14⁺HLA-DR^low DCs specifically within pDCs, cDC1s and cDC2s in septic patients, thereby providing additional validation to our initial observations ( Figures 6F, G ).