Two single-cell transcriptome datasets related to sepsis were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov), specifically GSE175453 (which includes 9 samples) and GSE167363 (which includes 12 samples). In total, 21 single-cell transcriptomic samples were utilized for various analyses, including cell annotation, identification of neutrophils and their subtypes, pseudo-time series analysis, and the identification of differentially expressed genes in neutrophils from healthy individuals compared to those from sepsis patients. These datasets, capturing peripheral blood mononuclear cells (PBMCs) and leukocytes from sepsis patients, reflect systemic immune responses potentially relevant to organ dysfunction, including SAKI, though they lack explicit kidney-specific clinical annotations based on KDIGO criteria (e.g., serum creatinine or urine output).

Cell tags, genes, and expression matrices were obtained directly from the GEO database, and the Seurat v4 package (10) was utilized to process the data. The filtering criteria applied to the pooled samples included: 300 < nFeature_RNA < 2000, nCount_RNA < 40000, percent.mt (ratio of mitochondrial genes) < 30, and log10GenesPerUMI > 0.85.

The Seurat pipeline was utilized to reduce the dimensionality of the combined cell data. First, principal component analysis (PCA) was conducted and visualized, followed by the creation of an elbow plot to identify the most significant principal components. The harmony algorithm (11) was then employed to address batch effects among various samples. After that, Dimensionality reduction was performed on the selected principal components (PCs) using both UMAP and t-SNE. Subsequently, cluster plots were generated to determine the optimal resolution parameter for partitioning cells into biologically distinct clusters.

The dataset was annotated into different cell types, such as T lymphocytes, B lymphocytes, neutrophils, and plasma cells, using the CellMarker2 database. Visualization was enhanced through bubble plots of marker genes and UMAP annotation plots.

Differential expression analysis was carried out using the FindAllMarkers function, aiming to identify differentially expressed genes within the same group across different clusters, as well as those within the same cluster but across different groups. Genes were selected as differentially expressed if they met the criteria of |logFC| > 0.25 and a P-value < 0.05.

Pseudo-time trajectory analysis was conducted using the Monocle package (12). Highly variable genes were selected based on an average expression threshold (>0.1) and ranked within the top 2000 by coefficient of variation. Dimensionality reduction was performed using the DDRTree algorithm (13), and cellular trajectories were visualized through scatter plots generated by the built-in plot_cell_trajectory function. Genes were partitioned into two distinct clusters based on their pseudo-temporal expression patterns. Differential gene expression along the pseudo-time axis was visualized using the plot_pseudotime_heatmap function. Gene Ontology (GO) enrichment analysis for temporally dynamic pathways was subsequently performed with the clusterProfiler package (14), employing a significance threshold of p < 0.05 and false discovery rate (FDR) < 0.25.

Cell-cell communication analysis was performed on annotated cell subpopulations using the CellChat package (15). Ligand-receptor interactions were predicted based on a predefined signaling pathway database. Key communication patterns were visualized through three types of diagrams: (1) Cell-cell interaction networks representing connections between cell clusters; (2) Incoming signaling heatmaps showing receptor activity in target cell types; (3) Outgoing signaling heatmaps displaying ligand expression in source cell types. All visualizations were generated using default settings of the CellChat package.

Metabolic pathway activities across cell types were analyzed using the scMetabolism package (16). Heatmaps were generated to visualize pathway activity differences.

Publicly available datasets were retrieved from the Gene Expression Omnibus (GEO) database The following datasets were included:

Raw transcriptomic data were downloaded from the GEO database using the GEOquery package (20). Platform-specific annotation files were utilized to map probes to gene symbols. Expression matrices were normalized via the RMA algorithm (21), including background correction, quantile normalization, and probe summarization. For genes with multiple probes, expression values were averaged to generate gene-level measurements.

Differentially expressed genes (DEGs) were identified in the E. coli-induced sepsis dataset using the limma package (21). Linear models comparing sepsis patients and controls were fitted after variance-stabilizing normalization. Statistical significance was defined as |log₂FC| > 0.5 and p < 0.05, with FDR correction applied to control false discoveries (FDR < 0.1).

GO and KEGG analyses were performed using the clusterProfiler package (14, 22) Gene sets included pseudotime-associated DEGs (Section 2.4) and E. coli-sepsis DEGs (Section 2.8). Significant terms (FDR < 0.05) were visualized as lollipop charts (top 9 KEGG pathways) and dot plots (top 3 GO terms per category).

The CNSKnowall Cloud Platform (https://www.cnsknowall.com/veen) was utilized to generate a Venn diagram illustrating the overlap between NETF genes, neutrophil differential genes, and differential genes identified through Bulk-RNAseq. This analysis enabled the identification of neutrophil-urosepsis differential genes. A linear SVM model was trained on 70% of GSE185263 and validated on 30%, with diagnostic performance assessed by ROC curves using the pROC package (23). Hub genes with AUC > 0.7 were prioritized.

To characterize the immune microenvironment in urosepsis, samples from the GSE185263 dataset were stratified into normal and sepsis groups. Immune cell infiltration profiles were quantified using the xCell algorithm (24), which estimates the abundance of 64 immune and stromal cell types based on gene expression signatures. The proportions of infiltrating immune cells were visualized through stacked bar plots generated with the ggplot2 package (version 3.4.0).

To identify differential immune infiltration between high-risk and low-risk subgroups (defined by the hub gene expression signature), the Wilcoxon rank-sum test was applied with a significance threshold of p < 0.05. Statistically significant cell types were further visualized using violin plots to display distribution differences. Associations between hub genes (PADI4, CR1, MAPK14, CASP4) and immune cell infiltration were analyzed via Spearman correlation and represented as lollipop plots, where the dot size and color intensity reflected correlation coefficients (ρ) and significance levels, respectively.

Additionally, the interplay between hub genes and immune checkpoint molecules (e.g., PD-L1, TIM-3, CTLA-4) was explored through heatmaps, with hierarchical clustering applied to both rows (genes) and columns (checkpoint genes). Immune cell types showing significant co-enrichment patterns were highlighted in enrichment bubble plots, where bubble size indicated the -log10(p-value) and color denoted enrichment scores.

To investigate the biological pathways associated with hub genes in the disease, Gene Set Enrichment Analysis (GSEA) was performed using the training dataset from GSE185263. The “c2.cp.kegg.v7.0.symbols.gmt” gene set from the Molecular Signatures Database (MSigDB) was used as the reference gene set. First, Spearman correlation coefficients between each hub gene and all other genes in the dataset were calculated using the R package “psych” (version 2.5.3) (25). For each hub gene, a ranked gene list was generated by sorting all genes in descending order based on their correlation coefficients. Subsequently, GSEA was conducted using the R package “clusterProfiler” (version 3.8) to explore the potential functional roles of the hub genes (14). Analysis parameters were set at p < 0.05, FDR < 0.25, and |NES| > 1, with hub gene enrichment results visualized and reported.

In this study, data preprocessing and consensus clustering were performed to identify the optimal clustering approach. Initially, a specific gene set, denoted as gene_Consensus, was selected from the raw dataset. Each gene (row) was then subjected to median centering by subtracting its row-specific median. For the consensus clustering analysis, Pearson correlation was employed as the distance metric, with partitioning around medoids (PAM) as the clustering algorithm. The analysis parameters were set as follows: a maximum cluster number of 6, 1000 iterations, a subsampling rate of 80% of the samples, inclusion of all features, use of complete observations for correlation calculations, and a fixed random seed of 814 to ensure reproducibility. The optimal number of clusters (k) was determined by incrementally increasing k from 2 to 6, with k = 2 identified as the optimal clustering solution based on clustering stability.

To investigate the subcellular distribution of hub genes, subcellular localization analysis was performed using the GeneCards database (https://www.genecards.org/). The localization of these genes to specific organelles, along with their confidence scores, was visualized as bar plots using the R package “ggplot2”.

To explore existing drugs associated with the identified hub genes, drug-gene interaction data were retrieved from the DGIdb database (https://www.dgidb.org/). The relationships between hub genes and related drugs, including their interaction scores, were visualized as bar plots using the R package “ggplot2”.

The experiment utilized male Wistar rats (8–10 weeks old, 160–200 g; SLAC Laboratory Animal Co., Ltd., Shanghai, China), housed at 22 °C with 50% humidity and a 12:12 h light–dark cycle, with ad libitum access to standard chow and water. Rats were acclimated to the environment for at least one week prior to experimentation. All procedures were approved by the Institutional Animal Care and Use Committee (IACUC) and complied with national and international guidelines for animal experimentation. Rats were randomly assigned to three groups (n=3 per group): Sham group (laparotomy without cecal ligation or puncture), CLP group (cecal ligation and double puncture with a 20-gauge needle to induce sepsis), and CLP + sh-PADI4 group (CLP surgery followed by tail vein injection of 1×10^9 TU PADI4-targeting shRNA adenovirus). The CLP model was performed under aseptic conditions. Rats were anesthetized with 2.5% isoflurane, underwent midline laparotomy, and the cecum was ligated below the ileocecal valve and punctured twice with a 20-gauge needle, with a small amount of fecal matter extruded before abdominal closure. The Sham group underwent cecal exposure without ligation or puncture. Adenovirus was administered via slow tail vein injection, with all interventions conducted under anesthesia. Rats were fasted for 8–12 hours pre-surgery, recovered in a warm environment post-surgery, and monitored closely, adhering to predefined humane endpoints to minimize suffering. At 48 hours post-CLP, rats were euthanized under deep anesthesia, blood was collected via cardiac puncture, and serum was stored at -80°C. Both kidneys were harvested—one fixed in 10% neutral-buffered formalin or 4% paraformaldehyde for histological and immunofluorescence analysis, the other snap-frozen in liquid nitrogen and stored at -80°C for molecular assays. Aseptic techniques were maintained throughout, and sepsis-related tissues were disposed of according to biosafety protocols. The animal study was approved by the Institutional Animal Care and Use Committee of The Second Xiangya Hospital, Central South University (Approval No: 2022711).

Kidney tissues from the CLP and Sham groups were homogenized in RIPA buffer containing protease inhibitors and centrifuged at 12,000 × g for 15 min at 4°C to extract proteins. Protein concentrations were determined using the BCA method. Thirty micrograms of protein were mixed with loading buffer, heated at 95°C for 5 min, separated by 10% SDS-PAGE, and transferred to a PVDF membrane. The membrane was blocked with 5% non-fat milk for 1 h, incubated overnight at 4°C with primary antibodies against PADI4 (1:1000, Abcam) and β-actin (1:5000, Sigma-Aldrich), washed with TBST, and incubated with HRP-conjugated secondary antibodies (anti-rabbit 1:5000, anti-mouse 1:10000) for 1h. After further TBST washes, bands were visualized using ECL detection. Band intensities were quantified with ImageJ, and PADI4 expression was normalized to β-actin.

Kidney tissues from experimental groups were fixed in 4% paraformaldehyde at 4°C for 24–48 h, dehydrated, and embedded in paraffin. Sections (4–5 µm) were cut, mounted on adhesive slides, and baked at 60°C for 30 min. After deparaffinization in xylene and rehydration through a graded ethanol series, antigen retrieval was performed in citrate buffer (pH 6.0) at 95°C for 15 min. Sections were cooled, washed with PBS, and blocked with 5% goat serum in PBS for 30 min at room temperature. Sections were then incubated with anti-citH3 antibody (1:200, Abcam) overnight at 4°C, washed with PBS, and incubated with Alexa Fluor 488-conjugated anti-rabbit secondary antibody (1:500, Thermo Fisher Scientific) for 1 h at room temperature in the dark. Nuclei were counterstained with DAPI for 5 min, and slides were mounted with anti-fade medium. Images were captured using a fluorescence microscope and quantified for citH3-positive cells using ImageJ software.

Kidney tissues from experimental groups were fixed in 10% neutral-buffered formalin at room temperature for 24–48 h, dehydrated through a graded ethanol series, cleared in xylene, and embedded in paraffin. Sections (4–5 µm) were cut using a microtome, mounted on adhesive slides, and baked at 60°C for 30 min. Sections were deparaffinized in xylene, rehydrated through decreasing ethanol concentrations, and stained with hematoxylin for 5 min. After rinsing in tap water, sections were differentiated in 1% acid alcohol, blued in ammonia water, and counterstained with eosin for 1 min. Sections were dehydrated, cleared in xylene, and mounted with a coverslip using mounting medium. Images were captured using a bright-field microscope, and histopathological changes were evaluated.

Peripheral blood samples were collected from 5 sepsis patients and 5 healthy controls. And the study of these samples was approved by the Ethics Committee of Zhuzhou Hospital Affiliated to Xiangya School of Medicine, Central South University (Approval No: 2020175-01). Total RNA was extracted from whole blood using TRIzol reagent (Invitrogen, USA) per the manufacturer’s instructions. cDNA was synthesized from 1 µg total RNA using the PrimeScript RT reagent kit (Takara, Japan). Quantitative real-time PCR was performed on a CFX96 system (Bio-Rad, USA) with SYBR Green premix (Takara, Japan). Expression of PADI4, CASP4, CR1, and MAPK14 was measured, normalized to β-actin (primer sequences in Table 1 ). Relative gene expression was calculated using the 2^−ΔΔCt method.

All analyses were executed utilizing R programming language (v 4.2.2). The Wilcoxon test was harnessed to contrast differences, employing a statistical threshold of p < 0.05.

¹ Material and Methods section can be placed in any of the following ways.

21 single-cell transcriptome samples were retrieved from the GEO database (GSE175453, 9 samples; GSE167363, 12 samples). Quality control was performed using Seurat v4 with the following criteria: 300 < nFeature_RNA < 2000, nCount_RNA < 40000, percent.mt < 30, and log10GenesPerUMI > 0.85. This confirmed a high-quality dataset ( Figures 1A–F ). Batch effects were corrected using the Harmony algorithm, followed by principal component analysis (PCA) to select significant components ( Figure 2A ). UMAP dimensionality reduction was applied with an optimal resolution determined by cluster tree analysis ( Figure 2B ), clustering cells into 16 groups with distinct marker gene expression profiles ( Figure 2D ). These were annotated as 9 distinct cell types, including B cells, T/NK cells, monocytes, and neutrophils ( Figures 2C, E ). In the sepsis group, the proportion of neutrophils was significantly elevated (10.63% vs. 8.94% in the healthy group), while T/NK cells showed a slight reduction (34% vs. 37%), and monocytes remained consistent (27%, Figure 2H ). Marker gene analysis validated cell identities, with neutrophils highly expressing S100A8 and CSF3R, and T/NK cells expressing CD3E and KLRD1 ( Figures 2F, G ), providing a foundation for subsequent neutrophil-focused analyses.

Neutrophils were subdivided into 8 clusters using UMAP, which were further categorized into four subtypes based on immune-related and metabolic gene expression profiles: mature neutrophils, immature neutrophils, anti-inflammatory neutrophils, and pro-inflammatory neutrophils ( Figures 3A, C ). Anti-inflammatory neutrophils, defined by elevated MT-ND1, MT-CO1, MALAT1, and CST3, were more prevalent in the control group (27.04% vs. 18.43% in sepsis, Figure 3D ), while pro-inflammatory neutrophils, distinguished by PF4, CXCR2, CXCL8, and CCL5, exhibited a slight increase in sepsis (6.72% vs. 5.33% in controls, Figure 3D ), suggesting an intensified inflammatory state. Anti-inflammatory neutrophils were identified by high expression of MT-ND1, MT-CO1, MALAT1, and CST3, which are associated with mitochondrial homeostasis and inflammation resolution (26–28). Pro-inflammatory neutrophils were characterized by elevated PF4, CXCR2, CXCL8, and CCL5, promoting chemotaxis and inflammatory amplification in sepsis (29–33). Minimal cross-expression of marker genes across subtypes validated the robustness of this classification ( Figures 3B, E ). The substantial increase in immature neutrophils in the sepsis group (40.53% vs. 4.19% in controls, Figure 3D ) underscored their mobilization during acute inflammation. The decrease in anti-inflammatory neutrophils (18.43% vs. 27.04%) and the slight rise in pro-inflammatory neutrophils (6.72% vs. 5.33%) further indicated an immune imbalance, potentially contributing to the exacerbated inflammation in sepsis.

Differential expression analysis (GSE237960, Escherichia coli-induced sepsis) identified 4646 differentially expressed genes (DEGs), with 2312 upregulated and 2334 downregulated ( Figure 4A ). To narrow down key neutrophil-specific genes, DEGs were intersected across neutrophil pseudotime, neutrophil extracellular trap formation (NETF) pathways, and E. coli-induced sepsis, resulting in 15 shared genes ( Figure 4B ). These DEGs were enriched in chemokine signaling, cytokine-cytokine receptor interaction, and Toll-like receptor pathways ( Figures 4C, D ), reflecting enhanced immune and inflammatory responses in sepsis. Cell communication analysis revealed increased IL-1 signaling (e.g., IL1A-FPR1, IL1B-FPR1) in the sepsis group, suggesting a stronger pro-inflammatory response and potential therapeutic relevance of IL-1 receptor antagonists (e.g., Anakinra). Suppression of HLA-DRB1-CD4 signaling ( Figures 5A–D ) indicated impaired antigen presentation, which may contribute to adaptive immune dysregulation. In the control group, metabolic activities such as glycolysis and fatty acid metabolism were predominant in mature neutrophils ( Figures 6A, B ). In the sepsis group, we observed a significant upregulation of biotin and lipoic acid metabolism in mature neutrophils ( Figure 6B )—a pattern also observed in immature and pro-inflammatory subtypes—suggesting metabolic reprogramming under inflammatory conditions. Pseudotime analysis revealed a continuous trajectory of neutrophil states ( Figure 7A ), with distinct state distributions reflecting developmental progression ( Figure 7B ). A shift toward an early pro-inflammatory state was observed in sepsis neutrophils compared to controls ( Figure 7C ), with immature subtypes enriched in early stages and mature and anti-inflammatory subtypes in later stages ( Figure 7D ). The expression patterns of the 15 shared genes along the pseudotime axis further highlighted their dynamic regulation, with genes like PADI4 showing elevated expression in early pro-inflammatory stages ( Figure 7E ). Gene clustering revealed two expression patterns: Cluster 1 (low early, high late) enriched in immune processes, and Cluster 2 (high early, low late) enriched in protein synthesis ( Figures 8A–C ), illustrating dynamic neutrophil regulation in sepsis.

To identify hub genes linked to neutrophil function in sepsis, differential gene expression data were integrated and validated externally using a support vector machine (SVM) model. An SVM model trained on 70% of the GSE185263 dataset identified PADI4, CASP4, CR1, and MAPK14 as significant feature genes ( Figure 9A ). Internal validation yielded an AUC of 0.9372, with only 2 misclassified samples ( Figures 9B, C ). External validation using GSE243217 confirmed high AUC values: PADI4 (92.4%), CASP4 (90.6%), CR1 (97.3%), and MAPK14 (93.9%) ( Figures 9D, E ). The consistent differential expression and high AUC values across datasets underscored their reliability as hub genes in sepsis.

Gene set enrichment analysis (GSEA) indicated that all four hub genes were enriched in immune and inflammation-related pathways, suggesting synergistic roles in sepsis immune regulation. Specifically, PADI4 was enriched in Pathogenic Escherichia coli Infection (NES = 1.62, p = 8.56e-07, FDR = 1.22e-05) and Fc Gamma R-Mediated Phagocytosis (NES = 1.62, p = 9.31e-10, FDR = 1.47e-07, Figure 10A ). CASP4 showed enrichment in Leishmania Infection (NES = 1.70, p = 1.03e-08, FDR = 2.80e-07), B Cell Receptor Signaling Pathway (NES = 1.69, p = 1.04e-09, FDR = 3.88e-08), and Proteasome (NES = 1.67, p = 1.10e-06, FDR = 1.14e-05, Figure 10B ). CR1 and MAPK14 were enriched in B Cell Receptor Signaling and Renal Cell Carcinoma pathways ( Figures 10C, D ). Subcellular localization analysis revealed PADI4 primarily in the cytosol and nucleus, CASP4 in the endoplasmic reticulum and cytosol, CR1 in the extracellular space and plasma membrane, and MAPK14 in the cytosol, extracellular space, and nucleus ( Figures 11A–D ), supporting their functional diversity. Drug-gene interaction analysis suggested potential therapeutic candidates, including uric acid for PADI4, emricasan for CASP4, and RO-3201195 for MAPK14 ( Figures 12A–D ), offering directions for targeted sepsis therapies. Consistent with this prediction, preclinical studies demonstrated that febuxostat treatment lowered tissue urate levels and significantly reduced citH3 expression in mice (34), indicating suppression of PAD4-dependent NETosis in vivo.

Immune microenvironment analysis revealed significantly increased neutrophil infiltration (P<0.001) and reduced T cell infiltration (e.g., CD4+ memory T cells, P<0.001) in the sepsis group ( Figures 13A, C ), with a strong negative correlation between neutrophils and T cells (P<0.001, Figure 13B ). PADI4, CASP4, CR1, and MAPK14 were positively correlated with neutrophil infiltration (P<0.001) and negatively correlated with T cell infiltration (P<0.001, Figure 13D ), and positively associated with immune checkpoint genes (CD24, CD47, P<0.01, Figure 13E ), suggesting immunosuppressive mechanisms. The increased neutrophil infiltration and decreased T cell presence indicated an imbalance between innate and adaptive immunity in sepsis, potentially linked to disease severity. To corroborate these xCell-derived profiles, CIBERSORT analysis was applied to the same bulk RNA-seq data, estimating proportions of 22 immune cell types. Consistent with xCell findings, it confirmed significantly elevated neutrophil fractions in the sepsis group (P < 0.001; Supplementary Figures 1A, C) and reduced T cell subtypes, such as CD4 memory activated T cells . Neutrophils and T cells showed a strong negative correlation (Supplementary Figure 1B), while the hub genes exhibited positive correlations with neutrophil infiltration (r = 0.3–0.7, Supplementary Figure 1C). These genes were also positively associated with immune checkpoint molecules (Supplementary Figure 1E), reinforcing the immunosuppressive microenvironment. Consensus clustering divided patients into two clusters: Cluster A exhibited lower PADI4 and CASP4 expression and higher CR1 and MAPK14 expression ( Figures 14A–D ), highlighting molecular heterogeneity. However, survival differences between clusters were not statistically significant (p=0.7594), possibly due to limited sample size, necessitating further clinical validation.

To investigate the role of hub genes in sepsis-associated kidney injury, we employed a cecal ligation and puncture (CLP) model to induce sepsis in mice, followed by extraction of kidney tissues for in vitro experimental analyses to assess the expression differences of PADI4. Western blot (WB) analysis revealed a significant upregulation of PADI4 expression in the kidney tissues of the CLP group compared to the sham-operated group, with quantitative analysis (normalized to β-actin) confirming a markedly higher relative expression level (p < 0.001, Figures 15A, B ). The full uncropped Western Blot images are provided in Supplementary Figure S2 . Immunofluorescence (IF) staining further demonstrated a pronounced increase in the red fluorescence signal of citrullinated histone H3 (citH3, a marker of neutrophil extracellular trap [NET] formation) in the kidney tissues of the CLP group, significantly exceeding that of the sham group (p < 0.001). Notably, knockdown of PADI4 resulted in reduced citH3 expression, suggesting a critical role for PADI4 in sepsis-induced NET formation (p < 0.001 compared to the CLP group, Figures 15C, D ). Hematoxylin and eosin (H&E) staining revealed significant pathological changes in the kidney tissues of the CLP group, characterized by tubular dilation, expansion of the tubulointerstitial area, and inflammatory cell infiltration, consistent with the features of sepsis-induced acute kidney injury. The injury score was significantly higher in the CLP group than in the sham group (p < 0.001). However, PADI4 knockdown markedly attenuated these pathological changes (p < 0.01 compared to the CLP group, Figures 15E, F ), indicating that PADI4 silencing significantly ameliorates sepsis-induced renal histopathological damage. Furthermore, we validated the elevated expression of the hub genes in the peripheral blood of patients with sepsis, consistent with our previous findings (compared to the control group, Figures 15G–J ).