Kidney formation requires the coordinated growth of multiple cell types including the collecting ducts, nephrons, vasculature and interstitium. There is a long-held belief that interactions between progenitors of the collecting ducts and nephrons are primarily responsible for kidney development. However, over the last several years, it has become increasingly clear that multiple aspects of kidney development require signaling from the interstitium. How the interstitium orchestrates these various roles is poorly understood. Here, we show that during development the interstitium is a highly heterogeneous patterned population of cells that occupies distinct positions correlated to the adjacent parenchyma. Our analysis indicates that the heterogeneity is not a mere reflection of different stages in a linear developmental trajectory but instead represents several novel differentiated cell states. Further, we find that β-catenin has a cell autonomous role in the development of a medullary subset of the interstitium and that this non-autonomously affects the development of the adjacent epithelia. These findings suggest the intriguing possibility that the different interstitial subtypes may create microenvironments that play unique roles in development of the adjacent epithelia and endothelia.
Development of the kidney relies on interactions between the metanephric mesenchyme (MM) and an epithelial structure known as the ureteric bud (UB) (Grobstein, 1955, 1953b). The MM is a heterogeneous population of cells containing at least two cell type-specific progenitor populations. A Six2/Cited1-positive nephron progenitor cell (NPC) population is located within the MM directly adjacent to the tips of the UB (Boyle et al., 2008; Kobayashi et al., 2008). NPCs undergo mesenchymal-to-epithelial transition (MET) to ultimately form the nephron, the functional unit of the kidney, which is patterned into functionally distinct segments: the glomerulus, proximal tubule, loop of Henle, distal tubule and connecting segment (Boyle et al., 2008; Kobayashi et al., 2008). There is a second, molecularly distinct progenitor population within the MM that surrounds the NPCs. These forkhead box D1 (Foxd1, formerly known as BF2)-expressing cells have been shown to give rise to a significant percentage of the stromal/interstitial cells and the renal capsule (Hatini et al., 1996; Humphreys et al., 2010; Kobayashi et al., 2014).
During development, the UB coordinates the proliferation and differentiation of nephron progenitors into the precursor of the nephron, the renal vesicle (Carroll et al., 2005; Karner et al., 2011). Reciprocally, the MM (both the NPC and interstitial populations) induces the outgrowth and branching of the UB until the UB has formed an arborized epithelial network of tubules referred to as the collecting ducts (Grobstein, 1955, 1953a). In addition to its role in regulating branching, the interstitium plays further roles in differentiation of the nephron progenitor population and patterning of the vasculature (Hatini et al., 1996; Das et al., 2013; Hum et al., 2014; Hurtado et al., 2015). Although frequently referred to in broad terms, the adult interstitial cell population includes renal fibroblasts and various smooth muscle cell types including vascular smooth muscle, pericytes, mesangial cells and the smooth muscle of the ureter and renal pelvis (Humphreys et al., 2010; Kobayashi et al., 2014). Moreover, much of the endocrine function of the kidney, including the production of renin and erythropoietin, is performed by interstitial cells (Maxwell et al., 1997; Berg et al., 2013).
Given that the interstitium has diverse roles in renal development, structure and function, it seems likely that this cell population is molecularly heterogeneous. However, extensive molecular characterization has been lacking. Here, we perform single cell RNA sequencing of Foxd1-derived interstitial cells from E18.5 mouse kidneys in order to define the heterogeneity and thus facilitate further inquiry into the development and function of these cells.
Our analysis revealed striking transcriptional heterogeneity in the renal interstitium, identifying 17 unique cellular clusters. Antisense mRNA in situ hybridization analysis demonstrated unexpected regionalization, uncovering at least 12 histologically similar but anatomically distinct domains along the cortical medullary axis. Importantly, comparison of mouse and human single cell data reveals that the interstitial heterogeneity is conserved. Analysis of transcription factor activity (regulons) showed cluster specificity, with the transcriptional regulator β-catenin active within a medullary sub-population of stroma. Genetic ablation studies showed β-catenin played a cell autonomous role in formation of this medullary interstitial region and a non-autonomous role in the development of the adjacent medullary epithelia. These findings stimulate multiple questions regarding the nature of interstitial-parenchymal crosstalk in development, physiology, homeostasis and disease of the embryonic and adult kidney.
Transcriptional analysis of the renal interstitium reveals molecular heterogeneity
Interstitial/stromal cells are present throughout the kidney, extending from the cortex to the most medullary regions of the renal papillae and surrounding the ureter. Within the interstitium, three molecularly and anatomically distinct regions of interstitial cells have previously been defined. These populations were annotated as nephrogenic, cortical and medullary interstitium based on their unique anatomical positions (Little et al., 2007). However, more detailed examination of kidneys stained with various antibodies to proteins expressed in interstitial cells suggests that the degree of heterogeneity may be underestimated (Kobayashi et al., 2014). For example, Foxd1 is expressed in a small subset of the interstitial cells cortical and lateral to the Six2-positive cap mesenchyme. Expression ceases adjacent to the renal vesicles (Fig. 1A,B). In comparison, tenascin-C is expressed in a subset of Foxd1-expressing cells lateral to the Six2-positive cap mesenchyme, but is not expressed in the cells cortical to the cap. Tenascin-C expression extends into the cells adjacent to the proximal tubule (Fig. 1C,D). Slug (also known as Snai2) expression appears to be largely, if not completely, non-overlapping with Foxd1, marking a subset of interstitial cells just medial to the cap mesenchyme and extending medially to the level of the proximal tubules (Fig. 1E,F). Acta2 is detectable in the interstitial cells adjacent to the proximal tubules but is not detectable in the medullary interstitium of the renal papillae (Fig. 1G,H). CDKN1c is expressed in the majority of, if not all, interstitial cells of the renal papilla. Lef1 is expressed in interstitial cells extending from just medullary to the renal vesicles through the entire papillae (Fig. 1I,J). However, in contrast to CDKN1c, Lef1 appears to only be expressed in a single layer of interstitial cells that lie directly adjacent to the collecting ducts (Fig. 1K,L). Tbx18 is only expressed in interstitial cells surrounding the renal pelvis and ureter (Fig. 1M). Interestingly, of all the proteins discussed, only Acta2 is expressed in the smooth muscle surrounding the vasculature and mesangial cells, both of which are Foxd1-derived interstitial cell types. Based on the spatial differences in expression of the proteins described above and recent studies, one would predict that there are may be as many as 10 distinct interstitial cell types within the kidney.
To gain a more complete understanding of interstitial heterogeneity within the developing kidney, we performed single cell RNA sequencing (scRNA-seq) on dissociated embryonic day (E) 18.5 wild-type mouse kidneys. Although previous single cell analyses have been performed on both adult and embryonic kidneys (Combes et al., 2019; Wu et al., 2019; Lindström et al., 2019; Park et al., 2018), the interstitium represents a relatively small percentage of the total number of cells and thus has been under-represented. To enrich for interstitium, we purified cells via fluorescence activated cell sorting (FACS) from Foxd1Cre;Rosa-Tomato kidneys, as Foxd1-positive cells have been shown to give rise to the majority of the renal interstitium (Kobayashi et al., 2014; Bohnenpoll et al., 2013). Further, because recent studies suggest that a sub-population of the interstitium (in particular the ureteric fibroblasts/smooth muscle) is derived from a distinct, Tbx18-positive progenitor, and thus might not be present in our isolated cells, we bioinformatically isolated interstitial cells from published datasets derived from whole kidneys (1482 total cells) and included them in our analysis (Combes et al., 2019). Unsupervised clustering was performed on all sequenced cells that met quality control standards. After pre-processing, quality control, normalization for cell cycle phase and lineage filtering, 8683 interstitial cells were selected for further analysis.
Shared nearest neighbor clustering based on similarity between gene expression profiles revealed 17 distinct clusters of cells that were identified as ‘interstitial’ (Fig. 2, Fig. S1A, Movies 1 and 2). Clusters ranged in size from 16 cells (cluster 13) to 1564 cells (cluster 10).
Importantly, the only cluster that was not derived from the Foxd1-Cre isolated cells was cluster 13, which represents the Tbx18-derived ureteric interstitium. These data suggest that our analysis includes most, if not all, renal interstitial cells, and that this population shows a high level of heterogeneity. We next sought to validate this heterogeneity in situ in order to gain insight into the significance and nature of the cellular diversity.
Renal interstitium shows spatial heterogeneity
To validate our results, we performed section in situ hybridization with over 50 differentially expressed genes (DEGs) from each cluster (Figs 3–5, Fig. S1C-K, Tables S1 and S2). Candidates were chosen based on relative abundance and any identified regionalized expression observed in publicly available databases (McMahon et al., 2008; Harding et al., 2011; Visel et al., 2004). Although all DEGs examined were expressed in the interstitium, some that were indicated as being differentially expressed between clusters did not show cell type-specific expression by in situ hybridization (Tables S1 and S2). Instead, these genes appeared to be ubiquitously expressed, suggesting that their being called a DEG was owing to differences in mRNA levels between different cell types rather than cell-type specificity. A subset of the queried genes were expressed broadly in epithelial and/or endothelial populations as well as the interstitium. Both of these classes of genes were largely excluded from further analysis. Several DEGs appeared to be enriched in the interstitium over other cell types and showed regionalized expression (see below). It is important to note that very few DEGs were detected in only a single cluster or a single cell type. Although individual examples for each cluster are presented in the main figures, anatomical assignment of clusters is based on in situ data from multiple DEGs (Table S2). All in situ data is available at the (Re)Building a Kidney website (https://www.rebuildingakidney.org).
Clusters 1-3 are highlighted by expression of Foxd1, a gene known to be expressed in the cortical interstitium (Fig. 4A,B, Fig. S1D, Table S1) and podocytes, a non-interstitial cell type. It is important to note that Foxd1-expressing cells that were disjointed from interstitial cells and expressed typical podocyte markers, e.g. Podxl and Nphs2, were excluded from analysis at the lineage-filtering stage. Thus clusters 1-3 represent Foxd1-expressing interstitial cells only.
To gain insight into the nature of clusters 1-3, we assayed the expression of several mRNAs that were differentially expressed between the clusters including netrin 1 (clusters 1 and 3), Fibin (clusters 2 and 3), Smoc2 (cluster 1) and Dlk1 (cluster 1) (Figs 3 and 4A-B′, Fig. S1C-G, Table S1). Netrin 1 (clusters 1 and 3) was uniformly expressed in the interstitial cells cortical but not lateral to the cap mesenchyme (Fig. 3C,D). Fibin (clusters 2 and 3) shows expression in a subset of the cortical interstitium as well as expression in cells that lie lateral to the cap mesenchyme (Fig. 3E,F). Smoc2 (clusters 1 and 9) and Dlk1 (cluster 1, 13 and 16) both showed mosaic expression in the cortical-most interstitium (Fig. 3G-J) as well as other, non-cortical cell types (described below). Owing to the limited resolution of in situ hybridization, we could not confirm that Fibin was expressed in distinct cell types from Smoc2 and Dlk1. However, this analysis suggests that the cortical Foxd1-expressing cells are not molecularly homogeneous and these observations suggest that clusters 1 and 3 are unique subsets within the cortical-most subset of Foxd1-expressing cells and cluster 2 represents a lateral sub-population of Foxd1-expressing cells. Thus, the region of stroma previously annotated as ‘cortical’ appears to have at least three molecularly distinct cell types.
Genes present in clusters 4-5 were expressed on the medullary side of the UB tips surrounding the newly forming renal vesicles and correlated with a region of interstitium previously referred to as the nephrogenic interstitium (Table S1). These clusters were highlighted by the expression of multiple genes involved in cell division (even after controlling for cell cycle phase) and computational assignment of cell cycle (Scialdone et al., 2015) suggests an enrichment of cells in G2/M transition (Fig. S2).
Lysyl oxidase (Lox) was enriched in clusters 6-8 and showed expression in the interstitium adjacent to the most-medullary proximal tubules (Fig. 4C-D′, Fig. S1H, Table S1). Clca3a1 (cluster 7) was also enriched in cells adjacent to the proximal tubules, although it appeared not to be expressed adjacent to the more medullary proximal tubules and its expression appeared more mosaic than Lox (Fig. 4E-F′).
Smoc2 (cluster 9) is enriched in a population of interstitium that lies just medullary to the proximal tubule in the outer medulla (Fig. 4G-H′, Fig. S1G, Table S1). Smoc2-expressing cells do not appear to expand cortically into the region adjacent to the proximal tubules.
The mRNA for proenkephalin (Penk; clusters 9-12) is most intense in a population of interstitium just medullary to Smoc2 in the outer stripe of the inner medulla (Fig. 4I-J′, Fig. S1I, Table S1). Although its expression does appear to expand into more cortical populations of cells, it does not appear to overlap with Smoc2, Lox or Clca3a1.
Wnt4 is enriched in clusters 11-12. In situ analysis confirms that, along with expression in the pre-tubular aggregates, it is highly expressed in the interstitium of the papillary region of the kidney, a region previously referred to as the medullary interstitium (Fig. 4K-L′, Fig. S1J, Table S1; Stark et al., 1994; Itäranta et al., 2006). Several Wnt/β-catenin target genes are also expressed in these clusters including Apcdd1, Axin2, Nkd1 and Nkd2, suggesting this is a region of high canonical Wnt signaling activity (Table S2).
Cluster 13 is the only population that was not derived from the Foxd1 lineage-sorted cells. All cluster 13 cells were derived from the whole kidney single cell data generated by Combes and colleagues (Combes et al., 2019). Of note, cluster 13 was composed of only 16 cells, most likely a reflection of the relatively small number of interstitial cells included in the whole organ studies. As expected, genes present in cluster 13 (e.g. Dlk1 and Tbx18) were predominantly expressed in the interstitium adjacent to the ureter (Figs 1M and 5A-B, Fig. S1C, Table S1).
Clusters 14-16 showed a high degree of overlap in gene expression. All three clusters expressed several genes that are accepted markers of perivascular cells. To determine whether the clusters represented unique perivascular cell types, we performed in situ hybridization with cluster-specific DEGs. Akr1b7 and Ren1 are both enriched in cluster 14 over other clusters and show expression in vascular smooth muscle as well as cells within the juxtaglomerular apparatus (Fig. 5C-D, Fig. S1B, Table S1 and data not shown). Dlk1 (clusters 15 and 16 along with clusters 1 and 13) was detectable in the glomerular mesangial cells (Fig. 5A-B, Fig. S1C, Table S1). Igf1 was enriched in cluster 15 and showed expression in a subset of interstitial cells in the papillary region of the kidney, similar to the anatomical location assigned to cluster 10 (Fig. 5E-F). Based on these observations, along with gene set enrichment analysis (Sergushichev, 2016 preprint) based on automated text-mining of protein-cell type associations from the biomedical literature (Fig. S3; Santos et al., 2015), we suggest that cluster 14 represents vascular smooth muscle, cluster 15 represents a papillary pericyte population and cluster 16 represents mesangial cells. However, given the molecular similarity between these three cell types, it is possible that greater sequencing depth/clustering will alter these annotations.
Cluster 17 contained many genes that were also DEGs in other clusters. Although we did identify a specific DEG list for cluster 17, the most differentially expressed genes encoded mitochondrial mRNAs and we have not been able to produce strong signal with any of these probes (Table S1).
In all, analysis of the transcriptome of single interstitial cells generated 17 clusters of cells that have been assigned to at least 12 anatomically distinct cell types.
Pseudotime analysis suggests multiple distinct developmental trajectories within the interstitium
Although interstitial heterogeneity is not completely unexpected, the heterogeneity in the histologically indistinguishable cells spanning the cortical-medullary axis was surprising. A simple explanation of the identity of these cells is that they merely represent different stages in a linear developmental trajectory from a cortical-most progenitor to a more medullary, differentiated cell type. To gain insight into the lineage relationship of the different clusters, we employed the complementary approaches of RNA velocity (La Manno et al., 2018) and pseudotemporal ordering (Haghverdi et al., 2016). Briefly, RNA velocity is determined by modeling the relationship between the unspliced and spliced forms of a transcript, with the reasonable assumption that a cell actively transcribing a gene will have a higher ratio of unspliced to spliced transcripts and so a higher ‘RNA velocity’. Conversely, cells that are not actively transcribing a gene (but still maintain its expression) will have a lower ratio of unspliced to spliced transcripts, as the majority of mRNA will represent processed transcripts (lower RNA velocity). Using single genes as representative examples, this analysis indicates that Lef1 is actively transcribed in cluster 10 and stabilizes in clusters 11 and 12, suggesting that clusters 11 and 12 represent cells derived from cluster 10 (Fig. S4A,B). In comparison, Smoc2 is actively transcribed and processed in cells present within cluster 9 but not in any other clusters (other than the Foxd1-positive progenitor cells), suggesting cluster 9 represents a terminal differentiation state (Fig. S4C,D). Extending this model to all genes and all cells within the purified interstitium, we created an RNA velocity field to predict transitions between clusters (Fig. S4E). It is important to note that the kidney capsule, which is derived from Foxd1-expressing cells (Kobayashi et al., 2014), was removed from our kidneys before cell sorting and thus is not included in our analysis. Therefore, we cannot comment on the derivation of this cell type at this time. We also ultimately excluded the ureteral smooth muscle cells (cluster 13) from this analysis as they did not show a relationship with any of the other clusters, as expected from previous work indicating that they are derived from a different lineage than the rest of the interstitium (Bohnenpoll et al., 2013).
Simulation of a reverse Markov process with RNA velocity-based transition probabilities (please refer to Materials and Methods for a more complete description) allowed us to identify a group of cells that likely represented the ontological parent of the analyzed interstitial cells (Fig. 6A). After selecting one such parent or ‘root’ cell, we were able to pseudotemporally order the cells. Annotating the diffusion map with pseudotime demonstrated that the parent/root population ramifies into multiple distinct branches that develop simultaneously rather than following a single linear differentiation trajectory (Fig. S4F).
The integration of the RNA-velocity and pseudotime analyses generates a model wherein at least two of the distinct clusters of Foxd1-positive ‘progenitor’ cells undergo independent branching events to give rise to the cycling/proliferating cells and at least six different trajectories that we refer to as the inner medullary fibroblast, outer medullary fibroblast, vascular smooth muscle, proximal tubule interstitium, pericyte and mesangial subtypes (Fig. 6B). The vascular smooth muscle cluster appears more distantly related than the other cell types, consistent with its distinct differentiation profile (Fig. 6B, Movie 3). Interestingly, cluster 6 appears to be located at a bifurcation or trifurcation point that gives rise to clusters 14, 15 and 16. These data suggest that the distinct interstitial subtypes do not represent transient stages in a linear developmental trajectory. Instead, these cells appear to represent previously undescribed, anatomically distinct interstitial cell types of distinct lineages.
Although frequently referred to as a homogeneous population of cells that functions primarily as a scaffold, several recent studies have shown that, in various systems, the interstitium sets up unique microenvironmental niches that direct tissue development and/or maintenance and can contribute to pathological conditions (Benias et al., 2018; Liotta and Kohn, 2001; Yu and Scadden, 2016; Greicius and Virshup, 2019; Driskell and Watt, 2015; Driskell et al., 2013; De Palma et al., 2017; Hepler et al., 2018). Indeed, within the kidney, studies have shown that the interstitium is important in numerous developmental processes affecting both the epithelia and endothelia (Hatini et al., 1996; Das et al., 2013; Hum et al., 2014; Hurtado et al., 2015). We next sought to re-visit previous work in light of our current findings.
Interstitial pattern affects development of the renal parenchyma
Several transcriptional regulators have been shown to be expressed in and necessary for the development of the stroma. To determine whether any of these factors were active in specific stromal subtypes, we employed the SCENIC package to reconstruct gene regulatory networks and measure their activity within cells and so within the cell parent clusters (Fig. 7; Aibar et al., 2017). Strikingly, certain transcriptional signatures (aka regulons) were active in distinct clusters or groups of clusters while inactive in others. Interestingly, we found that the Lef1 regulon was predominantly active in clusters 8-12 (Fig. 8A), which represents much of the medullary fibroblast population (Fig. 5). Previous studies have shown that inactivation of β-catenin in the interstitial progenitor population (using Foxd1Cre) leads to defects in the expression of several genes within the interstitium, especially β-catenin targets (Yu et al., 2009). To test whether β-catenin is necessary for the formation of specific interstitial subtypes (versus a general interstitial defect), we assessed the expression of regionalized genes in Foxd1Cre;catnb−/flox kidneys. Expression of genes normally expressed in clusters 1-7 (clusters in which β-catenin is not active) appeared unaffected or even expanded in Foxd1Cre;catnb−/flox mutants (Fig. 8B-D,H-J). In contrast, genes expressed in clusters 8-12 were markedly reduced or undetectable (Fig. 8E,F,K,L) in mutants. The absence of medullary interstitium (clusters 8-12) correlated with a severe deficit in the formation of the epithelia (including the loop of Henle) that lay adjacent to these zones, whereas the epithelia that lay adjacent to the unaffected cortical interstitium (including the proximal tubules) appeared to form normally, as previously reported by Yu et al. (2009) (Fig. 8G,M).
Human fetal kidney interstitium shows a similar degree of heterogeneity to the mouse
Recently, several groups have employed scRNA-seq on human fetal kidney at different stages (Lindström et al., 2018b). While analyzing the nephrogenic region of fetal kidneys, which contains cells from multiple lineages, only five unique interstitial sub-clusters were identified (Lindström et al., 2018b). Thus, we wondered whether extensive interstitial heterogeneity is unique to the mouse.
To understand whether interstitial heterogeneity is a generalizable phenomenon between these two species, we reanalyzed previously published (Lindström et al., 2018a; Tran et al., 2019) week 17 human fetal kidney scRNA-seq data. After identifying the major cell populations within the data (epithelia, endothelia, leukocytes, etc.), we bioinformatically isolated the cells defined as interstitium (see Materials and Methods). We then deployed the same clustering technique used to cluster mouse interstitium on the interstitium of the human fetal kidney. We found that the interstitium of the cortical region of human fetal kidney segments into 13 molecularly distinct clusters (Fig. 9A, Movie 4). Although the human dataset was limited to more cortical cell populations, we were able to identify nearly as many unique clusters as found in whole embryonic kidney. When the cluster assignments from E18.5 mouse interstitial cells were mapped onto the 17-week human interstitial cells, the majority were found to be represented in the human data, although they did not resolve as clearly (Fig. 9B). There are likely to be several factors underlying this imperfect alignment, including incomplete sampling in the human data, divergence between developmental time scales and time points between mouse and human, and fewer cells available for analysis in the human dataset resulting in less resolving power. Evidence of the sampling bias is evident in that we did not identify any human cells that were analogous to cluster 13, the most medullary interstitium. These data indicate that renal interstitial heterogeneity is a generalizable characteristic between mouse and human.
Although known to play a role in providing physical support, growing evidence from multiple systems indicates that interstitial cells play active roles in tissue development, maintenance and disease. Within the developing kidney, non-autonomous roles for the interstitium have been identified in UB branching, nephron differentiation and blood vessel formation (Das et al., 2013; Hatini et al., 1996; Hum et al., 2014; Hurtado et al., 2015). Mechanistic insight into how the interstitium has so many distinct functions has been hampered by a poorly defined transcriptome. Here, using scRNA-seq combined with in situ hybridization, we have generated a map of interstitial gene expression in the E18.5 kidney. Our analysis shows a previously unappreciated level of heterogeneity in the interstitium of the developing mouse kidney. Analysis of human embryonic interstitial cells show correspondence in heterogeneity found in the embryonic mouse.
Although previous work characterizing the heterogeneity of the developing mouse kidney uncovered four interstitial clusters (1482 cells) (Combes et al., 2019), our informed analysis identified 17 distinct clusters that we were able to spatially resolve into at least 12 anatomically distinct subtypes. Additional analysis will be required to determine whether the five remaining clusters represent additional unique cell types or whether the data is currently overclustered. Indeed, we were not able to obtain specific signal from in situ hybridization with any DEGs specific to cluster 17, a population enriched for mitochondrial genes. Thus, it is possible that this cluster represents an artifact of the dissociation protocol. Alternatively, although immunostaining identified an interstitial cell type associated with the collecting ducts, our clustering did not conclusively identify these cells. It is possible that they lie within one of the medullary populations (12 or 13) or alternatively, that rather than being overclustered, our data is underclustered. We think it is likely that higher resolution techniques (e.g. single cell resolution in situ hybridization, antibody staining or reporter gene generation from cluster-specific DEGs) will resolve additional unique cell types. For example, there are three clusters containing genes expressed in the interstitium surrounding the proximal tubules. Although there do appear to be spatial differences in the expression of some proximal tubule interstitial genes, without cellular resolution we cannot at this point determine whether these cells represent unique cell types. The fact that one of the populations, cluster 6, shows similarity to and is predicted to be the parent of clusters 14, 15 and 16, suggests that this cell type may represent a spatially distinct mural cell or perhaps even a mural cell progenitor.
Unexpectedly, we found at least three molecularly distinct clusters within the Foxd1 expression domain. This observation raises the question of the identity of the true interstitial progenitor cell. One possibility that is supported by our RNA velocity analysis, is that rather than containing a single multipotent progenitor cell, the Foxd1 domain is comprised of several lineage-restricted progenitors that give rise to the renal interstitium and capsule. Of note, we also identified two clusters, clusters 4 and 5, that appear to represent cells medullary to clusters 1-3, that differentially express many genes important for DNA-replication and are predicted to represent cells undergoing a G2/M transition. We speculate that these clusters represent a transient amplifying cell within the interstitium that correlates with the growth of multiple cell types in the kidney including the collecting duct and nephron progenitors. Testing these hypotheses will require more detailed molecular characterization, as well as lineage tracing with cell type-specific Cre drivers.
Interestingly, our analysis identified several zones of histologically indistinguishable but molecularly distinct fibroblast-like cell types that occupy unique spatial locations along the cortical medullary axis. These cell types map to distinct anatomical regions in the adjacent parenchyma, raising the intriguing possibility of stromal-epithelial cross talk that is involved in the patterning and/or differentiation of the kidney parenchyma. By analyzing our transcriptomic data, we were able to identify signaling pathways unique to the distinct interstitial clusters. Re-analysis of β-catenin mutant interstitium reveals a unique role for this factor in the development of the papillary stroma, which secondarily affected the development of the adjacent epithelia, supporting the hypothesis that the distinct interstitial zones represent unique microenvironments that influence development of the adjacent epithelia and endothelia. Although several studies have revealed cell autonomous mechanisms underlying nephron patterning, it is possible that distinct interstitial subpopulations produce factors that interact with intrinsic pathways to assure proper position and relative size and spacing of the nephron segments with the adjacent collecting ducts and associated renal vasculature and/or promote maturation. Such a phenomenon would not be unique, as previous work has shown that mesenchymal cells play instructive roles in patterning the adjacent epithelia during development of various organ systems including the vertebrate gut tube (Golosow and Grobstein, 1962; Kim et al., 2011; Le Guen et al., 2015; McCulley et al., 2015). RNA-seq data reveals multiple growth factors, small molecules, extracellular matrix components and metabolites that appear to be regionally produced. Further genetic analysis will be required to test specific roles. In addition, knockout data exist for several interstitial DEGs (Table S1) including Pbx1 (Hurtado et al., 2015; Schnabel et al., 2003), Foxd1 (Hatini et al., 1996; Levinson et al., 2005) and Wnt4 (Stark et al., 1994; Itäranta et al., 2006), and defects in development of renal parenchyma have been described. Re-characterization of these and other mutants taking interstitial heterogeneity into consideration will be informative.
Finally, as the adult kidney shows exquisite patterning along the cortical/medullary axis, it will be of great interest to determine whether a similar degree of heterogeneity and pattern exists in the adult interstitium and how this pattern correlates to normal anatomy, physiology, homeostasis, injury, regeneration and disease.
Given the growing evidence of the essential nature of the interstitium in multiple organ types and cellular processes (Benias et al., 2018; Liotta and Kohn, 2001; Yu and Scadden, 2016; Greicius and Virshup, 2019; Driskell and Watt, 2015; Driskell et al., 2013; De Palma et al., 2017; Hepler et al., 2018), a similar analysis of interstitial heterogeneity in different organ systems at different stages may reveal that the role of the interstitium in patterning and morphogenesis is a generalizable principle.
MATERIALS AND METHODS
All animals were housed, maintained and used according to protocols approved by the Institutional Animal Care and Use Committees at the University of Texas Southwestern Medical Center and following the guidelines from the NCI-Frederick Animal Care and Use Committee. For each experiment, female mice of 7-8 weeks of age were crossed with a male of 9-10 weeks of age. Plugs were checked and the embryos were collected at the desired time points for further analysis. Noon of the day on which the mating plug was observed was designated E0.5. The following mice were used in the studies described: Foxd1Cre (JAX Stock #012463), Rosa26Tomato (JAX Stock #007909), Rosa26DTA (JAX Stock #006331), Rosa26YFP (JAX Stock #006148), catnb null and catnb flox (Brault et al., 2001). Oligonucleotide primers used for genotyping can be found in Table S3.
In situ hybridization
For section in situ hybridization, kidneys isolated at specific stages were fixed overnight in 4% paraformaldehyde (in PBS) at 4°C and cryopreserved in 30% sucrose. Tissues were frozen in OCT (Tissue Tek) and sectioned at 10 μm. Sections were subjected to in situ hybridization as previously described (Karner et al., 2011). Antisense RNA probes against Foxd1, Lox, Smoc2, Penk, Wnt4, Lef1 and Tgfb1i1 were linearized and transcribed as previously described. Plasmids were unavailable for Dlk1 and Akr1b7; thus, single-stranded DNA for each gene was purchased with the T7 RNA polymerase binding site in the reverse orientation added to the 3′ end of the gene sequence. These probes were made through RNA transcription of these single-stranded DNA gblocks using T7 polymerase.
Histology, immunohistochemistry and immunocytochemistry
Kidneys isolated at birth were formaldehyde fixed and paraffin embedded. Sections (5 μm) from paraffin-embedded kidneys were subjected to Hematoxylin and Eosin staining. For immunohistochemistry, fixed kidneys were embedded in OCT and sectioned on a cryostat (10 μm). Frozen sections were washed with PBS and blocked with 5% serum for 1 h at room temperature and incubated with primary antibodies at 4°C overnight. After primary incubation, sections were washed and incubated with HRP-tagged secondary antibodies for 1 h at room temperature. Further, signal was detected with tyramide amplification. Slides were washed and re-stained with additional markers according to the above-mentioned immunohistochemistry protocol. Slides were then mounted with Vectashield and images were captured with a Zeiss LSM500, Zeiss LSM700 or a Nikon A1R confocal microscope. The details for the antibodies used can be found in Table S4. Figures provided are representative images from three kidneys from at least two different litters.
Preparation of single cells from mouse renal interstitium
E18.5 Foxd1Cre;Rosa26Tomato mouse kidneys were dissected in cold PBS without calcium or magnesium. Kidneys were cleaned and the adrenal gland, capsule and ureters were removed. Kidneys were washed in Hanks’ Balanced Salt Solution for 2 min at 37°C then minced using two razorblades on ice. The minced kidneys were digested for 8 min at 37°C in 2 ml of 0.25% w/v collagenase A/1% w/v pancreatin (Sigma-Aldrich; 101378001, P1750), with manual dissociation via pipetting through a P1000 tip every 2 min. No more than four pairs of kidneys were dissociated in 2 ml enzyme digest. Digestion was inactivated by adding 125 µl serum. After pelleting the cells at 400 g for 5 min, cells were resuspended in 1 ml of AutoMACS Running Buffer (Miltenyi Biotec, 130-091-221) and passed through a 30 µm pre-separation filter (Miltenyi Biotec, 130-041-407). Filters were immediately washed with 500 µl AutoMACS Running Buffer. Cells were resuspended in 500 µl AutoMACS Running Buffer and filtered at least two more times through a cell-strainer cap (Falcon, 352235) attached to a 5 ml polypropylene round bottom tube (Globe Scientific, 110428).
10x genomics single-cell library preparation
A suspension of 10,000 Foxd1Cre;Rosa26Tomato-positive cells, isolated via FACS, was run on a chromium 10× Single Cell Chip (10x Genomics). Libraries were prepared using Chromium Single Cell Library kit V2, and sequenced on an Illumina NextSeq using 75 bp paired-end sequencing. Sequencing resulted in an average of 12,000 reads/cell and 3000 genes/cell. The single cell data presented in this manuscript has been deposited in the GEO and (Re)Building a Kidney databases.
Single cell data analysis
When analyzing the single cell data collected from the experiments outlined above, we included, where possible, the dataset generated by Combes et al. (GSE108291; Combes et al., 2019). Each batch was processed independently using the scran Bioconductor package (Lun et al., 2016b). Unfiltered feature-barcode matrices were generating by running the CellRanger count pipeline. Cells were called from empty droplets by testing for deviation of the expression profile for each cell from the ambient RNA pool (Lun et al., 2018). Cells with large mitochondrial proportions, i.e. more than 3 mean-absolute deviations away from the median, were removed. Cells were pre-clustered, a deconvolution method was applied to compute size factors for all cells (Lun et al., 2016a) and normalized log-expression values were calculated. Variance was partitioned into technical and biological components by assuming technical noise was Poisson-distributed and attributing any estimated variance in excess of that accounted for by a fitted Poisson trend to biological variation. The dimensionality of the data set was reduced by performing principal component analysis and discarding the later principal components for which the variance explained was less than variance attributable to technical noise.
Masking of biological effects by expression changes due to cell cycle phase were mitigated by blocking on this covariate. The cell cycle phase was inferred using the pair-based classifier implemented in the cyclone function of scran. Corrected log-normalized expression counts were obtain by calling the removeBatchEffect from the limma (Ritchie et al., 2015) Bioconductor package with a design formula including G1 and G2M cell cycle phase scores as covariates.
A single set of features for batch correction was obtained by computing the average biological component of variation across batches and retaining those genes with a positive biological component. The batches were rescaled and log-normalized expression values recomputed after the size factors were adjusted for systemic differences in sequencing depth between batches. Batch effects were corrected by matching mutual nearest neighbors in the high-dimensional expression space (Haghverdi et al., 2018). The resulting reduced-dimensional representation of the data was used for all subsequent embeddings including t-SNE and UMAP.
Cells were clustered by building a shared nearest neighbor graph (Xu and Su, 2015) and executing the Walktrap algorithm (Pons and Latapy, 2006). Differential gene expression analysis was performed using the two-part generalized linear model that concurrently models expression rate above background and expression mean implemented in MAST (Finak et al., 2015). A one-versus-all strategy was employed comparing each cluster with all other identified interstitial clusters.
Gene sets for enrichment analysis were obtained from the TISSUES Text-mining Tissue Protein Expression Evidence Scores dataset (http://amp.pharm.mssm.edu/Harmonizome/; Santos et al., 2015). The gene sets were filtered to include only those genes with a standardized value greater than 1. Enrichment analysis was performed using the fgsea (Sergushichev, 2016 preprint) Bioconductor package.
Week 17 human fetal kidney scRNA-seq data were obtained from the Gene Expression Omnibus under series accession numbers GSE112570 and GSE124472 (Lindström et al., 2018a; Tran et al., 2019). Cluster assignments were transferred from the E18.5 mouse kidney dataset to the 17-week human fetal kidney dataset using a neural network classifier constructed using the Tensorflow system (Nagalakshmi and Yu, 2015). Graph regularization (Bui et al., 2017 preprint) considering eight shared nearest neighbors of each cell was used during training of a sequential network with two hidden layers, each containing 1024 hidden nodes, to classify the expression profiles of cells from the mouse dataset by cluster. Ortholog-mapped expression profiles of human cells were input to the classifier to assign each cell in the human dataset to the mouse cluster with greatest similarity. Expression profiles were cosine-normalized before training or prediction. Only genes with a biological-to-technical variance ratio greater than zero were utilized for classification.
RNA velocity was calculated to analyze the dynamic relationships between identified cell states (La Manno et al., 2018). The analysis was limited to the replicate single cell datasets produced for this paper as BAM files were not available for the Combes, et al. dataset (Combes et al., 2019). Read counts were partitioned between spliced, unspliced and ambiguous sets by calling velocyto run10x with an mm10 repeat mask obtained from the University of California, Santa Cruz genome browser and the genome annotation file that came prepackaged with CellRanger. The results of previously conducted cell filtering and feature selection were applied to these expression matrices. Normalization, principal component analysis, k-nearest neighbor smoothing, gamma fit, extrapolation, Markov process modeling and projection onto pre-defined low-dimensional embeddings were executed by calling the relevant functions from the velocyto.analysis module.
To reconcile clustering with trajectory inference, we performed partition-based graph abstraction (PAGA) to model the connectivity between clusters (Wolf et al., 2019). At cluster resolution, the edge-score threshold was varied until the graph was decomposed into connected components. Evaluation of marker gene expression within the components permitted assignment of components to the categories of epithelium, leukocyte, erythrocyte or interstitium. In many cases, RNA velocity allowed for assigning a direction to the edges between clusters, indicating a tendency of transition between the clusters.
Gene regulatory network reconstruction and measurement of regulon activity within each, was conducted using SCENIC (Aibar et al., 2017). All cells of acceptable quality were used for network inference. The regulon activity was binarized using the following strategy. An attempt was made to model the activity of the regulon as a mixture of two normal distributions using the mixtools R package. If such a model could be fit using expectation maximization, then those cells for which the regulon activity were assigned higher probability by the distribution with the greater mean were assigned a binarized regulon activity of one, and zero otherwise. If a two-component mixture of normal distributions could not be fit, then a β-distribution was fit to the activity of the regulon and cells for which the activity was >1 mean absolute deviation above the mean were assigned a binarized regulon activity of one and zero otherwise.
The authors would like to thank Ondine Cleaver, Denise Marciano and Phoebe Carter for their insight in preparing this manuscript.
Conceptualization: A.R.E., A.D., T.J.C.; Methodology: A.R.E., C.P.C., A.D., G.C.H., T.J.C., D.W.S., D.A.; Software: C.P.C.; Validation: A.R.E., C.P.C.; Formal analysis: A.R.E., C.P.C.; Investigation: A.R.E., C.P.C., A.D., M.P., A.M., D.A., K.A.D.; Resources: T.J.C.; Writing - original draft: A.R.E., T.J.C.; Writing - review & editing: A.R.E., C.P.C., T.J.C.; Visualization: A.R.E., C.P.C., A.D., M.P., K.A.D.; Supervision: G.C.H., D.W.S., K.A.D., T.J.C.; Project administration: A.R.E., G.C.H., D.W.S., K.A.D., T.J.C.; Funding acquisition: T.J.C.
This work was supported by the National Institutes of Health (R01DK095057, R01DK080004, R24DK106743, RO1DK090127 to T.J.C. and F31DK122670 to A.R.E.) and a fellowship from the University of Texas Southwestern Medical Center Hamon Center for Regenerative Science and Medicine (to A.R.E.). This work was also supported by the University of Texas Southwestern Medical Center George O'Brien Kidney Research Core (DK079328). Deposited in PMC for release after 12 months.
To increase rigor, reproducibility and transparency, raw image files and other data generated as part of this study were deposited into the (Re)Building a Kidney consortium database (https://doi.org/10.25548/16-WYG2). scRNA-seq data has been deposited in GEO under accession number GSE155794.
Peer review history
The peer review history is available online at https://dev.biologists.org/lookup/doi/10.1242/dev.190108.reviewer-comments.pdf
The authors declare no competing or financial interests.