Immune landscape in rejection of renal transplantation revealed by high-throughput single-cell RNA sequencing

Background: The role of the cellular level in kidney transplant rejection is unclear, and single-cell RNA sequencing (scRNA-seq) can reveal the single-cell landscape behind rejection of human kidney allografts at the single-cell level. Methods: High-quality transcriptomes were generated from scRNA-seq data from five human kidney transplantation biopsy cores. Cluster analysis was performed on the scRNA-seq data by known cell marker genes in order to identify different cell types. In addition, pathways, pseudotime developmental trajectories and transcriptional regulatory networks involved in different cell subpopulations were explored. Next, we systematically analyzed the scoring of gene sets regarding single-cell expression profiles based on biological processes associated with oxidative stress. Results: We obtained 81,139 single cells by scRNA-seq from kidney transplant tissue biopsies of three antibody-mediated rejection (ABMR) patients and two acute kidney injury (AKI) patients with non-rejection causes and identified 11 cell types, including immune cells, renal cells and several stromal cells. Immune cells such as macrophages showed inflammatory activation and antigen presentation and complement signaling, especially in rejection where some subpopulations of cells specifically expressed in rejection showed specific pro-inflammatory responses. In addition, patients with rejection are characterized by an increased number of fibroblasts, and further analysis of subpopulations of fibroblasts revealed their involvement in inflammatory and fibrosis-related pathways leading to increased renal rejection and fibrosis. Notably, the gene set score for response to oxidative stress was higher in patients with rejection. Conclusion: Insight into histological differences in kidney transplant patients with or without rejection was gained by assessing differences in cellular levels at single-cell resolution. In conclusion, we applied scRNA-seq to rejection after renal transplantation to deconstruct its heterogeneity and identify new targets for personalized therapeutic approaches.


Introduction
Kidney transplantation is a very effective treatment for endstage renal disease, but complications after kidney transplantation have been a major clinical problem and long-term survival rates have been suboptimal (Augustine, 2018). Rejection reactions remain a major threat and are the main independent risk factor for longterm survival of transplanted kidneys (Cooper, 2020;Wang et al., 2021). The two phases of immune response in allografts are different in early and late stages. In the past, graft acceptance without longterm immunosuppression was not uncommon in living transplants (Lim et al., 2017). Although modern immunosuppressive techniques are now available and the short-term results of kidney transplantation have improved to >90% within 1 year, the long-term results over 20 years have not improved significantly (Valenzuela and Reed, 2017;Eskandary et al., 2018). Fibrotic arterial thickening, interstitial fibrosis and tubular atrophy not only severely affect graft function, but also survival (Wekerle et al., 2017). A large amount of data has been established for bulk transcription analysis to predict graft outcome. However, these analyses are limited, and traditional bulk RNA-seq and kidney biopsy methods ignore transcriptome heterogeneity at single-cell resolution (Liao et al., 2020).
Single-cell RNA sequencing (scRNA-seq) plays a key role in identifying cellular subtypes and elucidating differential expression between genes (Saliba et al., 2014;Li et al., 2019;Zhao et al., 2020). Recently, there have been some new discoveries about scRNA-seq, revealing the mechanism of kidney disease (Wu et al., 2018). For example, scRNA-seq analysis of systemic lupus erythematosus nephritis identified high expression of interferon-inducible genes in renal tubular cells associated with disease severity (Der et al., 2017). Currently, our understanding of the single-cell level of rejection, especially antibody-mediated rejection (ABMR), is poor. Numerous analyses have shown that ABMR is the leading cause of late-stage allograft failure (Sellares et al., 2012). The balance of the inflammatory and anti-inflammatory responses of the immune system, and the cytokine storm caused by the immune response after transplantation may cause rejection, thereby affecting the survival rate of transplantation (Salehi et al., 2020). Therefore, analysis of cytokine production at the single cell level may be important to elucidate the fate of the graft.
In this study, we explored changes at the single-cell level in kidney transplanted living tissue from 5 patients with AKI (including 3 with ABMR and 2 nonrejection cause), which in turn determined the pro-inflammatory parenchymal response in the rejected kidney. Our data provide biological insights into renal transplant rejection by constructing a very comprehensive singlecell landscape of renal transplant rejection, describing the molecular function of the cells, which will aid reduce the adverse reactions after renal transplantation.

Data collection and processing
The GPL24676 platform-based renal transplant biopsy-related scRNA-seq dataset GSE145927 (Malone et al., 2020) was downloaded from Gene Expression Omnibus (GEO, https://www. ncbi.nlm.nih.gov/geo/). A total of 5 patients with AKI (including 3 with ABMR and 2 nonrejection cause) were included and their kidney allograft tissue was biopsied, with biopsies performed between 5 days and 7 years after transplantation. All patients diagnosed with ABMR had donor-specific antibodies (DSA) at the time of biopsy (Supplementary Table S1). In addition, wholeexome sequencing data for donor and recipient DNA in the dataset were not consistent with this study, and this data was not included in the study.

Construction of single cell atlas
Quality control of single-cell Seurat objects on double cell numbers, dead cells, and mitochondrial genes using Seurat package (Stuart et al., 2019) in R language. For all downstream analyses, we selected cells that have at least 1,000 unique molecular identifiers (UMIs) (indicating the number of captured transcripts) mapped to at least 200 unique genes and ensured that each gene is expressed in more than 3 cells. We excluded cells with poor viability and quality by removing more than 10% of the cells whose gene counts reflected mitochondrial genes or ribosomal RNA (Supplementary Figure S1). Single-cell data were merged using the IntegrateData function (Butler et al., 2018) of the Seurat package (Stuart et al., 2019) in R language to perform cell clustering analysis according to default parameters. The clustering results were downscaled and visualized (Becht et al., 2018) based on the uniform manifold approximation and projection (UMAP) technique, and projected onto a two-dimensional image defined as a single-cell atlas. In addition, the FindAllMarkers function of the Seurat package was used to identify genes specifically expressed in each cell cluster, with p-value < 0.05 and |log fold change (logFC)| > 0.5 being considered significant. In addition, cell types were annotated according to cellular markers known from previous studies (Malone et al., 2020).

Differential gene expression analysis
At the single-cell level, differentially expressed genes in each cluster between allograft biopsies from patients with and without rejection were determined using the "FindAllMarkers" function, and differences associated with adjusted p-value < 0.05 were considered significant.
Frontiers in Cell and Developmental Biology frontiersin.org

Functional enrichment and gene enrichment analysis
To explore the biological processes and pathways involved in the marker for each cell cluster, the R package clusterProfiler (Yu et al., 2012) was applied to the marker for enrichment analysis regarding the biological processes (BPs) of gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway, with p-value < 0.05 BPs and KEGG signaling pathways were considered significant.

Pseudotime analysis
Based on the changes in gene expression of different cell subpopulations over time, we reconstructed the developmental trajectory of dysregulated cell differentiation in renal transplant rejection using the Monocle 3 package (Trapnell et al., 2014) in R language to explore the changes in dysregulated immune cells during the development of renal transplant rejection.

Gene regulatory network (GRN) analysis
In addition, using the Python module tool pySCENIC ( Van de Sande et al., 2020), this study comprehensively reconstructed the transcription factor-centered gene regulatory network to further explore the regulatory mechanisms of dysregulated cells. The tool first uses a per-target regression approach (GRNBoost2) to infer co-expression modules, then uses cis-regulatory motif discovery (cisTarget) to prune indirect targets from these modules, and finally quantifies the activity of these regulators by enrichment scores of regulator target genes.

Cellular scoring of oxidative stressrelated gene sets
The AddModuleScore function (Tirosh et al., 2016) in the seruat package was used to score the specified gene sets. Oxidative stressrelated pathway dataset from the Gene Set Enrichment Analysis (GSEA) database (Supplementary Table S2). This function is to calculate module scores for feature expression programs in single cells in the oxidative stress-related pathway dataset.

Data Analysis and Statistics
All statistical analyses were performed in R package. Comparisons between the two groups were made using Student's t test and correlation coefficients were calculated using Spearman analysis. p-value < 0.05 was considered significant.

Global single-cell landscape of renal transplant patients
We obtained single-cell sequencing data from public databases of biopsied renal allograft tissue from three patients with ABMR and two patients with a nonrejection cause of AKI to further explore the potential ecological panoply of renal transplant patients. After standardized data processing and quality control, a total of 81,139 high quality single cells were obtained and clustered to generate 26 cell clusters ( Figure 1A), including 11 cell types such as endothelial cells (En), epithelial cells (Ep), fibroblasts, macrophages (Mac), loop of henle (LOH), stromal cells (Stroma), proximal tubule (PT), principal cells (PC), B cells, CD8 + T cells, and intercalated cells (IC) ( Figure 1B). The cell clusters positively expressed markers were consistent with the gene signatures published in recent scRNA-seq and laboratory studies, among others (Wu et al., 2018;Malone et al., 2020), in line with the phenotypic characteristics of the corresponding cells ( Figure 1C). Differential gene expression analysis revealed a general dysregulation of gene expression among different cell types in the ABMR group relative to the nonrejection patient group ( Figure 1D). Further comparison of the differences in cell composition between rejection-free, ABMR patients revealed a dramatic increase in fibroblast abundance in ABMR patients, while an increased abundance of macrophages was observed ( Figure 1E). In summary, we initially constructed a dynamic single-cell ecological global landscape of the microenvironment of biopsied renal allograft tissue in renal transplant patients by single-cell histology, and we found significant gene expression dysregulation accompanying between different cell types and explored the altered cellular ecology in nonrejection and ABMR patients.

Landscape of fibroblast subpopulations in renal transplant patients
Fibroblasts are thought to be closely associated with allograft fibrosis (Wang et al., 2017), and in our study, we found significant fiber abundance in ABMR patients. Based on a single-cell resolution cellular ecological atlas, we further explored the subpopulations of fibroblasts in renal transplant patients and identified a total of eight subpopulations of fibroblasts ( Figure 2A) and found that these subpopulations were heterogeneous among different subgroups ( Figure 2B). Among these fibroblast subpopulations, we found that Fibroblasts_PLVAP and Fibroblasts_CXCL9 had the highest proportions in ABMR patients, and Fibroblasts_MME and Fibroblasts_MT1H had the highest proportions in the group of patients with a nonrejection cause ( Figure 2C), while specific markers for these fibroblast subpopulations were mapped in the single cell atlas ( Figure 2D). By enrichment analysis, these specifically altered fibroblast subpopulations were found to be significantly enriched for Allograft rejection, Graft-versus-host disease and Renin secretion, which are relevant pathways regarding post-transplant rejection, and in addition, we observed that these subpopulations are also involved in the TGF-β signaling pathway, epithelial to mesenchymal transition (EMT) and angiogenesis-related pathways ( Figure 2E). These fibroblast subpopulations underwent a continuous developmental process, with Fibroblasts_PLVAP located at the beginning of development ( Figure 2F). Subsequent GRN analysis showed that this subpopulation of genes was organized into two modules ( Figure 2G), each subject to transcription factors (TFs) that regulate specific gene expression in kidney transplant patient-Frontiers in Cell and Developmental Biology frontiersin.org specific fibroblasts to guide cell fate selection ( Figure 2H). These results reflect that the Fibroblasts_PLVAP and Fibroblasts_ CXCL9 subpopulations were activated and associated with EMT of the renal tubules in patients with ABMR compared to patients with a nonrejection cause, which may have implications for progressive tubulointerstitial fibrosis.    Frontiers in Cell and Developmental Biology frontiersin.org 07

Landscape of epithelial cell subpopulations in renal transplant patients
Extensive immune system activation can cause necrosis of renal tubular epithelial cells, rupture of the basement membrane, and eventually fibrosis and loss of renal transplant function (Rogers et al., 2016). By subpopulation analysis, ten epithelial cell subpopulations were further identified in this study ( Figure 3A) and the distribution of these subpopulations in different patient groups was observed ( Figure 3B). Among them, Ep_S100A1 and Ep_EGR1 were highly represented in patients with ABMR compared to patients with a nonrejection cause ( Figure 3C), and the expression of specific marker genes of these subpopulations was further observed at single-cell resolution ( Figure 3D). Subsequently, enrichment analysis revealed that these specifically altered epithelial cell subpopulations were significantly enriched for signaling pathways associated with post-transplant rejection, such as Allograft rejection, Graft-versus-host disease, and Renin secretion, and were also involved in pan-death pathways, IL-17 signaling pathway, Cell adhesion molecules, Leukocyte transendothelial migration, and ECM-receptor interaction ( Figure 3E). In addition, pseudotime analysis revealed high mimetic time values for Ep_S100A1 and Ep_EGR1, indicating that they are at the end of the developmental trajectory ( Figure 3F). By GRN, we explored TFs that regulate epithelial cell subpopulations (Figures 3G, H).

Landscape of LOH cell subpopulations in renal transplant patients
By subpopulation analysis, we further identified six LOH cell subpopulations ( Figures 4A, B). In addition, we found that LOH_ CLCNKA and LOH_JUN were highly represented in patients with ABMR compared to patients with a nonrejection cause ( Figure 4C) and further explored the expression of specific markers for these subpopulations by single-cell resolution ( Figure 4D). Enrichment analysis revealed that cell subpopulations of LOH were significantly involved in signaling pathways related to post-transplant rejection such as Proximal tubule bicarbonate reclamation, Graft-versus-host disease and Allograft rejection, and also significantly involved in IL-17 signaling pathway and TNF signaling pathway of inflammation-related pathways ( Figure 4E). By pseudotime analysis, LOH_CLCNKA and LOH_JUN were found to be at the end of the developmental trajectory of LOH cells ( Figure 4F). In addition GRN showed that LOH cells were regulated by TFs such as EGR1, ZNF429, ZNF491, and HNF4A, respectively ( Figures 4G, H). In conclusion, this study identified LOH cell subpopulations in renal transplant patients and found that these subpopulations were significantly involved in a number of signaling pathways related to post-transplant rejection and inflammation-related pathways. In addition, the transcriptional regulation of cell subpopulations during their differentiation was further elucidated by inferring the developmental trajectories of LOH cell subpopulations.

Landscape of macrophage subpopulations in renal transplant patients
We also performed a more refined typing of macrophages in renal transplant patients. Ten macrophage subpopulations were identified by subpopulation analysis in this study, and again these subpopulations showed heterogeneity across patient groups ( Figures 5A, B). Notably, we found the highest proportion of Mac_FCN1 in patients with ABMR compared to patients with a nonrejection cause ( Figure 5C), and the specific marker of Mac_FCN1 was verified by single-cell resolution ( Figure 5D). We then further explored the pathways involved in the macrophage subpopulation in kidney transplant patients and showed that this subpopulation was significantly enriched in Chemokine signaling pathway, TNF signaling pathway and Viral protein interaction with cytokine and cytokine receptor, in addition to some post-transplant rejection-related pathways were also enriched ( Figure 5E). The pseudotime analysis showed that this subpopulation is early in the developmental trajectory of macrophages, indicating that this subpopulation also has a high developmental potential ( Figure 5F). Subsequently, by GRN, we identified a series of TFs regulating the macrophage subpopulation ( Figures 5G, H). In conclusion, by identifying macrophage subpopulations, we identified the Mac_ FCN1 subpopulation, that is, highly represented in the microenvironment of patients with ABMR after renal transplantation, which has a high differentiation potential and may play an important role in the inflammatory response after renal transplantation.

Landscape of stromal cell subpopulations in renal transplant patients
Stromal cells play a key role in pro-inflammatory and profibrotic activity in renal transplantation (Venner et al., 2016;Cippa et al., 2019). By subpopulation analysis, we identified a total of 9 subpopulations of stromal cells and observed the distribution of subpopulations in different patient groups ( Figures 6A, B). Among them, we also found that Stromal_ PDZK1IP1 and Stromal_MFAP5 had the highest proportion in patients with ABMR compared to patients with a nonrejection cause, and conversely, Stromal_FBLN1 and Stromal_FOSB had a lower proportion ( Figure 6C), and marker genes for each stromal cell subpopulation are shown in Figure 6D. Enrichment analysis revealed that these stromal cell subpopulations were significantly enriched for post-transplant rejection-related pathways and pro-inflammatory pathways, and specifically, cell signaling pathways for adhesion were also enriched ( Figure 6E). The pseudotime analysis showed that Stromal_ PDZK1IP1 and Stromal_MFAP5 were at the end of the developmental trajectory of stromal cells ( Figure 6F). In addition, GRN analysis revealed that stromal cell subpopulations were regulated by TFs such as TAL1, NR1H4, and FOXK1, respectively ( Figures 6G, H)

. In conclusion, by identifying subpopulations of stromal cells in renal transplant
Frontiers in Cell and Developmental Biology frontiersin.org

Oxidative stress in single cells in renal transplant rejection
Calculation of scores for oxidative stress-related biological processes in single-cell samples using AddModuleScore revealed much lower scores for GO_CELLULAR_OXIDANT_ DETOXIFICATION in patients with a nonrejection cause than in patients with ABMR ( Figure 7A). In addition, stromal cells ( Figure 7B), Ep ( Figure 7C), fibroblasts ( Figure 7D), LOH cells ( Figure 7E), and macrophage ( Figure 7F) were differentially active in response to oxidative stress, indicating a strong influence by reactive oxygen species in renal transplant rejection.

Discussion
In this study, we present a comprehensive overview of cell types and subpopulations in kidney biopsy samples from three patients with ABMR and two patients with a nonrejection cause of AKI at single-cell resolution. By assessing the ecological changes and different signaling profiles of subpopulations of different cell types and analyzing the developmental trajectories and transcriptional regulators of different cell subpopulations, our data will help advance our insights into transcriptional differences between cells of different origins in the posttransplant rejection response of kidneys.
First, our analysis identified 11 independent cell types, including a number of renal tubular cell types, epithelial cell types, immune cells, stromal cells, and endothelial cells. Second, each cell type exhibited different pathway characteristics and activity between Frontiers in Cell and Developmental Biology frontiersin.org patients with ABMR and patients with a nonrejection cause of AKI or between each other. Tissue remodeling is dependent on mesenchymal cells (fibroblasts and myofibroblasts), which is a prominent feature of chronic rejection in allograft kidneys (Grimm et al., 2001). In our study, the increased proportion of fibroblasts in the kidney biopsy tissue of patients with ABMR confirmed that loss of allograft function and the development of interstitial fibrosis in patients with allograft rejection are associated, consistent with previous observations (Toki et al., 2014). The fibroblast subpopulations Fibroblasts_PLVAP and Fibroblasts_CXCL9, plasma vesicle-associated protein (PLVAP) genes identified in ABMR patients are strongly associated with chronic allogeneic rejection and are associated with TG grading and proteinuria (Yamamoto et al., 2007). In contrast, earlier studies have demonstrated that urinary CXCL9 and CXCL10 levels may be predictive of T-cell-mediated rejection in the early posttransplant period and that measurement of urinary CXCL9 and CXCL10 levels may provide an additional tool for the diagnosis of rejection (Ciftci et al., 2019). In addition, it has been recognized that inflammation is a key process leading to progressive renal fibrosis (Meng et al., 2014). In an enrichment analysis of fibroblast subpopulations, these subpopulations were found to be significantly involved in a number of inflammationrelated pathways. Although fibroblasts may stimulate renal fibrosis in an indirect manner through the production of proinflammatory and pro-fibroblastic cytokines and growth factors (Nikolic-Paterson et al., 2011), a direct link between fibroblasts and myofibroblast accumulation during chronic allogeneic rejection remains undefined. In addition, TGF-β signaling pathway, epithelial-to-mesenchymal transition (EMT) and angiogenesis-related pathways are enriched by fibroblast subpopulations and TGF-β 1 is an important mediator of renal fibrosis after renal transplantation (Campistol et al., 2001). In the context of renal fibrosis, it is necessary for the transition of renal tubular epithelial and endothelial cells to mesenchymal cells (Sato et al., 2003;Yeh et al., 2010;Zhou et al., 2010). Exploring answers to the question of fibroblast subpopulations in allograft kidney biopsies in patients with ABMR and patients with a nonrejection cause of AKI has prompted research into the molecular clues that lead to greater recruitment of fibroblasts to kidneys with rejection, which may have enormous therapeutic implications.
ABMR is commonly associated with glomerular and interstitial inflammation (Magil and Tinckam, 2003;Sund et al., 2004). Activated resident cells, particularly renal tubular epithelial cells, represent a key feature of acute graft rejection and any inflammatory response in natural and transplanted kidneys (Woltman et al., 2000). Several recent evidences suggest that IL-17 can regulate the innate immune response and, to some extent, the adaptive immune response (Ishigame et al., 2009;Sutton et al., 2009). Our study further confirms that renal tubular epithelial cells may serve as immunomodulatory cells in renal allograft rejection.
There is evidence that intra-graft macrophages are associated with poor renal graft outcomes in human and animal models (Matheson et al., 2005;Toki et al., 2014;Ikezumi et al., 2015). In our study, we identified functional heterogeneity among macrophages. Macrophages can promote injury and repair of renal allografts depending on the underlying injury and the effectiveness of immunosuppressive therapy. A subpopulation of macrophages positive for FCN1 was identified in patients with ABMR, FCN1 is involved in the pathogenesis of autoimmune diseases and targeting FCN1 is a promising strategy for the treatment of these diseases (Katayama et al., 2019), however, no studies have yet confirmed the role of FCN1 in renal transplant rejection. In renal transplant rejection, pro-inflammatory macrophages contribute to graft loss and are characterized by the production of pro-inflammatory cytokines and chemokines (IL-1β, MCP-1 and TNF-α) as well as reactive oxygen species production (Qi et al., 2008;Salehi and Reed, 2015), which is consistent with our findings, which found that macrophages appear to be involved in a number of pathways including Chemokine signaling pathway, TNF signaling pathway and Viral protein interaction with cytokine and cytokine receptor in pro-inflammatory processes. In addition, one study found that macrophages can be converted into collagenproducing myofibroblasts, providing a novel mechanism for the direct involvement of macrophages in interstitial fibrosis during chronic renal allogeneic rejection (Wang et al., 2017). This led us to think about the link between macrophages and fibroblasts.
Oxidative stress is a major mediator of adverse outcomes throughout the transplantation process. Transplanted kidneys are susceptible to oxidative stress-mediated injury due to preand post-transplant conditions, leading to reperfusion injury or an imbalance between oxidants and antioxidants. In addition to adversely affecting allografts, oxidative stress and its persistent partners, inflammation can lead to cardiovascular disease, cancer, metabolic syndrome and other disorders in transplant recipients (Nafar et al., 2011). We observed significant oxidative stress in patients with rejection after renal transplantation in the present study, and oxidative stress may be associated with various acute injuries that occur after transplantation. However, the study of oxidative stress within single cells in the transplantation population remains elusive and needs to be explored in the future.
In conclusion, we have demonstrated that despite the apparent complexity and heterogeneity of rejection after kidney transplantation, we envision a future in which kidney biopsies subjected to scRNA-seq could be used as part of a molecular diagnostic test. This study reveals hidden transcriptional heterogeneity in population-average measures by analysing scRNA-seq data from kidney biopsies of patients with ABMR and patients with a non-rejection cause of AKI, providing ideas for identifying new targets for personalised therapeutic approaches. However, our analysis has several limitations. First, the sample size was small, and more scRNA-seq data from post-transplant kidney biopsy tissues should be collected for analysis and exploration in the future. Second, this study lacks the analysis of specific cell/gene expression more or less related to specific clinical events, and the single-cell profile behind human kidney transplant rejection obtained in this study can only serve as a reference for subsequent specific cellular/molecular studies, which require further analysis and validation. Finally, this study investigated rejection after kidney transplantation at the single-cell level by bioinformatics techniques, which should be validated by cellular experiments in the future.
Frontiers in Cell and Developmental Biology frontiersin.org