Single-cell RNA sequencing reveals PDGFRα⁺ stromal cell subpopulations that promote proacinar cell differentiation in embryonic salivary gland organoids

All developmental processes require precise cellular organization and coordination between many cell types. Branching morphogenesis in the parenchymal epithelium of salivary glands requires instructions from embryonic mesenchymal or adult stromal cells (Sakakura et al., 1976; Wei et al., 2007). These stromal-to-epithelial interactions shape epithelial differentiation, instructing either secretory acinar or ductal cell fates. Functional salivary acini produce saliva and promote good oral health (Pedersen et al., 2018). However, acinar cells accumulate damage from diseases, such as type 2 diabetes and Sjögren's syndrome, from oral drug use, and from radiotherapy for head and neck cancer (Jensen et al., 2014; Liu et al., 2012; Marmary et al., 2016; Plemons et al., 2014; von Bültzingslöwen et al., 2007; Weng et al., 2018; Yoo et al., 2014). Improved therapeutics for diseased and dysfunctional salivary glands require characterization of the molecular mechanisms involved in stromal-epithelial cell interactions that drive secretory acinar cell differentiation.

Tissue recombination experiments first demonstrated that mesenchymal cells in salivary glands impact epithelial differentiation (Hoffman, 2002; Kusakabe et al., 1985; Sakakura et al., 1976; Wei et al., 2007). Fibroblast growth factor (FGF) signaling mediates important stromal epithelial interactions, which facilitate developmental processes, including limb formation and wound healing (De Moerlooze et al., 2000; Martin, 1998). FGF signaling between the stroma and epithelium also affects how submandibular salivary glands (SMGs) develop (Chatzeli et al., 2017; Hoffman, 2002; Lombaert et al., 2013; Makarenkova et al., 2009; Steinberg et al., 2005; Wei et al., 2007). The FGF growth factor family is a large and complex family, with 23 known interacting members and four receptors with multiple isoforms (Yun et al., 2010). Mice with FGFR2IIIb knocked out exhibit abnormal salivary gland development (De Moerlooze et al., 2000; Ohuchi et al., 2000) because the epithelium ignores stromal signals. Similarly, when the FGFR2IIIb receptor ligand, FGF10, is knocked out, early gland development fails (Tina Jaskoll et al., 2005; Ohuchi et al., 2000; Steinberg et al., 2005). The salivary gland connective tissue includes diverse cell types, such as stromal fibroblasts, vascular endothelial cells and nerves, which synergize for instructing development (Kwon et al., 2017; Nedvetsky et al., 2014; Takebe et al., 2013, 2015). In previous work, we formed complex salivary organoids from embryonic day (E) 16 epithelial and stromal cells, which underwent robust FGF2-dependent proacinar cell differentiation (Hosseini et al., 2018, 2019). FGFR pharmacological inhibition or FGF2 knockdown specifically within the stroma inhibited epithelial proacinar cell differentiation. However, stromal factors that directly promote epithelial proacinar cell differentiation were not identified.

Using hierarchical clustering, single-cell RNA-sequencing (scRNA-Seq) has revealed high levels of cellular heterogeneity within organs (Enge et al., 2017; Hauser et al., 2020; Kwon et al., 2019; Lu et al., 2016; Macosko et al., 2015; McCarthy et al., 2020b; Sekiguchi et al., 2020; Song et al., 2018; Zepp et al., 2017, 2021). These large datasets and their computational models have identified diverse cell subsets within known cell types. This extra complexity in known cellular communications increases developmental nuances and reveals previously uncharacterized cellular specificities. In this study, we identified numerous stromal cell subpopulations in both E16 SMG and complex salivary organoids using scRNA-Seq. Moreover, using microarray and scRNA-Seq, we determined that FGF2 directly regulates the stromal cell phenotype. We identified effects of FGF2-dependent stromal BMP genes on epithelial proacinar cell differentiation. Pdgfra is a FGF2-controlled gene that is expressed in several stromal subpopulations in the salivary gland. PDGFRα⁺ stromal cells alone sufficiently promoted epithelial proacinar cell differentiation in salivary organoids. FGF2-dependent BMP genes synergized with FGF2, promoting proacinar cell differentiation. FGF2 sustained stromal competence that instructed epithelial proacinar cell differentiation, prevented stromal myofibroblast transition in culture and reversed TGFβ1-mediated proacinar cell dedifferentiation in organoid cultures.

Given that our previous research showed that proacinar cell differentiation in our complex salivary organoids required E16 SMG stromal cells (Hosseini et al., 2018), we hypothesized that specific stromal subpopulations may direct epithelial proacinar cell differentiation. We addressed this hypothesis by performing scRNA-Seq using enriched stromal cells isolated from E16 SMGs. The compiled scRNA-Seq data contained 9165 fully sequenced cells (Fig. S1). The transcriptional subpopulations were modeled using Seurat's principal component analysis (PCA) followed by k-nearest neighbor clustering and Louvain algorithms (Satija et al., 2015). We visualized the modeled subpopulations using graph-based dimensional Uniform Manifold Approximation and Projection (UMAP) reductions. Based on cells expressing the gene encoding epithelial cell adhesion molecule (Epcam), we sequenced 1195 epithelial cells and 7970 nonepithelial cells (Fig. 1A). Epithelial cell segregation and reanalysis (Fig. 1C) yielded nine epithelial subpopulations, including known epithelial subtypes, such as Krt5⁺ and Krt19⁺ ductal cells and Epcam⁺Acta2⁺ myoepithelial cells (Hauser et al., 2020; Song et al., 2018). Additionally, we found that the two distinct proacinar populations expressing either Smgc or Bpifa2, previously reported in postnatal day (P) 1 SMG (Hauser et al., 2020), were detectable at E16. We also identified many connective tissue cell types (Fig. 1A,B; Fig. S1). Based on Col1a1 expression, we sequenced 3790 stromal fibroblast cells. The remaining 4180 cells included endothelial cells, lymphatic endothelial cells, immune cells, neurons, Schwann cells and erythrocytes, as defined by marker gene expression.

The stromal fibroblasts were segregated and reclustered, identifying 13 previously undefined stromal cell subpopulations (Fig. 1D). Nine subpopulations are classifiable based on the known stromal identity markers Pdgfra, Pdgfrb, Acta2, Thy1 and Wnt2 (Contreras et al., 2019; Dominici et al., 2006; Karpus et al., 2019; Kishi et al., 2018; Li et al., 2018; Santini et al., 2020; Zepp et al., 2017, 2021), and proliferative subpopulations accounted for the other four stromal clusters. Both PDGFRα and PDGFRβ are especially interesting because cells expressing these receptors have been shown to play important roles in organ epithelial differentiation, wound healing and disease phenotypes (Kinchen et al., 2018; Kramann et al., 2015; Kuppe et al., 2021; Li et al., 2018; Liu et al., 2019; Ramachandran et al., 2019; Santini et al., 2020). PDGFRα is expressed in the neural crest-derived SMG mesenchyme, and PDGF signaling regulates early branching morphogenesis through FGFs (Hoffman, 2002; Morikawa et al., 2009; Sakakura et al., 1976; Soriano, 1997; Steinberg et al., 2005; Yamamoto et al., 2008). We annotated the transcriptomes of the Pdgfra⁺ subpopulations using Metascape (Zhou et al., 2019) and found transcriptional landscapes matching the traditionally understood functions of the stroma. These Pdgfra⁺ clusters showed differing Gene Ontology (GO) categories for extracellular matrix organization, vasculature development and supramolecular fiber organization (Fig. 1E). We evaluated subpopulations expressing two growth factors with known importance in proacinar cell differentiation: Fgf2 (Hosseini et al., 2018) and Fgf10 (Chatzeli et al., 2021; Patel et al., 2007; Wells et al., 2013) (Fig. 1F). Of the nine Pdgfra⁺ subclusters, seven expressed Fgf2 and Fgf10 at varying levels, suggesting that these stromal subpopulations may participate in proacinar cell differentiation instruction.

To explore the possible developmental relationship between heterogeneous stromal cell types, we modeled cell differentiation trajectories in pseudotime (Fig. 1G). Using Monocle DDRTree (Qiu et al., 2017), we inferred trajectories from the E16 scRNA-Seq transcriptomes. We set the root node as the noncell-cycling Pdgfra⁺ subpopulation because previous reports indicated that the stroma arises from the Pdgfra⁺ migratory cranial neural crest (Morikawa et al., 2009). The resulting trajectory analysis showed three branches: one branch directly linked the Pdgfra⁺ root subpopulation with the Wnt2⁺Thy1⁺ subpopulation; a second branch led to another Pdgfra⁺ subpopulation that expresses more ECM genes compared with the root population; and a third branch transitioned from a Pdgfra⁺ state through the cell cycle into a Pdgfrb⁺Acta2⁺ subpopulation. This analysis predicts that the Pdgfra⁺ progenitor population has multiple differentiation trajectories during SMG development.

Given that organoids can be used to evaluate stromal-epithelial cell signaling, we examined the changes in salivary gland organoid epithelial gene expression in response to FGF2. We produced organoids containing E16 SMG epithelial and stromal cells grown with or without FGF2, as performed previously (Hosseini et al., 2018, 2019). We separated the organoid epithelial cells from the stromal cells and identified epithelial differentiation transcriptomes using microarrays. We detected a clear separation between the control- and FGF2-treated epithelial cells in the first principal component (PC1) of the PCA (Fig. S2C). We generated a twofold change heatmap and noted specific acinar and ductal differentiation markers showing contrasting signal magnitudes between the FGF2 and control conditions (Fig. 2A). FGF2 increased proacinar gene expression, as expected based on our prior work (Hosseini et al., 2018, 2019), including known proacinar genes, Smgc [encoding mucin 19 (MUC19)] and Bpifa2 [encoding parotid secretory protein (PSP)], and known acinar genes that are not highly expressed at the protein level at E16, including Prol1 [encoding mucin 10 (MUC10)] and Bhlha15 [encoding muscle, intestine and stomach expression 1 (MIST1)]. Gene set enrichment analysis (GSEA) (Subramanian et al., 2005), using gene sets derived from P8 SMG epithelial scRNA-Seq data (Song et al., 2018), identified global differences showing acinar and ductal genes enriched under FGF2-cultured and control conditions, respectively (Fig. 2B). However, neither organoid condition showed myoepithelial gene enrichment. Using immunocytochemistry (ICC), we confirmed that, in organoids lacking FGF2, expression of the ductal marker KRT7 progressively increased and was inversely correlated with the acinar marker, AQP5 (Fig. S2A,B). Together, these results identify an FGF2-dependent proacinar transcriptome, and an emerging ductal transcriptome if FGF2 is absent.

Given that in vivo SMGs contain both proacinar cells and adjacent myoepithelial cells, forming bilayered acini, we looked for protein-level evidence that FGF2 regulates myoepithelial cell differentiation. Organoids were grown for 7 days with or without stroma and with or without FGF2, and ICC was performed with the myoepithelial marker, calponin 1 (CNN1). When stroma was present, CNN1 was detected in the organoids, whereas stroma-deficient organoids had little-to-no detectable CNN1. Stromal cell-containing organoids showed similar CNN1 levels, independent of FGF2 (Fig. 2C,D), consistent with transcriptional analysis. To confirm that CNN1 was expressed by myoepithelial cells, we grew epithelial-stromal organoids with and without FGF2 and performed ICC using CNN1 together with the epithelial marker, EPCAM, and the basement-membrane protein heparan sulfate proteoglycan 2 (HSPG2 or perlecan). With high-magnification confocal imaging, we confirmed that the EPCAM⁺ cells expressing CNN1 were enclosed within the HSPG2⁺ basement membrane (Fig. 2E), as expected in myoepithelial cells. These data demonstrate that, without FGF2 signaling, stromal cells can sufficiently stimulate myoepithelial cell differentiation.

Given that we had previously determined that stromal cells are the direct targets of FGF2 signaling in salivary gland organoids (Hosseini et al., 2018), we evaluated the effects of FGF2 on the stromal cell phenotype. Isolated E16 stromal cells were grown on glass coverslips or on porous polycarbonate filters in the presence or absence of FGF2. Stromal cells grown on glass with FGF2 showed a different cellular morphology compared with cells grown without FGF2. Although the PDGFRβ area was constant, there was a threefold increase in area in cultured cells expressing PDGFRα in the presence of FGF2 (Fig. 3A,B). Additionally, THY1 expression was generally elevated in stromal cells grown with FGF2 on glass (Fig. S3A,B). Thus, the organization of the stroma was dependent on the culture surface but not on FGF2. Although the stromal cells grew as a 2D monolayer on glass, they coalesced into 3D spherical aggregates when cultured on the filters, reminiscent of mesenchymal stem cells in a classical colony-forming assay (Friedenstein et al., 1970). Additionally, stromal cells grown on softer polycarbonate filters with FGF2 showed a threefold increase in relative cell area that was positive for THY1 (Fig. 3C,D). In organoids treated with or without FGF2, FGF2-dependent PDGFRα, and FGF2-independent PDGFRβ expression in stroma was consistent with that in the 2D cultures (Fig. 3E,F). Thus, stromal cell PDGFRα and THY1 protein levels are dependent on FGF2 signaling. To examine these cell populations in vivo, we performed immunohistochemistry on E16 and adult SMG tissue sections and identified key stromal markers: PDGFRα, PDGFRβ, vimentin (VIM), ENG and THY1 (Fig. 3G). We detected stromal cell subpopulations expressing markers generally comparable with the scRNA-Seq profile of the E16 gland (Fig. 1C). scRNA-Seq showed a small Thy1⁺ stromal subpopulation that was also Wnt2⁺ and Pdgfra⁺ (Fig. 1D). Interestingly, in tissue sections, these THY1⁺ stromal cells resided in a neurovasculature region near endothelial cells expressing platelet endothelial cell adhesion molecule (PECAM1/CD31) and neuronal cells expressing tubulin-β3 (TUBB3) in the adult SMG (Fig. 3G).

Although we previously reported that total E16 stromal cells can promote proacinar cell differentiation in adult epithelial cells in organoids (Koslow et al., 2019), we did not examine whether adult stromal cells can instruct proacinar cell differentiation in our complex salivary organoids. We tested this by isolating and expanding adult stromal cells in FGF2-containing media. Once enough stromal cells formed, they were combined into organoids using E16 SMG epithelial cells. Significantly, expanded adult SMG stromal cells retained PDGFRα expression and, based on AQP5 expression, promoted epithelial proacinar cell differentiation (Fig. 3H,I). These data indicate that adult stromal cells can instruct FGF2-dependent proacinar cell differentiation.

We next questioned whether PDGFRα⁺ stromal fibroblasts sufficiently promote proacinar cell differentiation in organoids. We isolated E16 PDGFRα⁺ stromal cells using magnetic-activated cell sorting (MACS). We compared organoids generated using MACS PDGFRα⁺-selected cells with organoids containing epithelial cells alone or epithelial cells with whole E16 stroma. Organoids containing only epithelial cells and treated with FGF2 demonstrated little to no PDGFRα⁺ stromal cell contamination (Fig. 4A,B), very low expression of AQP5 and no budding (Fig. 4C,D). By contrast, organoids formed using PDGFRα⁺ MACS-sorted stroma or total stroma showed a similar PDGFRα⁺ area (Fig. 4A,B). Significantly, these FGF2-treated organoids containing MACS-isolated PDGFRα⁺ stromal cells showed high AQP5 expression and branched morphology comparable with whole stroma-containing organoids (Fig. 4C,D). These data demonstrate that PDGFRα⁺ stromal cells are sufficient to stimulate the budding proacinar organoid phenotype.

To identify how FGF2 stimulation changes the transcriptome of SMG stromal cells and instructs epithelial proacinar cell differentiation, we performed scRNA-Seq analysis on enriched stromal cells derived from organoids in the presence or the absence of FGF2. After Seurat analysis, we detected multiple cell types, including the expected stromal enrichment in both datasets (Fig. 5A,B; Fig. S4). When both datasets were integrated, a significant shift was seen in the presence (yesFGF2) or absence (noFGF2) of FGF2 (Fig. 5C,D). The stroma of the FGF2 organoid contained a high percentage of cells expressing both Pdgfra and Thy1 (Fig. 5E,F). Both genes were predominantly in clusters 0, 1, 6 and 8, and showed GO categories associated with cell movement and vasculature development (Fig. 5G). The stromal subpopulations in control organoids lacking FGF2 contained higher numbers of cells expressing Acta2 and/or Pdgfrb (Fig. 5E,F). These clusters were labeled as 2, 3, 7 and 10, and were enriched for GO categories of actin cytoskeletal arrangement, cell adhesion, vascular development, response to wounding and ECM organization (Fig. 5G).

Using microarrays, we also examined the transcriptomic changes in organoid stromal cells depending on the presence or absence of FGF2. PC1 of the PCA accounted for over 50% of the variance and demonstrated significant differences between the two conditions (Fig. 6A). We examined comprehensive transcriptomic changes that were FGF2 dependent using GSEA (Subramanian et al., 2005) and observed that FGF2 treatment enriched many stromal gene sets associated with proliferation or targets of the proliferation-inducing E2F transcription factors (Fig. 6B). Gene enrichment of E2F targets was expected because FGF2 participates in driving the cell cycle (Krejci et al., 2004; Neary et al., 2005).

To identify FGF2-regulated paracrine factors that may instruct epithelial proacinar cell differentiation, we targeted the growth factors that were differentially expressed. From this list, we selected three growth factor genes with a greater than twofold difference in expression with and without FGF2: amphiregulin (Areg), epiregulin (Ereg), and bone morphogenetic protein 7 (Bmp7) (Fig. 6C). These three growth factors were also upregulated in scRNA-Seq-described subpopulations from the FGF2 organoid stroma (Fig. 6D). Prior research illustrated the mixed contributions of BMPs to promoting salivary gland branching morphogenesis (Hoffman, 2002; Steinberg et al., 2005), but epithelial cell differentiation was not addressed. However, in intestinal organoids, hair follicles and lungs, BMPs induce epithelial cell differentiation (Lu et al., 2016; McCarthy et al., 2020a; Zepp et al., 2017). Areg and Ereg are EGF homologs, and EGF signaling is a known SMG developmental promoter (Mizukoshi et al., 2016). We tested whether any of these three factors could substitute for the stroma. We created epithelial organoids lacking stroma but cultured with media containing FGF2 and one candidate growth factor. BMP7, but neither EGF family ligand (AREG or EREG) induced AQP5 in conjunction with FGF2 (Fig. 6E,F), although budding of the organoid decreased. We validated that BMP signaling in the epithelium is necessary for proacinar marker expression using a BMP inhibitor. Epithelial-only organoids were treated with FGF2, BMP7 and a selective BMP receptor inhibitor, LDN19318, which targets primarily ALK2 and ALK3. Combined treatment with FGF2, BMP7 and LDN19318 resulted in no proacinar marker expression in the epithelial organoids, which looked similar to FGF2-only cultured organoids (Fig. 6G,H). Given that BMP2 is another BMP family member that can signal through ALK2 and ALK3 (Sanchez-Duffhues et al., 2020) and was FGF2-dependently upregulated in both the microarray and scRNA-Seq (Fig. 6C,D; Fig. S4), we also tested BMP2 in FGF2-treated epithelial organoids. The response to BMP2 was similar to that to BMP7, synergizing with FGF2 and promoting AQP5 expression (Fig. 6I,J). Thus, these data identify BMP signaling as a FGF2-dependent, stromal-produced, paracrine-acting pathway that promotes proacinar cell differentiation in salivary organoids through FGF2 synergy.

Based on the increased Acta2 and Cnn1 expression in the stroma and reduced PDGFRα staining in organoids when FGF2 was absent, we questioned whether lacking FGF2 phenotypically changes stromal cells. We examined the microarray gene sets enriched in the stroma grown without FGF2 (Fig. 7A). Three of the most significantly enriched gene sets, identified by GSEA, were myogenesis, extracellular matrix and TGFβ signaling. All three gene sets align with myofibroblast characteristics. We hypothesized that, without FGF2, stromal cells differentiate into myofibroblasts. Myofibroblasts are known injury and fibrotic disease contributors that increase extracellular matrix deposition (Contreras et al., 2019; Kramann et al., 2015; Kuppe et al., 2021; Li et al., 2018; Nagaraju et al., 2019). We generated a heatmap for core myofibroblast genes that showed that all genes were more highly expressed in control organoid stroma compared with FGF2-treated organoid stroma (Fig. 7B). Based on the in vivo scRNA-Seq trajectory prediction (Fig. 1F, Fig. 7C) and organoid scRNA-Seq data (Fig. 5, Fig. 7D; Fig. S4), E16 stromal cells have the differentiation potential to express the classical myofibroblast gene smooth muscle alpha actin (αSMA) and TGFβ ligands. We used GSEA to determine which scRNA-Seq-described SMG subpopulation best represented the transcriptome of the total control organoid stroma (Fig. 7E), finding that the Pdgfrb⁺Acta2⁺, stromal cluster 10 (SC_10), genes were most representative. These data suggest that the SC_10 transcriptome dominates control organoid stroma and has myofibroblast-like characteristics.

To confirm whether SMG stromal cells differentiate into a myofibroblast-like state in the absence of FGF2, we cultured E16 stromal cells alone on coverslips with or without FGF2 and performed ICC using classical myofibroblast markers. Stromal cells grown without FGF2 showed a nearly fivefold elevation in the myofibroblast markers, αSMA and CNN1, relative to cells grown with FGF2 (Fig. 7F,G). To determine whether similar changes in gene expression occurred in the organoids, we examined αSMA and CNN1 expression in organoids grown with or without FGF2. We detected a similar myofibroblast-like transition in organoids in the absence of FGF2, but only with a threefold change (Fig. 7H,I). A time-course was performed with time points at 2, 4 and 7 days to examine explicitly myofibroblast conversion in organoids in the absence of FGF2. At day 2, there was no significant difference in CNN1 staining in the organoid stroma relative to FGF2 treatment (Fig. 7J,K). By day 4, there was noticeably more myofibroblast marker expression in the control organoid stroma, and this myofibroblast marker expression increased into day 7. These data support the hypothesis that SMG stromal cells progressively differentiate into a myofibroblast-like state in the absence of FGF2.

Two mRNAs that were increased in organoid stroma grown without FGF2 (no_FGF2) were Tgfb2 and Tgfb3 (Fig. 7B), which are TGFβ family profibrotic cytokines that can induce myofibroblast differentiation and fibrosis (Avery et al., 2018; Nagaraju et al., 2019; Ó Hainmhire et al., 2019; Woods et al., 2015). Notably, the Acta2⁺Pdgfrb⁺ stromal cluster (SC_10) showed appreciable Tgfb2 and Tgfb3 expression in the scRNA-Seq-defined stromal cell subsets (Fig. 7C). Both Tgfβs were also upregulated in specific organoid stromal subpopulations in control conditions (Fig. 7D). We tested TGFβ signaling in the organoids by treating them either for 4 days with TGFβ and subsequently with FGF2 for 3 days (Fig. S5A,B) or with FGF2 for 4 days and subsequently with TGFβ for 3 days (Fig. S5C,D). These experiments showed a correlation between proacinar inhibition and the myofibroblast phenotype, revealing an inherent plasticity in both the epithelial and stromal cell populations (Fig. S5). We also found that prolonged culturing, for up to seven passages, increased CNN1 expression in the stroma, indicative of a myofibroblast phenotype. These late-passage stromal cells did not promote proacinar cell differentiation when incorporated into FGF2-treated organoids (Fig. S6A,B). We further tested whether TGFβ signaling inhibits proacinar cell differentiation. First, organoid AQP5 expression was induced after culturing with FGF2 for 5 days. Then, FGF2 was washed out and replaced with TGFβ1-containing media. One day after TGFβ1 addition, the area positive for the proacinar marker, AQP5, had decreased threefold (Fig. 7L,M). Notably, this TGFβ1-induced AQP5 reduction was reversible after TGFβ1 media was washed away and FGF2 media was reapplied. AQP5 re-expression occurred 4 days after FGF2 re-introduction (Fig. 7L). These FGF2-TGFβ1-FGF2 sequentially cultured organoids had AQP5 staining similar to that in organoids grown with FGF2 for 10 consecutive days (Fig. 7M). Together, these data indicate that TGFβ treatment inhibits proacinar organoid differentiation in a reversible fashion and that FGF2 can reverse TGFβ1-mediated proacinar inhibition. This epithelial reversal occurs after the myofibroblast-like phenotype in the stroma reverses.

Even though it has been known for decades that stromal fibroblasts direct organ development, recent advances in scRNA-Seq are illuminating the specific contributions of this diverse population. Here, we performed scRNA-Seq on developing salivary glands and organoids and evaluated the contributions of stromal cell subpopulations to acinar cell differentiation. We defined subpopulations using specific markers describing functional stromal subsets in other organs, specifically Pdgfra, Pdgfrb, Acta2, Eng, Thy1 and Wnt2 (Karpus et al., 2019; Kramann et al., 2015; Kwon et al., 2019; Lee et al., 2017; Li et al., 2018; Santini et al., 2020; Stzepourginski et al., 2017; Yang et al., 2017; Zepp et al., 2021; Zhao et al., 2017). Our pseudotime analysis interpretation suggests that these burgeoning subpopulations arise from a Pdgfra⁺ population that differentiates into multiple subpopulations. Organoid stromal scRNA-Seq provides evidence that these divergent differentiating stromal cells are FGF2 dependent. UMAP regression and superimposed pseudotime models revealed that these stromal cells progress through a continuous transcriptional landscape rather than by discrete trajectories. By contrast, other developmental scRNA-Seq data using pseudotime analysis suggest that embryonic epithelial subpopulations in the SMG and lung become more transcriptomically isolated over time (Hauser et al., 2020; Zepp et al., 2021). Whether stromal cells differentiate into distinct subpopulations, as epithelial populations do, is not yet firmly established. However, our data show a dualistic response to the growth factors FGF2 and TGFβ1 and suggest that stromal cells remain more plastic during development. Recent studies reveal that, in other organs, specific stromal subpopulations help create microniches that localize distinct epithelial identities into specific regions and mediate disease responses (Contreras et al., 2019; Li et al., 2018; McCarthy et al., 2020b; Santini et al., 2020; Yang et al., 2017). Similarly, we identified stromal cell subpopulations expressing Pdgfra that support proacinar cell differentiation in embryonic organoids (Fig. 8).

Our developmental model suggests that E16 stromal cells have multiple transcriptomic trajectories. The model predicts that these PDGFRα⁺ cells can transition into two smaller subpopulations that increase either Thy1 or Acta2/Pdgfrb. Our in vitro studies, in which stromal cells were cultured either with or without FGF2, showed that PDGFRα⁺ cells can differentiate into at least two different subpopulations, which expressed markers analogous with the pseudotime subpopulations modeled as end nodes. This finding is consistent with that of others who have shown that PDGFRα⁺ stromal cells have a transcriptomic plasticity that is important for regulating both developmental and injury responses (Contreras et al., 2019; Li et al., 2018). Using organoids, we explored the control of a specific stromal population over the epithelium, finding that FGF2 maintains a PDGFRα⁺ stromal state resembling in vivo stromal cell gene expression. Cells in this stromal state can promote epithelial proacinar cell differentiation. We provide evidence that stromal cells expressing Pdgfrb and Acta2 resemble classically defined myofibroblasts that do not support the proacinar cell differentiation state. The current definition of a myofibroblast is controversial and myofibroblasts have different characteristics in vivo and in vitro (Bochaton-Piallat et al., 2016). Whether Pdgfra⁺ cells become myofibroblasts is not clear and, in addition, how PDGFRα and PDGFRβ signaling contributes to organoid differentiation also remains to be determined. Nevertheless, we have identified a clear relationship between the transcriptomic state of stromal cells and the differentiation state of salivary gland epithelial cells.

Others have suggested that PDGFRα stromal expression aligns with functional dexterity (Contreras et al., 2019; Li et al., 2018; Santini et al., 2020). Our organoid experiments showed that these PDGFRα⁺ cells can recapitulate the function of the whole stroma in promoting proacinar cell differentiation. Organoids containing embryonic or adult stroma both retained similar proacinar-promoting abilities, which suggests a functional preservation in some stromal cell populations. The continuing instructiveness but decreasing cellular proportions of the stroma during development may reflect a conversion of the epithelium from stromal paracrine signaling to homeostatic autocrine signaling. As epithelial cells reach maturity, they may exhibit reduced stromal dependence in homeostasis and the stroma may no longer be required at high density unless damage occurs. Organoids can model damage and repair mechanisms and, thus, under these conditions, we showed that adult stromal cells can direct epithelial differentiation. Interestingly, the adult stroma retains FGF2 sensitivity and the ability to promote proacinar cell differentiation when cultured with FGF2. Whether the mechanisms through which adult and E16 stromal cells induce proacinar cell differentiation are similar or whether adult cells that support proacinar cell differentiation are PDGFRα⁺ remain to be determined.

Our study additionally suggests that the stromal state affects the differentiation of epithelial progenitor cells into ducts or acini (Fig. 8). First, FGF2 stimulates stromal PDGFRα and BMP expression. BMP signaling directs the epithelium toward a proacinar rather than a ductal state. However, BMPs are insufficient on their own to directly induce proacinar cell differentiation; FGF2 signaling is also required but is also insufficient. FGF2 signaling may somehow ‘prime’ the epithelium to respond to the BMP signal. Results of our organoid experiments, using only epithelial cells and growth factors, also indicated that proacinar cell differentiation and branching morphogenesis are distinct processes. This idea is in line with studies showing that BMP7 does not complement branching morphogenesis (Hoffman, 2002; Steinberg et al., 2005), and studies demonstrating that cell-to-matrix adhesions direct branching morphogenesis (Sakai et al., 2003; Wang et al., 2021). Other FGF-dependent and FGF-independent signaling pathways contribute to morphogenesis. Stromal BMPs directing epithelial differentiation is a theme in many organs. A similar BMP gradient induced by PDGFRα⁺ cells also controls differentiation in the small intestinal crypt (McCarthy et al., 2020a). These intestinal PDGFRα^hi telocytes and our PDGFRα⁺ salivary stroma may be functionally and transcriptomically similar. This BMP shift contrasts with a second effect that FGF2 has in salivary organoids, whereby FGF2 signaling also prevents PDGFRα⁺ cells from differentiating into ACTA2⁺/CNN1⁺ myofibroblasts. This myofibroblast phenotype produces at least one known branching morphogenesis-inhibiting growth factor, TGFβ (Hall et al., 2010; Woods et al., 2015). TGFβ inhibits proacinar cell differentiation and morphogenesis, and likely promotes ductal differentiation, as previously reported (Iwano et al., 2002; Liu et al., 2016; Rastaldi et al., 2002; Rees et al., 2006; Termén et al., 2013). This BMP and TGFβ signaling axis may control the balance between developmental maturation and injury recovery. BMPs are known TGFβ signaling antagonizers, largely due to competition for shared intracellular signal transduction pathway mediators (Ning et al., 2019). TGFβ is important for early salivary gland development (Jaskoll and Melnick, 1999) and is co-opted during injury for repurposing cells. The BMPs that originally differentiate cells during development (Bandyopadhyay et al., 2006; Lu et al., 2016; McCarthy et al., 2020a,b) can be used to redifferentiate adult cells post-injury. However, further time-dependent stromal analyses during injury will be needed to confirm this stromal switch.

Evaluating stromal heterogeneity remains a difficult frontier. Organs with smaller stromal cell proportions, such as salivary glands (Hauser et al., 2020), require cell enrichment steps before informative scRNA-Seq and hierarchical clustering can be performed. Even if the adult stromal population is small, our findings suggest that the adult stroma in salivary glands retains a therapeutic potential. However, more studies evaluating the evolution of stromal cell subpopulations at later developmental stages and additionally in disease models will be needed to effectively manipulate stromal cells for regenerative medicine approaches. Our research demonstrates the advantage of the organoid model for evaluating predictions made with scRNA-Seq. Organoids facilitate modular testing of cellular and mechanistic hypotheses that enable big data approaches, such as scRNA-Seq. Overall, understanding the stromal complexity is necessary for understanding cellular interactions that could inspire useful regenerative therapies.

E16 timed-pregnant CD-1 female Mus musculus were ordered from Charles River Laboratories. First, the E16 embryos and then the SMGs from those embryos were dissected out following animal protocols approved by the University at Albany Institutional Animal Care and Use Committee (IACUC). SMG removal involved slicing the mandible with a sharp scalpel and then removing the glands using sterile forceps under a dissecting microscope (Nikon SNZ800). SMGs were microdissected in 2×collagenase/hyaluronidase (StemCell Technologies, 7912) diluted using 1×PBS (Life Technologies, 70011-044). The solution was further diluted, and SMGs were further dissected into lobules, when an additional volume of 1.6 U/ml of dispase II (Life Technologies, 17105041) was added, creating a final concentration of 0.8 U/ml dispase II. The lobules were further broken down after a 30-min incubation at 37°C followed by trituration, generating epithelial clusters and single stromal cells. Enzymatic activity was quenched with 10% fetal bovine serum (FBS) (Life Technologies, 10082147) in 1:1 DMEM/F12 (Life Technologies, 11039047) supplemented with 100 U/ml penicillin and 100 mg/ml streptomycin (Pen-Strep; Life Technologies, 15140122). Cell populations were separated by gravity sedimentation for ∼10 min until the large clusters formed a loose pellet. Any DNA-induced clumping was stopped using 0.05 mg/ml DNase I (STEMCELL Technologies, 07900) during the second gravity sedimentation. The stromal cells were enriched by filtration of the resulting supernatant through 70 µm (Falcon, 087712) and 40 µm cell strainers (Thermo Fisher Scientific, 22363547). Filtered stroma were pelleted at 300× g for 8 min and the buffer was replaced using DMEM/F12 containing 10% FBS and Pen-Strep.

Epithelial clusters were enriched, as previously described (Hosseini et al., 2018, 2019), from the gravity pellet with two additional gravity sedimentations and by plating the clusters into a 35 mm tissue culture dish, incubating at 37°C for 2 h and then collecting the floating cells. Any suspended stromal cells were removed from epithelial clusters using two washes with centrifugation at 10× g for 1 min, supernatant removal and resuspension in medium.

Stromal cells were enriched, as described above, from SMGs harvested from 14 E16 embryos. The embryos were excised from two CD-1 timed-pregnant female mice. Dead-cell depletion was performed using the Miltenyi dead-cell removal kit (Miltenyi Biotec, 130-090-101). The live-to-dead cell ratios and cell numbers were evaluated with Trypan Blue (Gibco, 15250-061) staining and a hemocytometer (Marienfeld Superior, Neubauer-improved). The preparation had >80% viable cells and was frozen down in cryovials using 90% FBS+10% dimethyl sulfoxide (Amresco, 200-664-3). Cells were frozen at −1°C/min in a 100% isopropanol freezing chamber overnight at −80°C. Frozen cells were stored in a liquid-nitrogen chamber prior to shipping overnight on dry ice. Cellular recovery goal for library preparation was 10,000 cells. scRNA library preps were created and sequenced by Singulomics or created in-house using the Chromium Next GEM Single Cell 3′ GEM, Library & Gel Bead Kit v3.1 (10X Genomics, 1000121) and sequenced at the University at Albany Center for Functional Genomics.

Bcl2fastq and CellRanger v4.0 were used to assemble and count the sequence data. Data files were imported using SEURAT v3.1.5 in R v3.6.3 (R Core Team, 2021) and R Studio v1.2.5042. Data clusters were calculated following the default pipeline (Satija et al., 2015; Stuart et al., 2019). Dead or apoptotic cells were removed if >7% of unique molecular identifiers (UMIs) mapped to mitochondrial genes. Any cells with <200 or >9000 genes were also removed before downstream analysis to remove noncells and doublets, respectively. The ‘ElbowPlot’ function and PC ‘Heatmaps’ helped determine how many PCs were used for unsupervised cluster modeling. Clustering calculations used linear dimensional reduction, and visualization was performed using nonlinear dimensional reductions and UMAP. In total, 40 PCs were used for E16 stroma calculations, 21 PCs and 40 PCs were used for the epithelial subset and the stroma subset calculations, respectively and 20 PCs were used for the integrated organoid datasets.

The genes defining the unsupervised clustering were determined using the function ‘FindAllMarkers’. Visualization in R required library packages cowplot v1, dplyr v0.8.5, ggplot2 v3.3.0, patchwork v1, magrittr v1.5 and stringr v1.4. Stromal subpopulation analysis in Metascape (Zhou et al., 2019) used all positively differentiated genes in the stromal clusters.

E16 and adult salivary glands were removed from CD-1 mice and preserved following a standard method (Shubin et al., 2017). Glands were fixed with 4% paraformaldehyde (PFA; Electron Microscopy Sciences, 15710) in 1×PBS for 2 h. The fixative was washed three times with 1×PBS. The glands were then dehydrated using three consecutive 1-h sucrose incubations at 4°C with 5%, 10% and 15% sucrose (Fisher, 57-50-1). Another two overnight incubations at 4°C were then performed first with 30% sucrose, then 15% sucrose and 50% Tissue Freezing Medium (TFM; Electron Microscopy Sciences, 72592). The fixed and dehydrated glands were frozen in TFM by floating the samples over liquid nitrogen. Cryosections were cut at 10 µM using a Leica CM1860 cryostat and adhered to charged histological slides (Electron Microscopy Sciences, 71869-10).

Enriched E16 SMG stromal cells were pelleted at 300× g for 8 min. The media was removed and the stroma was resuspended in 90 µl FACS buffer [1×Hanks' buffered saline solution (HBSS); Life Technologies, 14175095], 10% FBS and 5 mM EDTA (Falcon, S311-100), and positively selected using 10 µL PDGFRα antibody-coated microbeads (1:10; Miltenyi Biotec, 130-101-547). PDGFRα⁺ stoma was placed into MACS immunomagnetic cell separation columns (Miltenyi Biotec, 130-042-021) and isolated following the manufacturer's instructions. PDGFRα⁺ stomal cells were pelleted at 300× g for 8 min, resuspended in 1 ml of media, and counted using a hemocytometer.

Approximately 900 epithelial clusters (derived from one gland) grown with or without 20,000-50,000 stromal cells (derived from 0.2 glands) were embedded in Matrigel (Corning, 356234) at a 1:1 ratio to a total of 10 µl and seeded on 0.1 μm pore porous polycarbonate filters (Whatman Nuclepore, 0930051) floated on top of media in 50 mm glass-bottom dishes (MatTek P50G-1.5-14F). The cell preparations were incubated at 37°C in a tissue culture incubator (Thermo Fisher Scientific Forma Series II) with 5% CO₂ for 7 days to form organoids. Medium was replaced after incubating for 4 days. Media used in these experiments included DMEM/F12 with 10% FBS and Pen-Strep with growth factors or inhibitors added. Growth factors included: amphiregulin (AREG, 100 ng/ml; Peprotech, 315-36), BMP2 (100-200 ng/ml; Peprotech, 120-02), BMP7 (100-200 ng/ml; Peprotech, 120-03P), FGF2 (100 ng/ml; Peprotech, 450-33), EREG (100 ng/ml; Peprotech, 100-04) and TGF-β1 (TGFβ, 5 ng/ml; Peprotech, 100-12). Growth factors were solubilized in 0.2% bovine serum albumin (BSA; Sigma Millipore, A2934-100G) in water and stored at −20°C in single-use aliquots prior to use in epithelial cell cultures.

Primary E16 SMG stromal cells were cultured in 10 cm tissue culture dishes (Corning, 430167) containing 10 ml DMEM/F12+10% FBS+Pen-Strep with media replenished after plating and then every 2 days. Primary adult stromal cells were cultured in 35 mm tissue culture dishes (Corning, 353001) containing 2 ml DMEM/F12+10% FBS+Pen-Strep+50 ng/ml FGF2 with medium replaced every 3 days for 14-21 days until use in organoid preparations. All incubations were at 37°C with 5% CO₂ and cells were passaged when the density reached 90%. Cells were passaged and isolated using standard trypsinization (Gibco, 25200056) techniques. For growth factor treatments, 50,000 SMG stromal cell (∼0.2 glands’ worth of cells) at passage 1 were seeded on top of 50 mm glass bottom dishes (MatTek, P50G-1.5-14F) on day 0 and growth factors or pharmacological inhibitors were added from day 1 to 4.

For organoid cultures, the media under the filter was replaced with 4°C FACS buffer, which was then replaced with FACS buffer on top of the organoids, which were incubated at 4°C for 20 min, triturated again, and collected into a 15 ml conical tube containing an equal volume of 4°C FACS buffer. The cells were triturated at 4°C to separate the cells and dissolve the Matrigel. Epithelial clusters were collected by centrifuging at 10× g for 1 min and resuspended in 4°C FACS buffer. Total RNA was extracted from these enriched epithelial cells for microarray analysis. Single stromal cells were pelleted at 450× g for 5 min and then resuspended in 4°C FACS buffer. The stroma was further negatively selected against epithelium using EPCAM antibody-coated microbeads (1:10; Miltenyi Biotec, 130-061-101) and a Miltenyi MS column. Either total RNA was extracted from the stromal cells for microarray analysis or the cells were used for scRNA-Seq.

Total RNA was isolated using RNeasy mini kits (Qiagen, 74134) or RNAqueous micro kits (Thermo Fisher Scientific, AM1931). RNA quality was analyzed using a UV/Vis Spectrophotometer (NanoDrop ND-1000) at 260/280 nm and an Agilent 2100 Bioanalyzer. A 260-280 ratio of 1.7-2.2 and RIN score >7.0 were used as cut-offs, respectively. The RNA for transcriptome analyses were prepared with Clariom-S (Epithelium) (Thermo Fisher Scientific, 902930) or Clariom-D (Stroma) (Thermo Fisher Scientific, 902513) mouse microarrays. Data were analyzed using Transcriptome Analysis Console v.4.0.2.15, R functions and the Broads Institute's GSEA software v4.

SMG organoid and cell cultures grown on Whatman Nuclepore filters (Cytiva, 110405) or glass coverslips were fixed with 4% PFA in 1×PBS for 20 min or in −20°C methanol for 18 min. Downstream antibody staining determined the fixation method. ICC was performed as described previously (Daley et al., 2009; Sequeira et al., 2012), with the use of 0.4% Triton-X 100 (Sigma, T9284-100ML) for PFA-fixed samples. All primary antibody incubations were overnight at 4°C. Secondary antibody incubations were for 1-3 h at room temperature. The primary antibodies and dilutions used were: AQP5 (1:400; Alomone, AQP-005), CNN1 (1:600; Abcam, ab46794), KRT7 (1:200; Abcam, ab9021), EPCAM-647 (1:400; Biolegend, 118212), PDGFRα (1:200; Thermo Fisher Scientific, 14-1401-81), PDGFRα (1:100; R&D Systems, AF1062-SP), PDGFRβ (1:200; Abcam, ab32570), PDGFRβ (1:100; Thermo Fisher Scientific, ab32570), HSPG2 (1:200; SantaCruz, sc-33707), αSMA/ACTA2 (1:1000; Sigma, 12-1402-81), THY1 (1:200; Biolegend, 105301), ENG (1:100; Thermo Fisher Scientific, 14-1051-81) and VIM (1:2000; Sigma, V2258). Secondary antibodies were Cyanine and Alexa dye-conjugated AffiniPure F(ab″)₂ fragments purchased from Jackson ImmunoResearch Laboratories and used at a dilution of 1:250. DAPI (1 µg/ml: Life Technologies, D1306) was used for nuclei staining together with secondary antibodies. Organoids were mounted on Whatman Nuclepore filters using antifade mounting media containing 90% (vol/vol) glycerol (Sigma, G5516-1L) in 1×PBS with 1-4 diazobicyclo[2,2,2]octane (Sigma, D27802-100G) and n-propyl gallate (Sigma, P3130-100G) as antifade agents (Gerdes et al., 2013; Valnes and Brandtzaeg, 1985). Imaging was performed using Zeiss Z1 Cell Observer widefield (for most experiments) or a Zeiss LSM 710 confocal or a Leica TCS SP5 confocal microscopes (for high-resolution imaging) at 10×, 20×, 40× and 63× (oil immersion) magnifications, with the same microscope and laser configurations for all samples within an experiment, as appropriate to the imaging method.

SEURAT v3.1.5 (Satija et al., 2015; Stuart et al., 2019) uses nonparametric Wilcoxon rank sum tests for statistical significance. Microarrays used RMA normalization and log2 transform in Transcriptome Analysis Console v4.0.2 (Thermo Fisher Scientific) or in R (R Core Team, 2021). Adjusted P-values of <0.05 and log fold change >1.99 were parameters used to determine differentially expressed genes. For ICC images, FIJI v1.53c performed the rolling-ball background subtraction and then color thresholding was used to calculate the stained area. To determine the staining area for organoid stromal measurements, regions of interest drawn by hand were used to include only stromal cells. To determine the staining area from organoid epithelial measurements, color thresholding on DAPI was used to select the more intensely DAPI-stained epithelial regions. For image quantification, at least one representative image was captured from each sample at a 5× widefield magnification that framed 20-40% of the entire sample. Statistical significance (P<0.05) was calculated between stained areas using a Student's t-test for dual-sample comparison or single-factor ANOVA followed with post-hoc Tukey HSD test for multisample comparisons. Statistical tests were performed in Microsoft Excel (Microsoft Corporation) or R. Staining area results are expressed as mean±standard deviation (s.d.).

The authors thank the University at Albany RNA Institute and Life Science Core facilities and the Center for Functional Genomics for equipment and expertise in both microarray and sequencing data creation. Specifically, Dr Sridar Chittur and Marcy Kuentzel.

The peer review history is available online at https://journals.biologists.com/dev/article-lookup/doi/10.1242/dev.200167.