The orchestrated cellular and molecular responses of the kidney to endotoxin define the sepsis timeline

Clinical sepsis is a highly dynamic state that progresses at variable rates and has life-threatening consequences. Staging patients along the sepsis timeline requires a thorough knowledge of the evolution of cellular and molecular events at the tissue level. Here, we investigated the kidney, an organ central to the pathophysiology of sepsis. Single cell RNA sequencing revealed the involvement of various cell populations in injury and repair to be temporally organized and highly orchestrated. We identified key changes in gene expression that altered cellular functions and can explain features of clinical sepsis. These changes converged towards a remarkable global cell-cell communication failure and organ shutdown at a well-defined point in the sepsis timeline. Importantly, this time point was also a transition towards the emergence of recovery pathways. This rigorous spatial and temporal definition of murine sepsis will uncover precise biomarkers and targets that can help stage and treat human sepsis.


Introduction 29
Acute kidney injury (AKI) is a common complication of sepsis that doubles the mortality risk. In 30 addition to failed homeostasis, kidney injury can contribute to multi-organ dysfunction through 31 distant effects. Indeed, the injured kidney is a significant mediator of inflammatory chemokines, 32 cytokines, and reactive oxygen species that can have both local as well as remote deleterious 33 effects 1-4 . Therefore, understanding the complex pathophysiology of kidney injury is crucial for 34 the comprehensive treatment of sepsis and its complications. 35 We have recently shown that renal injury in sepsis progresses through multiple phases. These 36 include an early inflammatory burst followed by a broad antiviral response and culminating in 37 translation shutdown and organ failure 5 . In a non-lethal and reversible model of endotoxemia, 38 organ failure was followed by spontaneous recovery. The exact cellular and molecular 39 contributors to this multifaceted response remain unknown. Indeed, the kidney is architecturally 40 a highly complex organ in which epithelial, endothelial, immune and stromal cells are at 41 constant interplay. Therefore, we now examined the spatial and temporal progression of 42 endotoxin injury to the kidney using single cell RNA sequencing (scRNAseq). Our data revealed 43 that cell-cell communication failure is a major contributor to organ dysfunction in sepsis. 44 Remarkably, this phase of communication failure was also a transition point where recovery 45 pathways were activated. We believe this spatially and temporally anchored approach to sepsis 46 pathophysiology is crucial for identifying potential biomarkers and therapeutic targets. We harvested a cumulative amount of 63,287 renal cells obtained at 0, 1, 4, 16, 27, 36 and 48 55 hours after endotoxin (LPS) administration. The majority of renal epithelial, immune and 56 endothelial cell types were represented (Fig. 1a). Note the absence of podocyte and mesangial 57 cells, which can be a limitation of single cell RNAseq renal dissociation procedures 6 . Cluster 58 identities were assigned and grouped using known classical phenotypic markers (Fig. 1b,  59 Supplementary Fig. 1a) [7][8][9][10][11] . Interestingly, the UMAP-based computational layout of epithelial 60 clusters recapitulated the normal tubular segmental order in the nephron. This indicates that 61 gene expression gradually changes among neighboring tubular segments along the nephron. 62 Note that the expression of cluster-defining markers varied significantly during the injury and 63 recovery phases of sepsis ( Fig. S1b; Supplementary Table 1). Therefore, we also identified a 64 set of genes that are conserved across time for a given cell type (Fig. S1c). 65 In the integrated UMAP (Fig. 1a), we noted the presence of a proliferative cell cluster (Cdk1 and 66 Ki67 expression). By back mapping to time-specific unintegrated UMAPs, we determined that 67 these proliferating cells could be traced to specific cell types at various points along the sepsis 68 timeline (Fig. 1c). For example, within the first hour after LPS, these proliferative indices were 69 expressed primarily in S1 cells. These cells are the site of LPS uptake in the kidney as we have 70 previously shown 12-14 . At later time points, proliferative indices are seen in macrophages (4 71 hours) and S3 cells (36 hours) (Fig. 1c). These proliferative indices reflect cell cycle activity 72 which may be involved in injury, repair or recovery processes 15 . 73 We also noted the presence of a proximal tubular cluster expressing unique gene identifiers: 74 Agt, Rnf24, Slc22a7 and Slc22a13 (Fig. 2a). This is likely the proximal tubular S3-Type 2 75 (S3T2) reported by others 16 . This cluster maintained a separate and distinct identity throughout 76 the sepsis timeline (Fig. 1c). Because the location of S3T2 is currently unknown, we performed 77 in-situ spatial transcriptomics on septic mouse kidneys 17 . We then integrated our scRNAseq 78 with the in-situ RNAseq in order to map our scRNAseq clusters onto the tissue (Supplementary 79 Fig. 2a, S2b). We found that the classical S3 cluster localizes to the cortex while S3T2 is in the 80 outer stripe of the outer medulla (Fig. 2b, Supplementary Fig. 2b). We confirmed the location 81 of S3T2 to the OS-OM with single molecular FISH (Supplementary Fig. 2c). The differential 82 gene expression between S3 and S3T2 is likely dictated by regional differences in the 83 microenvironments of the cortex and the outer stripe. 84 Because angiotensinogen (Agt) was strongly expressed in S3T2, we examined the expression 85 of other components of the renin-angiotensin system (RAS). We first noted the absence of Ace 86 expression in S1 tubular cells (Fig. 2c, Supplementary Fig. 3). In contrast, Ace2 was strongly 87 expressed in S1, S3 and S3T2 cells. There is currently great interest in understanding the 88 biology of Ace2 because of its role in SARS-CoV-2 cellular invasion. Other essential 89 components of the SARS-CoV-2 entry mechanism include Tmprss2 and Slc6a19 18-22 . While 90 Tmprss2 was expressed in all proximal tubular segments, Slc6a19 was more strongly 91 expressed in S1 throughout the sepsis timeline. This may point to the S1 tubular segment as 92 one point of entry of SARS-CoV-2 into the kidney. 93

Cell trajectory and velocity field analyses of scRNAseq characterize subpopulations of immune 94 cells 95
The immune cell profile in the septic kidney was time-dependent and showed a five-fold 96 increase in immune cells, primarily macrophages (Fig. 3a, 3b). We noted two distinct 97 macrophage clusters denoted as Macrophage A and Macrophage B (Mϕ-A, Mϕ-B). Both of 98 these clusters expressed classical macrophage markers such as Cd11b (Itgam) (Fig. 3c). 99 accumulated macrophages were predominantly Mϕ-A. We noted the absence of proliferation 101 markers (Cdk1, Ki67) in this cluster, raising the possibility that this may be an infiltrative 102 macrophage type (Fig. 3d). The Mϕ-B cluster, located between Mϕ-A and conventional 103 dendritic cells (cDC) expressed also cDC markers such as MHC-II subunit genes (H2-Ab1) and 104 Cd11c (Itgax) indicating that it is an intermediary macrophage type (Fig. 3c). This continuum 105 between macrophages and dendritic cells in the kidney has been reported 23-26 . Interestingly, 106 Mϕ-B cells expressed proliferation markers (Cdk1, Ki67) and thus, may be differentiating 107 towards a Mϕ-A or cDC phenotype (Fig. 3c). Pseudotime and velocity field analysis suggested 108 that at earlier time points (1 hour) Mϕ-B was differentiating toward Mϕ-A phenotype. At later 109 time points (36 hours) the velocity field suggested that Mϕ-B was differentiating towards cDC 110 but pseudotime analysis was inconclusive (Fig. 3e). Similarly, the Mϕ-A cluster also showed two 111 subclusters on the RNA velocity map (Supplementary Fig. 4a). One of the subclusters showed 112 increased expression of alternatively activated macrophages (M2) markers such as Arg1 113 (Arginase 1) and Mrc1 (Cd206) 27 at later time points (36 hours, Supplementary Fig. 4b). 114 Therefore, RNA velocity analysis may be a useful tool in distinguishing macrophage subtypes in 115 scRNAseq data. 116 In T-cells, while Cd4 expression was minimal at all time points, the expression of Cd8 was 117 robust and relatively preserved over time (Fig. S4c). We also noted an increase of a distinct 118 plasmacytoid dendritic cell cluster at one hour (pDC). These pDCs, along with natural killer (NK) 119 cells, are known to signal through the interferon-gamma pathway and stimulate Cd8 expression 120 28,29 . This supports the early antiviral response we have previously reported in this sepsis model defined with pseudotime analysis. We note that at any given time point, directional progression 128 of states along pseudotime correlated well with real time state changes (Fig. 4a). Note that the 129 endothelium exhibited changes in states as early as 1 hour, while S1 showed changes at later 130 time points (4 hours). These sequential state changes may reflect the spatial and temporal 131 propagation of LPS signaling in the kidney. As sepsis progressed, many cell types lost function-132 defining markers while acquiring novel ones. For example, S1 and S3 lost classical markers like 133 Slc5a2 (SGLT2) and Aqp1 and expressed new genes involved in antigen presentation such as 134 H2-Ab1 (MHC-II) and Cd74 (Fig. 4b). Moreover, the highly distinct phenotypes that 135 differentiated S1 from S2/S3 at baseline merged into one phenotype for all three sub-segments 136 by 16 hours after LPS (Fig. 4c). However, despite the apparent convergent phenotype at 16 137 hours, additional analytical approaches such as RNA velocity revealed significant differences in 138 RNA splicing kinetics between S1 and S3 segments at this time point. In addition, RNA velocity 139 revealed the presence of two subclusters within the S3 segment at 16 hours (Fig. 4d). These 140 two velocity subclusters did not correlate with the two states seen in pseudotime analysis. This 141 indicates that multiple analytic approaches are needed to fully characterize cellular changes 142 along the sepsis timeline. 143

Sepsis induces time and cell-specific genes and pathways 144
We next show gene expression profiles in select cell types along the sepsis timeline. In this 145 analysis, we included endothelial cells, pericyte/stromal cells, macrophages and S1 tubular 146 cells. Within 1 hour of LPS exposure, most cell types showed decreased expression of select 147 genes involved in ribosomal function, translation and mitochondrial processes such as Eef2 and Rpl genes (Fig. 5a, Supplementary Fig. 5a). This reduction peaked at 16 hours and recovered 149 by 27 hours. Concomitantly, most cell types exhibited increased expression of several genes 150 involved in inflammatory and antiviral responses such as Tnfsf9, Cxcl1, Ifit1, and Irf7. However, 151 this increase was not synchronized among all cell populations. Indeed, it occurred as early as 1 152 hour in endothelial cells, macrophages and pericyte/stromal cells, all acting as first responders. 153 In contrast, epithelial cells were late responders, with increases in inflammatory and antiviral 154 responses occurring between 4 and 16 hours. In fact, four hours after LPS administration, 155 cluster-specific GO terms were indistinguishable among the majority of cell types with 156 enrichment in terms related to defense, immune and bacterium responses (Fig. 5b). One noted 157 exception was the S3T2 cells (outer stripe S3) which did not enrich as robustly as other cell 158 types in these terms. It mostly maintained an expression profile related to ribosomes, translation 159 and drug transport throughout the sepsis timeline (Supplementary Fig. 6). Other players of 160 interest in sepsis pathophysiology such as prostaglandin and coagulation factors are described 161 in Supplementary Figure 5b. 162 At the 48-hour time point, while S1 cells partially recovered to baseline, the macrophages 163 showed increased expression of genes involved in phagocytosis, cell motility and leukotrienes, 164 broadly representative of activated macrophages (e.g. Csf1r, Lst1, Capzb, S100a4, Cotl1, 165 Alox5ap, Fig. 5a). Intriguingly, at this late time point, the pericyte/stromal cells are enriched in 166 unique terms related to specific leukocyte and immune cell types such as lymphocyte-mediated 167 immunity, T cell mediated cytotoxicity and antigen processing and presentation. This suggests 168 that the pericyte may function as a transducer between epithelia and other immune cells. 169

Sepsis alters cell-cell communication in the murine kidney 170
Therefore, we next examined comprehensively cell-cell communication along the sepsis 171 timeline. We show select examples of cell type-specific receptor ligand pairs. For example, we 172 found that S1 and endothelial cells communicate with the Angpt1 (Angiopoetin 1) and Tek (Tie2) ligand-receptor pair at baseline and throughout the sepsis timeline ( Fig. 6a-b, Supplementary  174   Fig. 7a). In contrast, C3 was strongly expressed in pericyte/stromal cells, while its receptor 175 C3ar1 localized to macrophage/DCs. This communication, present at baseline, did increase 176 along the sepsis timeline with additional players such as S1 participating in the cross talk 177 (Supplementary Fig. 7). Another strong communication was noted between endothelial cells 178 and macrophage/lymphocytes using the Ccl2 and Ccr2 receptor-ligand pair. The architectural 179 layout of these four cell types, with pericytes and endothelial cells residing between proximal 180 tubule and macrophage/DCs may explain these complex communication patterns 6 . Such 181 communication patterns among these four cell types may also explain macrophage clustering 182 around S1 tubules at later time points in sepsis as we previously reported 13 . 183 When examined comprehensively, receptor-ligand signaling progressed from a broad pattern at 184 baseline into a more discrete and specialized one 4 hours after LPS (Fig. 6c, Supplementary  185   Fig. 7b-c) The murine sepsis timeline allows staging of human sepsis 201 Finally, we asked whether our mouse sepsis timeline could be used to stratify human sepsis 202 AKI. To this end, we selected the differentially expressed genes from all cells combined (pseudo 203 bulk) for each time point across the mouse sepsis timeline (Supplementary Table 4). We then 204 examined the orthologues of these defining genes in human kidney biopsies of patients with 205 sepsis and AKI. The clinical data associated with these human biopsies did not allow further 206 stratification or staging of the sepsis timeline (Supplementary Table 5). As shown in Figure 7d, 207 our approach using the mouse data succeeded in partially stratifying the human biopsies into 208 early, mid and late sepsis-related AKI. These findings suggest that underlying injury 209 mechanisms are conserved, and the mouse timeline may be valuable in staging and defining 210 biomarkers and therapeutics in human sepsis. 211

Discussion 212
In this work, we provide comprehensive transcriptomic profiling of the kidney in a murine sepsis 213 model. To our knowledge, this is the first description of spatial and temporal transcriptomic 214 changes in the septic kidney that extend from early injury well into the recovery phase. Our data 215 cover nearly all renal cell types and are time-anchored, thus providing a detailed and precise 216 view of the evolution of sepsis in the kidney at the cellular and molecular level. 217 Using a combination of analytical approaches, we identified marked phenotypic changes in 218 multiple cell populations along the sepsis timeline. The proximal tubular S1 segment exhibited 219 significant alterations consisting of early loss of traditional function-defining markers (e.g., 220 SGLT2). Similar losses of function-defining markers along the nephron may explain the 221 profound derangement in solute and fluid homeostasis seen in sepsis. Concomitantly, we 222 observed novel epithelial expression of immune-related genes such as those involved in antigen 223 presentation. This indicates a dramatic switch in epithelial function from transport and 224 homeostasis to immunity and defense. These phenotypic changes were reversible, thus 225 underscoring the remarkable resilience and plasticity of the renal epithelium. 226 In addition, our combined analytical tools clearly identified unique subclusters within each 227 epithelial cell population (e.g., cortical S3 and OS S3). These subclusters likely represent novel 228 populations that may be in part influenced by the complex microenvironments in the kidney. It is 229 likely that such microenvironments define unique features in epithelial subpopulations such as 230 the expression of complete SARS-CoV-2 machinery in S1. 231 Similarly, we also identified unique features in immune-cell populations. For example, the 232 combined use of RNA velocity field and pseudotime analyses uncovered differences in 233 macrophage subtypes relating to RNA kinetics and cell differentiation trajectories. Of note is that 234 these subtypes only partially matched the traditional flow cytometry-based classification of 235 macrophages (e.g., M1/M2). Therefore, the use of single-cell RNA seq is a powerful approach 236 that will add to and complement our current understanding of the immune cell repertoire in the 237

kidney. 238
Additional approaches such as receptor-ligand crosstalk and gene regulatory network analyses 239 identified unique cell-and time-dependent players involved in sepsis pathophysiology. Our work points to the urgent need for defining a more accurate and precise timeline for human 250 sepsis. Such definition will guide the development of biomarkers and therapies that are cell and 251 time specific. We show evidence supporting the relevance of murine models and their 252 usefulness in staging human sepsis. These precisely time-and space-anchored data will 253 provide the community with rich and comprehensive foundations that will propel further 254 investigations into human sepsis. 255

Experimental Model and Subject Details 257
Animal model: Male C57BL/6J mice were obtained from the Jackson Laboratory. Mice were 8- biopsies were used in this study, the Institutional Review Board determined that informed 269 consent was not required. 270

Isolation of single cell homogenate from murine kidneys 272
Murine kidneys were transported in RPMI1640 (Corning), on ice immediately after surgical 273 procurement. Kidneys were rinsed with PBS (ThermoFisher) and minced into eight sections. 274 Each sample was then enzymatically and mechanically digested with reagents from Multi-275 Tissue Dissociation Kit 2 and GentleMACS dissociator/tube rotator (Miltenyi Biotec). The 276 samples were prepared per protocol "Dissociation of mouse kidney using the Multi Tissue 277 Dissociation Kit 2" with the following modifications: After termination of the program "Multi_E_2", 278 we added 10 mL RPMI1640 (Corning) and 10% BSA (Sigma-Aldrich) to the mixture, filtered and 279 homogenate was centrifuged (300 g for 5 minutes at 4°C). Cell pellet was resuspended in 1 mL 280 of RBC lysis buffer (Sigma), incubated on ice for 3 minutes, and cell pellet washed three times 281 (300 g for 5 minutes at 4°C ). Annexin V dead cell removal was performed using magnetic bead 282 separation after final wash, and the pellet resuspended in RPMI1640/BSA 0.04%. Viability and 283 counts were assessed using Trypan blue (Gibco) and brought to a final concentration of 1 284 million cells/mL, exceeding 80% viability as specified by 10x Genomics processing platform. 285 286

Single cell library preparation 287
The sample was targeted to 10,000 cell recovery and applied to a single cell master mix with 288 lysis buffer and reverse transcription reagents, following the Chromium Single Cell 3' Reagent 289 Kits V3 User Guide, CG000183 Rev A (10X Genomics, Inc.). This was followed by cDNA 290 synthesis and library preparation. All libraries were sequenced in Illumina NovaSeq6000 291 platform in paired-end mode (28bp + 91bp). Fifty thousand reads per cell were generated and 292

Single cell data processing 299
The 10x Genomics Cellranger (v. 2.1.0) pipeline was utilized to demultiplex raw base call files to 300 FASTQ files and reads aligned to the mm10 murine genome using STAR 30 . Cellranger 301 computational output was then analyzed in R (v.3.5.0) using the Seurat package v. 3.0.0.9999, 302 31 . Seurat objects were created for non-integrated and integrated (inclusive of all time points) 303 using the following filtering metrics: gene counts were set between 200-3000 and mitochondrial 304 gene percentages less than 50 to exclude doublets and poor quality cells. Gene counts were log 305 transformed and scaled to 10 4 . The top 20 principle components were used to perform 306 unsupervised clustering analysis, and visualized using UMAP dimensionality reduction 307 (resolution 1.0). Using the Seurat package, annotation and grouping of clusters to cell type was 308 performed manually by inspection of differentially expressed genes (DEGs) for each cluster, 309 based on canonical marker genes in the literature 8-10,32,33 . In some experiments, we used 310 edgeR negative binomial regression to model gene counts and performed differential gene 311 expression and pathway enrichment analyses (topKEGG, topGO, Fig. 5, Supplementary Fig.  312 5a, Supplementary Fig. 6, and DAVID 6.8 Fig. 7b. 34,35 . 313 The immune cell subset was derived from the filtered, integrated Seurat object and included the 314 Macrophage/DC (cluster 10), neutrophil (cluster 19) and lymphocyte (cluster 13) cells. Gene 315 counts were log transformed, scaled and principle component analysis performed as for the 316 integrated object above. UMAP resolution was set to 0.4, which yielded 14 clusters. The 317 clusters were manually assigned based on inspection of DEGs for each cluster, and cells 318 grouped if canonical markers were biologically redundant. We confirmed manual labeling with 319 an automated labeling program in R, SingleR 36 . 320

Analysis of regulons and their activity in the integrated single cell dataset 321
SCENIC analysis 37 was performed using the default setting and mm9-500bp-upstream-322 7species.mc9nr.feather database was used for data display. 323

Pseudotemporal ordering of single cells 324
We performed pseudotime analysis on the integrated Seurat object containing all cell types as 325 well as the immune cell subset. Cells from each of the seven time points were included and 326 were split into individual gene expression data files organized by previously defined cell type.

Spatial Transcriptomics 375
A septic mouse kidney was immediately frozen in Optimal Cutting Temperature media. A 10 µm 376 frozen tissue section was cut and affixed to a Visium Spatial Gene Expression library 377 preparation slide (10X Genomics). The specimen was fixed in methanol and stained with 378 hematoxylin-eosin reagents. Images of hematoxylin-eosin-labeled tissues were collected as 379 mosaics of 10x fields using a Keyence BZ-X810 fluorescence microscope equipped with a 380 Nikon 10X CFI Plan Fluor objective. The tissue was then permeabilized for 12 minutes and RNA 381 was isolated. The cDNA libraries were prepared and then sequenced on an Illumina NovaSeq 382 6000. Using Seurat 3.1.4, we identified anchors between the integrated single cell object and 383 the spatial transcriptomics datasets and used those to transfer the cluster data from the single 384 cell to the spatial transcriptomics. For each spatial transcriptomics spot, this transfer assigns a 385 score to each single cell cluster. We selected the cluster with the highest score in each spot to 386 represent its single cell associated cluster. Using a Loupe Browser, expression data was 387 visualized overlying the hematoxylin-eosin image. 388

Single-molecule RNA in situ hybridization 389
Formalin-fixed paraffin-embedded cross sections were prepared with a thickness of 5µm. The 390 slides were baked for 60 minutes at 60 °C. Tissues were incubated with Xylene for 5 minutes 391 x2, 100% ETOH for 2 minutes x2, and dried at room temperature. RNA in situ hybridization was Fluorescein Plus Evaluation kit (PerkinElmer, Inc) was used as secondary probes for the 396 detection of RNA signals. All slides were counterstained with DAPI and coverslips were 397 mounted using fluorescent mounting media (ProLong Gold Antifade Reagent, Life 398 Technologies). The images were collected with a LSM800 confocal microscope (Carl Zeiss). 399

Quantification and Statistical Analysis 400
No blinding was used for animal experiments. All data were analyzed using R software 401 packages, with relevant statistics described in results, methods and Fig. legends. 402

Data availability 403
Data will be deposited to NCBI GEO. The authors declare that all relevant data supporting the 404 findings of this study are available on request. 405

Code availability 406
R scripts for performing the main steps of analysis are available from the Lead contact on 407 request. 408

Additional Information 409
Correspondence and requests for resources and reagents should be directed to and will be 410 fulfilled by the Lead Contact Takashi Hato (thato@iu.edu). 411 412