Transcription factor FOXF1 identifies compartmentally distinct mesenchymal cells with a role in lung allograft fibrogenesis

It is increasingly recognized that characterization of mesenchymal cell (MC) populations within distinct anatomic locations is key to understanding local homeostatic and regenerative processes as well as delineating unique pathogenic mechanisms of diverse fibrotic diseases. The lung is a vital organ with anatomically and functionally distinct conducting and gas-exchange units that harbor histologically identifiable, specialized cell populations. The presence of resident mesenchymal progenitors in human adult lungs was highlighted by our finding of donor-derived multipotent mesenchymal stromal cells (MSCs) in sex-mismatched allografts (1). Subsequent studies using transgenic mice and single-cell RNA-Seq (scRNA-Seq) have shed light on various MC subpopulations residing within the lung (2–4). However, despite the identification of a plethora of genes, the characterization of anatomically and functionally distinct collagen-expressing MC populations remains limited. Furthermore, the majority of studies have focused on interstitial pulmonary fibrosis using bleomycin-induced injury models, with a lack of investigation into other fibrotic processes such as the conducting airway fibrosis that marks chronically rejecting lung allografts.

Here, we used the expression of forkhead box F1 (FOXF1), a key transcription factor expressed in developing lung mesenchyme, to identify a unique subset of collagen 1⁺ (Col1⁺) MCs that resided within a specific peribronchial niche in an adult lung and contributed to allograft fibrogenesis. FOXF1 is a Sonic hedgehog/glioma-associated oncogene 1 (SHH/GLI1) target expressed in the lateral mesoderm during embryogenesis that is essential for branching morphogenesis and lung development (5, 6). FOXF1 is unique among embryonic transcription factors, in that its expression persists in adult lungs (5, 7). We have previously demonstrated that MSCs derived from human adult lung allografts express FOXF1 (1, 7, 8). However, whether persistent expression of FOXF1 in adulthood marks all or just a subset of collagen-expressing MCs is unknown, and the in situ niche of these cells remained to be determined. The present studies of Foxf1-tdTomato bacterial artificial chromosome (BAC) and bitransgenic Foxf1-tdTomato Col1a1-GFP mice reveal that FOXF1 expression demarcated a distinct collagen-expressing MC population with a unique transcriptional profile and anatomic localization. The demonstration of an in situ peribronchial subepithelial niche of FOXF1⁺Gli1⁺ cells, their in vitro ability to support epithelial cell organoids, and their in vivo expansion in response to airway-centric injury after lung transplantation reflects their participation in homeostatic/fibrotic responses within this specific compartment.

The embryonic transcription factor Foxf1 is retained in the Sca1⁺CD34⁺ mesenchymal stromal cell subset of Col1-expressing cells in the adult lung. To characterize the Foxf1-expressing, lung-resident MC population, we generated C57BL/6J-Tg(Foxf1-tdTomato) (referred to herein as Foxf1-tdTomato) mice using a modified BAC, in which exon 1 of the Foxf1 gene was replaced with the coding sequence of tdTomato (Supplemental Figure 1A; supplemental material available online with this article; https://doi.org/10.1172/JCI147343DS1). We identified a distinct tdTomato fluorescent population in the adult transgenic mouse by FACS analysis, with the viable sorted cells showing higher expression of Foxf1 and collagen 1α1 (Col1a1) (Figure 1A). The tdTomato expression faithfully recapitulated endogenous expression of Foxf1, with tdTomato⁺ cells demonstrating a lack of expression of the epithelial marker EpCAM and the immune marker CD45 (Figure 1B). As expected, the tdTomato⁺ population encompassed CD31⁺ endothelial cells, which are known to arise from embryonic lung mesenchyme and express FOXF1 in adulthood (9). Immunostaining revealed similar staining patterns for endogenous FOXF1 protein and tdTomato fluorescence driven by BAC-Foxf1 (Supplemental Figure 1B). The Lin(CD45/CD31/EpCAM)^–Foxf1-tdTomato⁺ population uniformly expressed the mesenchymal surface marker PDGFRα, and the progenitor cell markers stem cell antigen 1 (SCA1) and CD34 (Figure 1, C and D). A majority of Lin^–Foxf1-tdTomato⁺ cells also stained positive for CD44, whereas CD73 staining was limited to approximately one-third of the cells (Figure 1D). We further confirmed the mesenchymal phenotype of Lin^–Foxf1-tdTomato⁺ cells by staining for vimentin and prolyl-4-hydroxylase (Supplemental Figure 1C).

Figure 1

Foxf1-tdTomato labels Sca1⁺ MCs in the adult mouse lung. (A) Flow cytometric analysis of single-cell suspensions from lungs of adult C57BL/6J-Tg(Foxf1-tdTomato) BAC mice. Real-time PCR analysis of Foxf1-tdTomato⁺ (TD⁺) or Foxf1-tdTomato^– (TD^–) sorted cells. (B) TD⁺ populations from A were analyzed for the indicated markers. (C and D) Lin^–(CD45^–CD31^–EpCAM^–) Foxf1-tdTomato⁺ or Foxf1-tdTomato^– cells from Tg(Foxf1-tdTomato) BAC mice were analyzed for the indicated cell-surface markers. (E–H) Flow cytometric analysis of adult lungs from Foxf1-tdTomato Col1GFP mice. Foxf1-tdTomato Col1⁺ MCs were sorted and analyzed for FOXF1 expression by Western blotting (F) and CFU assay (G). n = 7. (H) SCA1 and CD34 expression in Lin^–Col1GFP⁺PDGFRα⁺Foxf1-tdTomato cells. (I) Characterization of the Lin^–SCA1⁺ cell population in adult Foxf1-tdTomato BAC mice by flow cytometry (n = 7). (J) SCA1 cell-surface expression in Foxf1-overexpressing (OE) murine lung fibroblasts. Data indicate the mean ± SD. *P < 0.05 and ***P < 0.001, by 2-tailed, unpaired t test (A and G). SSC-A, side scatter area.

To determine whether FOXF1 is expressed in all lung-resident Col1⁺ MCs or only in a MC subpopulation, Foxf1-tdTomato mice were crossed with Col1a1 promoter–driven GFP (referred to herein as Col1GFP) mice. FACS revealed 2 distinct Col1GFP populations demarcated by the presence or absence of Foxf1-tdTomato in Foxf1-tdTomato Col1a1-GFP–bitransgenic adult mouse lungs (Figure 1E). FOXF1 protein expression was restricted to the Foxf1-tdTomato⁺ subset of Col1–expressing MCs (Figure 1F). The Foxf1⁺Col1⁺ MCs formed significantly more CFU (Figure 1G) and had higher SCA1 and CD34 expression than that detected in the Foxf1^–Col1⁺ subpopulation (Figure 1H). The Lin^–SCA1⁺ cell population has been shown to mark multipotent mesenchymal progenitors in an adult mouse lung (10). Flow cytometric analysis of cells from adult Foxf1-tdTomato BAC mice showed that the Lin^–SCA1⁺ cells uniformly expressed FOXF1 (Figure 1I). Transient overexpression of Foxf1 in murine lung fibroblasts induced SCA1 expression, further suggesting a positive correlation between these markers in MCs (Figure 1J). Together, these data indicate that expression of the embryonic lung mesenchyme–associated transcription factor FOXF1 was restricted to a subset of collagen-expressing MCs in an adult mouse lung and overlapped with the progenitor cell marker SCA1.

Topographic confinement of FOXF1 to a 3D MC network along the conducting airways. Next, to investigate the localization of the FOXF1-expressing MC population in an adult lung, we performed FOXF1 staining of lung tissue from Col1GFP-transgenic mice. As can be appreciated by zooming in on the high-resolution wide-field image (Nikon Elements), colocalization of FOXF1 staining (red) with intrinsic collagen 1α1-GFP (green) was noted to be confined to spindle cells in the bronchovascular areas, in association with bronchi and vessels, and interspersed in the adjoining loose interstitium (Figure 2A). FOXF1⁺Col1GFP⁺ MCs lay in close apposition to club cell secretory protein–expressing (CCSP-expressing) bronchial epithelial cells (white), suggesting that these MCs are important components of the airway progenitor cell niche (Figure 2B and Supplemental Figure 2A). This network of cells extended in a 3D manner from the proximal to the distal airways (Figure 2B, insets 1–3). Subepithelial MCs in the lung airways have been identified histologically as distinct intermeshed, spindle-shaped cells; the unique appearance of this layer of thin MCs, with its large surface area and close proximity to the epithelial/environmental interface, led it to be termed the attenuated fibroblast sheath, with similar cells in the gut identified as the pericryptal fibroblast sheath (11). Interestingly, along with the finding that FOXF1 marked the MCs in this anatomically distinct location in the lung, our Foxf1-tdTomato BAC mice had expression of FOXF1 in the pericryptal MCs in the gut (data not shown).

Figure 2

In situ niche of FOXF1-expressing Col1⁺ MCs. (A–C) Adult lungs of Col1GFP mice were stained with FOXF1 (red) and CCSP (white). Representative high-resolution wide-field image of entire lung (A) and 4 locations along the bronchial tree (B). Scale bars: 50 μm (B). Arrows indicate FOXF1⁺Col1⁺ MCs along the BVB. (C) Quantification of nuclear FOXF1 intensity in Col1GFP⁺ cells in the bronchovascular space (Bronchi) and the alveolar interstitium (Alveoli). n = 3. (D) Confocal images and quantification of PDGFRα-GFP mouse lung tissue stained with anti-FOXF1 antibodies. Images show colocalization of FOXF1 (red) and PDGFRα (green) in the peribronchial region (cyan arrows), with no costaining noted in alveoli. n = 3. Scale bars: 20 μm. (E) Adult human lung tissue sections stained with FOXF1 (green) and the marker anti–smooth muscle actin (αSMA; red). Scale bars: 50 μm. AW, airway; BV, blood vessel. Data indicate the mean ± SD. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001, by 1-way ANOVA with Bonferroni’s test C) and unpaired t test (D).

Coexpression of FOXF1 was notably absent in the Col1GFP⁺ cells in the alveoli, where FOXF1 positivity was restricted to endothelial cells (Figure 2B, inset 4, and Supplemental Figure 2B). We performed confocal imaging and costaining of type 1 and type 2 pneumocytes with aquaporin 5 and thyroid transcription factor 1 (TTF-1) to confirm the absence of FOXF1 expression in MCs of the alveolar niche (Supplemental Figure 2C). On quantitation, FOXF1 nuclear staining was evident in 76.4% of Col1GFP⁺ cells within the bronchovascular unit, whereas 93.9% of MCs in the alveolar space were demonstrably negative for FOXF1 nuclear expression (Figure 2C).

We used the transgenic PDGFRα^EGFP mouse, in which the Pdgfra promoter drives expression of the H2B-eGFP fusion gene, to corroborate these findings. In the PDGFRα^EGFP adult mouse lung, FOXF1 staining colocalized with endogenous PDGFRα-GFP nuclear expression in MCs near the airways but was absent in PDGFRα⁺ alveolar MCs (Figure 2D). Immunofluorescence staining of adult human lung tissue sections confirmed the peribronchiolar niche of the FOXF1-expressing MCs, where FOXF1 expression was noted in MCs interspersed between the smooth muscle and bronchial epithelium (Figure 2E).

Bulk sequencing and scRNA-Seq confirm distinct transcriptional signatures of Foxf1⁺– and Foxf1^––resident Col1a1⁺ MC populations. The transcriptomic profiles of freshly sorted Col1GFP⁺Foxf1-tdTomato⁺ or Col1GFP⁺Foxf1-tdTomato^– MCs were compared by Affymetrix assay. Col1^–Foxf1⁺ endothelial cells were included as controls (Supplemental Figure 3A). Principal component analysis and unsupervised clustering indicated that the 2 MC populations have unique gene signatures (Figure 3A and Supplemental Figure 3B). Consistent with their unique anatomic niches, a comparison of biological processes gene ontology (GO) terms related to epithelium revealed statistically significant differences between the 2 populations (Supplemental Table 1). A previous study classified 2 distinct matrix fibroblast populations based on Col13a1 and Col14a1 expression (3). We found that Col1GFP⁺Foxf1-tdTomato⁺ cells expressed a Col14a1 matrix fibroblast signature (Col14a1, Pi16, Cygb, Rtp4, Rbp4, Meg3), whereas Col1GFP⁺Foxf1-tdTomato^– cells had a Col13a1 matrix fibroblast signature (Col13a1, integrin α 8 [Itga8], Npnt, Cxcl14, Cdh11, Tcf21) (Figure 3B). In addition to Foxf1, the SHH signaling gene Gli1, a transcriptional promoter of Foxf1 expression (12), was upregulated in Col1GFP⁺Foxf1-tdTomato⁺ compared with Col1GFP⁺Foxf1-tdTomato^– MCs (Figure 3B). This was further investigated by using tamoxifen-treated Gli1^CreERT2/WT Rosa^EYFP/WT adult mice. Flow cytometric analysis revealed that, similar to Foxf1-tdTomato⁺ MCs, greater than 80% of Gli1^CreERt2-EYFP–labeled cells were SCA1⁺ (Supplemental Figure 3C). Immunofluorescence staining showed nuclear FOXF1 expression in Gli1-GFP–expressing subepithelial peribronchial MCs; no significant GLI1 expression was noted in the alveoli (Figure 3C).

Figure 3

Bulk and single-cell transcriptomic analysis of Foxf1_MCs. (A) Principal component (PC) analysis plot from affymetrix analysis of sorted populations of MCs from adult Foxf1-tdTomato Col1GFP mice. (B) Heatmap comparing Foxf1-tdTomato⁺ (TD⁺) and Foxf1-tdTomato^– (TD^–) CD45^–Col1⁺ MCs for selected genes. (C and D) FOXF1 immunostaining of adult lung tissue from Gli1^CreERT2/WT Rosa26^mTmG/WT (C) and Col1GFP (D) mice. n = 3. Scale bars: 50 μm and 10 μm (insets). Arrowheads indicate FOXF1⁺Gli1GFP-expressing subepithelial peribronchial MCs (C) and ITGA8⁺FOXF1^–Col1GFP⁺ MCs in the alveoli (D). (E) Uniform manifold approximation and projection (UMAP) representation of 11,118 scRNA-Seq profiles; each dot represents a single cell. Cells are colored to indicate their membership in 1 of 14 clusters. (F) UMAP plots of specific genes. The coordinates are the same as in E, but the color indicates the expression level of the gene within each cell. (G) Correlation plot from scRNA-Seq analysis. Pearson’s correlation test.

Affymetrix analysis identified the cell-surface molecule Itga8 as being highly upregulated in Col1GFP⁺Foxf1-tdTomato^– cells (Figure 3B). Immunofluorescence staining of Col1GFP⁺ adult mouse lung tissue sections confirmed ITGA8 costaining in Col1GFP⁺FOXF1^– cells in the alveolar space and its absence in Col1GFP⁺FOXF1⁺ cells of the peribronchiolar space (Figure 3D). FACS sorting of lung MCs based on ITGA8 expression confirmed lower mRNA expression of Foxf1 and Gli1 in ITGA8⁺ MCs (Supplemental Figure 3, D and E). Freshly sorted Lin^–PDGFRα⁺ITGA8⁺ and Lin^–PDGFRα⁺ITGA8^– cells were studied by time-lapse microscopy and high-content analyses of live cell migration, proliferation, and morphology (Supplemental Figure 4). We observed very distinct phenotypes, with ITGA8^–(Foxf1⁺Gli1⁺) MCs demonstrating greater migratory and proliferative capacity. The shapes and areas of the whole cells and nuclei between the 2 cell populations were significantly different, as revealed by single-cell linear discriminant analysis (Wilcoxon rank-sum test; P < 0.0001).

We next performed scRNA-Seq on CD45^–CD31^– cell populations sorted from adult mouse lungs. Using non-negative matrix factorization and Louvain community detection, we identified 14 transcriptionally distinct cell populations (Figure 3E and Supplemental Figure 5). Most cells belonged to 2 large clusters characterized by high expression of Col1a1, with high expression of Foxf1 or Itga8, and we refer to them herein as Foxf1_MCs and Itga8_MCs, respectively (Figure 3, E and F). Consistent with our results above, the scRNA-Seq data indicated a significant correlation (Pearson’s P = 0.64, P < 2^–16) between the expression of Foxf1 and Gli1 within the same cell (Figure 3G). Foxf1_MC and Itga8_MC clusters also distinguished themselves by the expression of discriminatory Col13a1 and Col14a1 genes previously identified on scRNA-Seq analysis by Xie et al. (ref. 3 and Figure 3F). The Foxf1_MC cluster encompasses subbronchial and both subpopulations of cuff Col1⁺ MCs described by Tsukui et al. (2). The RNA velocity vectors for the MC clusters revealed an interesting trend: the arrows point from Foxf1_MC→Itga8_MC, suggesting that transcriptional changes proceed primarily in this direction, rather than from Itga8_MC→Foxf1_MC or in random directions (Supplemental Figure 6). The continuous nature of the Foxf1 and Itga8 MC populations, combined with the velocity analysis, further raises the intriguing possibility that there may be a spatial gradient of transcriptional variation among these cells.

Foxf1⁺Gli1⁺ MCs are significant contributors to airway fibrosis in an orthotopic mouse lung transplant model. Together, this work identified a Foxf1⁺Gli1⁺Itga8^–Col1⁺ MC population residing in a unique niche associated with BVBs (BVB MCs). BVBs were also the primary site of the allogeneic insult in a transplanted lung, with acute rejection marked by lymphocytic infiltration of blood vessels and airways and chronic rejection presenting as peribronchiolar and perivascular fibrosis. Our human studies revealed the mobilization of Foxf1-expressing lung MSCs during the post-transplant injury period, with an increase in their numbers early after transplantation and in chronically rejecting grafts (8). We have also recently demonstrated that Foxf1 acts as a suppressor gene in human lung MSCs, with its downregulation promoting cellular proliferation and migration (13). As FOXF1 is a SHH/GLI1 target (6, 12), we hypothesized that Shh expression in bronchial epithelial cells is key to maintaining Gli and Foxf1 expression in BVB MCs and that its dysregulation marks airway fibrosis.

We first studied the proposed Shh/Gli/Foxf1 paracrine signaling between epithelial and mesenchymal components in vitro using epithelial organoid formation. Col1GFP⁺Foxf1-tdTomato⁺ MCs promoted the formation of organoids when cocultured with freshly isolated Rosa26^mTmG mouse bronchioalveolar stem cells (BASCs) in Matrigel (47.5 ± 3.29/well, n = 6; Figure 4A). Gene expression analysis confirmed the exclusive expression of Shh and Foxf1 in the sorted epithelial and mesenchymal component of the organoids, respectively (Figure 4B). The organoid-forming capacity of Col1GFP⁺ Foxf1-tdTomato⁺ MCs was significantly higher than that of Col1GFP⁺Foxf1-tdTomato^– MCs (Supplemental Figure 7, A and B). Foxf1_MCs in organoid cocultures had higher Foxf1 and Gli1 mRNA expression than did those grown in Matrigel alone (Figure 4C). The presence of the Shh signaling and the Smoothened antagonist LDE225 inhibited Foxf1 and Gli1 mRNA expression in Foxf1_MCs grown in cocultures compared with those grown in Matrigel alone, demonstrating a contribution of epithelium-derived SHH in maintaining the expression of Foxf1 and Gli1 in Foxf1_MCs (Figure 4C).

Figure 4

FOXF1 and Gli1 expression and expansion of Gli1⁺ MCs in an orthotopic mouse lung transplant model. (A) Generation of 3D epithelial organoids in coculture of Foxf1-tdTomato⁺Col1GFP⁺ MCs and BASCs from Rosa26^mTmG mice. Scale bars: 500 μm. (B) Flow cytometry–sorted cell populations from organoids were analyzed by real-time PCR. (C) Foxf1⁺_MCs were cocultured with BASCs, with or without the Shh signaling inhibitor LDE225 (0.75 μM), and compared with Foxf1_MCs cultured in Matrigel alone. (D) Mesenchymal (Lin^–PDGFRα⁺) and epithelial (Lin^–EpCAM⁺) cells sorted from isografts and allografts in the orthotopic mouse lung transplant model were analyzed for the indicated genes by real-time PCR. D28, day 28. (E) Schema for the generation of Gli1^CreERT2/WT Rosa26^mTmG donor mice and lung transplant experiments. Tam, tamoxifen. (F) Representative immunofluorescence images of lineage-traced Gli1GFP⁺ MCs in 2 control and 3 transplanted lungs. Scale bars: 100 μm (zoom, cropped: 50 μm). Data indicate the mean ± SD. *P < 0.05, **P < 0.01, and ***P < 0.001, by 2-tailed unpaired t test (B and D) and 1-way ANOVA with Bonferroni’s test (C).

Next, we used a recently described F1→parent orthotopic murine lung transplant model of chronic lung allograft dysfunction to investigate whether disruption of Shh/Gli1/Foxf1 crosstalk accompanies airway remodeling. This B6D2F1/J→C57BL6/J model of lung transplantation demonstrates an evolution from marked lymphocytic infiltration with bronchial epithelium and endothelial injury to robust peribronchial and perivascular fibrosis (14). Affymetrix analysis of the whole-lung lysate at various time points after transplantation revealed a persistent decrease in Shh and Foxf1 expression in the allografts compared with expression in the isografts (Supplemental Figure 7C). To confirm their loss in specific cellular compartments, epithelial cells (CD45^–CD31^–PDGFRα^–EpCAM⁺) and MCs (Lin^–PDGFRα⁺) were sorted from graft single-cell suspensions on post-transplantation day 28. We observed a loss of Gli1 and Foxf1 mRNA expression in MCs and a concomitant decrease in Shh expression in epithelial cells in the allografts (Figure 4D). To track BVB MCs undergoing a loss of Foxf1 and Gli1 expression after lung transplantation, we elected to stably label these cells using the Gli1^CreERT2/WT Rosa26^mTmG/WT transgenic mouse (Figure 4E). By day 14, Gli1^CreERT2/WT Rosa26^mTmG/WT donor lung allografts showed a robust peribronchiolar and perivascular expansion of Gli1-labeled cells, with loss of CCSP-expressing epithelial cells and subbronchial fibrosis (Figure 4F). The day-28 post-lung transplantation allografts had a similar extensive expansion of Gli1⁺ MCs (Supplemental Figure 7, D and E). ITGA8 staining of these allografts demonstrated a continued lack of ITGA8 in Gli1⁺ airway MCs as well as no significant expansion or myofibroblast differentiation in ITGA8⁺ alveolar MCs, suggesting that fibrotic activation was limited to BVB MCs in the rejecting allografts (Supplemental Figure 7F). Comparison of MC subsets (Lin^–PDGFRα⁺ITGA8^– BVB MCs and Lin^–PDGFRα⁺ITGA8⁺ alveolar MCs) sorted from isografts and allografts confirmed the upregulation of Col1a1 mRNA in BVB MCs in allografts, with no change noted in alveolar MCs (Supplemental Figure 7G). Similarly, only BVB MCs showed upregulation of autotaxin (Atx), a lysophospholipase D enzyme that produces profibrotic mediator lysophosphatidic acid and has been shown to be key in autocrine fibrotic activation of lung MCs and allograft fibrogenesis (15). FOXF1 was identified in our recent study as a transcriptional repressor of ATX in human lung MSCs (13), suggesting a unique fibrotic differentiation program in rejecting allografts, with loss of FOXF1 in BVB MCs leading to ATX upregulation and fibrotic expansion.

Taken together, our results suggest that anatomically localized Foxf1⁺Gli1⁺Itga8^– BVB MCs are the key contributors to fibrogenesis in a rejecting allograft. Our data on the downregulation of Shh/Gli1/Foxf1 signaling during allograft fibrogenesis are consistent with prior studies, in which Shh was shown to regulate airway quiescence, with its inhibition leading to peribronchial fibrosis (16). The present study, combined with our recent demonstration of repressive functions of FOXF1 and its role as a brake on MC activation (13), suggests that FOXF1 is a likely mediator of SHH-induced quiescence of the airway compartment.

The identification in this study of simple cell-surface markers that isolate phenotypically, transcriptionally, anatomically, and functionally distinct BVB and alveolar MCs lays the foundation for future work investigating specific mechanisms by which these cells support their respective epithelial niches and compartmentalized fibrotic injury. Our data will also reinvigorate the interest in mesenchymal progenitors, as they identify the niche along the conducting airways of a FOXF1⁺SCA1⁺ MC population with CFU potential and an ability to sustain epithelial organoids. Our RNA velocity analysis of scRNA-Seq data identified the transcriptional changes proceeding in the direction from bronchovascular compartment–localizing Foxf1 MCs to alveolar compartment–residing Itga8 MCs. Further studies are needed to ascertain the role of these cells as potential progenitors for other MC populations and to characterize further specific subpopulations within this population.

Further information can be found in the Supplemental Methods.

Study approval. All animal experiments were performed in accordance with protocols approved by the IACUC of the University of Michigan.

Accession number. RNA-Seq and scRNA-Seq data were deposited in the NCBI’s Gene Expression Omnibus (GEO) database (GEO GSE176476, for the RNA-Seq data presented in Figure 3, A and B, and Supplemental Figure 3, A and B, and GSE179034, for the scRNA-Seq data presented in Figure 3, E–H, and Supplemental Figure 5).