All datasets were obtained from the GEO dataset (https://www.ncbi.nlm.nih.gov/geo/). The scRNA-seq dataset GSE167363 included 2 healthy subsets and 4 sepsis subsets. The transcriptomic data were downloaded from GEO datasets GSE57065 and GSE95233. The training dataset GSE57065 included 25 healthy subsets and 28 sepsis subsets. The validation dataset GSE95233 included 22 healthy subsets and 51 sepsis subsets. GSE57065 and GSE95233 were also used for meta-analysis. Information of GSE167363, GSE57065, and GSE95233 were showed in Supplementary Table S1 .

The R “Seurat” package (version 4.3.0.1) was used to process the scRNA-seq data. The expression profile was read in by Seurat package and screened out (nFeature_RNA>200 & nFeature_RNA<5000 & percent.mt < 20). 1Unsupervised clustering was conducted using Principal Component Analysis (PCA) and Uniform Manifold Approximation and Projection (UMAP) analysis after data normalization with the LogNormalize method, allowing us to see cell populations on a two-dimensional map. The Human Primary Cell Atlas was used as the source of reference material for the SingleR package’s cell annotation function. To identify each cluster’s marker gene, we used the “FindAllMarkers” tool and threshold values of fold change (FC).

High dimensional data, such as single cell RNA-seq, can be used to perform weighted gene co-expression network analysis (WGCNA) using hdWGCNA. hdWGCNA detects strong modules of linked genes and gives the biological context for these modules. It may be used to build co-expression networks in a cell type-specific manner. In this study, we investigated the hub genes of neutrophils under sepsis using hdWGCNA. And lastly, a subpopulation of hepatocytes undergoing regeneration was screened for distinctive genes. The “hdWGCNA” package version utilized in this research is 0.1.1.9010.

The trajectory analysis was performed via the R “Monocle” package (version 2.18.0). The cell data set was constructed using an integrated gene expression matrix, exported from the Seurat object into Monocle. SetOrderingFilter was used to sequence cells based on the variable genes defined by dispersionTable. Finally, we then reduced the dimensionality of the trajectory using the Darter method and estimated the arrangement of cells using the orderCells function. The differentiation time of cell subsets was then identified clustering characteristics and gene markers.

R “CellChat” package (version 1.0.0) was used to analyze cellular communications. In sepsis, we investigated the differences and connections between immune cells during the current regeneration process, as well as the interactions of regeneration-related subsets with other cell types via cell communication analysis. Finally, we selected 4000 immune cells for processes.

High-throughput molecular research results are often translated into biological applications by analyzing gene function. We analyzed gene function using the R “clusterProfile” package (version 3.17) and visualized hub genes containing Disease Ontology (DO), Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). The significance level was set at P<0.05.

Multiple machine learning algorithms were applied to this study to identify feature genes. The first step was to perform Least Absolute Shrinkage and Selection Operator (Lasso) to screen candidates by iteratively reweighting least squares. Feature variables were selected based on minimum criteria after running the algorithm for 1000 cycles. Next, seven machine learning algorithms (k-Nearest Neighbors, Linear Discriminant Analysis, Logistic Regression, Naive Bayes, Random Forest, Recursive Partitioning, and Support Vector Machine) were applied. Finally, the receiver operating characteristic (ROC) curve was used to assess the classification performance of the hub genes in both training and validation cohorts. The pROC software was used to create ROC curves and the area under the curve was calculated.

Single-cell regulatory network interference and clustering (SCENIC, version 1.2.4) was created specifically for single cell data. Following the introduction of the gene co-expression network derived from transcription factor (TF) motifs by SCENIC, high-reliability gene regulatory networks (GRNs) predominately composed of TFs were found. We further highlighted the transcription factors of hub genes.

CIBERSORT was used to determine immune cell proportions (13). Correlation analysis was then conducted to analyze the relationship between immune cells and eight hub genes (ACTB, ALPL, CD177, GAPDH, S100A8, S100A9, SLC25A37, and STXBP2).

Gene Set Enrichment Analysis (GSEA) (https://www.broadlnstitute.org/gsea/) and single sample Gene Set Enrichment Analysis (ssGSEA) were to explore the biological functions and gene expressions according to all gene sets in the Hallmark database.

The Gene Set Variation Analysis (GSVA) is an unsupervised, non-parametric method. Based on the R “GSVA” package (version 1.42.0) transformation of gene expression data, we analyzed the data as a feature derived from the expression matrix of each single gene. Rank statistics were used to calculate the gene sets for each feature, and expression matrices were converted into Enrichment Score (ES) matrixes. The GSVA enrichment score for each sample could be obtained, facilitating statistical analysis of the data.

Cecal Ligation and Puncture (CLP) mice sepsis models: After inhaling isoflurane, eight-week C57 mice were anesthetized, and a stump was punctured with a 22-gauge needle once to expel stool. Afterward, the cecum was repositioned intraabdominally, and the abdomen closed. For fluid resuscitation, 0.2 ml of saline was administered intraperitoneally. Sham-operated mice did not undergo puncture or ligation. Mice were returned to cages and sacrificed 4 hours after the operation. The mice were intraperitoneally injected with phosphate-buffered saline (PBS) or Lipopolysaccharide (LPS) (20 mg/kg) to establish LPS-induced sepsis model.

Mice of sepsis model were euthanized by cervical dislocation after injecting sodium pentobarbital intraperitoneally. Peripheral blood neutrophils were isolated from mice using the neutrophil isolation kit (130-097-658, Miltenyi).

Following the manufacturer’s instructions, TRIzol reagent (12183-555, Invitrogen) was used to extract RNA from neutrophils. Reverse transcription was carried out using Takara’s Prime-Script Rase. After the premix Ex-Taq (Takara), the gene expression level was determined by qRT-PCR and normalized to the beta-actin (β-Actin). The expression level was calculated using the 2^-ΔΔCt technique. Table 1 contains a list of the primer pairs utilized in the experiments. The study was approved by the Ethics Committee of General Hospital of Northern Theater.

To perform the statistical analysis, R software (version 4.2.3, available at https://www.r-project.org) and GraphPad Prism 8.0 software were used. The two-tailed Student t-test determined the statistical significance. Co-expression was adjusted using Pearson’s correlation. P<0.05 was deemed statistically significant.

The flowchart of this study was showed in Figure 1 .

To obtain high-quality cell samples from sepsis tissues, we preprocessed scRNA-seq data using stringent quality control metrics. We obtained 32,396 high-quality cell samples ( Supplementary Figure S1A ). The Pearson correlation coefficient of 0.82 showed a strong positive correlation between the number of genes detected and sequencing depth ( Supplementary Figure S1C ). The “Harmony” method and the “PCA” method was then used to reduce dimension and eliminate redundant samples ( Supplementary Figures S1B, D, E ). Unsupervised transcriptome analysis was performed using UMAP diagrams ( Figures 2A, B ). Cell abundance in healthy control and sepsis samples was showed ( Supplementary Figures S2A, B ). Seven separate clusters (B cells, Monocytes, T cells, Neutrophils, NK cells, Platelets, and CD34^- Pre-B cell) were showed on the UMAP diagram ( Figure 2C ). The ratio and proportion of those seven cell clusters in healthy control and sepsis samples were visualized ( Figures 2D–F ). In sepsis, Neutrophils, Platelets, and CD34-Pre-B cells tended to be increased, while other cell clusters (B cells, Monocytes, T cells, and NK cells) were on a decreasing trend. As neutrophils play a vital role in sepsis, our study then aimed to explore sepsis-related neutrophils. To explore the changes of neutrophils in sepsis, we first reduce dimension by the “PCA” method ( Figure 2G ). The ElbowPlot identified the optimal number of pcs as 11 ( Figure 2H ). UMAP diagram also showed abundances of neutrophils between sepsis and healthy controls ( Figure 2I ). Eight clusters of neutrophils were then identified. The ratio of eight separate clusters were also showed ( Figure 2J ). Cell numbers of each neutrophil cluster was showed in Table 2 . Cluster 1&2 was not found in healthy controls, however, cluster 1&2 significantly increased in sepsis. Cluster 3&4 also increased in sepsis and cluster 0 significantly decreased in sepsis. The result indicated cluster 1&2 of neutrophils may be the potential responder cells against infection. It was found that Sepsis_Neutrophils-Monocyte and Sepsis_Neutrophils-Common_Neutrophils had higher interactions ( Figure 2K ). Moreover, interactions of Sepsis_Neutrophils-Monocyte and Sepsis_Neutrophils-Common_Neutrophils may rely on RETN-CAP1 signal transduction pathway ( Figure 2L ). However, Sepsis_Neutrophils showed no signifcant correlation with CCL signaling pathway, which plays an important role in immune response ( Figure 2M ). Furthermore, incoming, and outgoing signaling patterns analysis demonstrated that Sepsis_Neutrophils correlated with Resistin, Annexin, Visfatin, and IL1 ( Figure 2N ). Moreover, Sepsis_Neutrophils and Common_Neutrophils showed a significant difference on both incoming interaction strength and outgoing interaction strength ( Figure 2O ).

After standardization of data, the trajectory analysis projects all neutrophils cells onto three states ( Figures 3A, B ). The result showed that eight neutrophils’ clusters appeared in the pseudo-timeline, and the cluster 1&2 were most abundant in the state 2&3 ( Figures 3C, D ). The hub genes in cluster 1&2 were showed. Interestingly, all 50 hub genes all showed an elevated expression in the late stage of neutrophil-development ( Figure 3E ). This result was in consist of the trajectory analysis. Furthermore, hdWGCNA analysis was used for exploring gene modules related to cluster 1&2 of neutrophils. The soft threshold was adjusted to 5 for scale-free network construction ( Figure 3F ). Then the adjacency matrix and the TOM was built. Finally, 3 modules (Turquoise, brown, and blue modules) were identified based on average hierarchical clustering and dynamic tree clipping ( Figure 3G ). K-means methods identified hub genes in each module ( Figure 3H ). RPS18, RPS23, RPL32, RPL22, RPL11, EEF1A1, RPL23A, RPS8, RPL13, and RPL18A in the turquoise module. PTPRC, C16orf27, TLE3, ROCK1, H2AFY, HIF1A, KLF6, ZEB2, CTNNB1, and MALAT1 in brown module. GADD45B, IFITM2, LCN2, CYSTIM1, S100A12, CST7, ALPL, MMP9, CD63, and CD177 in blue module. The correlation analysis showed that the blue module was positively correlated with the brown module and negatively correlated with the turquoise module ( Figure 3I ). Moreover, among three modules, the blue module expressed highest in cluster 1&2 of neutrophils ( Figure 3J ). As cluster 1&2 of neutrophils may be the potential responder cells against infection, the genes in the blue module may have potential value for research. DO, GO, and KEGG enrichment analysis were performed on the genes in cluster 1&2 neutrophils ( Figures 3K–M ). The results showed that most of genes were enriched in the neutrophil degranulation, neutrophil activation, and neutrophil mediated immunity. IL-17 signaling pathway may be the potential mechanism in the neutrophil-related function. The pathway-network and pathway-gene interaction analyses were then showed.

GSE57065 was regarded as training dataset and GSE95233 was validation dataset. Firstly, these two datasets were normalized before analysis ( Supplementary Figures S3A–D ). Secondly, datasets GSE57065 and GSE95233 were merged and normalized ( Supplementary Figures S3E, F ). Thirdly, meta-analysis was used to eliminate the batch effect ( Supplementary Figures S3G, H ). The differences in expression of 50 hub genes were showed ( Figures 4A–C ). Moreover, correlations between every two genes were also identified ( Figure 4D ). Hub genes of cluster 1&2 Sepsis_Neutrophils were subjected to LASSO regression and 8 features were selected (ALPL, ACTB, CD177, GAPDH, SLC25A37, S100A8, S100A9, and STXBP2) ( Figures 4E, F ). The coefficient of each feature was listed as Table 3 . Besides, seven machine learning algorithms (k-Nearest Neighbors, Linear Discriminant Analysis, Logistic Regression, Naive Bayes, Random Forest, Recursive Partitioning, and Support Vector Machine) were utilized and mean AUCs of seven algorithms were compared ( Figure 4G ). Only the Recursive Partitioning showed a lower AUC compared with other six algorithms ( Figure 4H ). Moreover, the AUC in the validation dataset was 0.988 and the model had a strong generalization ability ( Supplementary Figure S4 ). We demonstrated the expression of eight hub genes in the training dataset GSE57065. All eight hub genes showed statical higher expressions in neutrophils under sepsis ( Figure 4I ). ROC curves showed that five hub genes (except ACTB, S100A8, and SLC5A37) all showed higher AUCs ( Figure 4J ). In validation dataset GSE95233, six hub genes except ACTB and SLC5A37 expressed higher significantly in neutrophils under sepsis ( Supplementary Figure S5A ). ROC curves showed that six hub genes except ACTB and SLC5A37 all showed higher AUCs ( Supplementary Figure S5B ).

Significant differences were observed in B cell naïve, B cells memory, Plasma cells, T cells CD8⁺, T cells CD4⁺ naïve, T cells CD4⁺ memory resting, T cells follicular helper, T cells gamma delta, NK cells resting, Macrophages M0, Macrophages M1, Macrophages M2, Dendritic cells resting, Dendritic cells activated, Mast cells resting, Mast cells activated, Eosinophils, and Neutrophils ( Figures 5A, B ). Moreover, the correlation between hub genes and immune cells were identified and visualized ( Figure 5C ). Hub genes were positively associated with Macrophage M0 cells. Hub genes negatively correlated with T cells CD4⁺ memory resting, T cells CD8⁺, T cells gamma delta, and T cells CD4⁺ memory activated. To explore the correlation between each hub gene and immune cells, correlation between 20 types of immune cells and eight hub genes were showed ( Supplementary Figures S6 – S13 ). Among eight hub genes, only S100A9 and STXBP2 showed no significant correlation with neutrophils. While correlations with other immune cells were also showed.

ssGSEA and GSEA showed the potential biological pathway of eight hub genes in sepsis ( Figures 6A, B ). GSVA analysis showed that expression and correlation between eight hub genes and relative pathways in sepsis ( Figures 6C, D ). GSVA analysis showed that eight hub genes of neutrophils in sepsis were mainly enriched in Glycosaminoglycan_Degradation, Insulin_Signaling_Pathway, Prostate_Cancer, and Melanoma. Moreover, GSEA analysis reveals that related pathways of hub genes. Potential pathways that associated with up-regulation of eight hub genes(ALPL, ACTB, CD177, GAPDH, SLC25A37, S100A8, S100A9, and STXBP2)were showed in Supplementary Figure S14A . Results showed higher expression of eight hub genes mostly associated with starch and sucrose metabolism. Additionally, pathways of ribosome and graft-verus-host disease that associated with down-regulation of eight hub genes (ALPL, ACTB, CD177, GAPDH, SLC25A37, S100A8, S100A9, and STXBP2) were showed in Supplementary Figure S14B .

SCENIC was to identify TFs with gene expression activity among subclusters of neutrophils. We applied the SCENIC analysis pipeline to the neutrophils and resolved distinct regulons associated with each cluster in neutrophils. Regulons were robust in activity and specific for each cell type ( Figures 7A, B ). For example, FOS_extended_(36g) and FOS_(28g)were lower active in the sepsis-specific neutrophils. Moreover, SPI1_extended_(38g), SPI1_(37g), and HMGB1_(11g) were highly active in the sepsis-specific neutrophils ( Figure 7C ). Heatmap of 25 TFs between sepsis and healthy control were showed ( Figure 7D ). The top ten TFs with the most difference between sepsis and healthy control, including SPI1, HMGB1, JUN, NFIL3, RUNX1, CEBPD, FOS, REL, JUND, and FOSB ( Figure 7E ). The correlation analysis of the inferred TFs and 8 hub genes was performed. APLP, CD177, GAPDH, S100A9, and STXBP2 were significant associated with various TFs ( Figure 7F ).

Moreover, gene correlations were identified and visualized ( Figure 8A ). CD177 and S100A9 showed a highest correlation with significance. To evaluate the expression of eight hub genes in neutrophils under sepsis. We established two sepsis model: CLP model and LPS-induced model. qRT-PCR detected the relative expression of all eight genes in monocytes. Results showed that ALPL, CD177, S100A8, S100A9, and STXBP2 significantly up regulated in sepsis ( Figure 8B ). In conclusion, ALPL, CD177, S100A8, S100A9, and STXBP2 may be play a vital role in neutrophils function response to sepsis.