Single cell RNA sequencing reveals differential cell cycle activity in key cell populations during nephrogenesis

The kidney is a complex organ composed of more than 30 terminally differentiated cell types that all are required to perform its numerous homeostatic functions. Defects in kidney development are a significant cause of chronic kidney disease in children, which can lead to kidney failure that can only be treated by transplant or dialysis. A better understanding of molecular mechanisms that drive kidney development is important for designing strategies to enhance renal repair and regeneration. In this study, we profiled gene expression in the developing mouse kidney at embryonic day 14.5 at single cell resolution. Consistent with previous studies, clusters with distinct transcriptional signatures clearly identify major compartments and cell types of the developing kidney. Cell cycle activity distinguishes between the “primed” and “self-renewing” sub-populations of nephron progenitors, with increased expression of the cell cycle related genes Birc5, Cdca3, Smc2 and Smc4 in “primed” nephron progenitors. Augmented Birc5 expression was also detected in immature distal tubules and a sub-set of ureteric bud cells, suggesting that Birc5 might be a novel key molecule required for early events of nephron patterning and tubular fusion between the distal nephron and the collecting duct epithelia.

INTRODUCTION 1 associated with cell cycle/replication (see Supplemental Table 6). We then used SCE-1 NIC [49] to gain some insight into gene regulation driving the transcriptional changes we 2 observe across pseudotime between the two nephron progenitor cell types. Figure 3g 3 depicts the activity of inferred regulatory modules across pseudotime for 31 recovered 4 transcription factors. We observe three regulatory modules of down-regulated genes, at-5 tributed to the transcription factors Egr1, Maf and Fos. These findings corroborate previ-6 ous studies showing that these transcription factors are critical regulators of gene expres-7 sion, controlling transition from a pluripotent to differentiated state in nephron progenitor 8 and human embryonic stem cells [50]. For up-regulated genes, we observe modules as-9 sociated with cell cycle-related transcription factors like E2f8, Hcf1, Ezh2 and Mybl2, 10 which have previously been implicated in specific aspects of cell cycle progression and 11 cell fate decision in stem and progenitor cells [51][52][53][54][55]. Genes making up each module are 12 provided in Supplemental Table 8. 13 Birc5 expression in the tubular interconnection zone 14 The cell cycle-related genes Birc5, Ccnd1 and Tuba1a were up-regulated in im-15 mature distal tubules (Figures 2b, 4a). Immunostaining analysis confirmed augmented  Conserved features in mouse and human podocyte development 26 In the podocyte lineage, the genes most significantly defining the cluster are Pax8, 27 Podxl and Nphs1. Podxl and Nphs1 (in combination with Synpo, Nphs2 and VEGF-A) are 28 restricted to a sub-population of mature podocytes [28,31,57], which is consistent with 29 our observations (Figures 2b and 4b). In a sub-population of early podocytes, WT1 and 1 Mafb expression has been reported to overlap with the immature marker Pax8 [27,31], 2 and is expressed in parietal epithelial cells [58], also consistent with our findings ( Figures   3   2b and 4b). For a summary of gene expression changes during podocyte development 4 see Supplemental Figure 5a and Supplemental Table 9. Similar to previous scRNA-5 seq analysis in human fetal kidney [31], we observe enrichment in the PDZ domain pro-6 teins Magi2, Slc9a3r2 and Pard3b in mature podocytes (Figure 4b). We also observe 7 podocyte-specific activity of Cldn5 (while the claudins Cldn6 and Cldn7 are expressed in 8 tubular lineages, Figure 4a). Further on, the gene Sparc (a cystine-rich matrix-associated 9 protein) and the Tissue-Type Plasminogen Activator Plat are expressed specifically in the 10 podocyte lineage (as is Robo2, a gene known to be expressed and colocalized with 11 nephrin on the basal surface of mouse podocytes [59], while the cell-cycle regulator Gas1 12 (Growth Arrest Specific 1) is expressed in undifferentiated cells and mature podocytes, 13 but less so in mature tubular cells (Figure 4b). Together, these findings further define the 14 gene expression profile of the podocyte lineage, and they suggest substantial conserva-15 tion between mouse and human developing podocytes.

16
Gene expression differences between proximal and distal tubular cells 17 In addition to their respective marker genes Pdzk1, Slc34a1, Tmem52b and Shd 18 (see above), we observe that tubular lineages express the claudins Cldn6 and Cldn7, as   Table 10 show additional genes with dif-21 ferential expression between proximal and distal tubular cells). Ly6a is a member of the 22 murine L6 family and has been reported to mark cancer and tissue-resident stem cells in 23 mice [60]; however there is no known direct human ortholog for Ly6a, and also the func-24 tion of the LU domain, which characterizes Ly6a's superfamily of proteins, is currently 25 unknown in humans [60]. 26 We note that Mep1a, Aldob and Tmem174 mark proximal tubular cells in our data For distal tubular cells, we don't observe a selectively active kinase inhibitor but 6 note that Cyclin-Dependent Kinase-like 1 (Cdkl1) is specifically expressed in this cell type.

7
These findings pinpoint lineage-specific gene expression differences between the proxi-8 mal vs. distal tubular lineages, and they point towards lineage-specific control of the cell 9 cycle across nephron progenitor differentiation.

10
Recently published scRNA-seq papers have described differences in gene expres-11 sion across a variety of proximal tubule transcripts and lncRNAs in different sexes in the 12 adult mouse kidney [62,63]. We observe that female-enriched markers, including 13 Gm4450, Lrp2, Sultd1, Aadat, Hao2 were highly expressed in our proximal tubular cluster, 14 while most of the male-enriched markers (Slc22a12, Cndp2, Cesf1, etc.) were absent or 15 expressed at low levels. This data suggests that sexually dimorphic gene expression in 16 proximal tubule may occur at or before E14.5.

17
Expression of known lineage-marker genes in unexpected cell types 18 Expression of known lineage-marker genes in unexpected cell types has been re-19 ported based on the analysis of scRNA-seq data, for example that stromal cells express 20 Gdnf [27]. Consistent with this report, we found that nephron progenitor markers (Six2, 21 Cited1, Crym) are expressed in cells in the stromal cluster, and that stromal markers are 22 present in the nephron progenitor cluster (Meis1, Foxd1, Crabp1). We also confirm that 23 Gdnf is expressed in the stromal cluster (in addition to nephron progenitor cells), and that 24 Aldh1a2 RNA is present in stromal and nephron progenitor clusters ( Figure 5). 25 We note that nephron progenitor marker genes are not homogenously expressed 26 across different stromal cell types. For instance, Cited1 is detected (five or more reads) 27 in about 7% of cortical stromal cells, but in less than 1% of other stromal cell types. We 28 find similar enrichment (expression in ~7% vs less than 1% of cells) for Six2 and Crym in 29 cortical stroma, whereas Gdnf is more modestly enriched in cortical stroma (expressed in 30 ~3% vs less than 1% of cells, respectively). We next focused on cortical stroma and 1 looked at co-expression of nephron progenitor and stromal marker genes in the same 2 cells (binary expression "on" vs. "off") and find significant positive association between 3 the expression of stromal-and nephron-progenitor genes (Fisher exact test, Table 1). 4 This analysis demonstrates widespread co-expression of nephron-progenitor and stromal 5 markers in the same cortical/stromal cells, and on average we observe higher odds ratios 6 of association for Col1a1 expression with nephron progenitor lineage-markers, compared 7 with Meis1 (Table 1). 8 Further on, the cluster we identified as distal tubular cells contains cells with a  [25][26][27]29]. We find that cell cycle activity distinguishes between "primed" 28 and "self-renewing" sub-populations of nephron progenitors. Furthermore, augmented 29 BirC5 expression occurs in immature distal tubules and a sub-set of ureteric bud cells, 30 suggesting that Birc5 might be a novel key molecule required for early events of tubular 1 fusion between the distal nephron and the collecting duct epithelia.

2
All nephron segments derive from a multipotent self-renewing nephron progenitor 3 population, which co-expresses the transcription factor Six2 and the transcriptional acti-4 vator Cited1. Previous studies have identified two sub-types of nephron progenitors, with 5 Cited1 + / Six2 + progenitors transitioning to a Cited1 -/ Six2 + primed state as the nephrogen-6 esis proceeds [15,34,65]. Recent studies using time-lapse imaging and scRNA-seq anal-7 yses have indicated, however, that the nephron progenitor compartment is more hetero-8 geneous that initially supposed [26,27,29,[66][67][68]. Moreover, differences in cell cycle 9 length within progenitors appear to play a role in the sub-compartmentalization of the 10 progenitor population [45]. In agreement, our scRNA-seq analysis shows separation of 11 nephron progenitor cells into a "self-renewing" and "primed" sub-population, both co-ex-12 pressing Six2 and Cited1, but distinguished by higher cell cycle activity in the "primed" Pax8. These findings are consistent with previous scRNA-seq analyses of developing 26 human and mouse kidneys [28,69], but are in contrast to other studies on nephron pro-27 genitor subpopulations where Cited1 expression seems to be turned off prior to the acti-28 vation of pre-tubular aggregate genes [15,34,68]. Such discrepancies might be due to 29 differences in the technical sensitivity of the methods applied in each study (scRNA-seq 30 versus immunofluorescence or in situ hybridization). They also highlight the importance 31 of further analysis to confirm whether these nephron progenitor sub-populations coincide 1 with distinct spatial domains within the developing kidney.  Birc5 (also known as Survivin) has been implicated in a number of kidney condi-29 tions, including autosomal-dominant polycystic kidney disease, acute kidney injury and observed was considered embryonic day 0.5 (E0.5). All experimental procedures were 1 performed in accordance with the University of Pittsburgh Institutional Animal Care and 2 Use Committee guidelines (IACUC protocol #17091432), which adheres to the NIH Guide 3 for the Care and Use of Laboratory Animals. 4 We harvested two kidneys at E14.5 and generated a single cell suspension using 0.05% 5 trypsin at 37°C for 10 minutes. Kidneys were mechanically dissociated with pipetting at 5 6 and 10 min. 3% fetal calf serum in PBS was added to halt the trypsin. The cell suspension 7 was filtered using a 40 µm filter and pelleted. The cells were resuspended in 90% FCS in 8 DMSO and frozen, prior to shipment to GENEWIZ Ing. Single-cell library preparation and  Briefly, we fit a mean-variance trend to the gene variances using the trendVar function 9 and identified the biological component of the total variance with decomposeVar. All 10 genes with an FDR < 0.01 and proportion of biological variance of at least 25% are con-11 sidered as highly variable genes (HVG). Principal component analysis (PCA) was then 12 performed using denoisePCA and two-dimensional representation was then derived us-13 ing runTSNE.

15
Cells were grouped into clusters using the scran R package by building a shared k-near-16 est-neighbors graph using buildSNNGraph (with use.dimred=PCA and k=25 options), fol-17 lowed by clustering with the Walktrap community finding algorithm as implemented in the 18 igraph package (https://igraph.org), cutting the graph at 10 clusters. We used the expres-19 sion of a curated list of marker genes for major components of the developing kidney (see 20 Figure 1c) to assign cluster labels. Cluster-specific markers were derived using the find-21 Markers function (Supplemental Table 1). We note that at this resolution tubular distal  Focusing on nephron progenitor and descendant cell types (termed "nephron-progenitor", 28 "mixed/differentiating" (at that point containing "distal_tubular" cells as well) , "podocytes and "proximal_tubular" in Figure 1) and requiring expression of each gene with at least 1 three counts in three cells yielded a gene expression matrix of 9,611 genes across 1,273 2 cells. Following the same procedure as before we derived a low-dimensional representa-3 tion and identified six clusters of cells, corresponding to two types of nephron progenitor 4 cells ("self-renew" and "primed"), "mixed/differentiating" cells as well as distal tubular cells 5 and proximal tubular cells (Figure 2). Cluster-specific marker genes were derived as 6 before and are reported in Supplemental Table 2. Enrichment analysis for Gene Ontol-7 ogy terms enriched between "self-renew" and "primed" and between "primed" and 8 "mixed/differentiating" (Figures 3b and 3c

18
Next, we fitted a multinomial log-linear model (using the nnet package [106]) relating 19 pseudotime with the annotated clusters. For cells with more than one annotated lineage 20 (in the "self-renew", "primed" and early "differentiating" clusters) lineage-pseudotimes 21 from slingshot were averaged. This enabled us to define NP-cells as cells with annotated 22 pseudotime less than the (pseudo)timepoint between "primed" and "differentiating" where 23 the probability of the "primed" cluster has declined to 50% (i.e., 50% probability for "dif-24 ferentiating". These cells contain all "self-renew" cells, 28 "differentiating" cells and all but 25 15 "primed" cells and were used in subsequent pseudotime analyses comparing self-26 renewing and differentiating cells. Ontology terms for genes differentially expressed between self-renewing and primed NP         and mature distal and proximal tubular clusters, and also from the ub/cd cluster. We see 4 that in the tubular distal lineage, in contrast to the proximal lineage, a significant fraction 5 of cells expresses these ub/cd marker genes.  Table 1: Marker genes that distinguish between kidney cell types. 21 For each cluster marker genes with an FDR < 0.05 are reported (see Figure 1). (see Figure 2).