PPIF neutrophils promote mtROS driven NETosis mediated progression of colorectal cancer.

Introduction
Colorectal cancer (CRC) is the third most common malignant tumor and the second leading cause of cancer death, with an estimated 1.9 million new cases and 900,000 deaths globally in 2020. It is worth noting that the incidence of CRC is still rapidly increasing, and even worse, the prognosis for postoperative recurrence and metastasis in CRC patients remains poor, with a 5-year survival rate of only 10% for stage IV CRC patients [1]. Therefore, researching the pathogenesis of CRC and urgently seeking new drug targets and treatment methods is crucial to improve the prognosis of CRC patients.
The tumor microenvironment (TME) is comprised of cancer cells and a diverse cohort of stromal cell populations, including infiltrating immune cells like lymphocytes and macrophages, neutrophils, fibroblasts, and adipocytes. It also incorporates the extracellular matrix (ECM) and an array of soluble factors and signaling molecules generated by these cells, forming an intricate cellular network [2, 3]. Various types of immune cells can impact the progression and metastasis of cancer cells [4], and Immunotherapy works by altering TME to inhibit tumor progression, but its efficacy may be influenced by TP53 status, as P53 not only regulates DNA repair but also modulates immune responses within the tumor microenvironment [5, 6]. In the realm of CRC immunotherapy, immune checkpoint inhibitors are a prevalent treatment modality. For instance, programmed cell death protein 1 (PD-1) plays a pivotal role in modulating T cell function and preserving immune system homeostasis, thus establishing itself as one of the most extensively researched regulatory agents [7]. Additional immune checkpoints encompass cytotoxic T-lymphocyte-associated protein 4 (CTLA-4), T cell immunoglobulin 3 (Tim-3), and lymphocyte activation gene 3 (Lag3) [8]. These immunotherapy targets engender advantageous conditions for the immune-based management of CRC [9]. Therefore, a deep understanding of stromal cells and immune cells can provide new insights for developing and utilizing novel therapeutic targets to control the progression of CRC.
It is worth noting that neutrophils account for 50–70% of circulating white blood cells in human peripheral blood. As short-lived cells, neutrophils play an indispensable role in both healthy and tumor tissues [10]. During innate immune response, neutrophils can form sticky web-like structures known as NETs [11, 12]. NETs are network structures released by neutrophils and primarily composed of decondensed chromatin and granule proteins. However recent investigations have uncovered the widespread presence of NETs within the tumor microenvironment of CRC and shown that NETs play a pro-cancer role in CRC progression [13]. They can promote tumor cell proliferation, invasion, and metastasis, providing tumor cells with a microenvironment conducive to growth and survival [14]. Meanwhile, the activation of epithelial mesenchymal transition like processes in CRC cells by NETs may also contribute to the progression of CRC [15]. Therefore, in-depth exploration of the specific mechanisms of action of NETs in CRC is of great significance for revealing the pathogenesis of CRC and identifying new therapeutic targets.
Single-cell RNA sequencing (scRNA-seq) is an innovative technology for profiling the transcriptome of individual cells. It empowers the examination of gene expression alterations at the single-cell level, providing novel insights into the interplay among diverse cell subgroups in tumor tissues [16]. In this study, we obtained scRNA-seq data from public databases, re-identified and annotated cell populations, determined specific subgroups of malignant cells and immune cells involved in immune responses, constructed cell differentiation trajectories, and explored the complex mechanisms of cell-cell interactions. We focus on how tumor associated NETs play a role in the TME, which may become a potential therapeutic target for improving treatment strategies and clinical outcomes.

Materials and methods

Data collection
scRNA-Seq data related to CRC were obtained from the Gene Expression Omnibus (GEO) database [17]. This study included 5 datasets, namely GSE161277 [18], GSE178318 [19], GSE188711 [20], GSE200997 [21], and GSE201348 [22]. A total of 37 CRC samples, 32 adjacent cancer normal samples, and peripheral blood mononuclear cell (PBMC) samples from 4 CRC patients were included. The detailed clinical characteristics of the patients were obtained from the original text and are listed in Table S1. Additionally, transcriptome data related to CRC (TCGA-COAD) and corresponding clinical information were downloaded from The Cancer Genome Atlas (TCGA) [23], encompassing 471 CRC samples and 41 adjacent noncancerous tissue samples for subsequent bulk analysis. Subsequently, data preprocessing was carried out and the Seurat package [24] was used for data standardization. GEO and TCGA are both open-source databases, and the patients involved in the databases have obtained ethical approval.

Processing scRNA seq data
Integrate the downloaded and preprocessed sample data to generate a feature count and digital gene expression matrix for a single cell with default parameters, forming an expression profile for a single cell. Subsequently, the scRNA-Seq data were integrated and converted into Seurat objects. Further quality control processing was conducted for nFeature (total number of genes), nCount (total gene expression), and the proportion of mitochondrial genes. Cells with the highest and lowest 1% feature counts and mitochondrial gene proportions exceeding 10% were filtered out.

Dimensionality reduction, clustering, and differential expression analysis
Standardize the data using the SC Transform function, and then perform unsupervised clustering and differential gene expression analysis on the preprocessed data using the Seurat software package. And use the FindClusters function to obtain cell clusters. FindAllMarkers function results Subsequently, uniform manifold approximation and projection (UMAP) algorithm [25] was utilized to visualize the single-cell landscape in low-dimensional space. Furthermore, the FindAllMarkers function in the Seurat package was employed to identify differentially expressed genes between specific cell clusters and other clusters, with FDR correction and a significance threshold of |log2FC| > 1. The results of the FindAllMarkers function are listed in Table S2. Based on the identified marker genes, as well as markers validated in previous single-cell studies and in the laboratory, different cell clusters were further categorized into known cell types.

Analysis of copy number variation in chromosomes
Copy number variation (CNVs) between individuals is a major source of genetic variation in cancer individuals. InferCNV is used to explore tumor single-cell RNA-Seq data (https://github.com/broadinstitute/infercnv) to determine whether somatic cells exhibit large-scale chromosomal CNVs. To verify whether epithelial cells and malignant cells are correctly distinguished, this study used the InferCNV method to estimate the chromosome copy number of epithelial cells and malignant cells and visualized the relative expression intensity of genes on each chromosome through a heatmap. This analysis selected hg19 as the reference genome.

Functional enrichment analysis
To explore the potential biological roles of identified cell subgroups in the CRC microenvironment, gene ontology biological process (GO BP) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed using the clusterProfiler package in R language [26]. Enrichment results with p < 0.05 were considered significant. Additionally, the AddModuleScore function [27] in the Seurat package was used to calculate the scores of the interested pathways or gene sets in different cell subgroups, and the results were visualized as violin plots.

Gene set enrichment analysis
Using the clusterProfiler package, gene set enrichment analysis (GSEA) was performed with the molecular signature database (MSigDB) [28] c5.bp.v7.4.symbols.gmt and c2.cp.kegg.v7.4.symbols.gmt gene sets as the background set to explore the activation or inhibition of the interested gene sets.

Pseudotime analysis
To explore the differentiation trajectory of cell subgroups in CRC development, pseudotime analysis of the interested cell populations was performed using the Monocle3 method [29] in this study. The developmental trajectories of cell subgroups were reconstructed, and the results were visualized using the UMAP algorithm.

Gene set variation analysis and survival analysis
Gene set variation analysis (GSVA) is a non-parametric unsupervised analysis method mainly used to evaluate gene set enrichment results from microarray and transcriptome data. In this study, the genes of interested cell subgroups were selected (LogFC > 1.5), and the GSVA scores of different cell subgroups in the TCGA-COAD clinical cohort were computed using the GSVA package [30]. Divide the samples into high and low score groups based on the median score. Then, use the R software package to analyze the survival rate and surfminer, and conduct an analysis of overall survival rate (OS) and recurrence-free survival rate (RFS).

Gene regulatory network analysis
In this study, gene regulatory networks (GRN) were analyzed using the Python module tool pySCENIC [31]. This workflow begins with a counting matrix that describes the abundance of all cell genes. Firstly, use the GRNBoost2 method to infer co expression modules. Subsequently, indirect targets were removed from these modules using cis regulatory motif discovery (cisTarget). Finally, enrichment scoring was performed on the target genes (AUCells) of regulatory factors to quantify their activity, and the regulatory factor specificity score (RSS) was calculated. The final result is visualized using the ComplexHeatmap software package [32].

Cell communication analysis
Receptor-ligand mediated cell-cell communication is crucial for coordinating various biological processes such as development, differentiation, and disease. In this study, cell communication analysis was conducted using the iTalk package (Doi: https://doi.org/10.1101/507871) to infer the interactions between different cell subgroups. iTalk treats cell populations as interacting objects, calculates the expression of receptors and ligands in each cell subgroup, and uses this as an indicator of interaction to study the cell communication relationships between subgroups.

Construction of protein-protein interaction (PPI) network
In order to elucidate the targets of PPIF+ neutrophils enriched in the tumor environment promoting CRC progression, the PPI network was used to obtain the protein interaction targets of the receptor-ligand gene list in cell-cell communication based on the STRING database [33](https://string-db.org/). Subsequently, Cytoscape [34] was employed to visualize the PPI network.

Clinical sample
Fresh CRC and surrounding tissue specimens were obtained from patients who had undergone surgery for CRC at the Department of General Surgery at Zhujiang Hospital, Southern Medical University. Tumor samples were obtained from the central region of the tumor, while paired normal samples were collected from macroscopically normal tissue 8–10 cm away from the tumor margin. The Zhujiang Hospital’s Protection of Human Subjects Committee approved the study, and all participants provided informed consent.

Cell lines
Human CRC cell lines, including SW480, LOVO,HCT-116 were purchased from Procell Company (Wuhan, China). The cell line tested negative for mycoplasma contamination. SW480,LOVO and HCT-116 cells were maintained in DMEM medium (C11995500BT, Gibco) supplemented with 1% penicillin-streptomycin (P/S, 15,140,122, Gibco), and 10% FBS (164210-50, Procell Company Wuhan, China) in the incubator with 5% CO2 at 37 °C.

Western blotting
Cells and tissue were lysed in RIPA buffer containing protease and phosphatase inhibitors (P1049, Beyotime). Protein concentration was measured using the BCA protein assay kit (P0010, Beyotime). Equal amounts of protein (10–20 μg) were separated by SDS-PAGE and transferred onto PVDF membranes (IPVH00010, Merck millipore). Membranes were blocked with 5% skimmed milk for 1 h and incubated overnight at 4 °C with primary antibodies against DACH1, NKD1, APC, с-Мус, Wnt-1, Catenin-βand CyclinD1 (1:2,000). After washing, membranes were incubated with HRP- conjugated secondary antibody (1:10,000, SA00001-1, Proteintech) for 1 h at room temperature. Bands were visualized using an ECL reagent (BL520B, Biosharp) and quantified with ImageJ software. Primary antibodies used for Western blotting were listed in Table S3.

Real-time quantitative PCR (RT-qPCR) analysis
Total RNA was extracted using the FastPure Cell/Tissue Total RNA Isolation Kit V2 (Vazyme, RC112) following the manufacturer’s instructions. cDNA was synthesized using the HiScript III RT SuperMix for qPCR (+gDNA wiper) (R323, Vazyme). qPCR was performed using the ChamQ Universal SYBR qPCR Master Mix (Vazyme, Q711) on a QuantStudio6 Real-Time PCR system (Thermo Fisher, Waltham, MA, USA). Each sample was run in three technical replicates. The mRNA expression levels of target genes were normalized to GAPDH and calculated using the ΔΔCt method. The qPCR primers used are listed in Table S4.

Flow cytometry
Collecting blood samples from colorectal cancer patients to extract PBMC (Peripheral Blood Mononuclear Cells) and neutrophils. Fluorescently labeled antibodies specific for macrophages (CD68), Neutrophils(CD66B) were used to stain the cells. After washing and resuspending, Cells were then stained using fluorescently labeled S100A8,Cit-H3 and MPO antibodies. After washing and resuspending again，immune cell subsets were analyzed for their proportions and surface marker expression profiles using flow cytometry. The results were analyzed using a CytoFlow cell analyzer (Beckman). Primary antibodies used for Flow cytometry were listed in Table S3.

Immunofluorescence staining
Clinical samples were fixed for 24 hours, then dehydrated and embedded in paraffin. Sections were cut into 4 μm slices and baked at 60 °C for 1 hour to ensure adhesion. Staining was performed using the Treble-Fluorescence Immunohistochemical Mouse/Rabbit Kit (pH 9.0) (Immunoway Biotechnology, Plano, TX, USA). Following the manufacturer’s protocol, slides were dewaxed, subjected to antigen retrieval, and blocked. The primary antibody (diluted) was applied and incubated overnight at 4 °C in a humidified chamber. After incubation, slides were warmed at 37 °C for 30 minutes, followed by incubation with polymer and fluorescent secondary antibodies. For sequential staining, sections were fully immersed in antibody stripping buffer and microwaved for 15 minutes. After three PBS washes, a second primary antibody was applied, and the staining steps were repeated until all targets were labeled. Primary antibodies used for Immunofluorescence staining were listed in Table S3.

Neutrophil isolation
Neutrophils were isolated using Solarbio Peripheral Blood Neutrophil Separation Kit (P9040). Fresh anticoagulated whole blood (within 2 h of collection) was layered onto a pre-established density gradient (4 mL Reagent A + 2 mL Reagent C) and centrifuged at 1000 g for 30 minutes at room temperature. The granulocyte layer and approximately 1 mL of supernatant separation solution were aspirated, washed with 10 mL cell wash buffer (250 g, 10 min), and treated with red blood cell lysis buffer if erythrocyte contamination occurred. After repeating the wash step (250 g, 10 min), cells were resuspended for further use.

Co-immunoprecipitation (co-lP)
To analyze the endogenous binding of DACH1 and NKD1, HCT-116 cells cultured in 100 × 100 mm dishes were processed according to the rProtein A/G Magnetic IP/Co-IP Kit protocol (ACE Biotechnology, BK0004-01). Cells were lysed with 1000 μL ice-cold 1× Lysis/Wash Buffer (Enhanced) supplemented with protease inhibitor cocktail. Lysates were centrifuged at 13,000 × g for 10 min to collect the supernatant; 20 μL rProtein A/G MagPoly Beads were magnetically washed twice with 180 μL 1× Lysis/Wash Buffer (Enhanced) and then incubated with 5 μg of either anti-DACH1, anti-NKD1, or anti-IgG antibodies in 500 μL buffer at room temperature for 30 min. Subsequently, 500 μL of the HCT-116 cell lysate supernatant was added to the bead-antibody complexes and incubated for another 30 min at room temperature. After incubation, the complexes were washed twice with 500 μL buffer using a magnetic stand. The bound complexes were then eluted and denatured in 50 μL 1× SDS loading buffer by boiling at 100 °C for 10 min. Finally, the eluates were analyzed by Western blot. Primary antibodies used for Co-Immunoprecipitation and Western blotting were listed in Table S3.

Data statistics and analysis
All the bioinformatics analyses and research endeavors encompassed in this study are conducted on a specialized bioinformatics cloud platform (http://www.bioinforcloud.com). All single-cell data analysis was performed using Seurat v4.0. When the p value is less than 0.05, it is deemed to possess remarkable statistical significance. All experiments were performed in triplicate and repeated at least three times independently. For comparisons across multiple groups, statistical significance was assessed using ordinary one-way ANOVA. Two-group comparisons were analyzed by two-tailed Student’s t-test. Data are presented as mean ± SD, with significance denoted as follows: *p < 0.05, **p < 0.005, ***p < 0.0005, ns (not significant) All analyses were conducted using GraphPad Prism 9.0 (GraphPad Software, USA).

Results

The global single-cell microenvironment landscape of CRC
This study included the analysis of 37 CRC samples, 32 adjacent cancer normal samples and 4 PBMC samples from CRC patients, as shown in (Fig. 1A). After rigorous quality control processing, a total of 268,931 high-quality transcripts were obtained, clustered into 89 cell clusters. Based on known cell lineage marker genes, these cell clusters were further annotated into different cell types (Fig. 1B–C), including epithelial cells (EPCAM), endothelial cells (VWF), fibroblasts (COL1A2), macrophages (CD68), neutrophils (FCGR3B), T cells (CD3D, CD3E), B cells (CD79A), and plasma cells (XBP1) (Fig. 1D). Chromosomal heatmaps confirmed the presence of copy number alterations in the identified malignant CRC cells (Fig. 1E). Subsequently, to describe the differences between TME and normal control samples, we compared the cellular ecosystem composition in the three types of samples. Subsequently, to describe the differences between TME and normal control samples, we compared the composition of the cellular ecosystem. Compared with the control group, neutrophils were significantly enriched in CRC samples, and the infiltration of T cells, macrophages, and plasma cells was significantly increased in CRC samples, with the most significant increase in T cell infiltration. T cells dominated the cellular ecosystem in PBMC samples (Fig. 1F).

PPIF+ neutrophil subgroup promotes NETosis through mtROS in TME
The neutrophil population in TME exhibits heterogeneous phenotypes and functional diversity, playing a dual role in promoting or suppressing tumors [35]. Through clustering analysis of 1,799 neutrophils, this study identified 5 neutrophil subgroups (Fig. 2A), with distinct expression of different markers (Fig. 2B). Compared to controls, the PPIF+ Neutrophil subgroup was specifically enriched in CRC, the same conclusion was supported by multiplex immunofluorescence results, while the RPS15A+ Neutrophil subgroup was specifically enriched in PBMC (Figs. 2C-D, S1A). To explore the biological roles of each neutrophil subgroup in the microenvironment, enrichment analysis of marker genes revealed that AKAP9+ neutrophil, DCAH1+ neutrophil, and NR3C2+ neutrophil subgroups exhibited similar expression patterns and enriched signaling pathways related to cell proliferation and angiogenesis (p < 0.05) (Figs. 2E-F). Additionally, PPIF+ Neutrophil and RPS15A+ Neutrophil subgroups also showed similar expression patterns, significantly enriching immune-related biological processes and chemotaxis-related KEGG pathways (p < 0.05) (Figs. 2E-F).
Due to the important role of NETs in promoting tumor progression and metastasis, and their involvement in causing autoimmune inflammation and bystander tissue damage [36, 37], this study also investigated the expression of NETs scores in each subgroup. Remarkably, we found that the PPIF+ Neutrophil subgroup not only had the highest NETs score but also significantly activated the NETs gene set (Figs. 2G-H). Further research revealed that the specific marker of the PPIF+ Neutrophil subgroup, PPIF, encodes a major component of the mitochondrial permeability transition pore (mPTP) [38], and its overexpression may lead to inducing mPTP opening, resulting in increased levels of reactive oxygen species (ROS). The generation of mtROS and the opening of the mPTP both promote the induction of neutrophil extracellular traps formation (NETosis) [39]. Therefore, the high expression of PPIF in the PPIF+ Neutrophil subgroup may promote NETosis by inducing mtROS. Multiplex immunofluorescence analysis revealed significantly elevated expression of neutrophil extracellular trap (NET)-associated biomarkers, including citrullinated histone H3 (Cit-H3) and myeloperoxidase (MPO), in tumor tissues compared to normal controls (Fig. 2I). Meanwhile, we treated peripheral blood neutrophils from tumor patients with the PPIF inhibitor Cyclosporin A and subsequently assessed the proportion of Cit-H3+ and MPO+ neutrophils using flow cytometry. Results showed that PPIF inhibition significantly reduced NETs expression (Fig. 2J). Pseudotime analysis further reconstructed the differentiation trajectory of neutrophils, revealing the continuity of similar expression patterns among subgroups (Fig. 2K). In addition, GRN analysis identified the essential regulators and most specific regulators within each subgroup (Fig. 2L), with the PPIF+ Neutrophil subgroup primarily regulated by IRF9 (Fig. S1B). Survival analysis indicated that a higher gene set score of the PPIF+ Neutrophil subgroup was associated with poor prognosis in CRC (Fig. 2M), suggesting that this subgroup may promote tumor progression and deterioration.
In conclusion, by identifying neutrophil subgroups in the CRC microenvironment, we discovered a subgroup, PPIF+ Neutrophil subgroup with the potential to promote tumorigenesis. The high expression of PPIF in this subgroup may promote NETosis through mtROS, thereby facilitating tumor progression and metastasis. In addition, NETs may trigger pro-inflammatory effects through various pathways such as releasing inflammatory mediators and forming immune complexes, leading to autoimmune inflammation and bystander tissue damage, thereby promoting the progression of CRC.

The DACH1+ NKD1+ CRC subgroups have more significant malignant characteristics
In order to determine the subgroups of malignant cells in CRC, clustering was performed on the corresponding 27,948 cells in the dataset and 9 major CRC subgroups were identified (Fig. 3A). These CRC subgroups expressed specific markers (Fig. 3B), among which, CRC subgroups for SLC26A3+, DACH1+ NKD1+, and CXCL8+ occupied the main components of the CRC cell ecology (Fig. 3C). Subsequently, this study further analyzed the BP and pathways enriched in the CRC subgroups to reveal the expression patterns of specific genes among different subgroups. According to BP, CRC subgroups can be divided into two major categories, one significantly activating functions related to energy metabolism and substance transport, and the other significantly activating functions related to epithelial differentiation and histone modification (p < 0.05), suggesting that these subgroups may have a stronger invasive potential (Fig. 3D). In addition, we noted that the DACH1+ NKD1+ CRC subgroups significantly activated tumor proliferation and metastasis-related signaling pathways such as MAPK, ErbB, and Wnt (p < 0.05) (Fig. 3E), indicating that this subgroup may be more prone to metastasis. By calculating the tumor stemness score (Fig. 3F), the SYME2+ CRC subgroup had the highest tumor stemness score. Therefore, we determined it as the starting point for subgroup differentiation and reconstructed the differentiation trajectory of the CRC subgroup (Fig. 3G).
As key markers of the DACH1+ NKD1+ CRC subgroups, NKD1 was found to be a novel CRC-related gene and its binding to the Wnt signaling body and β-catenin was found to be crucial [40]. In addition, DACH1, as a transcription factor [41], has been confirmed to promote the growth and self-renewal characteristics of CRC cells [42]. These previous findings underscore the possible close relationship between the expression of DACH1 and NKD1, which was also confirmed in this study. In the DACH1+ NKD1+ CRC subgroups, DACH1 and NKD1 were co-expressed and significantly positively correlated (Fig. 3H–I), and the gene set scores from this subgroup also demonstrated a correlation with poor survival in CRC, suggesting that DACH1 may promote the activation of the Wnt signaling pathway by regulating the expression of NKD1, thereby mediating the malignant progression of CRC (Fig. 3J). Multicolor immunofluorescence assays further confirmed the specific expression of DACH1 and NKD1 in tumor tissues (Fig. 3K). Next, the interaction of endogenous DACH1 and NKD1 in HCT-116 colorectal cancer cell lines was confirmed by Co-IP assay (Fig. 3L). To further investigate the relationship among DACH1, NKD1 and the Wnt signaling pathway, we performed Western blot and quantitative polymerase chain reaction (qPCR) analyses on clinical specimens. The experimental results were consistent with our bioinformatics predictions (Figs. 3M,3P, S2A). To investigate DACH1 and NKD1-mediated Wnt signaling pathway activation, we conducted gene perturbation experiments in SW480 and LOVO colorectal cancer cell lines. Following the knockdown and overexpression of DACH1 and NKD1 the expression of Wnt pathway-related molecules (APC, с-Мус, Wnt-1, Catenin-β, CyclinD1) was evaluated. The results of qPCR and Western blotting analyses demonstrated a significant correlation between DACH1and NKD1 expression and the activation of the Wnt signaling pathway (Figs. 3N-O, 3Q-R, S2B-C). In addition, the CRC subgroups for IGKC+ and CXCL8+ were also associated with poorer OS in CRC (Figure S2D). Subsequently, we reconstructed the GRN of the CRC subgroups, clustering cell type-specific regulators into 3 major modules and predicting the essential regulators in each CRC subgroup (Figure S2E). By ranking the RSS, we further identified the most specific regulators for each subgroup (Figure S2F).
In conclusion, we identified CRC subgroups, identified the subgroups with higher malignancy, and provided specific regulators for regulating different CRC subgroups, which will provide new insights into strategies for inhibiting tumor metastasis and proliferation.

C1QC+ Mac subgroup induces M2 polarization to promote tumor cell growth
High density tumor associated macrophages play an important role in malignant progression [43]. Through cluster analysis of macrophages, 9 different subgroups were identified, each expressing different specific markers (Figs. 4A-B). Compared with the control group, the CIQC+ Mac subgroup and AQP9+ Mac subgroup were specifically enriched in CRC, while the S100A8+ Mac subgroup was specifically enriched in PBMC (Fig. 4C). Notably, the C1QC+ Mac subgroup and the AQP9+ Mac subgroup are specifically enriched in CRC. Multiplex immunofluorescence results were consistent with the cellular ecological analysis (Fig. 4D-F). Flow cytometry analysis demonstrated a significant increase in the proportion of S100A8+ macrophages in the CRC patients compared to the normal controls (Fig. 4G). To explore the potential role of macrophage heterogeneity in regulating the tumor immune microenvironment, enrichment analysis of subgroup marker genes was performed. The enrichment analysis shows that the C1QC+ Mac subgroup is enriched in negative regulation of immune effector processes, cell chemotaxis, and biological functions of the ECM tissue, as well as chemokine signaling pathways, cell adhesion molecules, and NOD−like receptor signaling pathway (p < 0.05) (Fig. 4H–I). This study also explored the expression of M1 and M2 scores in each subgroup. Importantly, we found that the C1QC+ Mac subgroup had stronger M2 signal expression compared to M1 signal, and the angiogenic signal was also stronger (Fig. 4J), suggesting that this subgroup may be in a polarized state between M1 and M2, potentially promoting the growth and invasion of the tumor. Proposed pseudotime analysis indicates that the developmental trajectory of Mac evolves from the NOP53+ Mac subgroup as a starting point to other subgroups (Fig. 4K). Survival analysis confirms that lower abundance of the NOP53+ Mac subgroup and the RACK1+ Mac subgroup, and higher abundance of the IFI27+ CRC subgroup, are associated with worse OS and RFS in the TCGA clinical cohort of CRC patients (Fig. 4L). Furthermore, we also constructed a GRN for Mac and analyzed and identified the necessary and most specific regulators in each subgroup, including transcription factors (TFs), such as TCF7 and NR2F1 (Figure S3A).
In summary, we identified that the C1QC+ Mac subgroup is enriched in negative regulation of immune effector processes, cell chemotaxis, and biological functions of the ECM tissue, as well as chemokine signaling pathways, cell adhesion molecules, and NOD−like receptor signaling pathway (p < 0.05), which may weaken the host immune response, affect the migration and invasion capabilities of tumor cells. It is worth noting that it exhibits M2 polarization, and through M2 polarization of the C1QC+ Mac subgroup, it may secrete multiple growth factors and cytokines to promote the growth of tumor cells, thereby promoting the progression of CRC.

Oxidative phosphorylation of RACK1+ tem and Tcm promotes tumor immune escape
In order to identify the functional states of different T cell subgroups, we conducted a deep analysis of T cells using unsupervised clustering based on UMAP, which identified 10 T cell subgroups (Fig. 5A). To assess the functional states of different T cell subtypes, specific markers were identified for 6 subgroups: CD8+ T cells, natural killer T cells (NKT), central memory T cells (Tcm), effector memory T cells (Tem), T helper 17 cells (Th17), regulatory T cells (Treg), and naive T cells (Fig. 5B), each with distinct expression of specific markers (Fig. 5C). Notably, the RACK1+ Tem、RACK1+ Tcm and JUND+ Treg subgroups showed a significantly higher specific enrichment in CRC compared to other subgroups (Fig. 5D). Multiplex immunofluorescence staining of clinical samples from colorectal cancer patients also confirmed the RACK1+ Tem cells were significant enrichment in tumor tissues (Figs. 5E-F). This subgroup may be closely related to the occurrence, development, or prognosis of CRC. Enrichment analysis showed that the RACK1+ Tem and Tcm subgroups not only enriched biological functions such as protein targeting, cell trend, and T cell activation, but also enriched oxidative phosphorylation signaling pathways (p < 0.05) (Figs. 5G–H), suggesting that the RACK1+ Tem and Tcm subgroups may interact with immune cells and exert immune efficacy, but their participation in oxidative phosphorylation, thereby affecting the metabolism and function of subgroups themselves, leading to weakened immune killing function. To better understand the dynamic transformation of T cells in the biological timeline, a pseudotime analysis was performed to reconstruct the differentiation trajectory of T cell subgroups, revealing that the ZBTB20+ Naive subgroup serves as the starting point for the differentiation of other subgroups (Fig. 5I). In addition, GRN analysis showed that these cell subgroups are divided into four modules regulated by transcription factors such as NKX3-1, MTA3, and IRF9 (Figs. 5J-K).
These results indicate that the infiltration of RACK1+ Tem and Tcm subgroups in CRC is increased, and both significantly activate the oxidative phosphorylation pathway (p < 0.05), indicating that this subgroup undergoes impaired metabolism and function, and immune function is suppressed, leading to immune escape.

Lineage differentiation of B-plasma cell subgroups in the CRC TME
For a long time, B cells that secrete antibodies have been considered to be the core element of intestinal homeostasis, and plasma cells and B lymphocyte differentiation maintain antibody production. However, tumor-associated B cells in CRC have not been well characterized [44, 45]. Therefore, we conducted a joint study on B cells and plasma cells, identifying 10 B subgroups and 7 plasma subgroups, respectively. Specifically expressing MS4A1, CD22, and CD19 were identified as B cells, while those specifically expressing TNFRSF17, CD38, JCHAIN, and SDC1 were identified as plasma cells (Figs. 6A-B). Each subgroup specifically expressed different biomarkers. Among them, CELF2+ positive and IGHG2+ positive were highly expressed in the B cell subgroup and plasma cell subgroup, respectively (Figs. 6C-D).
It is worth noting that the RACK1+ B subgroup and RACK1+ Plasma subgroup are specifically enriched in CRC (Fig. 6E). Multiplex immunofluorescence results showed higher RACK1 expression in B cells and plasma cells within tumor tissues compared to the normal controls (Figs. 6F-G, S4A-B). Enrichment analysis showed that the RACK1+ B subgroup was enriched in pathways such as protein processing and focal adhesion in the endoplasmic reticulum (p < 0.05) (Fig. 6H). The RACK1+ Plasma subgroup is rich in the RIG−I−like receptor signaling pathway (p < 0.05) (Fig. 6I). The pseudotime analysis shows that the developmental trajectory of B cells evolves from the CELF2+ B subgroup as a starting point to other subgroups, and the plasma’s trajectory evolves from the ARID1B+ Plasma subgroup as a starting point to other subgroups (Figs. 6J–K). We reconstructed the gene regulatory network of B cell and plasma cell subgroups, these cell subgroups are divided into two modules, regulated by transcription factors such as TGIF1 and FOXO1, and specific scores of the regulatory elements corresponding to the subgroups were calculated and ranked (Figures S4C-D).
In summary, we identified the B subgroup and Plasma subgroup, finding that they both exhibit specific high expression of RACK1+ in the subgroup and may promote the progression and metastasis of tumors by affecting signaling pathways and immune response processes.

PPIF+ neutrophils regulate TME through intercellular communication to support CRC progression
To evaluate the impact of PPIF+ neutrophil enrichment in the tumor environment on cancer progression, we explored the mechanisms of global cell-cell interactions in CRC, including immune checkpoints, cytokines, growth factors, and other factors. In the context of immune checkpoint, neutrophils are closely associated with Mac, T cells, B cells, and plasma cells through the TNFSF14–TNFRSF14 axis. Interaction between TNFSF14-LTBR axis and CRC cells. (Figure S5A). TNFSF14, TNFRSF14 and LTBR are expressed in different cell subgroups (Fig. 7A). Furthermore, in terms of cytokines, we found that PPIF+ neutrophils are closely linked to NKT/macrophages through the CXCL8–CXCR2 axis (Figure S5B). Among growth factors, we observed that TGFB1 is present in Mac, NKT, Tem, Treg, plasma, and B-cell subgroups, and communication between neutrophils and these cells is achieved through the VEGFA-ITGB1 axis (Figure S5C). In the extracellular matrix, we found that the PPIF+ neutrophil subgroup is mainly closely related to other cell subgroups through B2M (Figure S5D).
Subsequently, we obtained the protein - protein interaction targets of TNFSF14, TNFRSF14, and LTBR based on the STRING database, and performed enrichment analysis based on these target genes. We found that they are significantly enriched in the TGF-β signaling pathway and the colorectal cancer signaling pathway (Figs. 7B–D). Studies have shown that the TGF-β signaling pathway plays important roles in many biological processes, including cell growth, differentiation, apoptosis, migration, and cancer occurrence and progression [46]. To support our predictions regarding PPIF+ neutrophils, we also investigated the regulatory interactions inferred by GRN (Fig. 7E).
In conclusion, our analysis preliminarily indicates that neutrophils may regulate the functions of other cells through various ligand - receptor interactions, regulate other immune cell subtypes, suppress the immune system, promote immune evasion of tumor cells after metastasis, and reprogram the pre - metastatic microenvironment.

Discussion
This study analyzed CRC scRNA-seq data from GEO database to present the spatial transcription and regulation of CRC at single-cell resolution, revealing the mechanism by which PPIF+ neutrophil subgroups promote mtROS-driven NETosis. We further investigated six cell types closely associated with CRC, including neutrophils, CRC cells, Mac, T cells, B cells, and plasma cells, as well as the interactions between neutrophils within the TME and cells affected by NETs.
Neutrophils are crucial innate immune cells that execute a range of defense mechanisms [47]—including phagocytosis, degranulation, and the formation of neutrophil extracellular traps (NETs). Increasing evidence suggests that, beyond their antimicrobial roles, NETs contribute to the pathogenesis of various non-infectious conditions, including inflammation, allergic disease [48], autoimmune disorders, and cancer. In the context of tumors, Tumor-associated neutrophils can also form and release NETs into TME [49, 50], potentially influencing disease progression. In this study, we identified a specific PPIF+ neutrophil subgroup enriched in CRC tissues, which may act as a key driver of NETosis via mtROS production. These NETs potentially contribute to tumor progression by recruiting and activating other cell types within the tumor microenvironment. Inhibiting PPIF expression or its function in NET formation could offer a novel therapeutic strategy in CRC [47].
Our research has confirmed the close intercellular communication between neutrophils and cells affected by NETs in the TME. Specifically, we found that neutrophils are closely related to CRC cells through the TNFSF14-LTBR axis. Notably, CRC malignant cells co-expressing DACH1 and NKD1 display more pronounced malignant characteristics. Additionally, this subgroup was found to be enriched for Wnt signaling(p < 0.05) and exhibited more aggressive malignant features [51]. Furthermore, neutrophils interact with immune cells through the TNFSF14–TNFRSF14 axis. NETs can activate macrophages, leading to alterations in their phenotypes and functions [52]. Macrophages can polarize into two functionally distinct phenotypes: M1 and M2 [53]. We identified an increased abundance of C1QC+ macrophages, which appear to be transitioning from M1 to M2 states, suggesting an immunosuppressive shift [54]. Interestingly, in the CRC group, the abundances of the RACK1+ subgroups in Tem, B, and plasma cells are significantly higher compared with other subgroups. Notably, both the RACK1+ Tem and Tcm subgroups significantly activated the oxidative phosphorylation pathway (p < 0.05), which may contribute to impaired functions, thereby promoting the immune escape of tumor cells. However, this association remains hypothetical. Further mechanistic studies are needed to validate these findings.
Our findings highlight the activation of multiple keys signaling pathways in the CRC microenvironment, including the Wnt signaling, TGF-β signaling, and TNFSF14-related immune checkpoint pathways. These pathways are tightly linked to tumor proliferation, immune evasion, and metastasis. The TNFSF14-LTBR axis emerges as a novel immunotherapeutic target in CRC [55]. Although this pathway has been studied in other cancers, its specific role in the CRC tumor microenvironment remains poorly characterized. Given its involvement in tumor-immune interactions, further exploration of this pathway may offer promising avenues for the development of targeted immunotherapies.
This study explored the global single-cell landscape of CRC by analyzing scRNA-seq data. Despite revealing the cellular interactions and regulatory mechanisms of CRC at the single-cell level, several limitations of this study should be acknowledged. First, the relatively small sample size may restrict the representativeness and statistical power of the findings. Second, several conclusions are based on computational analyses, lacking direct experimental confirmation. In vivo validation remains particularly challenging. In our preliminary experiments, we were unable to identify a suitable in vitro model to support this study, largely because neutrophil isolation still relies on primary patient samples, and allogeneic neutrophils are difficult to maintain in a quiescent, non-activated state. Therefore, humanized mouse models with functional neutrophil systems are expected to serve as indispensable tools for future studies. As an initial step, we plan to generate genetically engineered mice targeting PPIF, thereby providing a platform for in-depth investigation of the PPIF+ neutrophil subpopulation.
Future studies will aim to expanded the sample size to control heterogeneity-related bias, following methodologies established in cardiovascular biomarker research [56]. We also plan to investigate the distribution of PPIF+ neutrophils across clinical subtypes, disease stages, and tumor localization. Meanwhile, functional experiments will be conducted to validate the proposed mechanisms and therapeutic targets, providing a solid foundation for precision therapy in colorectal cancer.

Conclusions
In summary, our study revealed the characteristics and heterogeneity of CRC in the immune microenvironment, providing a new perspective for a deeper understanding of the immunotherapy of CRC. At the same time, it also laid the foundation for further exploration of its biological mechanisms and clinical applications in related fields.

Electronic supplementary material
Below is the link to the electronic supplementary material.