Cellular heterogeneity in the ureteric progenitor niche and distinct profiles of branching morphogenesis in organ development

Epithelial branching morphogenesis shapes the development of a wide variety of organ systems in invertebrate and vertebrate embryos (Ochoa-Espinosa and Affolter, 2012; Iber and Menshykau, 2013). The airways of the Drosophila tracheal system and mammalian lung, and the urine-transporting network of the collecting duct system of the mammalian kidney, are particularly well-studied examples (Affolter and Caussinus, 2008; Costantini, 2012; Little and McMahon, 2012). In each, branching growth generates a complex tubular network for fluid transport laying down an organic-specific pattern of cellular diversity with a characteristic system-specific morphology. Branching within each system is stereotypical to varying degrees (Metzger et al., 2008; Short et al., 2014). In the salivary gland, a distinct process of epithelial clefting generates branch tips (Patel and Hoffman, 2014). Despite differences in the cellular responses and the interplay with adjacent cell populations, common regulatory themes have been identified in multi-organ control of branching growth, most notably fibroblast growth factor (FGF) pathway signaling (Trueb et al., 2013; Affolter and Caussinus, 2008).

In the mammalian kidney, active development initiates with the swelling and outgrowth of a subset of cells from the nephric duct, the ureteric bud (UB), into a pre-specified population of metanephric mesenchyme (MM) (McMahon, 2016). The cells of the UB generate the highly branched ureteric epithelium of the collecting duct network to which MM-derived nephrons attach to form the complete tubular network of the kidney's filtration and transport system. Outgrowth of the UB, and subsequent branching morphogenesis of the ureteric epithelium, involves a complex regulatory interplay among the collecting duct progenitors within the ureteric branch tips (UBTs) and the adjacent mesenchymal progenitors for nephron and interstitial cell lineages within the MM (Little and McMahon, 2012). UBTs lay down a complexity of epithelial cell types over the course of development through the differentiation of cells emerging from the UBT progenitor pool (Costantini, 2012). Together, the interactions within the UBT niche coordinate growth and differentiation of distinct cell lineages to generate a functional organ of appropriate size, shape and cellular diversity.

In the mammalian lung, mesenchymal-epithelial interactions drive an initially invariant pattern of branching growth at lung branch tips (LBTs) controlling airway growth, differentiation and morphogenesis (Herriges and Morrisey, 2014; Rock and Hogan, 2011; Metzger et al., 2008). Unlike the mammalian kidney where UBT-derived Wnt signals trigger the transformation of mesenchymal nephron progenitors to epithelial nephrons responsible for much of the mass and functional role of the kidney, the lung mesenchyme remains mesenchymal and shows a more limited differentiation (McCulley et al., 2015). Thus, each organ shows common developmental themes, such as the mesenchyme-mediated control of epithelial growth and branching, and the maintenance and controlled commitment of UBT- and LBT-localized progenitors, that might be underpinned by shared regulatory processes. However, the different morphology of UBT- and LBT-derived epithelial cells, the distinct cell types generated by these populations, and the different properties and responses of their adjacent mesenchyme populations predict organ-specific mechanisms within each tip population. Active branch outgrowth is not the only mechanism of branching morphogenesis. In the salivary gland, salivary branch tips (SBTs) arise as epithelial clefts that subdivide the epithelium into end buds and ducts (Patel and Hoffman, 2014).

UBT, LBT and SBT development is stimulated by mesenchyme-derived signals that enhance and maintain high proliferative activity within tip progenitor populations. In the kidney, MM production of glial-derived neurotrophic growth factor (GDNF) initiates RET/GFRA1 receptor pathway signaling in the ureteric bud (Costantini and Shakya, 2006). RET encodes a receptor tyrosine kinase and a variety of studies have shown a larger number of receptor tyrosine kinases, including FGFR, MET, EGFR and VEGFR2 (KDR), contribute to varying degrees at different stages of ureteric epithelial branching (Walker et al., 2016; Tufro et al., 2007; Ishibe et al., 2009). In the lung and salivary gland, mesenchyme-derived FGF signaling, predominantly mesenchyme-derived FGF10 acting through FGFR2, plays the central role in branching morphogenesis (Min et al., 1998; Sekine et al., 1999; Arman et al., 1999; De Moerlooze et al., 2000; Ohuchi et al., 2000). GDNF and FGF pathways converge on ETV- and SOX-regulated transcriptional pathways in both UBT and LBT progenitors though organ-specific differences in the use of ETV and SOX family members are observed (Reginensi et al., 2011; Lu et al., 2009; Herriges et al., 2015; Zhu et al., 2012). UBTs and LBTs also differ in the signals they produce: for example, Wnt11 in UBTs and sonic hedgehog in LBTs; both signals target the underlying mesenchyme in organ-specific reciprocal signaling networks between epithelium and mesenchyme to regulate branching processes (Majumdar et al., 2003; Pepicelli et al., 1998).

Importantly, genetic analyses in mouse and human have highlighted conserved roles within each organ for specific regulatory actions. In particular, congenital anomalies of the kidney and urinary tract (CAKUT) comprise a relatively common syndrome (20-30% of all congenital malformations) (Brown et al., 1987; Queisser-Luft et al., 2002). CAKUT is associated with mutations in RET, BMP4 and ROBO2, and a number of other genes that are known to alter branching growth of the ureteric network in the mouse kidney (Nicolaou et al., 2015).

To gain a more comprehensive understanding of epithelial branch tip regulation, we have extended previous studies using micro-dissection and microarray techniques to identify tip-specific regulatory factors in the developing mammalian kidney (Schwab et al., 2003; Stuart et al., 2003; Schmidt-Ott et al., 2005). The screen performed here identified a large number of previously unknown UBT-enriched genes, as well as a large gene set enriched in maturing regions of the non-branching ureteric epithelial network. Temporal and spatial analysis of their expression in the mammalian kidney confirmed UBT enrichment and defined spatial subdomains within UBTs and conservation of UBT expression between mouse and human. Comparative studies with the lung and salivary gland highlight gene sets that are UBT specific and those enriched within epithelial tips of multiple organ systems. The stratified gene sets identified here provide a rich resource for mechanistic investigation of branching programs within the developing kidney and across other mammalian organs. Furthermore, the gene sets highlight novel candidates to explore in the context of human congenital anomalies and diseases. Mechanistic analysis of the novel gene sets is likely to shed new light on how branch tip programs regulate the size, cellular diversity and distinct morphologies of mammalian organ systems.

To identify genes that are differentially expressed in tip and stalk regions of the developing ureteric epithelium of the mouse kidney, we utilized two fluorescently labeled transgenic mouse strains: Wnt11-RFP (Harding et al., 2011) and Hoxb7-GFP (Srinivas et al., 1999). Wnt11, a target of RET/GDNF signaling in kidney development (Pepicelli et al., 1997; Majumdar et al., 2003), is expressed within the tip population whereas Hoxb7 is expressed throughout the entire UB epithelium (Fig. 1A-C). Kidneys from a single litter were pooled together for each biological replicate from embryonic day (E) 16.5 day Wnt11-RFP; Hoxb7-GFP embryos and dissociated to single cells; tip (RFP⁺, GFP⁺) and stalk (RFP⁻, GFP⁺) enriched cell fractions were separated by fluorescence-activated cell sorting (FACS), and RNA sequencing (RNA-seq) was performed on three biological replicates of each cell population. Primary analysis of RNA-seq data mapped 37,149 transcripts. To identify more specifically RNAs strongly enriched in each population, we focused on only those genes with a reads per kilobases per million (RPKM) value greater than ten in at least one sample, and at least a 5-fold expression differential between the UBT and stalk fractions (Fig. 1D).

This analysis generated a stringent list of 116 unique tip-enriched and 393 unique stalk-enriched genes (Tables S1 and S2, Fig. S1). This tip-enriched gene list includes many genes previously characterized for UBT activity, including Ret, Gfra1, Wnt11, Etv4, Etv5 and Sox8 (for a review, see O'Brien and McMahon, 2014). Gene ontology (GO) analysis of the UBT-enriched set recovered terms associated with signaling pathways and structural aspects known to be important in UBTs, and more generally in branching morphogenesis (Fig. 1E). For instance, the top GO terms for biological processes are associated with well-established signaling pathways in the embryonic kidney, such as Wnt and MAP kinase (MAPK) pathways. Moreover, the GO terms for cellular components are linked to the extracellular matrix. During branching morphogenesis, matrix remodeling and degradation is likely to facilitate branching events (Kim and Nelson, 2012). For molecular functions, key actions such as heparin binding and tyrosine kinase activity relate to FGF signaling and downstream signal transduction pathways.

As expected, given the FACS approach and high levels of expression in UBTs, the most highly enriched UBT gene was Wnt11 (Table S1). Furthermore, the list includes many well-characterized genes encoding receptor complexes for GDNF and transcriptional mediators of UBT signaling pathways that are known to be UBT enriched and have been shown to play central roles in UBT regulation, including Ret, Gfra1, Etv4 and Etv5 (Pachnis et al., 1993; Schuchardt et al., 1994; Cacalano et al., 1998; Enomoto et al., 1998; Keefe et al., 2013; Lu et al., 2009). Sox8 co-regulates UBT development together with a UBT-enriched partner Sox9 but Sox9 does not appear in this high stringency list (Table S1). However, Sox9 expression was enriched in the UBT fraction below the 5-fold threshold cutoff [fold change (FC)=4.15, RPKM=23.92]. This result indicates that among the set of genes with <5-fold enrichment there are likely to be additional UBT-enriched genes of biological significance to UBT-directed developmental processes.

The most highly stalk-enriched gene sets includes a number of genes associated with differentiation of mature ureteric epithelial cell types, including Aqp2 (encoding a principal cell-specific water channel protein) (Kwon et al., 2013); Foxq1 (encoding a transcriptional determinant of pH-regulating intercalated cells) (Blomqvist et al., 2004); several uroplakin (Upk) genes, which encode components forming an impermeable epithelial barrier within the urinary tract that extends from the medullary collecting duct to the bladder (Wu et al., 2009); and signal-encoding genes such as Wnt7b and Shh that drive morphogenesis of the deep medullary ureteric epithelium and development of smooth muscle associated with this epithelium, respectively (Yu et al., 2002, 2009; Table S2). In summary, UBT and stalk epithelial populations reflect distinct cell types with expected differences in gene expression and predicted biological action giving confidence that the differential expression screen achieved our primary goal of generating robust datasets for regional gene activity in the developing ureteric system.

To compare expression of all the UBT-enriched genes from the stringent dataset with a subset of the most stalk-enriched gene set (FC≤−50), whole-mount in situ hybridization (WISH) was performed on E15.5 urogenital systems (UGSs) using Wnt11 (UBT specific, high expression) and Shh (stalk specific, low expression) as positive controls for sensitivity, probe penetration and regional expression. We chose E15.5 for follow up in situ analysis because E15.5 samples enable more effective complete penetration in WISH procedures; this stage has been the standard for whole-mount screens in large scale projects within the Genito-Urinary Molecular Anatomy Project (GUDMAP; www.GUDMAP.org) (McMahon et al., 2008; Yu et al., 2012). Furthermore, the expectation is that branching programs are similar at E15.5 and E16.5; to our knowledge no differences in tip gene expression have been described between any active stage of branching tip outgrowth in the mammalian kidney. We also compared our data with data in GUDMAP and Eurexpress (a transcriptome atlas database for mouse embryo; www.eurexpress.org) to evaluate and corroborate our new datasets (McMahon et al., 2008; Yu et al., 2012; Harding et al., 2011; Diez-Roux et al., 2011).

Predicted UBT-enriched genes were stratified into four kidney expression categories: (1) specifically expressed in UBTs, (2) UBT-enriched expression and additional non-ureteric epithelial population(s) expression, (3) expression only detectable outside of the ureteric epithelium, and (4) no expression detectable (Fig. 2). To ensure that none of the genes in categories 3 and 4 scored as negative from a mismatch between the E16.5 screen and the E15.5 in situ analysis, all genes in these categories (32 total) were secondarily screened for UBT expression at E16.5; as expected, none showed detectable UBT expression at E16.5 (data not shown). Furthermore, online (www.gudmap.org and www.eurexpress.org) section in situ hybridization (SISH) and microarray profiling of many of these genes verifies the absence of tip expression. These genes likely reflect true false positives from the primary screen or potentially low levels of expression below the threshold for robust detection.

Eighty-four of the predicted UBT-enriched genes (72%) showed UBT-specific expression patterns: 19 in category 1 expression and 65 in category 2 expression. Genes in these categories were assessed further by SISH on E15.5 UGSs to establish a more detailed view of which specific cell types show expression (Fig. 2). The SISH results confirmed the general WISH annotation but also provide more insights into other expression domains outside of the ureteric epithelium, which include the cap mesenchyme, interstitial cells, and segments of the developing nephron, and give more detailed spatial information on the domain of expression within the UBT.

A close examination of UBT expression identified three distinct categories of genes based on their expression domain within the UB: Class I (20 genes) with expression restricted to the distal tip such as Kcnn4, which encodes a calcium-activated potassium channel protein (Ishii et al., 1997; Joiner et al., 1997); Class II (34 genes) with broader tip expression such as Sema6a, which encodes a transmembrane protein linked to axon guidance (Leighton et al., 2001; Suto et al., 2007; Matsuoka et al., 2011; Belle et al., 2016); and Class III (11 genes), for which gene expression extends beyond the UBT into the near adjacent stalk epithelium, such as Vldlr, which encodes a reelin-binding lipoprotein receptor (Bock and May, 2016; Fig. 3A-F). Class I genes, with the most distal expression, have the highest mean fold change (Fig. 3G). The mean fold change decreases as the expression zones widen in Class II and III. Thus, there is a relationship between the spatial distribution of gene activity within the UBTs and UBT enrichment in the datasets. To confirm the differing zones of UBT expression observed, the RNAscope Duplex Detection kit was utilized to show the expression of two genes on a single tissue section. Wnt11 (Class I) is restricted to the distal tip whereas Ret (Class II) has a broader expression pattern (Fig. 3H). The expression of Vldlr (Class III) extends further than Ret into the neighboring stalk cells (Fig. 3I). Collectively, these data provide strong evidence for molecular heterogeneity among cells within the branch tips.

In the human kidney, UB outgrowth and branching initiates at gestation week 5 and is thought to continue until around week 14-15 when UBTs continue to extend cortically, undergoing growth and differentiation, but in the absence of repeated branching (Osathanondh and Potter, 1963). We used SISH to examine the expression of a selected group of UBT-restricted genes in 15-16 week human fetal kidney samples (Fig. 4A-R). Of these, WNT11, ETV4, ETV5, HS3ST3A1, HS3ST3B1, CRLF1 and SOX8 showed UBT-restricted expression at UBTs in the human kidney (Fig. 4A-P). In contrast, KDM2B, which encodes a lysine demethylase linked to Wnt pathway regulation, was the only gene we observed from this subset that was not detected in the human UBTs at this time, though strong conserved expression was observed within early forming nephron structures (Fig. 4Q,R).

To obtain a broader view of the expression of UBT-specific genes at different stages of the branching process, we examined the expression of category 1 and 2 genes at E12.5 by WISH, shortly after the initiation of UB branching, when a maximum of six branching generations have taken place (Short et al., 2014). Sixty genes (86%) had similar expression at both embryonic stages (Fig. 5): these included genes known to be associated with receptor tyrosine kinase (RTK) (Ret, Gfra1, Etv4, Etv5, Sox8, Spred1, Spred2) and Wnt signaling (Axin2, Ccnd1, Ctnnd2, Sfrp2, Tcf7, Wnt11). For the subset of genes not detected at E12.5, this could reflect dynamic temporal expression such that the gene is off early, or below the limits of detection, at this earlier time point.

To determine which UBT-specific genes in the kidney may have broader roles in growth and branching in other organs, we extended WISH expression studies for this set to E12.5 and E15.5 lungs. Comparative analysis highlighted three distinct categories of expression: (A) genes with branch tip-specific expression in both the kidney and lung, (B) genes with branch tip-specific expression in the kidney and expression in a non-tip population in lung, and (C) genes with branch tip-specific expression in the kidney and no detectable expression in the lung (Fig. 6). Forty genes (48%) of the E15.5 UBT-enriched gene set were expressed also within LBTs, 22 genes (26%) showed expression outside of the LBTs, and 22 genes (26%) were not detectably expressed within lung samples. Category A genes are candidates for playing conserved, generic roles in epithelial branch-tip regulation across organ systems, whereas those genes with kidney-restricted expression are candidates for kidney-specific actions that could link to different branching events driving the distinct pattern of kidney morphogenesis, or distinct tip-niche interactions that distinguish the developing kidney from the lung. Interestingly, a subset of the UBT-enriched genes was detected in other cell type(s) within the lung close to the LBT epithelium (category B). Seventy-six (90%) genes had a similar expression pattern between E12.5 and E15.5 lungs, demonstrating that the developmental process is largely consistent over time. Fig. 7 shows a comparative summary of kidney and lung expression of a selected set of genes at each developmental stage.

The mammalian salivary gland is another branched epithelial organ system that shares some regulatory mechanisms with the lung and kidney (Varner and Nelson, 2014). However, a major distinction is that SBTs arise by a process of clefting, where invaginations of the epithelial basement membrane establish new branching events. This highly dynamic process requires extracellular matrix (ECM) remodeling, including fibronectin assembly and rearrangements in cell-cell and cell-matrix adhesions (Patel and Hoffman, 2014). We investigated the expression of the UBT-enriched genes from categories 1 and 2 within the salivary gland utilizing two online resources: Eurexpress and the Salivary Gland Molecular Anatomy Project (http://sgmap.nidcr.nih.gov) (Diez-Roux et al., 2011; Musselmann et al., 2011) (Table S3). Fifty-five percent (46 genes) of the UBT-enriched genes were also enriched within the developing SBTs, and 31% (26 genes) showed conserved epithelial branch tip-enriched expression in the kidney, lung and salivary gland highlighting potential core branch program regulatory factors (Table S3). Notably, 25 of the 84 (30%) UBT-enriched genes were only detected in UBTs, which could indicate kidney-specific roles for this cohort.

Prior microarray analysis studies examining gene expression in subpopulations of the embryonic kidney offer informative insights into the regulatory networks involved in kidney development (Schwab et al., 2003; Stuart et al., 2003; Schmidt-Ott et al., 2005). Our UBT-enriched list expands these previous studies using RNA-seq technology. This gene set includes well-verified UBT-specific genes, such as Wnt11, Etv4, Etv5, Gfra1, Ret and Sox8/9. Although the RNA-sequencing screen identified established UBT-enriched genes, the majority of the genes identified here have not been documented previously in kidney development and/or branching morphogenesis. The expanded insight into gene expression linked to UBT cells reflects both the improved sensitivity of the RNA-seq-based approach and the absence of inherent limitations using defined probe sets on microarrays.

Whole-mount in situ hybridization on E15.5 mouse kidneys revealed that 84/116 (72%) of the primary UBT-enriched gene set from RNA-seq analysis display tip-restricted expression patterns within the ureteric epithelium. A number of these genes are associated with signaling pathways important in branching morphogenesis, notably RTK, Wnt and retinoic acid pathways (Table S4). Several genes have been implicated in the regulation of the cell cycle and cell proliferation; for example, Ccnd1, Mycn and Cdca7 (Ma et al., 2015; Eilers and Eisenman, 2008; Gill et al., 2013). This observation is consistent with the enhanced rates of cell proliferation seen in UBTs, which is likely to be important for growth and morphogenesis of the ureteric epithelium (Michael and Davies, 2004). The ECM is also of special interest given the requirement of branch tips to move through the extracellular space; several genes in the UBT-enriched gene set are either components of the matrix (Col9a3, Flbn2, Frem2) (Brachvogel et al., 2013; Olijnyk et al., 2014; Timmer et al., 2005) or predicted to remodel matrix components (Adamts18, Cadm1) (Kelwick et al., 2015; Moiseeva et al., 2014) (Table S4).

A subset of UBT-specific genes has been suggested to play a role in kidney development (Table S5). Several of these genes are well-established in kidney development, such as Wnt11, Ret, Gfra1, Etv4 and Etv5. The remaining genes have been linked through individual patient case studies and/or sparse research (Arnt2, Crispld2, Crlf1, Epha4, Mdk, Mycn, Ros1, Sfrp2) (Table S5); their presence in the UBT during kidney development is interesting. Lastly, most genes with UBT-enriched expression have not been implicated in kidney development previously. These novel genes expand our understanding of the molecular features of kidney development identifying new targets for mechanistic follow-up studies.

Distinct zones of gene expression within the UBT are highlighted by distinct subsets of UBT-enriched genes. Class I and II genes demarcate the UBT population discretely and therefore are potentially part of the ‘tip identity’. By contrast, Class III genes extend beyond the UBT into neighboring stalk cells. These differences could reflect distinct regulatory actions in the dynamic events of branching morphogenesis, or in the maintenance and commitment of progenitor cells within the UBT niche. For example, a subset of tip cells lose ‘tip identity’ generating stalk cells in a process that is poorly understood (Shakya et al., 2005; Riccio et al., 2016). Interestingly, Spred2 and Axin2, both Class III genes and negative regulators of MAPK and Wnt signaling, respectively (Wakioka et al., 2001; Yamamoto et al., 1998), show persistent expression in branch tips that could link to an ongoing suppression of UBT signaling promoting branch tip progenitors.

On the other hand, Class I gene expression, exemplified by Wnt11, Tcf7 and Sox8, is tightly restricted to the most distal ends of the UBT. Notably, these genes are all positive activators of branching morphogenesis (Pepicelli et al., 1997; van de Wetering et al., 1997; Reginensi et al., 2011). The broader epithelial expression of negative regulators and tight, tip-most expression domain for positive regulators of branch tip-associated signaling activity could ensure appropriate spatial regulation of these signaling pathways.

Though the temporal analysis was limited, the finding that 86% of genes showed UBT-enriched expression at both E12.5 and E15.5 kidneys argues against a dramatic change in UBT regulatory programs during early stages of mammalian kidney development. Branching growth of the UBT continues until approximately the day of birth (E19.5) in the mouse (Cebrián et al., 2004; Short et al., 2014). However, the UBT niche persists for another 2-3 days of postnatal development, a critical period as over half the nephron precursors (the renal vesicles) form at this time.

Comparing the gene expression pattern of UBT-enriched genes between kidneys and lungs has illuminated potential kidney-specific regulators as well as general regulators of branching morphogenesis. Of the 84 genes that display UBT-specific expression in kidneys, 40 have LBT-specific expression. Furthermore, of these genes with tip expression in both kidney and lung, 67% also show expression in the embryonic mouse salivary gland epithelial tip (according to publicly available online SISH and/or microarray data; www.eurexpress.com; sgmap.nidcr.nih.gov) (Table S3) though the salivary gland forms branches through a distinct process of epithelial cleft formation. This dataset is likely to highlight core branching components. Of note, many of these shared tip-specific genes are known to be associated with ECM properties (Adamts18, Cadm1, Col9a3, Hs3st3a1), and cell cycle and proliferation (Ccnd1, Cdca7, Cdca7l, Mdk, Nasp, Uhrf1) highlighting the importance of specific matrix components to branching and enhanced growth (Michael and Davies, 2004).

A few of these genes have been connected to lung and/or salivary gland branching morphogenesis. For example, Etv5 has been previously implicated downstream of FGF signaling in the regulation of lung branching (Liu et al., 2003), and Vldlr, Hs3st3a1 and Hs3st3b1 have been shown to increase FGFR2b-mediated signaling and proliferation during salivary gland branching morphogenesis (Rebustini et al., 2012; Patel et al., 2014). Furthermore, conditional ablation of Ctnnd2 in the salivary gland resulted in decreased levels of E-cadherin (cadherin 1), and defects in cell-cell adhesion and resulting salivary gland morphology (Davis and Reynolds, 2006).

Interestingly, 25 (30%) UBT-enriched genes appear to be absent from LBTs and SBTs (Table S3). These genes include Gfra1, Sox8 and Wnt11, which have been characterized for their roles in branching morphogenesis in the kidney (Cacalano et al., 1998; Enomoto et al., 1998; Reginensi et al., 2011; Pepicelli et al., 1997). However, no functional analysis has been reported for most of this gene set. As RET/GDNF signaling is a central, kidney-specific pathway stimulating UBT branching and progenitor maintenance, differences in the RTK signaling profile are expected and these are reflected by the kidney-specific expression of Dok6, Gfra1, Sox8, Wnt11, and probably several other genes.

Surprisingly, several genes associated with calcium-related biology appear to be kidney specific. These include Cib2, encoding a calcium- and integrin-binding protein that is associated with photoreceptor and hair cell maintenance (Riazuddin et al., 2012), and Capn6, which belongs to the calpain gene family and encodes a calcium-dependent cysteine protease linked to microtubule stabilization (Tonami et al., 2007).

The majority of developmental kidney studies have used the mouse as the standard model system for human kidney development. However, there are important differences that may reflect different regulatory actions in mouse and human UBTs (Osathanondh and Potter, 1963; Saxen and Sariola, 1987). Whereas the mouse kidney is unilobular, the human kidney is a fusion of multiple independently developing lobes. The differences here may indicate variability in the earliest branching events. Furthermore, much of the cortical expansion of the human kidney beyond week 15 (the human kidney initiates UB outgrowth at around 5 weeks and continues active nephrogenesis until 36 weeks; Osathanondh and Potter, 1963), occurs without further branching.

Despite these marked differences, our initial studies revealed conservation in the expression of a mouse UB-enriched gene set in the human kidney at 15-16 weeks. Thus, human and mouse kidneys likely utilize conserved regulatory mechanisms within and outside of UBTs. Interestingly, Kdm2b, which encodes a lysine demethylase linked to regulation of the Wnt signaling pathway, was not detected in human UBTs (Lu et al., 2015; Boulard et al., 2016). Given our rudimentary knowledge of the actual connection between discrete spatial signaling mediated by WNTs, GDNF, FGFs and retinoic acid, and the distinct cellular outcomes of branching and proliferative cell growth in the mouse kidney, subtle modulation in signaling responses through regulatory factors such as Kdm2b could play a role in the differential dynamics of mouse and human UBTs.

Our study has uncovered a large number of genes of potential interest in kidney-specific or organ-wide programs of epithelial branching growth in the mammalian embryo. For example, Adamts18 is highly enriched in the UBT population (FC=47.47, RPKM=15.25), and also shows elevated expression in LBTs and the epithelial tips in the developing salivary gland. Given that Adamts18 encodes a protease in a family that targets ECM components (Kelwick et al., 2015), it is an attractive candidate for a conserved role across different organs to facilitate matrix remodeling in branching growth.

Hs3st3a1 and Hs3st3b1 encode closely related heparan sulfate 3-O-sulfotransferases; the former is enriched in the tips of the kidney, lung and salivary gland whereas the latter is only enriched in SBTs and UBTs. These genes encode enzymes that catalyze the addition of a sulfate group to heparan sulfates. Given the requirement for heparan sulfate as a co-receptor for FGF receptor activation, and the shared requirement for FGF signaling in branching growth across multiple organs, these genes are predicted to play a conserved role in branching morphogenesis. Indeed, initial evidence points to a role for these enzymes in salivary gland development (Patel et al., 2014).

The next phase of analysis requires a detailed systematic functional dissection of the gene sets identified here in the context of mammalian organogenesis. Although much progress has been made in mutational studies in the mouse, this remains a considerable challenge in studying in vivo development, especially where conditional mutagenesis is required to reveal organ-specific roles. Kidney and lung organoid models could facilitate mutational studies (Takasato et al., 2015; Dye et al., 2015) if it can be shown that these organoid models recapitulate normal epithelial morphogenesis. Connecting these gene sets with genome-wide association studies suggesting gene-disease associations within specific organs systems when taken together with insights from the analyses here may help in prioritizing specific targets of interest for functional analysis.

UBT cells, and their committed epithelial stalk descendants, were isolated from E16.5 kidneys carrying Wnt11-RFP [Harding et al., 2011; B6;D-Tg(Wnt11-TagRFP/cre/ERT2)28Amc/J] and Hoxb7-GFP [Srinivas et al., 1999; 129S.Cg-Tg(Hoxb7-EGFP)33Cos/J] transgenes by FACS. All mouse handling, husbandry and procedures were performed in compliance with the guidelines set forth by the Institutional Animal Care and Use Committees (IACUC) at the University of Southern California. In the ureteric epithelium at this time, Wnt11 is expressed specifically in UBTs, whereas Hoxb7 is present throughout the entire ureteric epithelium [Wnt11 is also expressed within the interstitial mesenchyme cells starting from E15.5 (Yu et al., 2009)]. RNA was isolated from UBT-enriched (RFP⁺, GFP⁺) and stalk-enriched (RFP⁻, GFP⁺) fractions from three biological replicates, and mRNA samples subjected to Next Generation RNA sequencing: 150 bases, paired-end reads were obtained on a HiSeq 2000. Raw sequencing reads were trimmed with the Quality Score method. TopHat2 was utilized to map these trimmed reads to the mouse genome (mm10 annotation). The Partek E/M method was applied using Partek suite bioinformatics tool to quantify the aligned reads to genes/transcripts with Ensembl 72 annotation. Upper quartile normalization was used to normalize the read counts per gene/transcript in all samples (Bullard et al., 2010). A minimum expression filter was used to discard any transcripts if total raw reads across all samples was less than ten. This data was analyzed for differential expression using the PartekGene Specific Analysis method. The resulting RNA-seq data set comprised 37,149 mapped transcripts. We narrowed down these transcripts to focus on those genes that most likely modulate kidney tip morphogenesis by applying a reasonably stringent cutoff with a minimum RPKM value of ten in at least one of the two cell populations and a 5-fold or greater difference in tip or stalk expression to identify robust tip- and stalk-enriched gene sets. The RNA-seq data are available in Gene Expression Omnibus (accession number GSE93267).

Total RNA was extracted from Swiss Webster E15.5 kidneys or human fetal kidney tissues (de-identified human tissue was donated anonymously from elective terminations with informed consent and internal review board approval obtained for the experimental use of the tissue samples) using the RNeasy micro kit (Qiagen). RNA integrity was confirmed using a NanoDrop 2000c (Thermo Scientific) and the SuperScript VILO cDNA synthesis kit (Invitrogen) was used to generate cDNA from the total RNA. Primers for each gene of interest in the RNA expression screen were designed using Primer3, adding a T7 promoter sequence to the 5′ end of the reverse primer. A cDNA probe template was amplified using these primers, the appropriate PCR product confirmed by gel electrophoresis, and digoxigenin-labeled RNA synthesized was from each DNA template using T7 RNA polymerase. Digoxigenin-labeled RNA probes were purified on Micro Bio-Spin columns (Bio-Rad) and the transcript integrity confirmed by gel electrophoresis. Probes were diluted to a concentration of 10 ng/µl in pre-hybridization buffer (50% formamide, 5× SSC, 1% SDS, 50 μg/ml heparin, 50 μg/ml yeast tRNA) and stored long-term at −80°C.

In situ hybridization was performed using our previously reported procedure (Yu et al., 2012). Briefly, for WISH, urogenital systems and lungs of E12.5 and E15.5 Swiss Webster mice were harvested and fixed overnight in 4% paraformaldehyde (PFA), then dehydrated through a methanol series and stored in methanol at −20°C before use. Samples were rehydrated and bleached with 6% hydrogen peroxide, incubated in proteinase K (10 μg/ml), fixed in 4% PFA, and then pre-hybridized in hybridization buffer for several hours at 70°C before hybridization with each RNA probe. Samples were transferred to a BioLane HTI machine for formamide washes, antibody incubation, and MBST [100 mM maleic acid, 150 mM NaCl, 0.1% Tween-20 (pH 7.5)] washes. To reveal in situ hybridization, samples were incubated with BM Purple for up to 48 h, post-fixed in 4% PFA, then transferred to 80% glycerol for imaging on an AxioZoom.V16 stereozoom microscope (Zeiss).

For SISH, mouse and human fetal tissues (anonymously donated from elective terminations and obtained in accordance with institutional guidelines) were harvested and briefly fixed in 4% PFA. Tissues were placed in 30% sucrose overnight at 4°C, embedded in OCT, and sectioned at 12 μm on a Zeiss Microm HM550 cryostat. Slides were fixed in 4% PFA overnight, treated with proteinase K, followed by a brief 4% PFA fixation. To decrease background staining, tissues were acetylated in an acetylation solution (1 M triethanolamine, 0.65% HCl and 0.375% acetic anhydride) and then dehydrated in 95% ethanol. Probes were applied to the slides and incubated overnight at 70°C, then washed with 50% formamide, TNE [10 mM Tris (pH 7.5), 500 mM NaCl, 1 mM EDTA], 2× SSC, 0.2× SSC and MBST [100 mM maleic acid, 150 mM NaCl, 0.1× Tween-20 (pH 7.5)] solutions. Slides were incubated with blocking solution [2% blocking reagent (Roche), 20% heat-inactivated sheep serum] for 1 h, then overnight at 4°C in antibody solution. Samples were stained with BM Purple for up to 14 days, fixed in 4% PFA, and mounted with Glycergel (Dako). SISH samples were imaged on an AxioScan.Z1 (Zeiss).

To compare expression of two genes on the same section, we utilized an RNAscope 2.5 Duplex Detection Kit (Advanced Cell Diagnostics). Following routine dissection, embedding and sectioning as for other SISH experiments, tissues were serially incubated in hydrogen peroxide and protease in the kit according to manufacturer's instructions, then hybridized with RNAscope's probes at 40°C in the HybEZ oven (Advanced Cell Diagnostics). Slides were briefly stained as appropriate for each probe labeling following manufacturer's recommendations, and then mounted with VectaMount (Vector Laboratories). Tissues were imaged on an AxioImager.Z1 (Zeiss).

Two publicly available online resources were utilized to analyze SBT expression. Eurexpress (www.eurexpress.org) provides E14.5 whole embryo SISH data. Salivary Gland Molecular Anatomy Project (http://sgmap.nidcr.nih.gov) showcases microarray data from a range of developmental stages from microdissected regions of the submandibular gland (Diez-Roux et al., 2011; Musselmann et al., 2011). E13.0 epithelial bud to duct expression ratio was calculated, and ratios above 1.15 were considered SBT-enriched.

We thank Yibu Chen and Meng Li for their help with the RNA sequencing, and Gohar Saribekyan for sectioning the human tissue.