Corrigendum: Integrated temporal transcriptional and epigenetic single-cell analysis reveals the intrarenal immune characteristics in an early-stage model of IgA nephropathy during its acute injury

Chen Xu ¹ Yiwei Zhang ¹ Jian Zhou ¹ Jiangnan Zhang ² Hui Dong ¹ Xiangmei Chen ³ Yi Tian ^1* Yuzhang Wu ^1*

• ¹ Institute of Immunology, Third Military Medical University, Chongqing, China
• ² Xinqiao Hospital, Shapingba, Chongqing, China
• ³ People's Liberation Army General Hospital, Beijing, Beijing Municipality, China

1IntroductionIgAN is the most common primary glomerular disease worldwide characterized by the deposition of IgA-containing immune complexes in the mesangium, followed by mesangial hypercellularity and matrix expansion(1-4). Its glomerular histopathology also includes infiltration of inflammatory cells, proliferation of endocapillary, and formation of crescents. IgAN patients show great heterogeneity in their clinical manifestation, ranging from asymptomatic microscopic hematuria to rapidly progressing glomerulonephritis(3, 4). At present, effective targeted therapies are very limited for this tissue-specific autoimmune disease, and 30% to 40% of cases progress to uremia within 20 to 30 years(5).The infiltration of macrophages and T cells in the glomerulus and/or interstitial compartments has been corroborated to play an initiative role in regulating renal inflammation and contribute to IgAN progression(6, 7). Recently, scRNAseq analyses have unveiled the complexity of cellular phenotypes in both normal and injured human and mouse kidneys, thereby revolutionizing our understanding of kidney disease pathologies from diverse perspectives(8-19). However, the phenotypes of specific immune cell types in the IgAN kidney are not well understood. Integrated analysis of multimodal single-cell omics data is an emerging tool to describe an unbiased and comprehensive view of cell atlas. For example, Giles et al demonstrated that scRNAseq data obtained from the LCMV model of CD8+ T cell differentiation had less resolution in defining cell identity whereas the paired scATACseq data outperformed in determining cell ‘fates’(20). Joint profiling by longitudinal scRNAseq and scATACseq can provide a framework for understanding how development of cell subset diversity as well as how chromatin accessibility regulates transcription(16, 21). And immune cells in different pathology contexts may have distinct chromatin accessibility profiles that change as they differentiate. In the current study we used an early IgAN model with an extra LPS injection to mimic disease exacerbation resulting from mucosal infection and/or a reversible form of AKI, which is frequently observed in clinical practice among IgAN patients(5). LPS, a classical endotoxin derived from the outer membrane of Gram-negative bacteria, serves as an indispensable virulence factor that triggers Toll-like receptor-4 (TLR-4) inflammatory signaling, eliciting potent innate and adaptive immune responses(22). Therefore, we designated the established IgAN model with and without systemic administration of LPS as “immune activation (IA)/ immune remission (IR)” and “immune stability (IS)”, respectively (Methods). Accordingly, we generated temporal scRNAseq data and scATACseq data for CD45-positive immune cells isolated from kidneys in different immune responses to delineate population heterogeneity and identify gene expression as well as accessible chromatin patterns associated with major branches of immune cell phenotype. To date, this study represents the first multiomics analysis that characterizes the dynamic immune cellular landscape at a high resolution in an early IgAN mouse model.2Results 2.1 A murine model with LPS injection to phenocopy 3 states in early human IgAN-associated glomerulopathyEstablished IgAN model in this study showed prominent glomerular immune deposits (Fig. 1A) with histopathologic changes (Fig. 1B) as observed in early IgAN patients. Despite the LPS injection in the IA or IR group not inducing more severe detectable glomerular abnormalities, physiological measurements confirmed that the IA group (short for IA) exhibited more severe proteinuria, compared to the other groups (Fig. 1C). Therefore, these three phases simulate to a certain extent the onset and resolution of self-limited AKI in early IgAN.2.2 CD45+ cell transcriptional atlas in mouse kidneysDissociated kidney cells from CTRL, IS, IA and IR (4 or 5 mice each) were processed and enriched for CD45-positive immune cells (Methods and Fig. S1). Isolated CD45+ cells were then pooled into a single sample per experimental condition for scRNAseq and scATACseq analysis respectively using the 10× Genomics platform. After quality control and filtering (Methods and Fig. S2), we obtained 72304 CD45+ single-cell transcriptomes from vehicle and IgAN mice.Seven major cell types were identified on the basis of the canonical marker expressions as B cells, plasma cells, T and NK cells (TandNK), neutrophils, basophils, mononuclear phagocytes (MPs), and plasmacytoid dendritic cells (pDCs) (Methods and Fig. 1D). MPs could be further divided into macrophages, monocytes, conventional type 1 dendritic cells (cDC1), and conventional type 2 dendritic cells (cDC2) (Fig. S4). The proportion of immune cell types showed difference in disease phases, and we detected a prominent dynamic change of the myeloid cells and the lymphoid cells (Fig. 1E). To investigate the distinctive immune profiles of different disease phases, we performed Gene Ontology (GO) pathway analysis of DEGs for IA, IR and IS. We found increased expression of genes that regulate leukocyte migration and chemotaxis, especially mononuclear cells and granulocytes in IA compared to IR (Fig. 1F), indicating that circulatory myeloid activation altered the kidney myeloid compartments during the early onset of AKI. While in IR and IS, DEGs were more enriched in leukocyte differentiation, cell-cell adhesion and T cell activation than in IA (Fig. 1G and 1I). When compared with IS, IA was characterized by elevated levels of cellular metabolism, including generation of metabolites and cellular respiration (Fig. 1H), consistent with the activated state of immune cells. 2.3 Temporal scRNAseq reveals transcriptional heterogeneity in macrophage subsetsTo examine the kidney macrophage heterogeneity, we performed the unsupervised clustering of macrophages from vehicle and IgAN immune cells, which identified 8 subclusters, including infiltrating macrophages (infiltrating Mac, highly expressing infiltrating macrophages genes e.g., Plac8, Msrb1, Anxa2, Gngt2, Ear2), resident Mac, inflammatory Mac (highly expressing inflammatory macrophages genes e.g., Ccl3, Ccl4, Il1b), 2 subsets of interferon (IFN) gene signature high Mac (IFN Mac), proliferating Mac (highly expressing proliferating macrophages genes e.g., Stmn1, Hist1h1b, Hist1h2ae, Birc5, and Mki67), NLR family, pyrin domain containing 1B (Nlrp1b)-high expressing Mac (Nlrp1b Mac) and Mediterranean fever (Mefv) -high expressing Mac (Mefv Mac) (Fig. 2A-B) (23). The annotations were identified based on previously defined marker genes by Fu et al.(23). The cell types in the current study were largely similar with that in the mouse model of early diabetic kidney disease (DKD), implicating that IgAN and DKD might share conserved phenotypic spectrum in local macrophage transcriptome, despite of different mouse model from different background(23). Mannose receptor C-type 1 (Mrc1)-high expressing (Mrc1hi) Mac, one of resident macrophage subsets as Fu et al. have reported, showed relatively higher expression of resident marker genes (e.g., C1q, Cd81 and Mgl2) than other subclusters in our study, thus we annotated them as resident Mac (Fig. 2B). Although we did not detect any populations with a relatively high expression of Trem2 (Fig. S4), we distinguished 2 subsets Nlrp1b Mac and Mefv Mac that have not been previously described and might be specific to IgAN. And we observed notable changes in the proportions of Nlrp1b Mac and Mefv Mac in IA (Fig. 2C), suggesting their potential role in the acute phase of IgAN or kidney injury. Besides, we split IFN Mac into 2 subsets as it clustered “apart” in the scRNAseq UMAP space and they displayed different IFN-stimulated gene (ISG) expression patterns that IFN Mac Cxcl9 with a high expression of Gbp2, Gbp2b and Cxcl9, while IFN Mac Ifit3 with a high expression of Isg15, Ccl12, Ifit3 and Ifit2 (Fig. 2A-B). We first examined the macrophage subsets for the expression of canonical M1 markers and M2 markers, which showed an increasing trend for both M1- and M2-like macrophage subtypes, rather than having discrete M1 or M2 phenotypes (Fig. S5). Recent evidence from scRNAseq analysis has pointed to a more dynamic and continuous spectrum of macrophage polarization phenotypes of tissue macrophages under various conditions(23, 24).To understand the potential developmental transitions of macrophage clusters, we applied RNA velocity analysis to construct the developmental trajectories of 8 macrophage clusters. Two major trajectories were observed that started from infiltrating Mac, then bifurcated and ended up as Nlrp1b Mac and resident Mac respectively, demonstrating that there might be a constant flow of infiltrated macrophages entering these clusters and terminally differentiating as them (Fig. 2D). In line with previous publications, resident macrophages or other subclusters are both locally proliferating and partially replenished by the circulation especially when space in the resident macrophage pool is created under inflammatory conditions(24, 25). Figure 2D shows examples of genes that are highly expressed along the trajectories (e.g., pro-inflammatory gene S100a6 and ISG ifi205 in infiltrating Mac, complement activation gene Cd5l and phagocytosis-associated gene Myo1g in Mefv Mac).By examining signature genes defined by prior data (Methods and Table S1), we observed distinct functional status for each macrophage subset as follows: infiltrating Mac with the highest lipid metabolism score and a high monocyte-derived macrophage (MoMF) score (Fig. 2E). And GO analysis revealed that many pathways specifically enriched in infiltrating Mac were related to energy metabolic process and respiration (Fig. 2F). Macrophages augmenting lipogenesis and the utilization of glucose upon the stimulation of TLR4 by LPS has been linked to enhanced phagocytosis and cytokine production(26-28). We also observed high pro-inflammatory scores in resident Mac, inflammatory Mac and IFN Mac Ifit3, the highest resident score for resident Mac, as well as the highest interferon-responsed score for IFN Mac (Fig. 2E). Consistent with its pro-inflammatory potential, IFN Mac Ifit3 DEGs were enriched in pathways related to defense response to virus and cytokine production, whereas IFN Mac Cxcl9 in pathways associated with nucleotide metabolic process (Fig. 2F). Next, we asked whether IFN Mac Cxcl9 had anti-inflammatory potential by examining the Cd274 (encoding PD-L1) expression of macrophages. IFN Mac Cxcl9 was found to show a marked expression (Fig. 2G). Together, these results further confirmed the bifurcation of IFN Mac. Of interest, research has shown that tumor T cells cocultured with the PD-L1+ macrophages exhibited an impaired production of IFN-γ(29). In other words, IFN Mac Cxcl9 might inhibit T cell activity and suppress the overexpression of IFN. Furthermore, IFN Mac Cxcl9 displayed a relatively high score of immune regulation and high expression of immuneregulatory gene Sod1 (Fig. 2E and S5). Altogether, these observations indicated a key role for this subset in the ability to attenuate and resolve the inflammation. There was also a proliferating cluster characterized by the expression of genes consistent with DNA replication and repair, as well as cell division, such as Hist1h1b, Hist1h2ae, Birc5, and Mki67 (Fig. 2B and 2F). Because cell cycle genes can obscure underlying transcriptional identity, we projected these cells back onto the remaining clusters as described in previous publications (Methods)(20). Most proliferating cells belonged to the IA Mefv cluster, with part of cells derived from inflammatory Mac and Nlrp1b Mac present at IS (Fig. 2H). Thus, these clusters shared transcriptional features of proliferative activity that might drive colocalization in scRNAseq space. GO pathway analysis of Nlrp1b Mac DEGs disclosed their various regulatory effects on immune response and apoptotic signaling pathway, together with the upregulation of genes linked to leukocyte migration, chemotaxis and cell-cell adhesion, whereas Mefv Mac DEGs were positively correlated with the regulation of leukocyte differentiation and GTPase activity (Fig. S6). However, these “regulatory-alike” clusters did not belong to the canonical immunoregulatory cells due to their low immune-regulatory score as shown in Figure S6. They did not acquire anti-inflammatory functions that contributed to tissue fibrosis, angiogenesis or inhibition of T cell responses either (Fig. S6). Further looking at other immunoregulatory genes defined by Zhang et al revealed a significant upregulation of Lyn, Mtss1 and Pecam1(30). These analyses unveiled a previously unappreciated macrophage population with potent immunoregulatory effects. To gain more molecular insights into these clusters, next we used gene activity, a metric of local gene accessibility, to approximate gene expression and avoid the interferences caused by overexpression of mitochondrial genes or ribosomal genes. First we compared Nlrp1b Mac and Mefv Mac to resident Mac, the representative populations of the two trajectories. Nlrp1b Mac and Mefv Mac were distinguished by another immunoregulatory gene Ptprd, genes associated with cell adhesion (Nrxn3, Tenm2 and Grid2) and transcription factor (TF) ESRRG (Fig. S7). Of note, the largest number of genes with differential activity found in this comparison were down-regulated genes of Nlrp1b Mac and Mefv Mac, including Rab7b (expression involved in response to IFN), Cd52 (expression involved in response to bacterium), Cd83 (expression involved in positive regulation of CD4+ T cell differentiation), C1qc and H2-Eb1 (expression involved in phagocytosis and antigen presentation), and migration-related genes Ccl3 and Ccl4 (Fig. S7). Then based on comparison with each other, Nlrp1b Mac had higher expression of Ptprd, whereas more immunoregulatory genes Tnfrsf1b and Pik3ap1 were higher in Mefv Mac (Fig. S7). These data further confirmed the distinct immunoregulatory effects of Nlrp1b Mac and Mefv Mac.In summary, these non-canonical immunoregulatory cells with high expression of multiple related genes and various biological effects did not demonstrate a typical pro- or anti-inflammatory phenotype. Next, we performed ligand-receptor (L-R) analysis to explore molecular crosstalks between these cells and the rest macrophages or CD4+ T cells. We found that APP_CD74 showed the highest interaction potential (Fig. S8). Amyloid beta precursor protein, encoded by the App gene, is involved in negative regulation of blood circulation; positive regulation of endothelin production; and positive regulation of tumor necrosis factor (TNF) production, and CD74 is mainly involved in macrophage migration inhibitory factor signaling pathway. In addition, it has been recently reported that the APP_CD74 axis contributes to Treg-exhaustion CD8+ T cells (Tex) interaction in HBV-infected patients(30). Together, these data may provide valuable information for functional verification of these immunoregulatory cells in future studies regarding immune pathogenesis as well as therapeutic attempts for IgAN.2.4 Dynamic variances in neutrophil profile across different phasesNeutrophils were further divided into Ccrl2, Csf3r, Ly6g and Mmp8 subclusters based on their high expression of the marker gene Ccrl2, Csf3r, Ly6g and Mmp8 respectively(Fig. 3A). Supplementary Figure S9-10 show enrichment scores of neutrophil subsets for specific gene signatures and GO pathways that are specifically enriched in each neutrophil subset. Ccrl2 neutrophils were characterized by high pro-inflammatory score and involved in multiple signaling pathways such as NF-κB signaling, pattern recognition receptor (PRR) signaling pathway, and extrinsic apoptotic signaling pathway (Fig. S9-10). Ly6g neutrophils with a high ISG score and the highest pre-neutrophil score (Fig. S9), serving as the starting root of the pseudo-time trajectories (Fig. 3B), exhibited characteristics indicative of both early differentiated and viral response. Mmp8 neutrophils were distinguished by chemotaxis and migration-related genes and had more pro-angiogenic potential (Fig. S9-10).The proportion of neutrophil cell types was largely similar between CTRL and IS (Fig. 3C). However, in IA there was skewing away from Csf3r neutrophils and toward the Ly6g subset (Fig. 3C). While in IR, there was a substantial increase in Csf3r neutrophils with concomitant loss of Ly6g subset (Fig. 3C), suggesting the specific neutrophil response might contribute to control the acute exacerbation of IgAN. Based on gene activity, neutrophils in IA highly expressed Il10, Adora2a, Tnip1, Pou2f2, Clec4a1, Lta, etc., genes mainly involved in positive regulation of inflammatory response through multiple signaling pathway (Fig. 3D), including MYD88-dependent TLR signaling pathway, NF-κB signaling, G protein-coupled receptor signaling pathway, IL17-mediated signaling pathway, type I IFN-mediated signaling pathway. While looking at neutrophils from IR, representative genes Il1r2, Il15, Zfp36l2 (negative regulation of cell differentiation) and Wbp1l (CXCR4 signaling pathway) were relatively up-regulated (Fig. 3D). Taken together, these findings highlighted the potential roles of neutrophils in the pathogenesis of IgAN that have rarely been explored.2.5 The integrated immune landscape distinguishes Nlrp1b Mac as a distinct macrophage subsetDistinct cell type is the result of an epigenetically distinct developmental path driven in part by specific TF(31). Next, we proceeded to explore the epigenetic and transcriptomic landscape of immune cells to identify the regulatory elements that define the biological states among the four groups based on parallel scRNAseq and scATACseq on the same samples. As shown in the UMAP, following the integration of scATAC- and scRNA-seq datasets, 2 major clusters and 7 celltypes were identified (Fig. 4A and S11), validating that scATACseq reveals fewer cell “fates” underlying multiple transcriptional states(20). Unexpectedly, Nlrp1b Mac stood out as a distinct population outside of other clusters (Fig. 4A). Given the epigenetic divergence of Nlrp1b Mac versus other macrophage clusters, we next compared chromatin accessibility changes between them. Among regions with increased or lost accessibility, only a finite number of them were shared (Fig.S11). Next, we investigated epigenetic programs used by different macrophage subsets as described by Giles previously(20). We visualized all differentially accessible chromatin regions (DACRs) (Fig. 4B), then assessed the number of DACRs in each gene locus (Fig. 4C). This approach revealed four representative global patterns of distinct ACRs among macrophage clusters, namely Nlrp1b Mac, infiltrating Mac, inflammatory Mac and resident Mac (Fig. 4B-C). Accordingly, Nlrp1b Mac had a unique ACR profile and had the most accessible DACRs among macrophage clusters (Fig. 4B-C).Nlrp1b Mac have a unique TF motif profile characterized by enrichment in ESRRG, NR4A2, RBPJL, WT1, ZIC1::ZIC2, etc. (Fig. 4D and S12). In this case, we observed a specific TF regulation activity pattern in these clusters with low deviation score of NFE2L2 in Nlrp1b Mac and proliferating Mac, and high score in IFN Mac Cxcl9, infiltrating Mac and resident Mac (Fig. 4E). Since most DEGs of Nlrp1b Mac were shared with other Mac (Fig. S13), these underscore the potential of chromatin accessibility to provide additional information beyond transcriptional data(20). 2.6 Tex/CX3CR1-transTeff is the major CD8+ effector during the acute phaseFurther clustering of scRNAseq T and NK cells yielded 12 clusters, mainly including type 2 innate lymphoid cells (ILC2), γδ T cells (GDT), naïve T cells, CD4+ T-follicular helper cells (Tfh), regulatory CD4+ T cells (Treg), effector memory CD8+ T cells (Tem), 2 effector CD8+ T cells (Teff_Jun and Teff_Stat4 with their high expression of the marker gene Jun and Stat4 respectively) and exhaustion CD8+ T cells (Tex) (Fig. 5A-C). These clusters were similar to those resolved and annotated from the integrated data (Fig. S14). This further validated the accuracy of the cell type annotation, and demonstrated that unlike the Mac subpopulation, which has a common lineage "fate", the heterogeneity of T cell subpopulations is relatively stable and conserved from both transcriptomic and epigenetic perspective. Therefore, subsequent analyses are mostly based on the integrated scRNAseq and scATACseq data to reduce potential errors and improve the reliability of the results.Although Tex herein displayed a typical exhausted state, which was featured by high expression of inhibitory markers (Tox, Havcr2, Lag3, Pdcd1 and Tigit) and low expression of effector markers (Tbx21 and Ifng), their high expression of Gzmk and intermediate expression of Lag3 among other genes could term them in a different way as predysfunctional cells identified in a melanoma cohort(32) (Fig. 5C). Distinguished from the “authentic” progenitor Tex or Tex precursor (Tpex) with naïve-associated genes (Lef1, Sell, Ccr7 and Il7r) expression, Tex or predysfunctional cells herein with expression of Teff genes Tbx21 and Zeb2 more resemble transitory Teff (transTeff) (Tpex→transTeff→termTex), the latter were previously described as Tpex-derived effector-like CX3CR1+ T cells during mouse lymphocytic choriomeningitis virus (LCMV) infection(20, 33, 34) (Fig. 5C and S15). Despite transTeff in our study with loss of CX3CR1 (Fig. S15), their effector pattern was defined by high cytotoxic score and effector score as shown in figure 5D. Given the observed increase in the proportion and exhaustion score of Tex in IR (Fig. 5B and 5D), we have deduced that CX3CR1-transTeff represents the predominant CD8+ effector population during IA, which subsequently evolves towards termTex from IA to IR. We next asked how key changes in chromatin accessibility identified by scATACseq were associated with developmental trajectories. We identified distinct epigenetic patterns associated with expression of key genes (Gzmk, Pdcd1 and Tigit) across time points using scATACseq (Fig. 5E). The increased accessibility of the Gzmk locus in IA and the Pdcd1 locus in IR was observed as anticipated (Fig. 5E). Interestingly, we observed increased accessibility at the promoter and enhancers regions of Gzmk in IS (Fig. 5E). Lastly, we noted a unique feature of heat-shock protein (HSP) genes (Hspa8, Hspa1b, Hspa1a, Dnajb1, Hspd1, Hspe1) overexpression in Teff_Jun (Fig. 5F and S14). A CD8+ T cell population with HSP genes upregulated was previously reported in Tex isolated from the mouse spleen during LCMV infection, but their biological function remained unelucidated(20). Further studies are needed to clarify the significance of Teff_Jun in IgAN.Similarly, to observe distinct functional status for each T-cell and NK-cell subset, we examined signature genes defined previously on the integrated scRNAseq data (Methods and Table S1). As expected, the highest regulatory score for Treg and the highest exhaustion score for Tex confirmed their transcriptomic signature (Fig. 5D). We also observed the highest cytotoxic score for Tem (Fig. 5D). As for NK subsets, NK_Gzma (with high expression of marker gene Gzma) had a significantly higher cytotoxic score, while NK_Xcl (with high expression of marker gene Xcl) had a relatively higher expression of general stress-related genes, a pattern reminiscent of CD8+ Tem and Teff_Jun (Fig. 5G and S14). Tem and NK_Gzma might belong to one epigenetic group that featured by accessibility at As3mt, Dagla, NK-associated gene Fcgr2b and cytotoxic gene Gzma (Fig. 5F). Similar to recent work, our analyses also identified specific T-cell cluster expressing genes associated with NK cells (Klr genes or Fcgr2b, for example)(20, 35, 36) (Fig. S1). 2.7 Cell-cell interaction between immune cells and endothelial cells To dissect the potential crosstalk of immune cells and renal resident cells, we performed scRNAseq on kidney cells obtained from an independent cohort of IS and IA mice without immune-cell isolation (Fig. 6A). Several clusters or subsets were undetected because of their relatively low overall proportions (Fig. 6A). Then we used Cellphone DB to infer cell–cell interaction between macrophage or CD8+ T-cell subsets and endothelial cells (ECs), considering prior findings that ECs have a key role in the activation/recruitment of leukocytes at the initial stages of IgAN(37) (Fig. 6). We detected a wide range of interaction events among ECs, Teff, Tex, resident Mac, inflammatory Mac and proliferating Mac, and the most prevalent ones occurred in ECs, resident Mac and inflammatory Mac (Fig. 6B). By examining L-R pair interactions, several L-R pairs between vascular endothelial growth factor A (VEGFA) and its receptors (KDR and FLT1) exhibited strong potential interaction among ECs themselves and between ECs and proliferating Mac (Fig. 6B). This growth factor induces proliferation of vascular ECs and angiogenesis. While FASLG_FAS and integrin_a4b1_complex_PLAUR showed higher interaction potentials between CD8+ T-cell subsets and ECs (Fig. 6B), which are involved in the apoptotic signaling pathway. These results highlighted the different roles for macrophages and CD8+ T cells in maintaining homeostasis of ECs under stress.3Discussion In our study, an IgA nephropathy (IgAN) model with an extra injection of LPS was established to imitate the onset and resolution of self-limited AKI typically seen in early stages of IgAN. Through longitudinal use of scATAC- and scRNA-seq, we generated a comprehensive cellular landscape of renal immune microenvironment. We obtained a total of 72,304 single-cell transcriptomes from CTRL, IS, IA and IR mice, enabling us to perform a high-resolution mapping of all major immune cell types with an additional detailed analysis of macrophage, neutrophil, CD8+ T-cell and NK-cell subsets. Our data provide an unbiased illustration of the immunological hallmarks as the proportions of immune cell types varied during disease phases. IS is largely similar to the indolent clinical course seen in most IgAN patients, exhibiting low-grade inflammation in the kidney. Although there was a greater shift in the proportion of MPs in comparison to the control kidneys, we nevertheless detected subtle alterations in the proportion of cell subtypes of MPs (Fig. 1E). In IA, the initial phase of AKI, there was an increased proportion of infiltrating Mac, Mefv Mac, Ly6g neutrophils, NK_Gzma and Tem, together with a relatively decreased proportion of inflammatory Mac, Nlrp1b Mac, Csf3r Neutrophils and Tex (Fig. 2B, 3B and 5B). While in the later stage, IR was characterized by an increased proportion of Csf3r Neutrophils and Tex, along with a relatively decreased proportion of Ly6g neutrophils, contrary to IA (Fig. 3B and 5B). These changes collectively contribute to their distinct immune profiles. As neutrophils was the ones with the most dramatic alteration among these cell types, our findings underscore the potential significance of neutrophils in the pathogenesis of IgAN, which has been rarely explored. These analyses also uncovered new subpopulations, including an epigenetically distinct macrophage subset with potential immunoregulatory effects (Nlrp1b Mac), and an early subset of Tex distinguished by effector and cytolytic potential (CX3CR1-transTeff). Our insights into these cell identities may help identify specific targets or pathways for future therapeutic manipulation. However, further investigation incorporating higher-dimensional profiles in the context of IgAN or kidney disease is imperative to resolve the different levels of dysfunctionality within the CD8+ T cell compartment and create a consensus nomenclature(38). Lastly, our analysis detected a wide range of interaction events among endothelial cells, Teff, Tex, resident Mac, inflammatory Mac, and proliferating Mac, emphasizing the different roles for macrophages and CD8+ T cells in maintaining homeostasis of endothelial cells under stress. In conclusion, the study provides a comprehensive analysis of the dynamic changes in immune cell profiles in a model of IgAN, identifying key cell types and their roles and interactions. While these findings contribute to the understanding of IgAN pathogenesis and may provide potential targets for therapeutic intervention, it is important to acknowledge the limitations of the study. The injection of LPS alone cannot replicate the complex renal changes seen in various clinical AKI cases, therefore further validation of these results in human samples or in vivo is necessary once the experimental conditions are established.4Materials and Methods4.1. Mouse modelMale C57BL/6J mice, aged 11±1 weeks and weighting 25 ± 2 g were obtained from SpePharm (Beijing) Biotechnology Co., Ltd. (Beijing, China) and maintained in a clean-grade room at controlled stable temperature and humidity. Our study exclusively examined male mice. It is unknown whether the findings are relevant for female mice. The IgAN mouse model was induced by “BSA + CCl4 + LPS” method as previously described(39). BSA (Sigma) dissolving water (200 mg/kg body weight) was gavaged every other day for 10 weeks, combined with subcutaneous injection weekly and intraperitoneal injection biweekly of CCl4 dissolved in castor oil (1:5; 0.1 ml). Then LPS (Sigma) (50 μg) was injected into tail vein at week six and eight. For the male weight-matched littermate controls served as the normal controls (designated as “CTRL” group), saline was used instead of the above reagents. The model was established at the end of the 10th week. IgA deposits in the glomeruli was observed by direct immunofluorescence, and transmission electron microscopy was also utilized to evaluate model establishment.4.2. Experiment designTwenty-four-hour urine samples were collected the day before renal tissues collected for subsequent experiments. For the IS group, IgAN mice were kept for at least 1 month without treatment after model establishment. For the IA and the IR groups, IgAN mice were subjected to tail intravenous injection of 1 mg/kg LPS the day and 5~7 days before samples and tissues collection, respectively (Fig. S1). 4.3. Kidney CD45-positive immune single-cell isolation and processingKidneys from 4 or 5 mice per experimental group (CTRL, IS, IA, IR) were processed together. Following exsanguination by perfusion of phosphate-buffered saline (PBS) containing 1 mM EDTA, mouse kidneys were cut into small pieces and digested on a rotor at 37 ℃ in Tissue Dissociation Mix (Singleron), according to the manufacturer’s instructions. Dissociated cell suspensions were filtered through a 40-μm cell strainer and further washed with calcium- and magnesium-free PBS. Leukocytes were enriched through 36%~72% Percoll (GE Healthcare) density gradient centrifugation before microbeads (Stem Cell Technologies) isolation, according to the manufacturer’s recommendations. 4.4. Histology and immunofluorescenceFresh renal tissues were fixed in 4% paraformaldehyde and embedded in OCT separately. Four-μm thick formalin-fixed and paraffin-embedded sections were cut and subsequently deparaffinized for periodic acid-Schiff (PAS) staining. Images were acquired by PreciPoint M8 dual digital microscope & scanner. 10 μm frozen sections were made and slides were washed in 0.01M PBS for 20 min. Then the slides were further blocked with blocking buffer (5% bovine serum albumin with 0.1% Triton X-100) for 1 h. Fluorescein isothiocyanate (FITC)-labeled goat anti-mouse IgA (1:400, Abcam) was diluted by blocking buffer and incubated at 4 °C overnight. The next day, slides were washed by 0.01M PBS for three times, 10 min each. Fluorescence images were acquired by Olympus VS120 Virtual Slide Microscope. 4.5. scRNAseq library generationscRNAseq library construction was performed using the GEXSCOPE® Single Cell RNA Library Kit Tissue V2 (Singleron Biotechnologies) as per the manufacturer’s protocol. Purified libraries were sequenced on an Illumina Hiseq X sequencer with 150 bp paired-end reads.4.6. scATACseq library generationscATACseq library construction was performed using the Nuclei Isolation for Single Cell ATAC Sequencing (10× Genomics; CG000169 Rev. E) and Chromium Next GEM Single Cell ATAC Reagent Kits v2 User Guide (10x Genomics; CG000496 Rev. A) as per the manufacturer’s protocol. Libraries were sequenced on an Illumina HiSeq X with 50 bp paired-end reads. 4.7. scRNAseq data processing The gene expression profiles were generated from the raw reads using CeleScope (v1.5.2, Singleron Biotechnologies) with default parameters. Briefly, the barcodes and unique molecular identifiers (UMIs) for each gene-cell combination were extracted and corrected from R1 reads. Adapter sequences and poly A tails were trimmed from R2 reads, which were then aligned against the GRCm38 (mm10) transcriptome using STAR (v2.6.1b). Uniquely mapped reads were assigned to genes using Feature-Counts (v2.0.1). Reads with identical cell barcode, UMI, and gene information were grouped together to create a gene expression matrix for further analysis.Quality control (QC) and clustering analyses were performed using Scanpy (v1.8.1) under Python 3.7(40). Cells that expressed less than 200 genes or had fewer than 5 cells expressing a particular gene were excluded in order to filter out low-quality cells based on three metrics: 1) genes expressed in less than 5 cells; 2) cells expressing less than 200 genes; 3) cells with more than 50% mitochondrial gene expression. A total of 72,304 high-quality cells with an average of 1085 genes per cell remained for downstream analyses after filtering steps.Raw counts data was normalized and transformed into logarithmic scales before selecting the top variable genes by setting flavor = ‘seurat’ among the top-ranked differentially expressed genes. The first twenty principal components identified through Principle Component Analysis (PCA), based on scaled expression profiles, were used for unsupervised clustering utilizing the Louvain algorithm with a resolution parameter set at 1.2. Cell clusters obtained from clustering analysis were visualized using Uniform Manifold Approximation and Projection (UMAP).4.8. Differentially expressed genes analysis The differentially expressed genes (DEGs) were identified using the Scanpy function “sc.tl.rank_genes_group” based on the Wilcoxon rank sum test with default parameters. Genes with a p-value of ≤0.05 and a log2 fold change (log2FC) ≥1 or ≤ -1 were considered significantly up- or down-regulated. 4.9. Cell type annotation The major cell types were identified based on the expression of canonical markers obtained from the reference database SynEcoSysTM (Singleron Biotechnology). SynEcoSysTM encompasses a comprehensive collection of canonical cell type markers derived from CellMakerDB, PanglaoDB, and recently published literature for single-cell sequencing data(41). An additional round of clustering was conducted within major cell types, followed by a comparative analysis of their global gene expression patterns with previously identified murine macrophage and CD8+ T cell gene signatures from scRNAseq studies(23,42).4.10. Pathway enrichment analysisThe potential functions of each subcluster were investigated using Gene Ontology (GO) analysis, implemented with the R package “clusterProfiler” (v 3.16.1)(43). The significantly enriched pathways were determined based on a significance level of p<0.05. Bar plots were generated to visualize the selected significant pathways. GO gene sets representing biological processes (BPs) were used as the reference for pathway analysis.4.11. UCell gene set scoring To elucidate the functional properties of each subcluster of macrophages and neutrophils, gene sets curated from relevant literature were compiled, followed by computation of signature scores for specific gene sets(30,32,44-68). The genes corresponding to each score are listed in supplemental table 1. Gene set scoring was performed using the R package UCell (v 1.1.0)(69). The UCell method is a rank-based scoring approach that is well-suited for analyzing large datasets with multiple samples and batches.4.12. Cell-cell interaction analysisCell-cell interactions (CCIs) involving macrophages and CD4+ T cells, macrophages and endothelial cells, as well as CD8+ T cells and endothelial cells were predicted based on known ligand-receptor pairs using Cellphone DB (v2.1.0)(70). The permutation number for calculating the null distribution of average ligand-receptor pair expression in randomized cell identities was set to 1000. Individual ligand or receptor expression was thresholded using a cutoff based on the average log gene expression distribution across each cell type. Predicted interaction pairs with a p-value < 0.05 and an average log expression > 0.1 were considered statistically significant and visualized using heatmap_plot and dot_plot in CellphoneDB.4.13. RNA velocityThe analysis of RNA velocity was performed using velocyto (v 0.2.3) and scVelo (v0.17.17) in python with default parameters, utilizing BAM files containing macrophages and neutrophils respectively, along with the reference genome GRCm38 (mm10)(71,72). The resulting data was projected onto the UMAP plot derived from Seurat clustering analysis to ensure visual consistency in visualization.4.14. scATACseq data processing and clusteringCells with low quality were filtered out based on TSS enrichment < 4 and nFrag < 1,000, while bin regions overlapping with ENCODE Blacklist regions were excluded from downstream analysis. Subsequently, dimensionality reduction was performed using the iterative latent semantic indexing (LSI) approach through the ArchR (v1.0.1.) addIterativeLSI function(73,74). The ArchR “addClusters()” function was utilized for performing cell clustering. UMAP projection was conducted following the aforementioned procedure.4.15. Integration of scRNA and scATAC dataThe integration of scRNAseq data was performed using Seurat’s integration framework to identify corresponding cell pairs between two modalities of data. In this step, the scRNAseq data served as the reference dataset for training the classifier, and each scATAC cell was assigned a cell type based on its similarity to scRNA cells. Specifically, the “FindTransferAnchors” function (reduction = ‘pcaproject’) was employed to identify shared correlation patterns between scATACseq gene activity and scRNAseq gene expression. Subsequently, the “TransferData” function predicted the cell type label for each cell in the scATACseq data. After filtering, a total of 30,139 cells were retained. The filtered scATACseq objects underwent reprocessing with LSI and clustering using SLM algorithm. To assess consistency between predicted cell identities through label transfer and curated annotations based on known marker gene activities, we utilized Jaccard index.4.16. Peaks CallingThe peaks were identified by MACS2 based on the aggregated insertion sites from all cells of each cell type. A consensus set of uniform-length non-overlapping peaks was obtained by selecting the peak with the highest score from each set of overlapping peaks. In brief, the significance of peaks was ranked and only the most significant peak was retained for further analysis, while any directly overlapping peak was excluded. This process was repeated iteratively until no more peaks remained.4.17. Motif enrichment analysisThe ArchR “getMarkerFeatures()” function was utilized to obtain differential peaks between two clusters. Subsequently, the “addMotifAnnotations()” function followed by the “peakAnnoEnrichment()” function was employed to calculate transcription factor (TF) motif enrichment in these differential peaks. The JASPAR2020 motif dataset was applied in the “addMotifAnnotations()” function to determine the presence of motifs in the peak set. Finally, the “addDeviationsMatrix()” function was used to compute TF activity enrichment on a per-cell basis across all motif annotations based on chromVAR.4.18. Venn diagram Differential peaks between Nlrp1b macrophages (Nlrp1b Mac) and other macrophages were filtered based on a log2FC threshold of 0.125 and an adjusted P value ﹤ 0.05, while overlapping peaks between the two clusters were calculated using Bedtools (v2.31.0) and visualized in a Venn diagram as previously described.4.19. scATAC gene score analysisThe gene scores, which estimate the gene expression and TF motif activity based on scATACseq data, were computed using the “addGeneScoreMatrix()” function with gene score models implemented in ArchR. Subsequently, the gene scores matrix was employed to identify differentially expressed genes, following a similar approach utilized in scRNAseq data analysis.4.20. StatisticsStatistical tests were described in the Methods section. Nonparametric tests in bar graphs were analyzed using a two-sided Student’s t-test with Benjamini-Hochberg correction at P < 0.05. Differential gene expression analysis used an adjusted P < 0.05 cutoff and absolute (average log2FC) > 0.25 threshold. Data distribution was assumed to be normal but not formally tested, and no data were excluded from analyses. Blinding was not performed due to cage labeling requirements; data analysis focused on quantitative rather than qualitative measures. For scRNAseq and scATACseq, 19,000-26,000 cells per sample were collected from pools of 4~5 mice as biological replicates.4.21. Study approvalAll mouse experiments were conducted in accordance with the guidelines of the Institutional Animal Care and Use Committees of the Third Military Medical University.5Conflict of InterestThe authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.6Author ContributionsConceptualization, C. X.; Formal analysis, H. D.; Funding acquisition, X. C.; Investigation, C. X., Y. Z., J. Z. (Jian Zhou), and J. Z. (Jiangnan Zhang); Methodology, Y. Z., J. Z. (Jian Zhou), J. Z. (Jiangnan Zhang), and H. D.; Project administration, Y. W.; Software, C. X., and H. D; Supervision, C. X. and Y. T.; Writing – original draft, C. X.. All authors have read and agreed to the published version of the manuscript.7FundingThis work was supported by grants from the Special Program of the National Natural Science Foundation of China (No. 32141005), General Program of the National Natural Science Foundation of China (No. 32170887 and 31770986), and the Major Research Plan of the National Natural Science Foundation of China (No. 92269110), as well as the Chongqing Talent Program of China (cstc2022ycjh-bgzxm0020), Special fund for performance incentive guidance of research institutions in Chongqing (cstc2021jxjl130036) and Chongqing International Institute for Immunology Project (2022YJC01). The funders had no role in the study design, data analyses, or the decision to publish.8AcknowledgmentsNot applicable. 9References1.Jarrick, S.; Lundberg, S.; Welander, A.; Carrero, J.J.; Höijer, J.; Bottai, M.; Ludvigsson, J.F. Mortality in IgA Nephropathy: A Nationwide Population-Based Cohort Study. J. Am. Soc. Nephrol.. 2019, 30, 866-876.2.Jarrick, S.; Lundberg, S.; Welander, A.; Fored, C.M.; Ludvigsson, J.F. Clinical validation of immunoglobulin A nephropathy diagnosis in Swedish biopsy registers. Clin Epidemiol. 2017, 9, 67-73.3.Magistroni, R.; D'Agati, V.D.; Appel, G.B.; Kiryluk, K. New developments in the genetics, pathogenesis, and therapy of IgA nephropathy. Kidney Int.. 2015, 88, 974-989.4.Wyatt, R.J.; Julian, B.A. IgA nephropathy. N. Engl. J. Med.. 2013, 368, 2402-2414.5.Lai, K.N.; Tang, S.C.; Schena, F.P.; Novak, J.; Tomino, Y.; Fogo, A.B.; Glassock, R.J. IgA nephropathy. Nat Rev Dis Primers. 2016, 2, 16001.6.Heymann, F.; Meyer-Schwesinger, C.; Hamilton-Williams, E.E.; Hammerich, L.; Panzer, U.; Kaden, S.; Quaggin, S.E.; Floege, J.; Gröne, H.J.; Kurts, C. Kidney dendritic cell activation is required for progression of renal disease in a mouse model of glomerular injury. J. Clin. Invest.. 2009, 119, 1286-1297.7.Pei, G.; Zeng, R.; Han, M.; Liao, P.; Zhou, X.; Li, Y.; Zhang, Y.; Liu, P.; Zhang, C.; Liu, X.; et al. Renal interstitial infiltration and tertiary lymphoid organ neogenesis in IgA nephropathy. Clin J Am Soc Nephrol. 2014, 9, 255-264.8.Wilson, P.C.; Wu, H.; Kirita, Y.; Uchimura, K.; Ledru, N.; Rennke, H.G.; Welling, P.A.; Waikar, S.S.; Humphreys, B.D. The single-cell transcriptomic landscape of early human diabetic nephropathy. Proc. Natl. Acad. Sci. U.S.A.. 2019, 116, 19619-19625.9.Park, J.; Shrestha, R.; Qiu, C.; Kondo, A.; Huang, S.; Werth, M.; Li, M.; Barasch, J.; Suszták, K. Single-cell transcriptomics of the mouse kidney reveals potential cellular targets of kidney disease. Science. 2018, 360, 758-763.10.Wu, H.; Malone, A.F.; Donnelly, E.L.; Kirita, Y.; Uchimura, K.; Ramakrishnan, S.M.; Gaut, J.P.; Humphreys, B.D. Single-Cell Transcriptomics of a Human Kidney Allograft Biopsy Specimen Defines a Diverse Inflammatory Response. J. Am. Soc. Nephrol.. 2018, 29, 2069-2080.11.Ransick, A.; Lindström, N.O.; Liu, J.; Zhu, Q.; Guo, J.J.; Alvarado, G.F.; Kim, A.D.; Black, H.G.; Kim, J.; McMahon, A.P. Single-Cell Profiling Reveals Sex, Lineage, and Regional Diversity in the Mouse Kidney. Dev. Cell. 2019, 51, 399-413.e7.12.Rudman-Melnick, V.; Adam, M.; Potter, A.; Chokshi, S.M.; Ma, Q.; Drake, K.A.; Schuh, M.P.; Kofron, J.M.; Devarajan, P.; Potter, S.S. Single-Cell Profiling of AKI in a Murine Model Reveals Novel Transcriptional Signatures, Profibrotic Phenotype, and Epithelial-to-Stromal Crosstalk. J. Am. Soc. Nephrol.. 2020, 31, 2793-2814.13.do Valle Duraes, F.; Lafont, A.; Beibel, M.; Martin, K.; Darribat, K.; Cuttat, R.; Waldt, A.; Naumann, U.; Wieczorek, G.; Gaulis, S.; et al. Immune cell landscaping reveals a protective role for regulatory T cells during kidney injury and fibrosis. JCI Insight. 2020, 5.14.Chen, L.; Chou, C.L.; Knepper, M.A. Targeted Single-Cell RNA-seq Identifies Minority Cell Types of Kidney Distal Nephron. J. Am. Soc. Nephrol.. 2021, 32, 886-896.15.Kuppe, C.; Ibrahim, M.M.; Kranz, J.; Zhang, X.; Ziegler, S.; Perales-Patón, J.; Jansen, J.; Reimer, K.C.; Smith, J.R.; Dobie, R.; et al. Decoding myofibroblast origins in human kidney fibrosis. Nature. 2021, 589, 281-286.16.Muto, Y.; Wilson, P.C.; Ledru, N.; Wu, H.; Dimke, H.; Waikar, S.S.; Humphreys, B.D. Single cell transcriptional and chromatin accessibility profiling redefine cellular heterogeneity in the adult human kidney. Nat Commun. 2021, 12, 2190.17.Zimmerman, K.A.; Bentley, M.R.; Lever, J.M.; Li, Z.; Crossman, D.K.; Song, C.J.; Liu, S.; Crowley, M.R.; George, J.F.; Mrug, M.; et al. Single-Cell RNA Sequencing Identifies Candidate Renal Resident Macrophage Gene Expression Signatures across Species. J. Am. Soc. Nephrol.. 2019, 30, 767-781.18.Conway, B.R.; O'Sullivan, E.D.; Cairns, C.; O'Sullivan, J.; Simpson, D.J.; Salzano, A.; Connor, K.; Ding, P.; Humphries, D.; Stewart, K.; et al. Kidney Single-Cell Atlas Reveals Myeloid Heterogeneity in Progression and Regression of Kidney Disease. J. Am. Soc. Nephrol.. 2020, 31, 2833-2854.19.Zheng, Y.; Lu, P.; Deng, Y.; Wen, L.; Wang, Y.; Ma, X.; Wang, Z.; Wu, L.; Hong, Q.; Duan, S.; et al. Single-Cell Transcriptomics Reveal Immune Mechanisms of the Onset and Progression of IgA Nephropathy. Cell Rep. 2020, 33, 108525.20.Giles, J.R.; Ngiow, S.F.; Manne, S.; Baxter, A.E.; Khan, O.; Wang, P.; Staupe, R.; Abdel-Hakeem, M.S.; Huang, H.; Mathew, D.; et al. Shared and distinct biological circuits in effector, memory and exhausted CD8+ T cells revealed by temporal single-cell transcriptomics and epigenetics. Nat. Immunol.. 2022, 23, 1600-1613.21.Leblay, N.; Ahn, S.; Tilmont, R.; Poorebrahim, M.; Maity, R.; Lee, H.; Barakat, E.; Alberge, J.B.; Sinha, S.; Jaffer, A.; et al. Integrated Epigenetic and Transcriptional Single-Cell Analysis of t(11;14) Multiple Myeloma and Its BCL2 Dependency. Blood. 2023 .22.Meijnikman, A.S.; Gerdes, V.E.; Nieuwdorp, M.; Herrema, H. Evaluating Causality of Gut Microbiota in Obesity and Diabetes in Humans. Endocr. Rev.. 2018, 39, 133-153.23.Fu, J.; Sun, Z.; Wang, X.; Zhang, T.; Yuan, W.; Salem, F.; Yu, S.M.; Zhang, W.; Lee, K.; He, J.C. The single-cell landscape of kidney immune cells reveals transcriptional heterogeneity in early diabetic kidney disease. Kidney Int.. 2022, 102, 1291-1304.24.Loke, P.; Lin, J.D. Redefining inflammatory macrophage phenotypes across stages and tissues by single-cell transcriptomics. Sci Immunol. 2022, 7, eabo4652.25.Wood, M.K.; Daoud, A.; Talor, M.V.; Kalinoski, H.M.; Hughes, D.M.; Jaime, C.M.; Hooper, J.E.; Won, T.; Čiháková, D. PD-L1-expressing macrophages play a protective role in the joint during arthritis. Arthritis & rheumatology. 2023 .26.Everts, B.; Amiel, E.; Huang, S.C.; Smith, A.M.; Chang, C.H.; Lam, W.Y.; Redmann, V.; Freitas, T.C.; Blagih, J.; van der Windt, G.J.; et al. TLR-driven early glycolytic reprogramming via the kinases TBK1-IKKɛ supports the anabolic demands of dendritic cell activation. Nat. Immunol.. 2014, 15, 323-332.27.Cader, M.Z.; Boroviak, K.; Zhang, Q.; Assadi, G.; Kempster, S.L.; Sewell, G.W.; Saveljeva, S.; Ashcroft, J.W.; Clare, S.; Mukhopadhyay, S.; et al. C13orf31 (FAMIN) is a central regulator of immunometabolic function. Nat. Immunol.. 2016, 17, 1046-1056.28.Lee, J.H.; Phelan, P.; Shin, M.; Oh, B.C.; Han, X.; Im, S.S.; Osborne, T.F. SREBP-1a-stimulated lipid synthesis is required for macrophage phagocytosis downstream of TLR4-directed mTORC1. Proc. Natl. Acad. Sci. U.S.A.. 2018, 115, E12228-12228E12234.29.Lu, L.G.; Zhou, Z.L.; Wang, X.Y.; Liu, B.Y.; Lu, J.Y.; Liu, S.; Zhang, G.B.; Zhan, M.X.; Chen, Y. PD-L1 blockade liberates intrinsic antitumourigenic properties of glycolytic macrophages in hepatocellular carcinoma. Gut. 2022, 71, 2551-2560.30.Zhang, C.; Li, J.; Cheng, Y.; Meng, F.; Song, J.W.; Fan, X.; Fan, H.; Li, J.; Fu, Y.L.; Zhou, M.J.; et al. Single-cell RNA sequencing reveals intrahepatic and peripheral immune characteristics related to disease phases in HBV-infected patients. Gut. 2023, 72, 153-167.31.Xu, Y.; He, Z.; Du, J.; Chen, Z.; Creemers, J.; Wang, B.; Li, F.; Wang, Y. Epigenetic modulations of immune cells: from normal development to tumor progression. Int. J. Biol. Sci.. 2023, 19, 5120-5144.32.Li, H.; van der Leun, A.M.; Yofe, I.; Lubling, Y.; Gelbard-Solodkin, D.; van Akkooi, A.C.J.; van den Braber, M.; Rozeman, E.A.; Haanen, J.B.A.G.; Blank, C.U.; et al. Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma. Cell. 2019, 176, 775.33.Lan, X.; Zebley, C.C.; Youngblood, B. Cellular and molecular waypoints along the path of T cell exhaustion. Sci Immunol. 2023, 8, eadg3868.34.Zander, R.; Schauder, D.; Xin, G.; Nguyen, C.; Wu, X.; Zajac, A.; Cui, W. CD4+ T Cell Help Is Required for the Formation of a Cytolytic CD8+ T Cell Subset that Protects against Chronic Infection and Cancer. Immunity. 2019, 51, 1028-1042.e4.35.van Montfoort, N.; Borst, L.; Korrer, M.J.; Sluijter, M.; Marijt, K.A.; Santegoets, S.J.; van Ham, V.J.; Ehsan, I.; Charoentong, P.; André, P.; et al. NKG2A Blockade Potentiates CD8 T Cell Immunity Induced by Cancer Vaccines. Cell. 2018, 175, 1744-1755.e15.36.Good, C.R.; Aznar, M.A.; Kuramitsu, S.; Samareh, P.; Agarwal, S.; Donahue, G.; Ishiyama, K.; Wellhausen, N.; Rennels, A.K.; Ma, Y.; et al. An NK-like CAR T cell transition in CAR T cell dysfunction. Cell. 2021, 184, 6081-6100.e26.37.Zambrano, S.; He, L.; Kano, T.; Sun, Y.; Charrin, E.; Lal, M.; Betsholtz, C.; Suzuki, Y.; Patrakka, J. Molecular insights into the early stage of glomerular injury in IgA nephropathy using single-cell RNA sequencing. Kidney Int.. 2022, 101, 752-765.38.van der Leun, A.M.; Thommen, D.S.; Schumacher, T.N. CD8+ T cell states in human cancer: insights from single-cell analysis. Nat. Rev. Cancer. 2020, 20, 218-232.39.Meng, T.; Li, X.; Ao, X.; Zhong, Y.; Tang, R.; Peng, W.; Yang, J.; Zou, M.; Zhou, Q. Hemolytic Streptococcus may exacerbate kidney damage in IgA nephropathy through CCL20 response to the effect of Th17 cells. PLoS ONE. 2014, 9, e108723.40.Wolf, F.A.; Angerer, P.; Theis, F.J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol.. 2018, 19, 15.41.Zhang, Y.; Li, B.; Duan, J.; Chen, X.; Zhang, X.; Ye, J.; Veloso, A.; Fan, J.; Fang, N. SynEcoSys: a multifunctional platform of large-scale single-cell omics data analysis ,2023.42.Daniel, B.; Yost, K.E.; Hsiung, S.; Sandor, K.; Xia, Y.; Qi, Y.; Hiam-Galvez, K.J.; Black, M.; J Raposo, C.; Shi, Q.; et al. Divergent clonal differentiation trajectories of T cell exhaustion. Nat. Immunol.. 2022, 23, 1614-1627.43.Yu, G.; Wang, L.G.; Han, Y.; He, Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012, 16, 284-287.44.Guo, X.; Zhang, Y.; Zheng, L.; Zheng, C.; Song, J.; Zhang, Q.; Kang, B.; Liu, Z.; Jin, L.; Xing, R.; et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nat. Med.. 2018, 24, 978-985.45.Spranger, S.; Bao, R.; Gajewski, T.F. Melanoma-intrinsic β-catenin signalling prevents anti-tumour immunity. Nature. 2015, 523, 231-235.46.Herbst, R.S.; Soria, J.C.; Kowanetz, M.; Fine, G.D.; Hamid, O.; Gordon, M.S.; Sosman, J.A.; McDermott, D.F.; Powderly, J.D.; Gettinger, S.N.; et al. Predictive correlates of response to the anti-PD-L1 antibody MPDL3280A in cancer patients. Nature. 2014, 515, 563-567.47.Zemmour, D.; Zilionis, R.; Kiner, E.; Klein, A.M.; Mathis, D.; Benoist, C. Single-cell gene expression reveals a landscape of regulatory T cell phenotypes shaped by the TCR. Nat. Immunol.. 2018, 19, 291-301.48.Prager, I.; Watzl, C. Mechanisms of natural killer cell-mediated cellular cytotoxicity. J. Leukoc. Biol.. 2019, 105, 1319-1329.49.Huntington, N.D.; Cursons, J.; Rautela, J. The cancer-natural killer cell immunity cycle. Nat. Rev. Cancer. 2020, 20, 437-454.50.Sangro, B.; Melero, I.; Wadhawan, S.; Finn, R.S.; Abou-Alfa, G.K.; Cheng, A.L.; Yau, T.; Furuse, J.; Park, J.W.; Boyd, Z.; et al. Association of inflammatory biomarkers with clinical outcomes in nivolumab-treated patients with advanced hepatocellular carcinoma. J. Hepatol.. 2020, 73, 1460-1469.51.Guillot, A.; Tacke, F. Liver Macrophages: Old Dogmas and New Insights. Hepatol Commun. 2019, 3, 730-743.52.Yan, J.; Horng, T. Lipid Metabolism in Regulation of Macrophage Functions. Trends Cell Biol.. 2020, 30, 979-989.53.Hao, J.; Yan, F.; Zhang, Y.; Triplett, A.; Zhang, Y.; Schultz, D.A.; Sun, Y.; Zeng, J.; Silverstein, K.; Zheng, Q.; et al. Expression of Adipocyte/Macrophage Fatty Acid-Binding Protein in Tumor-Associated Macrophages Promotes Breast Cancer Progression. Cancer Res.. 2018, 78, 2343-2355.54.Aizarani, N.; Saviano, A.; Sagar; Mailly, L.; Durand, S.; Herman, J.S.; Pessaux, P.; Baumert, T.F.; Grün, D. A human liver cell atlas reveals heterogeneity and epithelial progenitors. Nature. 2019, 572, 199-204.55.Ma, R.Y.; Black, A.; Qian, B.Z. Macrophage diversity in cancer revisited in the era of single-cell omics. Trends Immunol.. 2022, 43, 546-563.56.Tokunaga, R.; Zhang, W.; Naseem, M.; Puccini, A.; Berger, M.D.; Soni, S.; McSkane, M.; Baba, H.; Lenz, H.J. CXCL9, CXCL10, CXCL11/CXCR3 axis for immune activation - A target for novel cancer therapy. Cancer Treat. Rev.. 2018, 63, 40-47.57.Zilionis, R.; Engblom, C.; Pfirschke, C.; Savova, V.; Zemmour, D.; Saatcioglu, H.D.; Krishnan, I.; Maroni, G.; Meyerovitz, C.V.; Kerwin, C.M.; et al. Single-Cell Transcriptomics of Human and Mouse Lung Cancers Reveals Conserved Myeloid Populations across Individuals and Species. Immunity. 2019, 50, 1317-1334.e10.58.Pittet, M.J.; Michielin, O.; Migliorini, D. Clinical relevance of tumour-associated macrophages. Nat Rev Clin Oncol. 2022, 19, 402-421.59.Winkler, J.; Abisoye-Ogunniyan, A.; Metcalf, K.J.; Werb, Z. Concepts of extracellular matrix remodelling in tumour progression and metastasis. Nat Commun. 2020, 11, 5120.60.Martin, P.; Gurevich, D.B. Macrophage regulation of angiogenesis in health and disease. Semin. Cell Dev. Biol.. 2021, 119, 101-110.61.Cheng, S.; Li, Z.; Gao, R.; Xing, B.; Gao, Y.; Yang, Y.; Qin, S.; Zhang, L.; Ouyang, H.; Du, P.; et al. A pan-cancer single-cell transcriptional atlas of tumor infiltrating myeloid cells. Cell. 2021, 184, 792-809.e23.62.Nakamura, R.; Sene, A.; Santeford, A.; Gdoura, A.; Kubota, S.; Zapata, N.; Apte, R.S. IL10-driven STAT3 signalling in senescent macrophages promotes pathological eye angiogenesis. Nat Commun. 2015, 6, 7847.63.Struckmeier, A.K.; Radermacher, A.; Fehrenz, M.; Bellin, T.; Alansary, D.; Wartenberg, P.; Boehm, U.; Wagner, M.; Scheller, A.; Hess, J.; et al. IDO1 is highly expressed in macrophages of patients in advanced tumour stages of oral squamous cell carcinoma. J. Cancer Res. Clin. Oncol.. 2023, 149, 3623-3635.64.Salcher, S.; Sturm, G.; Horvath, L.; Untergasser, G.; Kuempers, C.; Fotakis, G.; Panizzolo, E.; Martowicz, A.; Trebo, M.; Pall, G.; et al. High-resolution single-cell atlas reveals diversity and plasticity of tissue-resident neutrophils in non-small cell lung cancer. Cancer Cell. 2022, 40, 1503-1520.e8.65.Vono, M.; Lin, A.; Norrby-Teglund, A.; Koup, R.A.; Liang, F.; Loré, K. Neutrophils acquire the capacity for antigen presentation to memory CD4+ T cells in vitro and ex vivo. Blood. 2017, 129, 1991-2001.66.Xie, X.; Shi, Q.; Wu, P.; Zhang, X.; Kambara, H.; Su, J.; Yu, H.; Park, S.Y.; Guo, R.; Ren, Q.; et al. Single-cell transcriptome profiling reveals neutrophil heterogeneity in homeostasis and infection. Nat. Immunol.. 2020, 21, 1119-1133.67.Xu, J.; He, B.; Carver, K.; Vanheyningen, D.; Parkin, B.; Garmire, L.X.; Olszewski, M.A.; Deng, J.C. Heterogeneity of neutrophils and inflammatory responses in patients with COVID-19 and healthy controls. Front Immunol. 2022, 13, 970287.68.Mizuno, R.; Kawada, K.; Itatani, Y.; Ogawa, R.; Kiyasu, Y.; Sakai, Y. The Role of Tumor-Associated Neutrophils in Colorectal Cancer. Int J Mol Sci. 2019, 20.69.Andreatta, M.; Carmona, S.J. UCell: Robust and scalable single-cell gene signature scoring. Comput Struct Biotechnol J. 2021, 19, 3796-3798.70.Efremova, M.; Vento-Tormo, M.; Teichmann, S.A.; Vento-Tormo, R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. 2020, 15, 1484-1506.71.La Manno, G.; Soldatov, R.; Zeisel, A.; Braun, E.; Hochgerner, H.; Petukhov, V.; Lidschreiber, K.; Kastriti, M.E.; Lönnerberg, P.; Furlan, A.; et al. RNA velocity of single cells. Nature. 2018, 560, 494-498.72.Bergen, V.; Lange, M.; Peidli, S.; Wolf, F.A.; Theis, F.J. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat. Biotechnol.. 2020, 38, 1408-1414.73.Granja, J.M.; Klemm, S.; McGinnis, L.M.; Kathiria, A.S.; Mezger, A.; Corces, M.R.; Parks, B.; Gars, E.; Liedtke, M.; Zheng, G.; et al. Single-cell multiomic analysis identifies regulatory programs in mixed-phenotype acute leukemia. Nat. Biotechnol.. 2019, 37, 1458-1465.74.Granja, J.M.; Corces, M.R.; Pierce, S.E.; Bagdatli, S.T.; Choudhry, H.; Chang, H.Y.; Greenleaf, W.J. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat. Genet.. 2021, 53, 403-411.10Data Availability StatementThe accession number for the scRNAseq and scATACseq data generated in this study is GSA: CRA019499.