Authors: Sophie E. Edelstein, Satoshi Mizoguchi, Maria Tomàs Gracia, Nuoya Wang, Vi Lee, Tomohiro Obata, Hahram Kim, Connor Haynes, Colten Danelski, Tomoshi Tsuchiya, Maor Sauler, Micha Sam Brickman Raredon
Categories: Article, Immunology, Bioengineering, Tissue engineering, Developmental biology
Source: iScience
Authors: Sophie E. Edelstein, Satoshi Mizoguchi, Maria Tomàs Gracia, Nuoya Wang, Vi Lee, Tomohiro Obata, Hahram Kim, Connor Haynes, Colten Danelski, Tomoshi Tsuchiya, Maor Sauler, Micha Sam Brickman Raredon
Aberrant epithelial regeneration and immune remodeling are hallmarks of chronic lung diseases such as idiopathic pulmonary fibrosis and chronic obstructive pulmonary disease. How cellular context shapes these trajectories remains unresolved. We present a lung organoid model that varies immune, epithelial, and mesenchymal inputs to reveal how composition dictates epithelial plasticity and macrophage polarization. We observed condition-dependent emergence of transitional cell states, including Sox9^+^ stressed progenitors, RAS-like intermediates, and hillock-like cells, alongside macrophage activation profiles. In mesenchyme-rich contexts, epithelial-immune-mesenchymal crosstalk reinforced inflammatory signaling and stabilized transitional cells, while immune-dominant inputs favored ATI-like repair and squamous remodeling. Hillock-like cells displayed context-dependent activation and expressed immune-regulatory genes, suggesting a role as epithelial orchestrators calibrating inflammatory response during regeneration. Regenerative outcomes were associated with multicellular signaling networks integrating stress sensing, immune coordination, and epithelial resilience. This platform facilitates modeling of milieu-specific regenerative mechanisms and informs strategies to redirect epithelial fate in chronic lung disease.
Lung homeostasis, injury repair, and disease progression are all governed by epithelial plasticity.^1^ Regeneration of the lung is regulated by a finely tuned system of interactions between epithelial, immune, and mesenchymal cells. Based on cues from their microenvironment, epithelial cells transition through multiple phenotypic and transcriptomic states, giving rise to transitional populations that include, but are not limited hillock, activated respiratory airway secretory (aRAS), alveolar type 0 (AT0), Sox9^+^ progenitors, bronchoalveolar stem cells (BASCs), and damage-associated transient progenitor (DATP) cells.^2^^,^^3^^,^^4^^,^^5^^,^^6^ These transitional cell types have become of significant interest to the pulmonary biology community in recent years due to their flexible lineage potential, transient activation following injury, and capacity to mediate divergent outcomes ranging from effective regeneration to maladaptive remodeling.^6^^,^^7^ Importantly, many of these same transitional programs are activated in chronic lung diseases such as idiopathic pulmonary fibrosis (IPF), chronic obstructive pulmonary disease (COPD), and asthma, where they have been implicated in aberrant repair and progressive tissue remodeling.^2^^,^^3^^,^^4^^,^^6^^,^^8^
Traditional 2D epithelial culture systems fail to recapitulate native tissue architecture and lack immune and mesenchymal components, which are essential for regulating epithelial steady state and fate.^9^^,^^10^ Conversely, while in vivo animal models provide a physiologically relevant context, they do not allow for the systematic dissection of how specific cellular lineages, and the specific ratios of these lineages, shape epithelial behavior under controlled conditions.^11^
Mesenchymal populations are central regulators of epithelial and immune cell behavior, mediating a vast array of context-dependent behaviors.^12^^,^^13^ Epithelial-mesenchymal signaling coordinates spatial patterning and cell fate decisions in the developing lung, primarily through canonical pathways such as WNT and BMP.^14^^,^^15^^,^^16^ Moreover, mesenchymal cells are also capable of modulating immune cell polarization and cytokine production, thus modifying the inflammatory milieu and calibrating repair pathways.^17^^,^^18^ These regulatory axes come together during injury and involve epithelial and immune coordination through IL-33, IL-1β, and TNF-α to orchestrate tissue remodeling.^19^^,^^20^^,^^21^ That being said, the specific mechanisms by which mesenchymal presence or absence dictates epithelial-immune coordination in the rat lung are difficult to identify in vivo.
Pulmonary organoids are promising models for studying epithelial plasticity and multicellular dynamics in defined systems.^22^^,^^23^ However, most depend on genetic manipulations, exogenous injury cues, or pre-programmed differentiation protocols to induce epithelial transitions.^24^^,^^25^^,^^26^ Here, we present a lung organoid model that leverages primary, rat-derived populations and allows for the study of how epithelial plasticity is linked to immune cell behavior.
To investigate the influence of the relative lineage ratios at seeding on epithelial plasticity under uniform culture conditions, we isolated and characterized two distinct populations from the adult rat (1) a pulmonary dissociation (PD), generated by the enzymatic digestion of whole lung tissue and containing a heterogeneous mix of epithelial, mesenchymal, endothelial, and immune cells; and (2) a bronchoalveolar lavage (BAL) isolate collected via perfusion of the lungs, through the trachea, with cold saline to collect non-adherent immune cells from the surfaces of the airways. These populations were combined in different ratios to generate three organoid 100% PD (PD_3D), 100% BAL (BAL_3D), and a 1 mix of BAL to PD (Mixed_3D).
The spontaneous formation of reparative-like conditions occurred in all cultures, consistent with previous work demonstrating that cellular systems often remodel under stress when taken out of native physiological regulation.^27^^,^^28^^,^^29^ Despite the same media, matrix, and environmental parameters being used for all cultures, similar as well as distinct epithelial and immune cell states reproducibly emerged in the three conditions (PD_3D, Mixed_3D, BAL_3D). Among these derivatives were epithelial sub-types associated with injury, such as Krt13^+^ hillock-like cells and Scgb3a2^+^/Sftpc^+^/Scgb1a1^+^ activated RAS cells, along with macrophage polarization states that were condition-specific.^2^^,^^4^ We propose that the emergence of these cell states is an intrinsic feature of organoid culture, but that the relative weighting of epithelial and immune subsets is influenced by the initial cellular composition. Our system does not deliberately impose damage using chemical, mechanical, or physical stimuli such as bleomycin or cigarette smoke treatment.^30^^,^^31^^,^^32^ Organoid cultures, however, inherently mimic some stress-inducing characteristics such as local hypoxia, limitations in nutrient diffusion, and lack of physiological feedback, leading to the priming of cells toward a transitional or reparative state.^33^^,^^34^
In the results reported herein, we present a model system that implies that the ratios of epithelium, immune, endothelium, and mesenchyme relative to one another at the start of culture dictate the trajectory of multicellular interactions in a stress-permissive milieu, laying a foundation for understanding how lung tissues transition between homeostatic, reparative, and pathological states. Differences in the starting cell suspension composition influenced the emergence of transitional epithelial states and polarized macrophage patterns. Our findings suggest that cellular ratios can direct divergent regenerative trajectories, with immune enrichment favoring certain reparative or inflammatory cell states and mesenchymal presence facilitating epithelial state transitions. By enabling the systematic variation in discrete lineage inputs, our model provides an experimental framework to interrogate how compositional balance influences developmental and injury-associated programs—a principle with broad implications for regenerative medicine and therapeutic reprogramming.
The generation of our model started with two biologically accessible inputs derived from the native rat lung, an organ whose developmental trajectory, anatomical architecture, and cellular composition closely resemble those of the human lung and thus serves as an experimentally informative system for studying epithelial dynamics in health and disease.^35^^,^^36^^,^^37^^,^^38^ The two populations utilized were as (1) a pulmonary dissociation (PD) obtained by the enzymatic digestion of whole lung tissue resulting in a heterogeneous mixture of epithelial, mesenchymal, endothelial, and immune cells; and (2) a bronchoalveolar lavage (BAL) isolate enriched for non-adherent immune cells such as macrophages, neutrophils, B cells, and T cells, as well as a small, rare population of miscellaneous epithelium (Figure 1A).Figure 1A tunable 3D lung organoid system enables controlled modeling of epithelial-immune-mesenchymal interactions(A) Schematic overview of experimental design. Cells were isolated from rat lungs via either pulmonary dissociation (PD) or bronchoalveolar lavage (BAL) at 8 weeks of age, then combined to generate three experimental PD_3D (native lung dissociation), BAL_3D (immune-enriched), and Mixed_3D (1:1 mixture of PD and BAL cells). Organoids were maintained in LPM-3D progenitor medium through Day 10. Organoids were harvested for single-cell RNA sequencing (10× Genomics Chromium Next GEM, n = 3 replicates per condition) and processed for downstream transcriptomic and ligand-receptor signaling analysis (NICHES). See also Table S1.(B–F) Quantification of major cell class distributions across starting populations using scRNA-seq data. Bar graphs show proportions (B) major cell lineages, including a “General” category representing cycling cells that could not be definitively assigned to a single lineage, (C) epithelial subsets, (D) immune subsets, (E) mesenchymal populations, and (F) endothelial cells. See also Figures S1 and S2.(G and H) Quantification of morphological metrics from 3D organoid cultures over time. Organoid number (G) and mean organoid area (H) were quantified using QuPath on stitched brightfield images (n = 3 wells per condition per timepoint) (see organoid image analysis and quantification). ∗∗∗, p ≤ 0.001; ∗∗, p ≤ 0.01; ∗p ≤ 0.05; ns, not significant; statistical comparisons performed using Welch’s t-tests with Benjamini-Hochberg correction for multiple comparisons. See also Figure S3.(I) Brightfield images of 3D organoid cultures across time (days 1, 3, 5, 7, 9). Representative wells are shown for each condition. Insets on Day 7 show zoomed views of organoid morphology. Scale bars, 1 mm (main panels) and 200 μm (insets).
scRNA-seq of these starting populations confirmed that PD (n = 3, 9,730 cells) contained a diverse set of epithelial subsets, namely alveolar type 1 (Pdpn^+^/Ager^+^)(ATI)(3.2 ± 1.0%), alveolar type 2 (Sftpc^+^/Napsa^+^/Lamp3^+^)(ATII)(5.6 ± 1.4%), basal(Krt5^+^/Tp63^+^)(0.8 ± 1.2%), secretory (Scgb3a1^+^/Scgb3a2^+^)(11.2 ± 2.9%), ciliated (Ccdc153^+^/Pifo^+^)(9.2 ± 3.5%), and tuft cells (Dclk1^+^/Trpm5^+^)(2.3 ± 0.6%), as well as mural (Acta2^+^/Actc1^+^) mesenchyme (1.4 ± 0.3%), Msln^+^ mesothelium (4.4 ± 2.4%), arterial (Gja5^+^/Dll4^+^)(5.5 ± 1.6%), venous (Slc6a2^+^/Tmem252^+^), and capillary (Wif1^+^/Aplnr^+^) endothelial cells (1.1 ± 1.2%), and a variety of innate and adaptive immune components (Figures 1B–1F and S2D–S2X). In contrast, BAL (n = 3, 22,009 cells) was largely restricted to myeloid populations (Pparg^+^/Prodh2^+^) and lymphoid populations (Cd3e^+^/Cd79b^+^), with only a small population of epithelial cells and no mesenchyme (Figures 1B–1D and S1A–S1D).
Although the BAL starting population includes a subset of the immune populations also captured in PD, it is the proportional composition of the individual immune cell types relative to one another that makes the two starting populations most distinct from one another. The BAL isolate was dominated by alveolar macrophages (83.6 ± 1.3%), with only small fractions of B cells (0.31 ± 0.24%), T cells (2.0 ± 2.0%), and monocytes (1.8 ± 0.6%) (Figure 1D). In contrast, PD immune populations were more evenly distributed, with alveolar macrophages accounting for 33.7 ± 3.9% and B cells and monocytes together representing ∼45% of the immune cells (30.2 ± 2.6% and 14.6 ± 5.2%, respectively) (Figure 1D). PD also contained immune cell types not appreciably captured in BAL, including interstitial macrophages (5.7 ± 1.0%), NK cells (3.3 ± 2.9%), and pDCs (4.0 ± 2.7%) (Figure 1D).
All cell suspension mixtures were resuspended in Matrigel and maintained under identical culture embryonic progenitor media supplemented with broadly supportive growth factors and small molecule inhibitors (Fgf10, Egf, TGFβ inhibitor, p38-MAPK inhibitor, and GSK-3 inhibitor) in a 5% CO2, 37°C environment (Figure 1A and Table S1).^39^ To confirm that the media supported the survival of all four major lineages, we performed 2D mono-culture experiments for 7 days with isolated epithelial, immune, endothelial, and mesenchymal populations (Figures S10A and S10B). After validating lineage survival, organoid experiments were conducted using three unique starting populations of varying compositions. PD and BAL isolates were used individually or in PD_3D (100% PD), BAL_3D (100% BAL), and Mixed_3D (a 50 mix of BAL and PD), and three biological replicates were generated per condition (Figure 1A).
Organoids self-assembled in all three systems over the 10-day culture period, but each showed distinct growth kinetics and morphology (Figures 1G and 1H). By day 3, PD_3D organoids had expanded rapidly, in both number of organoids (count) and average surface area, showing early branching with increased structural complexity. By day 7, PD_3D organoids expanded to the point of physical contact, thereby forming a dense, interconnected web-like structure. Even though the organoid area continued to increase (Figures 1H and 2A), the number of discrete organoids in this condition decreased over time, demonstrating that individual organoids fused as they grew larger. Hematoxylin and eosin (H&E) staining demonstrated that the largest structures in this condition were enriched for ATI-like cells—the thin, flattened epithelial cells responsible for gas exchange in the body (Figures 2K and 2N).Figure 2Day 10 organoid systems exhibit condition-specific morphologies, epithelial architectures, and cellular compositions(A–F, H–J) Brightfield images of representative organoid cultures on Day 10 for PD_3D (A–C), BAL_3D (D–F), and Mixed_3D (H–J). (A, D, and H) Show stitched whole-well images (scale bars, 1 mm); insets show zoomed organoids highlighting distinct morphologies across conditions (scale bars provided on individual images due to variation).(K–W) Hematoxylin and eosin (H&E) staining reveals condition-specific epithelial architectures, including variations in epithelial layering, lumen formation, and structural complexity. Scale bars are shown on each image. Arrows in K and N highlight thin, flattened epithelial cells morphologically consistent with alveolar type I (ATI)-like cells.(X–Z) UMAP projections of condition-level single-cell RNA-seq data for PD_3D (X), BAL_3D (Y), and Mixed_3D (Z), with cells colored by annotated identity (see legend in figure). Circled populations mark epithelial, immune, and mesenchymal subtypes of relevance, including Hillock cells, secretory-alveolar intermediates (AT0/RAS), polarized macrophages, Rspo3^+^ mesenchyme, and Pdgfrb^+^ pericytes. Embeddings shown using node-aligned global object cell type annotations mapped back onto the condition level embeddings. (aa–dd) Within-class breakdowns of cell type composition across conditions on day 10. Stacked bar plots show the relative proportions of specific cell types within each major cell classes (aa), epithelial subsets (bb), immune subsets (cc), and mesenchymal populations (dd), normalized to 1 within each class per condition. PD_3D and Mixed_3D cultures maintained epithelial, immune, and mesenchymal populations, while BAL_3D showed dominance by immune and epithelial populations, as is consistent with the input populations. (ee) Dot plot shows the expression of canonical epithelial lineage markers across all epithelial subtypes, split by condition. Each row represents a condition-specific epithelial cluster, and each column corresponds to a curated marker gene associated with the respective epithelial program. Dot size reflects the percentage of cells expressing each marker within a cluster, and color indicates scaled average expression. See also Figure S7.
In both growth kinetics and structural complexity, Mixed_3D organoids displayed an intermediate phenotype between that of PD_3D and BAL_3D structures. Organoid counts over time were greater than in BAL_3D, but less than in PD_3D (Figure 1G). Similarly, the average area of the structures increased over the culture period, though they remained smaller on average than BAL_3D by day 10 (Figure 1H). Structurally, the Mixed_3D cultures contained a blend of dense, compact structures analogous to spheroids as well as larger, branched structures with early complexity comparable to PD_3D, but without the extensive fusion of structures (Figures 2H–2J).
Organoids in the BAL_3D system self-assembled at a slower rate, and these cultures produced fewer structures compared to Mixed_3D and PD_3D (Figures 1G and 1H). However, the structures that did form were consistently larger (Figure 2H). A subset of the structures in this condition displayed a floret-like morphology with branching extensions indicative of lung-budding architecture (Figures 2E–2G).
To further test whether the observed BAL_3D morphologies could be explained simply by a lower seeding density of epithelial cells at day 0 (see Figure S1 for BAL starting population breakdown), we quantified epithelial fold-expansion (EFE), defined as the epithelial fraction at Day 10 relative to Day 0, for each condition (see method details: epithelial fold expansion (EFE) analysis). Despite starting with the epithelial fraction at seeding (4.8%), BAL_3D exhibited the highest EFE (mean: 48.0, 95% CI: 7.7–138.1), compared to 3.05 (95% CI: 2.7–3.7) in PD_3D and 4.23 (95% CI: 3.7–5.0) in Mixed_3D (Figures S9A and S9B). These results indicate that BAL_3D cultures most efficiently expand epithelial cells, even though their initial epithelial representation was the smallest.
To relate the observed morphologies to underlying cell states, we performed scRNA-seq on each system (condition) at Day 10 (see method details: organoid dissociation for scRNAseq and other downstream applications). Lineage outputs were observably different across PD_3D (n = 3, 8,724 cells) and Mixed_3D (n = 3, 6,344 cells) contained epithelial, immune, and mesenchymal cell populations, whereas BAL_3D contained only epithelium and immune cells (Figures 2X–2Z and 2aa–2cc). Endothelial populations were lost by Day 10 in both PD_3D and Mixed_3D (Figure 1F). Pdgfrb^+^ pericytes and Rspo3^+^ fibroblasts emerged in Mixed_3D and PD_3D, in varying proportions (Figures 2bb and 2dd).
BAL_3D (n = 3, 9,676 cells), which lacked mesenchymal cells, had the second largest percentage of ATI-like cells (5.16 ± 3.88%) after Mixed_3D (5.60 ± 1.70%), and more than PD_3D (2.68 ± 0.54%). This finding is inconsistent with previous reports demonstrating the reliance of alveolar epithelial differentiation on mesenchyme-derived signaling pathways, such as Wnt, BMP, FGF, TGFβ, and YAP/TAZ.^40^^,^^41^^,^^42^ Although our data do not rule out a supportive role for mesenchyme, they indicate other cues, such as immune-derived signals, may be more associated with the increased ATI-differentiation observed in BAL_3D. Specifically, the expression of insulin-like growth factor 1 (Igf1) and interleukin-1 beta (Il1b) was highest among the immune cells in BAL_3D (Figures S8C and S8D). These findings are in line with previous work demonstrating that Igf1 can stimulate ATII-to-ATI differentiation via the activation of Wnt5a signaling.^43^ They are also consistent with results that Il1b primes ATII cells for alveolar regeneration by inducing a damage-associated state necessary for their eventual differentiation into ATI cells.^6^
Epithelial composition and relative fractions of subtypes varied by input. PD_3D cultures generated a full spectrum of epithelial states, including ATII-like, ATI-like, basal, secretory, hillock, and RAS-like subtypes (Figure 2X). BAL_3D cultures, despite beginning with only 4.8 ± 7.3% epithelial cells, gave rise to a comparable range of epithelial types (Figures 1C, 2bb, and 2ee). Basal cells were most abundant in BAL_3D (9.3 ± 7.7%), compared to PD_3D (7.2 ± 5.7%) and Mixed_3D (5.9 ± 5%) (Figures 2X–2Z, 2bb, and 2ee). Transitional remodeling populations such as RAS-like cells and Sox9^+^ cells with stress-associated transcription profiles were identified in all three systems (Figures 2X–2Z, 2bb, and 2ee). BAL_3D, however, uniquely gave rise to a basal-polarized hillock-like state, whereas the mesenchyme-containing systems, Mixed_3D and PD_3D, produced a luminal-polarized hillock-like state (Figures 2bb, 2ee, and 2X–2Z). Taken together, these observations suggest that differences in starting composition do not create unique epithelial lineages de novo; rather, they bias the distribution and combinations of common transitional and differentiated cell states.
Global analysis of the epithelial compartment at day 10 (organoids) revealed variation across conditions in the feature-level characterization of a Krt13^+^ lineage, prompting us to examine this population in greater depth. Hillock cells were first described in the mouse and human proximal airway epithelium, specifically over the junctions between cartilage rings of the trachea and along the posterior membrane.^2^^,^^44^ Importantly, these cells have been reported to exist in healthy humans and mice, but following denuded injury, they expand rapidly and display increased turnover, contributing to the regeneration of the airway epithelium.^2^^,^^44^^,^^45^ To ground our data in native biology, we first tested whether or not these cells are present in the native rat trachea, a finding previously not reported. We performed immunohistochemical (IHC) staining for Krt13 on longitudinally sectioned native rat trachea, where hillock cells have not previously been described. We identified organized clusters of Krt13^+^ cells at both the basement membrane (basal) and atop the cells at the basement membrane (luminal) (Figures 3B and 3C). These structures mimicked the characteristic anatomy of hillock domains seen in mouse and human a layer of stratified squamous Krt13^+^ cells overlaying a basal stem cell pool expressing Krt5.^2^^,^^45^Figure 3Hillock cells emerge as a polarized epithelial state with distinct differentiation trajectories and immune-regulatory features(A–C) Immunofluorescence staining of native rat trachea shows clusters of Krt13^+^ squamous cells (cyan) with Tuba1a (green) and Foxj1 ciliated cells (red). These stratified structures resemble hillock domains described in mouse and human airways. Scale bars, 50 μm.(D) UMAP of the integrated Day 10 single-cell dataset, combining all organoid sequencing data. The right panel shows the subsetted proximal epithelial population, used for downstream trajectory inference and marker analysis. Hillock cells are subdivided into Hillock_Basal (Krt13^+^/Krt5^+^) and Hillock_Luminal (Krt13^+^/Krt5^-^) populations.(E) Heatmap of marker gene expression across hillock subtypes, including squamous markers (Krt13, Sprr1a), basal markers (Krt5, Krt14), and immune-associated genes highlighted in pink.(F–G) Slingshot pseudotime analysis of the proximal epithelial subset reveals a trajectory from cycling progenitors through Basal_Like, Hillock_Basal, and ultimately to Hillock_Luminal cells.(H) Gene expression trends across pseudotime show stepwise activation of differentiation and immune programs.(I) MSigDB Hallmark enrichment analysis compares Hillock_Luminal and Hillock_Basal populations.(J–O) Immunofluorescence staining of day 10 organoids shows condition-specific spatial organization of hillock cells. In Mixed_3D (J and M) and PD_3D (K and N), Krt13^+^/Krt14^-^ luminal cells form flattened, apical epithelial layers. In BAL_3D (L and O), Krt13^+^/Krt14^+^ basal-like cells assemble into multilayered squamous structures (white arrowheads). A subset of hillock luminal cells co-express Scgb1a1 (white arrowheads), supporting a luminal hillock identity. Scale bars, 100 μm.(P–Q) CellChat analysis chord plots show ligand-receptor interactions between hillock-like cells and other system populations. (P) Outgoing signals from hillock-like cells to proximal and distal epithelial, mesenchymal, and immune compartments, highlighting regenerative/developmental and inflammatory axes, including EGFR, Notch, Wnt/Fzd, chemokines, and Fgf. (Q) Incoming communication to hillock-like cells from epithelial, immune, and mesenchymal cell types, including regenerative, regulatory, structural, and inflammatory signaling axes such as Areg—Egfr, Sema4a—Plxnb2, Fn1—Sdc1, and Spp1—Cd44, respectively.
Through the analysis of scRNAseq data, we found that all hillock-like cells in 3D culture simultaneously expressed Sprr1a, a member of the small proline-rich (Sprr) protein family, alongside Krt13. We highlight this finding as Sprr1a has been implicated in terminal squamous differentiation.^46^^,^^47^^,^^48^ Therefore, the co-expression of Krt13 and Sprr1a by these cells supports a model in which they are engaging in a mature squamous differentiation program.
Further analysis of scRNAseq data allowed us to elucidate the condition-specific polarization event that presented in our data. Sub-clustering analysis revealed that these Krt13^+^ cells clustered distinctly within the proximal epithelial subset (Figure 3D) and could be further subdivided into two transcriptionally polarized hillock basal (Krt13^+^/Krt5^+^/Sprr1a^+^) and hillock luminal (Krt13^+^/Sprr1a^+^/Krt5^-^), in line with previously characterized hillock subtypes in the murine and human lung.^2^^,^^44^^,^^45^ Interestingly, we found that the relative proportions of each subtype were condition-specifically BAL_3D cultures derived from an immune-rich, mesenchyme-poor starting population exhibited an abundance of hillock basal cells, whereas Mixed_3D and PD_3D conditions, derived from a more compositionally diverse starting population, favored the hillock luminal phenotype (Figure 3D). These findings indicate that the immune-rich context of BAL_3D favors a basal-skewed hillock state, whereas epithelial-mesenchymal heterogeneity promotes luminal differentiation (Figure 3E). Given that the initial epithelial compartment in BAL largely mirrors that of PD, albeit at a much smaller scale (Figures S1B, S1D, S2A, and S2D), these results suggest that microenvironmental cues, rather than epithelial-intrinsic differences, are major drivers of hillock subtype polarization.
To understand how these cells differentiate, we used Slingshot pseudotime analysis. From this analysis, we found these cells progress along a well-defined trajectory from proximal cycling epithelial cells through canonical basal cells (Krt5^+^/Sprr1a^+^/Krt13^-^) to basal-hillock-like (Krt5^+^/Sprr1a^+^/Krt13^+^), ultimately terminating in the luminal hillock state (Krt5^-^/Sprr1a^+^/Krt13^+^) (Figures 3F and 3G). Earlier findings confirm this trajectory as they suggest that the more squamous, hillock luminal cells may be derived from hillock basal cells and thus, the luminal cells represent a more differentiated or specialized form of the hillock lineage.^2^
Given evidence that epithelial cells can adopt immunological roles in repair contexts, examining the immune-associated genes expressed by the hillock-like cells provides insight into their potential contributions in our systems and disease contexts.^1^^,^^2^^,^^49^ scRNAseq data revealed that luminal hillock cells upregulate a distinct suite of immune-associated genes, including Slpi, Mal, S100a7, Il1a, and S100a8 (Figure 3E). Slpi, enriched in the Mixed_3D and PD_3D systems, encodes a serine-protease inhibitor that protects airway epithelium from neutrophil-derived proteases while also exerting anti-inflammatory and antimicrobial properties.^50^^,^^51^^,^^52^ This defensive profile was complemented by Il1a, a pro-inflammatory cytokine that functions as an alarmin and rapidly recruits immune cells in response to stress, necrosis, or damage.^53^ Additionally, Mal, a proteolipid involved in mediating apical sorting among polarized epithelial cells, ensures the proper delivery of proteins to the apical membrane, as well as facilitating trafficking in T cells for TCR-immune activation.^54^^,^^55^ The upregulation of calcium-binding proteins S100a7 (psoriasin) and S100a8 underscores the inflammatory potential of these cells*.* Whereas S100a7 is reported to be constitutively expressed by the bronchial epithelium and is heightened during bacterial infection, S100a8 induces Muc5ac production by airway epithelial cells in COPD and activates alveolar epithelial cells via TLR-4, triggering chemokine and cytokine release.^56^^,^^57^ Collectively, the expression of these genes suggests that luminal hillock cells may act not only as transitional intermediates, but as immunomodulatory effectors within the remodeling epithelium.^21^
We further explored the functional specialization within the hillock population by comparing luminal hillock cells to basal hillock cells using MSigDB Hallmark pathway enrichment analysis.^58^ Luminal hillock cells were significantly enriched for pathways related to inflammatory signaling and immune communication (TNFα signaling via NF-κB, IL2/STAT5 signaling), as well as stress adaptation (complement, TGFβ signaling, and hypoxia) (Figure 3I). In contrast, basal hillock cells were enriched for processes related to proliferation (MYC targets, E2F targets, and G2M checkpoint), consistent with a more canonical basal epithelial state (Figure 3I).^59^
To validate these transcriptional and pathway-based findings, we also more deeply studied the spatial distribution of hillocks using immunofluorescence staining of day 10 organoids across BAL_3D, PD_3D, and Mixed_3D conditions. We found that in BAL_3D, Krt13^+^ cells were of a multilayered squamous appearance and were localized throughout the apical face (Figure 3I). Alternatively, Krt13^+^ cells in Mixed_3D and PD_3D exhibited a flattened, apically oriented morphology with less extensive coverage of the epithelial surface (Figures 3J, 3K and 3N). Of note, in both PD_3D and Mixed_3D organoids, some luminal cells also co-expressed Scgb1a1, consistent with their transcriptional signature and suggesting that Krt13^+^/Scgb1a1^+^ cells represent a differentiated luminal hillock state (Figures 3J and 3N).
Previous studies have suggested that progenitor epithelial states, such as hillock cells, can act as both sensors and effectors in multicellular systems to facilitate repair.^60^ We therefore sought to connectomically determine whether hillock-like cells function as signaling hubs within our three organoid systems. To do this, we applied CellChat to model the probability of ligand-receptor communication between hillock cells and all other cells, as well as using the package’s network centrality to determine which cell types assumed roles of sender versus receiver.^61^^,^^62^ Chord plots of statistically enriched ligand-receptor interactions (p < 0.05) revealed hillock cells serving dual roles, broadcasting regenerative and inflammatory cues (Egfr, Fgf, Notch, Wnt/Fzd, and chemokines), while simultaneously integrating adhesive, growth factor, and immune-regulatory inputs (Figures 3P and 3Q). Together, these results provide another line of evidence that hillock cells serve as dynamic hubs of multicellular communication that coordinate epithelial, immune, and mesenchymal responses during remodeling.
To better understand how the immune compartment changed over the course of culture, we next examined the persistence and spatial dynamics of immune cells across conditions. While most immune populations were completely absent in our organoid cultures by day 10 of culture, a population of activated macrophages represented an exception. These macrophages were maintained in all three engineered conditions, with the most consistent maintenance in the BAL_3D and Mixed_3D cultures, where the initial (day 0) ratio of immune cells to all other cells was greatest (Figures 2aa, S1A, and S2A). Spatially, the macrophages remained dispersed in the extracellular matrix rather than embedded within organoid structures or contained within luminal spaces (Figure 4E). They were spatially segregated from epithelial cells, indicating that their effects likely were facilitated via paracrine signaling rather than cell-cell contact (Figure 4E).Figure 4Macrophage polarization trajectories reveal condition-specific immune modulation(A) UMAP of the integrated day 10 organoid dataset highlighting macrophage populations, which cluster as a transcriptionally distinct immune compartment within the broader dataset.(B) Unsupervised Seurat graph-based clustering (Louvain; see method details) of the macrophage subset identifies three clusters due to the presence of cycling genes.(C) UMAP of macrophage subset, annotated into two transcriptionally distinct anti-inflammatory-like (salmon) and pro-inflammatory-like (lavender).(D) Cell cycle phase scoring of macrophages reveals that while most cells are in G1, a subset of pro-inflammatory-like macrophages are cycling, with enrichment in S and G2/M phases.(E) Stacked bar plot shows condition-specific contributions to the macrophage population. BAL_3D cultures are enriched for anti-inflammatory-like macrophages, whereas PD_3D and Mixed_3D conditions contain a higher proportion of pro-inflammatory-like macrophages.(F) Immunofluorescence staining of day 10 BAL_3D cultures shows Ptprc^+^ alveolar macrophages (cyan) dispersed in the surrounding matrix, often adjacent to but spatially distinct from epithelial organoids (outlined by dashed lines).(G) Pseudotime trajectory inferred using Slingshot reveals a transcriptional continuum from anti-inflammatory-like to pro-inflammatory-like macrophage states.(H) Volcano plot shows differentially expressed genes between pro- and anti-inflammatory-like macrophages. Pro-inflammatory-like macrophages express elevated levels of inflammatory and alarmin genes (e.g., Cxcl2, Trem1, Il1a, and S100a8), while anti-inflammatory-like macrophages upregulate genes associated with tissue remodeling and resolution (e.g., Pparg, Ms4a4a, and Cela1).(I) Changes in average gene expression (log-normalized) between alveolar macrophages in the starting populations and polarized macrophages in the three organoid systems.(J) Inferred mesenchyme-immune signaling axes at day 10. Directed edges denote ligand-receptor interactions, line width scales with inferred connectivity, and node size reflects the fraction of each cell type within the subset system. Statistically significant (p < 0.05) interactions and biologically relevant interactions are Rspo3—Sdc4, Cxcl12—Itgb1, Cxcl12—Sdc4, and Vcam1—Itgb2.
To quantify whether these separations fall within a physiologically plausible signaling radius, we measured the linear distance between each Ptprc^+^ macrophage and the nearest organoid boundary in the immunofluorescence images displayed in Figure 4F. Using ImageJ, we extracted X-Y coordinates for organoid epithelial edges and individual macrophages and used these to compute the shortest Euclidean distance for each macrophage-organoid edge pair (Figure S12). These measurements showed that macrophages generally resided ∼100–400 μm from the epithelial surface (Figure S12), placing them within a range that is compatible with paracrine communication rather than direct contact-dependent (juxtracrine) signaling.
Contextualizing these distances, we considered representative ligands produced by the macrophage population. Cxcl1, upregulated in the pro-inflammatory subset, has a reported diffusion coefficient of ∼39 μm^2^/min^63^ Using standard diffusion distance approximations, as put forth by Fick’s second law of diffusion (L≈Dt), such a coefficient supports diffusion over tens to hundreds of microns which is well within the macrophage-organoid separation range we observed (Figure S12).^64^ As an additional comparison, Il1b exhibits a substantially higher reported diffusion coefficient (∼4.44 × 10^4^ μm^2^/min) implying an even larger characteristic diffusion length.^65^ Larger ligands such as Serpina1, however, are less likely to have a significant effect on the epithelial-immune signaling circuit given their restricted diffusion capabilities. See Tables S8 and S9 for more information on upregulated ligands and receptors, specific to the polarized macrophage populations.
To further characterize macrophage heterogeneity, we subclustered the immune compartment and identified two transcriptionally distinct macrophage clusters (Figures 4B and 4C). To investigate whether the macrophages that persisted across conditions adopted distinct activation states, we performed differential gene expression analysis on this subset, comparing cluster 0 (enriched in BAL_3D) to clusters 1 and 2 (enriched in Mixed_3D and PD_3D). The resulting transcriptional programs did not fully align with classical M1/M2 macrophage paradigms but instead resembled anti-inflammatory-like and pro-inflammatory-like states (Figures 4B and 4G).^66^
Anti-inflammatory-like macrophages, which predominated in BAL_3D, expressed a coordinated program of immune regulatory and tissue-adaptive genes, including Pparg, Cd84, Cadm1, and Cela1 (Figure 4G). Pparg, a nuclear receptor known to repress inflammatory gene expression, likely acts upstream to promote alternative activation and lipid handling.^67^^,^^68^ Cd84, an SLAM family receptor, may reinforce this anti-inflammatory response by modulating NF-κB signaling during inflammation.^69^ Cadm1 and Cela1 further point to tissue Cadm1 has been implicated in macrophage-epithelial interactions, while Cela1, although canonically epithelial, contributes to extracellular matrix remodeling and has been detected in macrophages during lung development.^70^^,^^71^ Taken together, these data suggest that macrophages in BAL_3D not only play immunoregulatory roles but may be involved in restructuring the physical properties of the engineered niche.
Pro-inflammatory macrophages, which were enriched in Mixed_3D and, to a lesser extent, PD_3D, upregulated genes associated with acute inflammatory signaling, leukocyte recruitment, and matrix remodeling (Figure 4G). These cells showed a marked increase in the expression of several potent neutrophil chemoattractants such as Cxcl1, Cxcl2, and Cxcl3 (Figure 4H). These chemokines are known to signal through Cxcr1 and Cxcr2. Cxcr1 was expressed by a subset within the activated macrophage population, suggesting that the pro-inflammatory macrophage upregulating Cxcl1/2/3 were operating via an autocrine loop to elicit their pro-inflammatory effects.
Concurrent expression of Il1a and Nos2 by the pro-inflammatory macrophage subset reflects classical M1-like Il1a being a master cytokine for amplifying local inflammation, and Nos2 catalyzing nitric oxide production, contributing to microbicidal activity and tissue damage (Figure 4G).^72^^,^^73^ Supporting our findings that the macrophages in PD_3D and Mixed_3D assumed a pro-inflammatory signature, a subset of these macrophages were inferred to be in the G2M and S phases, according to cell cycle scoring performed in Seurat (Figure 4C). This observation may indicate that these cells were undergoing continuous cell cycle progression, a feature that has previously been associated with the activation-induced expansion of inflammatory macrophage populations.^74^ This pro-inflammatory program was more common in the systems that maintained mesenchymal cells (Figures 2aa and 4D), suggesting that mesenchymal-derived cues could drive macrophage inflammatory polarization, even under otherwise similar culture conditions.
We performed pseudotime analysis on the immune subset to further characterize the relationship between these states and their system-level specificity (Figure 4F). This analysis revealed a continuous trajectory from the anti-inflammatory cluster to the pro-inflammatory cluster. Again, given our findings that the anti-inflammatory macrophages were unique to the BAL_3D system that lacked mesenchyme, pseudotime analysis reinforces our interpretation that mesenchyme-derived systems drive differentiation (Figure 4D). The continuous trajectory in pseudotime, combined with the system-specific distribution of cell states, suggests that these populations may not represent fixed or terminal identities but instead exist along an activation spectrum and that mesenchymal signals are critical in driving such activation.
Beyond the polarization in macrophage state across conditions, it is notable that while all three systems contained some proportion of polarized macrophage by the end of culture, the Mixed_3D condition maintained the greatest amount. These findings prompted us to consider how the presence of mesenchymal cells may play a critical role in the maintenance of immune cells. To test this, we created a subset object containing only the immune and mesenchymal populations present in Mixed_3D and PD_3D and ran cell-to-cell connectomic analysis on the subset using NICHES.^75^ From this analysis, we identified several statistically enriched ligand-receptor signaling axes that were upregulated between mesenchymal and immune compartments. These signaling axes included Rspo3—Sdc4, Cxcl12—Itgb1, Cxcl12—Sdc4, and Vcam1—Itgb2, which together support a Wnt-modulatory mechanism and chemotactic support of the immune cells from the mesenchyme (Figure 4J). Of particular interest, the Rspo3—Sdc4 mechanism has been reported to potentiate Wnt signaling through the linkage of R-spondin ligands to syndecans, prompting downstream activation of planar cell polarity.^76^^,^^77^ Additionally, the production of the chemokine Cxcl12 by mesenchymal cells, which was predicted to interact with Itgb1 and Sdc4 on immune cells, suggests that mesenchymal cells are supporting the survival of immune cells and calibrating their inflammatory response.^20^^,^^78^^,^^79^
The transcription factor Sox9 is best known as a marker of multipotent distal epithelial tip progenitors during branching morphogenesis in the developing lung, but recent studies have also identified Sox9^+^ transitional progenitors in adult models.^80^^,^^81^ These specific progenitors have been found to become activated in response to injury and subsequently participate in distal lung regeneration. To test whether a similar population could be identified in the rat lung and in our engineered systems, we mapped the distribution and transcriptional features of Sox9^+^ cells under homeostatic and culture conditions. We identified Sox9^+^ cells in the native adult rat lung near the bronchoalveolar duct junction (BADJ), as well as in organoid cultures (Figures 5A and 5H–5J). In the adult rat lung, these cells occasionally co-expressed Scgb1a1, consistent with their role as a bipotent progenitor population (Figure 5A). We also observed Scgb3a2^+^ secretory cells concentrated in the conducting airways, mirroring the spatial distribution reported in human and ferret respiratory bronchioles (Figure 5B).^4^ In native adult rat lungs not exposed to injury cues, Scgb3a2^+^ cells did not co-express Sftpc (Figure 5B). Rather, Sftpc expression was confined to canonical ATII cells within the alveolar region (Figure 5B).Figure 5Mesenchyme-rich conditions promote the emergence of secretory-alveolar intermediates, while Sox9^+^ stressed progenitors emerge across all engineered systems(A and B) Immunofluorescence staining of native adult rat lung reveals rare epithelial progenitor and secretory populations. (A) Sox9^+^ cells (cyan) are observed near the bronchoalveolar duct junction (BADJ), occasionally co-expressing Scgb1a1 (red). Pdpn (green) marks alveolar type I (ATI) cells. (B) Scgb3a2^+^ secretory cells are localized to the conducting airways. Arrowheads highlight individual positive cells.(C) UMAP of the integrated day 10 organoid dataset highlights transcriptionally distinct transitional populations, including aRAS/AT0-like, stressed progenitor (Sox9^+^), and canonical secretory subsets.(D) Boxplot quantifying the percentage of Sox9^+^ cells (threshold: Sox9 expression >0.1) across five conditions (PD, BAL, PD_3D, Mixed_3D, BAL_3D). Sox9^+^ cells are rare in the starting populations but expand markedly after 10 days of culture in all conditions. Each point represents an individual biological replicate. Statistical comparisons were performed using two-proportion Z-tests on raw Sox9^+^ counts, corrected using the Benjamini-Hochberg method (∗∗∗p < 0.001).(E) Scaled heatmap shows differential gene expression across transitional epithelial subtypes. aRAS/AT0-like cells (Scgb3a2^+^/Scgb1a1^+^/Sftpc^+^) are distinct from stressed progenitors, which express stress-associated markers (Cox4i2, Ptges, Stc1) alongside the progenitor marker Sox9.(F) Violin plots show the expression of signature markers distinguishing canonical secretory and activated RAS-like epithelial populations, split by condition.(G) Boxplot shows the percentage of epithelial cells per replicate positive for a RAS-like signature (Scgb1a1^+^/Scgb3a2^+^/Sftpc^+^) across engineered conditions. Signature positivity was determined using scaled expression Scgb1a1 > 2, Scgb3a2 > −2.5, and Sftpc > −1.5. Triple-positive cells were observed across all three conditions. Statistical comparison via a Kruskal-Wallis test revealed no significant differences in overall prominence; however, the PD_3D system exhibited greater replicate variability, as demonstrated by the wider interquartile range.(H–J) Immunofluorescence staining of day 10 organoids reveals Sox9^+^ epithelial cells (cyan) variably distributed across all conditions. Scale bars, 100 μm.(K and L) Immunofluorescence staining of day 10 organoids (PD_3D and Mixed_3D) shows co-expression of Scgb3a2 (green) and Sftpc (magenta) (RAS-like cells). Scale bars, 100 μm (K) and 50 μm (L).
Drawing on the findings that Sox9^+^ progenitors proliferate in the adult mouse lungs in response to chronic injury cues (e.g., Il4 and bleomycin-induced lung fibrosis), we sought to determine if similar secretory alveolar intermediates were present in our organoid cultures.^81^ Leveraging scRNAseq data, we identified an Scgb1a1^+^/Sftpc^+^/Sox9^+^ expressing progenitor-like population in all three organoid systems (Figures 5D and 5H–5J). While these cells were present in our PD and BAL starting populations, they were not nearly as abundant (relative fraction of sample) (Figures 1B and 1C), nor did they express markers associated with cell stress.
To further explore the transcriptomic shift in the Sox9^+^ cells present in the starting populations versus those that expanded in organoid culture, we ran differential gene expression analysis between not only the “Stressed_Progenitors” and Sox9^+^ cells in the PD and BAL starting populations, but also between the three engineered conditions themselves. This analysis revealed “Stressed_Progenitors” to be enriched for the expression of stress-associated genes, including Cox4i2, Ptges, and Stc1, in addition to the already present expression of BASC signature genes (Figure 5E). The co-expression of BASC signature markers and stress-associated genes suggests that these Sox9^+^ transitional cells may be an injury-responsive intermediate state reflecting the intrinsic plasticity and extrinsic cues within our controlled culture conditions (Figures 5C and E). Cox4i2, a cytochrome c oxidase subunit, is normally expressed in healthy lung epithelium, but is upregulated in response to hypoxia, supporting the idea that these progenitors underwent metabolic adaptation to sustain survival and expansion in organoid conditions.^82^^,^^83^ Upregulation of Ptges and Stc1 further highlights the reparative potential of these cells, whereby Ptges contributes to PGE2-mediated inflammatory modulation, while Stc1 promotes anti-apoptotic signaling during injury.^84^^,^^85^
Alongside the enrichment for Sox9^+^ stressed progenitors in 3D culture, we also identified an activated RAS-like epithelial state, as defined by the co-expression of Scgb1a1, Scgb3a2, and Sftpc, across all three organoid systems. Of note, the fraction of the epithelium that they represented varied significantly across systems. These cells were most abundant in the PD_3D and Mixed_3D conditions where mesenchyme was preserved (as demonstrated by both RAS signature scoring and expression of individual genes, Figures 5F, 5G, S7D, S7H, and S7I). In addition to the RAS-like transitional state, we also identified canonical secretory cells (Scgb1a1^+^/Scgb3a2^+^/Sftpc^-^) across all three systems, but unlike the RAS-like cells, these cells did not express alveolar-associated genes (e.g., Sftpc, Napsa, and Lamp3). Organoid-inherent cues aside, the RAS-like intermediate population was not starting population-specific, nor did it require exogenous injury to emerge, suggesting that epithelial-mesenchymal-immune interactions affect the balance between transitional and differentiated epithelial states under homeostatic culture conditions.
Understanding how diverse cellular inputs shape epithelial and immune behavior requires not only classifying cell states but also mapping the signaling cues that guide them. Having established condition-dependent emergence of specific cellular polarities, the mechanisms coordinating these shifts remained difficult to elucidate at the individual-lineage level. With an understanding that mesenchymal cells facilitate paracrine signaling, producing signals that affect both epithelial differentiation and immune polarization, we hypothesized that variation in the presence of mesenchyme modulates these outcomes by regulating context-specific signaling programs.^86^^,^^87^^,^^88^ To test this, we ran cell-to-cell connectomic analysis using NICHES on a subset object including only cell populations of “Hillock_Like,” “Polarized_Mac,” “Rspo3+_Mes,” and “Pdgfrb+_Pericytes” (Figure 6B). We then performed differential expression analysis on the subset connectomic object (log2 fold-change >0.25, expressed in ≥25% of cells) and calculated a power metric (expression ratio ∗ avg. log2FC) to identify upregulated ligand-receptor pairs that were both differentially expressed and population-specific. Through this analysis, we were able to resolve how communication within and across epithelial, mesenchymal, and immune compartments changes under different engineered conditions (Figure 6A). We identified signaling axes upregulated in mesenchyme-rich (Mixed_3D, PD_3D) versus mesenchyme-poor (BAL_3D) contexts, including many involving pro-regenerative, pro-inflammatory, and immunomodulatory ligands (Figure 6C) (see Tables S10 and S11 for a comprehensive list of ligands and receptors upregulated by the Hillock-like population).Figure 6Hillock cells participate in reciprocal signaling circuits that coordinate mesenchymal and immune inputs(A) Schematic summary of ligand–receptor interactions between cell populations of interest, inferred using NICHES. Arrows represent directional communication among Polarized_Mac, Pdgfrb^+^_Pericyte, Rspo3^+^_Mes, and Hillock_Like epithelial cells. Hillock cells act as central orchestrators in these networks, orchestrating and responding to diverse mesenchymal and immune-derived signals.(B) UMAPs showing sending (left) and receiving (right) lineages used for analysis of cell-cell connectivity (NICHES output).(C) Violin plots shows the expression of candidate ligands and receptors enriched in hillock cells.(D) Scaled heatmap of ligand and receptor expression across cell types and conditions, highlighting transcriptional diversity of connectomic inputs.(E–H) Circuit plots illustrate inferred cell-cell signaling across selected ligand-receptor axes. Each plot shows directional connectivity between sending and receiving populations based on NICHES scoring. Line thickness represents scaled communication strength; node size reflects the relative abundance of each cell population. (E) Circuits for Serpina1—Lrp1 and Cd24—Selp axes. (F) Circuits for Il1a—Il1r2, Il1rn—Il1r1, and Il1rn—Il1r2, highlighting inflammatory and regulatory signaling pathways. (G) Circuits for developmental remodeling signals, including Mdk—Tspan1, Efnb2—Ephb3, and Ereg—Erbb3. (H) Circuits for matrix-associated and EGFR signaling, including Fn1—Itgav and Tgfa—Egfr.
This analysis validated our hypothesis, demonstrating a clear signaling axis between hillock-like cells and activated macrophages. Broadly, this network analysis revealed that all lineages included were communicating via both autocrine and paracrine signaling (Figure 6B). The ligands Serpina1, Cd24, Ill1rn, and Il1a were expressed by both hillock and macrophage populations and targeted receptors on mesenchymal and immune cells (Figures 6C and 6D). For example, the upregulation of the Serpina1—Lrp1 and Cd24—Selp signaling axes point to the possible role hillock cells play in regulating macrophage adhesion and activation (Figure 6E).^89^^,^^90^^,^^91^ By comparison, Il1rn—Il1r1, Il1rn—Il1r2, and Il1a—Il1r2 signaling converged on the Rspo3^+^ mesenchyme node, indicating that epithelial and immune sources may work in tandem to influence mesenchymal behavior (Figure 6F). The frequency of the IL-1 family interactions, as well as their specificity, made them stand out as a potential core organizing principle within the hillock-mesenchyme-macrophage circuit established. The juxtaposition of Il1a, a pro-inflammatory cytokine, with Il1rn, a natural antagonist that dampens IL-1 signaling, suggests not only the delivery of signal, but also a mechanism to balance activation and resolution, akin to a built-in regulatory logic gate.^91^^,^^92^^,^^93^ These findings suggest that epithelial and myeloid cells do not directly activate mesenchyme engagement in a unidirectional, “fibroinflammatory” response, but instead initiate and modulate mesenchymal engagement to ensure the amplitude and duration of signaling are appropriate for maintaining tissue structure and avoiding excessive remodeling.^93^
We also identified reverse signaling from mesenchymal populations to hillock cells, including Mdk—Tspan1, Efnb2—Ephb3, and Ereg—Erbb3 (Figure 6G). These axes may imply that mesenchymal cells contribute actively to epithelial regulation, delivering cues that support survival, polarity, and lineage plasticity. Midkine (Mdk) signaling through Tspan1 may help maintain epithelial cells in a reparative and plastic state capable of responding to environmental cues when necessary.^94^ Tspan1 has been identified as a regulator of epithelial-to-mesenchymal transition (EMT) in alveolar epithelial cells in the context of IPF, whereby its expression downregulates Smad2/3 and beta-catenin signaling.^95^ This activity helps preserve epithelial identity and prevents the transition toward a mesenchymal state.^95^
Efnb2—Ephb3 interactions may contribute to the spatial containment of the hillock-like epithelium. Here, Efnb2, secreted by Rspo3^+^ mesenchyme, acts on the Ephb3 receptor of the hillock-like population. The upregulation of this signaling axis between these two cell populations likely reinforces compartment boundaries, thereby preventing hillock cell expansion into surrounding alveolar or transitional territories.^96^^,^^97^ This spatial guidance allows hillock-like cells to maintain their identity by preserving their stacked-squamous-like architecture and required spatial insulation. Ereg—Erbb3, meanwhile, may promote hillock survival and stress adaptation, consistent with the findings that Erbb3 protects against injury-induced cell death and facilitates epithelial repair via PI3K/AKT activation.^98^^,^^99^^,^^100^
Simultaneously, we observed broader multi-nodal pathways such as the Tgfa—Egfr axis, which exhibited robust signaling from hillock cells to multiple compartments, including Rspo3^+^ mesenchyme, Pdgfrb^+^ pericytes, and themselves (Figure 6H). Tgfa, a potent Egfr ligand involved in epithelial proliferation and morphogenesis under stress, was broadly expressed in hillock cells, while Egfr was distributed across hillock-like and mesenchymal populations.^101^^,^^102^ This specific network implies that Tgfa-producing hillock cells may serve as central broadcasting effector cells to coordinate multicellular responses (sensor cells) via autocrine and paracrine loops. While Tgfa—Egfr signaling reinforces epithelial identity, it may also modulate mesenchymal responsiveness or act as a feedback hub.^102^^,^^103^ Other cross-compartmental axes, such as Fn1—Itgav, flowed predominantly from mesenchyme to epithelium, indicating the support of matrix remodeling and epithelial repair (Figure 6H).^103^ Collectively, these patterns revealed a multi-lineage signaling circuit that includes epithelial, immune, and mesenchymal populations. In this circuit, hillock cells act as central nodes capable of broadcasting signals to modulate immune and mesenchyme-driven processes while also responding to instructive cues from themselves.
Here, we report on a lung organoid model defined by compositional inputs that enables multicellular self-organization, lineage-specific fate emergence, and microenvironmental signaling analysis. By varying the relative contribution of immune- and mesenchyme-containing inputs under uniform matrix, media, and environmental conditions, we established a platform to investigate how initial cellular composition shapes epithelial and immune dynamics. The ability to maintain cell-type diversity and track context-dependent behavior under standardized conditions enables a degree of experimental control rarely achievable in vivo and allows for the interrogation of regenerative signaling at both cellular and network levels. Moreover, because epithelial, mesenchymal, endothelial, and immune lineages can be isolated from mouse and human lungs using established methods,^104^^,^^105^^,^^106^^,^^107^^,^^108^^,^^109^^,^^110^^,^^111^ the framework described here is directly translatable to murine and human systems.
While recent work has used fluorescence-activated cell sorting (FACS) to expand BAL-derived progenitors,^112^ our BAL_3D organoid model exhibits substantial epithelial expansion from rat BAL samples that had initially comprised only a minor population of epithelium and no mesenchyme. These cultures, which were predominantly immune in composition and devoid of any deliberate injury cues or progenitor enrichment, reliably produced a wide range of epithelial cell types, including the transitional and secretory subtypes of increasing interest to the pulmonary biology community. Additionally, this system, as well as the PD_3D and Mixed_3D systems, allowed for the expansion of both proximal and distal populations without requiring specialized “proximal” or “distal” growth media. The ability to generate complex epithelial structures from a cell suspension with trace epithelial input, under defined culture conditions, poises our BAL_3D system as having significant translational potential for regenerative medicine and cell therapy development.
Transitional epithelial states (i.e., Sox9^+^ stressed progenitors, hillock cells, and secretory-alveolar intermediates) emerged without deliberate injury cues such as exposure to noxious agents or mechanical agitation, suggesting that matrix and media alone promote epithelial plasticity. Nevertheless, we do recognize that the culture medium utilized was originally designed to expand embryonic lung progenitors.^39^ Despite this, the relative abundance and character of transitional populations in our systems varied between conditions, indicating that local community composition, especially regarding immune and mesenchymal cells, can influence epithelial remodeling. Among the transitional states observed, hillock-like cells expanded with striking consistency, demonstrating the utility of this system to expand these cells in vitro, thus affording researchers the ability probe their molecular identity, functional role, and potential contributions to tissue homeostasis in greater depth. Hillock-like cells were most prominent in conditions retaining mesenchymal populations, particularly Mixed_3D and PD_3D. These hillock-like cells exhibited stratified organization, context-specific polarization, and expression of immune-modulatory genes, features consistent with barrier-associated remodeling.^113^^,^^114^^,^^115^ To our knowledge, this is the first study to begin delineating the cell-cell signaling architecture associated with hillock cells. We identified signaling axes linking hillock cells to mesenchymal and immune partners, including Il1a—Il1r2 and Serpina1—Lrp1, that position these cells as both sensors and broadcasters of inflammatory and regenerative cues. Rather than passive responders or indicators of dysfunction, hillock cells may serve as active modulators of tissue recovery, helping to calibrate local inflammation and stabilize epithelial integrity.^2^
As with the more differentiated luminal hillock-like cells generated by our models, the RAS-like intermediates expanded most robustly in systems containing mesenchyme (e.g., Mixed_3D and PD_3D). Therefore, we propose that the presence of mesenchyme, and subsequent presence of mesenchymal-derived ligands such as WNT and BMP, critical for epithelial patterning and resolution, are supportive of transitional epithelial states, providing them with stability to exist at the “ridges” of Waddington’s landscape.^16^^,^^116^^,^^117^^,^^118^^,^^119^ This observation is consistent with the growing body of evidence that epithelial fate is not strictly lineage-encoded but rather dynamically tuned by environmental signals.^120^^,^^121^ Additionally, it lends further credence to the hypothesis that the maladaptive stabilization of transitional states underlies disease pathogenesis.^10^^,^^48^ Our model provides a tractable system to study the regulation and persistence of these alveolar-secretory intermediates, just as it does for the study of hillock-like cells.
Like the polarization displayed by several epithelial intermediates in our systems and their implications regarding plasticity in response to local cues, the macrophages that were sustained in each 3D system can be thought of in a similar way. The pro-versus anti-inflammatory polarization exhibited by the activated macrophages aligns with frameworks that poise macrophage plasticity as environmentally regulated rather than preprogrammed, especially in the context of chronic lung diseases.^122^^,^^123^ Macrophages that were maintained in our 3D culture platform in the presence of mesenchymal cues (Mixed_3D and PD_3D) were biased toward a pro-inflammatory state and thus, upregulated canonical inflammatory mediators such as Il1a, Il1b, Cxcl1/2/3, and Nos2, many of which are well within a physiologically plausible paracrine range. These findings are in line with previous reports of mesenchyme guiding a pro-inflammatory state in development and injury processes.^92^^,^^103^^,^^124^ When mesenchyme was not present to provide secreted ligands such as Il6, Ccl2 and Cxcl12 that have been demonstrated to activate macrophage cells to a pro-inflammatory state,^125^^,^^126^^,^^127^ macrophages adopted an anti-inflammatory phenotype, amplifying programs associated with immune resolution and epithelial support via the expression of genes such as Pparg, Cela1, and Lpl.^70^^,^^128^^,^^129^ While our findings are unable to confidently and accurately define the exact molecular mechanisms driving the shifts observed in macrophage polarization in chronic lung disease, they do suggest that the crosstalk driving macrophage polarization is not mediated by a singular axis between mesenchyme and immune or epithelium and immune. Rather, the three lineages (epithelium, immune, and mesenchyme) are in dynamic conversation with one another to define regenerative versus maladaptive trajectories.
Across the various cell types described and our study of their intra-lineage communication, we draw on the framework set forth by Medzhitov^130^ whereby inflammation is not considered a byproduct of damage but a regulatory process that maintains tissue function through adaptive coordination. From this perspective, the emergence of Sox9^+^ stressed progenitors, hillock-like cells, and aRAS-like intermediates can be thought of as an epithelial response to environmental cues to ensure regenerative preparedness. These transitional states may serve as context-sensitive intermediates ready to modulate proliferation, stress response, and signaling plasticity.^131^ The reproducibility of these intermediates, which are of increasing interest to the lung biology community, speaks to the value of our platform in studying successful regeneration as well as its derailment in chronic disease.
First, while we performed single-cell RNA sequencing on the initial PD and BAL isolates that were seeded in various ratios for our three systems, we did not prospectively enrich, deplete, or barcode individual lineages before the start of culture. Therefore, while we were able to characterize the lineages and their respective cell types present at the start of culture in each system, we cannot reliably determine which exact cells differentiated into the resulting output populations in culture. That being said, our approach did allow us to identify which cell types, more broadly, persisted and which ones emerged spontaneously in 3D culture. If a particular cell type with a specific transcriptional program was not present in either of our starting samples (PD or BAL), but was present in the sequencing data generated from our engineered samples (PD_3D, BAL_3D, and Mixed_3D), it is guaranteed that these transitional cell types emerged as a product of the system in which their origin cell was grown in, highlighting the plasticity and context-dependent emergence of these intermediates. Future experiments employing prospective lineage tracing, genetic barcoding, or FACS-based enrichment protocols would be useful to resolve the origin and fate of particular cell types within the proposed regenerative circuit.
Second, our analysis was conducted at a single, terminal time point (day 10). Without sampling over time, it remains unclear whether the epithelial populations we observed, such as the stressed Sox9^+^ progenitors, hillock cells, and secretory-alveolar intermediates, are stable endpoints or transient endpoints along the path to terminal differentiation. Defining such snapshots in time for these transitional cells warrants time-course transcriptomic profiling in combination with the already utilized trajectory inference and fate-mapping tools.
Third, although we confirmed the presence and localization of the various highlighted epithelial populations of interest via immunohistochemistry, we were not able to effectively validate the spatial relationships between all signaling partners (e.g., macrophage, mesenchymal, and hillock cells). The reasoning behind these limitations are 2-fold: first, we found that our organoids exhibited such significant structural complexity that when it came time to dissociate them (i.e., digesting the Matrigel), they had remodeled the culture environment to such a large degree that dissociation was quite difficult. Additionally, organoids in the BAL_3D system, specifically, frequently lost structural integrity quickly upon dissociation, making them difficult to manipulate for whole-mount staining. Furthermore, there remains a limited availability of high-quality rat-reactive antibodies for the detection of the specific molecules of interest in our systems. Future application of spatial transcriptomics or high-resolution in situ hybridization approaches will be crucial to determine whether observed ligand-receptor circuits correspond to physically adjacent cell types and to define how niche topology may constrain or facilitate regenerative signaling.
Finally, this model captures many aspects of regenerative signaling and tissue self-organization, but it does not incorporate biomechanical inputs (e.g., stretch, airflow, and perfusion), which are known to affect epithelial differentiation, cytoskeletal dynamics, and extracellular matrix remodeling in vivo. For example, cyclic stretch and airflow have been shown to regulate alveolar type I/II fate balance, in addition to modulating surfactant production and progenitor cell activation in the native lung.^132^^,^^133^ Mechanical strain has also been implicated in the mediation of mesenchymal-epithelial interactions during branching morphogenesis^134^^,^^135^^,^^136^ and reported to bias immune cells toward polarized states via mechanotransduction pathways such as YAP/TAZ and PI3K/AKT/mTOR.^137^^,^^138^^,^^139^ Incorporating biomechanical forces into future versions of this model, whether via microfluidic perfusion, cyclic stretch chambers, or bioreactors, may enhance physiological fidelity and yield new insight into how physical stimuli intersect with lineage composition and intercellular communication.
Our findings suggest that transitional cells are part of a larger regenerative plastic intermediates that integrate local cues to orchestrate repair, calibrate inflammation, and potentially resist progression to disease.
Based on this work, we propose several questions for future study.•What determines whether transitional states resolve or persist? Can Sox9^+^ or hillock-like cells revert, differentiate, or become pathologically stabilized depending on microenvironmental context?•Is the hillock niche an inducible, facultative module for immune-epithelial coordination in the proximal lung or a universal regenerative structure across airway compartments? How does its formation relate to chronic remodeling in disease?•How tightly coupled are macrophage polarization and epithelial transitions? What additional signaling axes beyond the IL-1 family and Serpina1 axes may govern their coupling or uncoupling?•Can lineage manipulation at the start of culture redirect the subsequent regenerative trajectories we observed, and how do immune or mesenchymal biases (imposed at seeding) influence the fate spectrum or propensity for repair versus maladaptive remodeling?
We present these experiments and findings to the community to advance understanding of how compositional and signaling cues guide epithelial state transitions. This knowledge may open new therapeutic avenues for restoring tissue integrity in chronic lung disease.
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Micha Sam Brickman Raredon (michasam.raredon@yale.edu).
This study did not generate new unique reagents or devices.
•Raw FASTQ files, extracted digital gene expression matrices, and all R objects, containing relevant metadata, used in this article are available at Gene Expression Omnibus GEO accession [GSE299479](GSE299479) and are publicly available as of the date of publication.•Original immunofluorescence and brightfield imaging data have been deposited at Mendeley at https://doi.org/10.17632/rkcvf4299m.2 and are publicly available as of the date of publication.•All original code used for the analysis and figure craft has been made publicly available on GitHub at https://github.com/RaredonLab/Edelstein2025.•Any additional requests for data or software may be directed to and will be fulfilled by the lead contact, Micha Sam Brickman Raredon (michasam.raredon@yale.edu).
We would like to thank Yale Center for Genome Analysis (YCGA), the Yale Pathology Tissue Services (YPTS), and Yale Center for Research Computing (YCRC) for their efforts in making this work possible. We would also like to acknowledge Dr. Allison Greaney and Dr. Themis Kyriakides for their feedback and support in story craft, single-cell analysis methodology, and overall mentorship. We would like to acknowledge Ako Ndefo-Haven for his assistance with line-editing.
Conceptualization, S.E.E., M.S.B.R, S.M., and M.S.; methodology, S.E.E., M.S.B.R., S.M., and N.W.; software, S.E.E., M.S.B.R., N.W., and H.K.; validation, N.W., H.K., S.M., and C.H.; formal analysis, S.E.E., N.W., and M.S.B.R; investigation, S.E.E., S.M., M.T.G, and C.D.; resources, M.S.B.R; writing – original draft, S.E.E. and M.S.B.R; writing – review and editing, S.E.E., M.S.B.R., M.S., V.L., C.H., and C.D.
M.S.B.R. holds stock in and consults for Humacyte Inc., a regenerative medicine company. Humacyte did not influence the conduct, description, or interpretation of the findings in this report.
REAGENT or RESOURCESOURCEIDENTIFIERAntibodiesRabbit anti-CD45 (Ptprc)AbcamCat #ab10558; RRID: AB_442810Mouse anti-CD45 (Ptprc)BD BiosciencesCat #554875; RRID: AB_395568Goat anti-CCSP (Uteroglobin/Scgb1a1)EMD MilliporeCat #ABS1673; RRID: AB_2910611Goat anti-UGRP1/SCGB3A2R&D SystemsCat #AF3545; RRID: AB_2183543Mouse anti-RTI-40 (Podoplanin)Terrace BiotechCat #TB-11ART1-40; RRID: AB_2892067Mouse anti-Acetylated TubulinSigmaCat #T7451; RRID: AB_609894Rabbit anti-Krt13AbcamCat #ab92551; RRID: AB_2134681Mouse Anti-Cytokeratin 14 (clone RCK107)AbcamCat #ab9220; RRID: AB_307087Rabbit Anti-Prosurfactant Protein CMilliporeCat #AB3786; RRID: AB_91588Rabbit Anti-FOXJ1SigmaCat #AV38038; RRID: AB_1849069Rabbit Anti-Aquaporin 5MilliporeCat #AB3559-50UL; RRID: AB_2141915Rabbit anti-Sox9AbcamCat #ab185966; RRID: AB_2728660Rabbit anti-Krt8AbcamCat #ab59400; RRID: AB_942041Rabbit Anti-DCAMKL1AbcamCat #ab31704; RRID: AB_873537Donkey anti-Rabbit IgG (H+L) Alexa Fluor 647InvitrogenCat #A31573; RRID: AB_2536183Donkey Anti-Goat IgG (H+L), Alexa Fluor 555InvitrogenCat #A21432; RRID: AB_2535853Donkey Anti-Mouse IgG (H&L), DyLight 488AbcamCat #ab96875; RRID: AB_10698084Donkey Anti-Mouse IgG (H+L), Alexa Fluor 488InvitrogenCat #A21202; RRID: AB_ 141607Chemicals, peptides, and recombinant proteinsAdvanced DMEM/F12 Basal MediumThermoFisherCat #11320033Collagenase/DispaseRocheCat #10269638001ElastaseWorthingtonCat #LS002292RNase-Free DNaseQiagenCat #792540.4% Trypan BlueGibcoCat #15250061Amphotericin BHyCloneCat #SH30071.03Penn/Strep (P/S)GibcoCat #15140-122HeparinSigma AldrichCat #H3149-10KUTransferrinSigma AldrichCat #10652202001Bovine Serum AlbuminGemini Bio-ProductsCat #700-100PACK Lysing BufferGibcoCat #A1049201Fetal Bovine Serum (FBS)GibcoCat #A5256701Dulbecco's Phosphate-Buffered Saline (DPBS)GibcoCat #14190250M-280 Sheep anti-Rabbit IgG DynabeadsThermoFisherCat #11203DFibroblast Growth Factor 10 (FGF10)PeprotechCat #400-42Fibroblast Growth Factor 9 (FGF9)Novus BiologicsCat #NBP-35196Epidermal Growth Factor (EGF)PeprotechCat #400-25CHIR99021 (GSK-3 Inhibitor)Cayman ChemicalsCat #13122-10BIRB796 (p38-MAPK inhibitor)Cayman ChemicalsCat #10460-10Y27632 (ROCK Inhibitor)Cayman ChemicalsCat #10005583-10A8301(ALK-5/TGF-beta/NODAL Inhibitor)Cayman ChemicalsCat #9001799-10InsulinSigma-AldrichCat #910077C-100GlutamineThermoFisherCat #25030081Anti-AntiThermoFisherCat #15240062G418 (Genticin)ThermoFisherCat #10131035Donkey SerumMillipore SigmaCat #D9663DAPIThermoFisherCat #62247PVA-DABCOThermoFisherCat #10981Dimethyl Sulfoxide (DMSO)Fisher ScientificCat #AAA132800EGrowth Factor-Reduced Matrigel, Phenol-Red FreeCorningCat #3562310.05% Trypsin-EDTAGibcoCat #25300-054HistoGelThermoFisherCat #HG-4000-012Critical commercial assaysChromium Next GEM Single Cell 3’ Reagent Kit (v3.1)10x GenomicsCat #PN-1000268Deposited datascRNAseq dataThis manuscriptGEO: GSE299479Brightfield and immunohistochemistry microscopy dataThis manuscriptdoi: https://doi.org/10.17632/rkcvf4299m.2Experimental Cell linesPrimary Rat Lung Microvascular Endothelial Cells (RLMVECs)VEC Technologies–Experimental Organisms/strainsSAS Sprague-Dawley Rattus norvegicusCharles RiverStrain Code: 400Software and algorithmsSeurat (v5.2.0)Satija et al., 2015https://satijalab.org/seurat/Monocle3 (v1.3.7)Cao et al., 2019https://cole-trapnell-lab.github.io/monocle3/Slingshot (v2.10.0)Street et al., 2018https://github.com/kstreet13/slingshotTradeSeq (v1.16.0)Van den Berge et al., 2020https://github.com/statOmics/tradeSeqNICHES (v0.2.3)Raredon et al., 2023https://msraredon.github.io/NICHES/SeuratWrappers (v.0.4.0)–https://github.com/satijalab/seurat-wrappers/CellChat (v1.6.1)Jin et al., 2021; Jin et al., 2025https://github.com/sqjin/CellChatclusterProfiler (v4.16.0)Yu et al., 2012; Yu, 2024https://guangchuangyu.github.io/software/clusterProfiler/ComplexHeatmap (v2.24.0)Gu et al., 2016https://github.com/jokergoo/ComplexHeatmapEnhancedVolcano (v1.26.0)Blighe et al., 2025https://github.com/kevinblighe/EnhancedVolcanoImageJNIHhttps://imagej.nih.gov/ij/QuPathBankhead et al., 2017https://qupath.github.io/Infinity Analyze (v7.1.1.66)–https://www.teledynevisionsolutions.com/products/infinity-analyze/?model=infinityanalyze&vertical=tvs-lumenera&segment=tvsOther100 μm Cell StrainerCorningCat #43175270 μm Cell StrainerCorningCat #43175140 μm Cell StrainerCorningCat #431750ReadyProbes™ Hydrophobic Barrier Pap PenThermoFisherCat #R3777DynaMag™-5 MagnetThermoFisherCat #12303DFalcon Permeable Support for 24-well plate (0.4 mm)CorningCat #353095
All animal procedures (lung dissociation and bronchoalveolar lavage (BAL)) were approved by the Yale University Institutional Animal Care and Use Committee (IACUC protocol #11190) and performed in accordance with NIH Guidelines for the Care and Use of Laboratory Animals. Wild-type Sprague Dawley rats (Rattus norvegicus) were used as source animals for primary dissociation (PD) and BAL isolations. Adult male rats aged 8-10 weeks (n = 9; 250 g ± 25 g) were used for PD-derived starting populations and downstream organoid experiments. BAL collections were performed from male rats aged 8-10 weeks (n = 9; 250 g ± 25 g). Only male animals were used. Sex was not evaluated as a biological variable in this study and therefore, potential sex-associated differences cannot be excluded and represent a limitation in terms of study generalizability.
All animal procedures were approved by the Yale University Institutional Animal Care and Use Committee (IACUC protocol #11190) and performed in accordance with NIH Guidelines for the Care and Use of Laboratory Animals. Primary cells used in this study included rat lung-derived epithelial, immune, mesenchymal, and endothelial populations, sourced either directly from lung tissue or from commercial vendors. Primary pulmonary dissociation-derived cells and bronchoalveolar lavage-derived cells were isolated from adult male Sprague Dawley rats, as described in method details. For two-dimensional lineage culture survival studies, epithelial populations were derived from pulmonary dissociation isolates via Dclk1-based enrichment (see method details) and subsequent 2D culture expansion (multiple passages). Immune populations were obtained similarly, but from BAL isolates. Primary rat lung microvascular endothelial cells (RLMVECs) were purchased from VEC technologies (passage 1, isolated from 4-6-week-old rats). Primary neonatal rat lung fibroblasts (passage 1) were isolated from male Sprague Dawley pups (n = 5, 7 g ± 5) aged 7-9 days, as previously described^140^ and used as the mesenchymal population. No immortalized cell lines were used in this study. Culture conditions and handling procedures for all primary cells are described in greater depth throughout method details. Cell lines used were not tested for mycoplasma contamination.
8-10-week-old male Sprague Dawley rats (n = 9; 250 g ± 25) were sacrificed for cell isolation. Rats were sedated in an induction chamber containing Isoflurane-soaked (20% w/v) gauze, followed by intraperitoneal administration of 0.25 mL Ketamine-Xylazine solution (K: 75 mg/mL; X: 5 mg/mL) and 0.15 mL heparin (1,000 U/mL). Once fully sedated, animals were sterilized with 70% ethanol and povidone-iodine prep pads (Dynarex). Through the thoracic cavity, the lungs were allowed to deflate, and the thymus was removed. The clavicle was dissected to expose the trachea, which was cannulated using a barbed Y 1/16” connector. A second barbed Y 1/16” connector was used to cannulate the pulmonary artery. The cardiac apex was excised to allow fluid flow, and the lung tissue was perfused at 50 mL/min with heparin (100 U/mL) and sodium nitroprusside (SNP; 0.1 mg/mL). The lungs were then excised and placed in a petri dish for perfusion with dissociation buffer (DMEM HG (Gibco), 1 mg/mL Collagenase/Dispase (Roche), 3 U/mL Elastase (Worthington), 20 U/mL DNase (Qiagen)) through both the airway and vasculature via gravity. 10 mL of dissociation buffer was used to inflate the lungs three times. The lungs were transferred to conical tubes containing dissociation buffer and incubated on a rocker at 37°C for 25 minutes. All animal procedures were conducted in accordance with Yale IACUC, as outlined in Experimental Model Details.
Cells were dissociated from native rat lung tissue following enzymatic digestion (described in Animal Handling and Lung Dissociation) and then mechanically disrupted using a spatula using a protocol previously reported by our lab.^140^^,^^141^^,^^142^ Large collagenous structures were allowed to pass through a strainer. The strainer was then washed with 20 mL of quenching medium (DMEM HG (Gibco, 11965092), 10% FBS, 1% P/S (Gibco), 0.1% gentamycin (Worthington), 1% Amphotericin B (HyClone)), and the suspension was collected. The suspension was centrifuged at 300 × g for 5 minutes at 4°C to pellet cells and supernatant was then discarded. The pellet was resuspended in an equal volume of ACK Lysing Buffer (Gibco, A1049201) and gently agitated by tapping for 2 minutes at room temperature to lyse red blood cells. This suspension was then diluted with 10 mL of MACS buffer (0.1% BSA (Gemini) in PBS (Gibco), filtered), and then centrifugation at 300 × g for 5 minutes at 4°C to pellet. Supernatant was aspirated and the pellet resuspended in 5 mL of MACS buffer before being passed through a 70 μm strainer (Falcon, 352350) and then a 40 μm strainer (Falcon, 352340). The filtered cell suspension was then centrifuged again at 300 × g for 5 minutes at 4°C to pellet cells, supernatant was aspirated, and the pellet resuspend in 5 mL of MACS buffer. To ensure a single-cell suspension, the suspension was passed through a 40 μm strainer twice more. Cells were then counted and viability evaluated using 0.4% Trypan Blue solution (Gibco, 15250061). Cells were then resuspended in MACS buffer at 1 mL per 10^7^ total cells before proceeding to cell selection using Dclk1^+^ tagged beads.
Following enzymatic digestion, cells were dissociated from the lung tissue using a spatula and passed through a strainer, leaving behind larger collagenous structures. The strainer was rinsed with 20 mL of quenching medium (DMEM HG (Gibco, 11965092), 10% FBS, 1% P/S (Gibco), 0.1% gentamycin (Worthington), 1% Amphotericin B (HyClone)), and the resulting cell suspension was collected. The suspension was centrifuged at 300 × g for 5 minutes at 4°C, and the supernatant was discarded. The cell pellet was resuspended in an equal volume of ACK Lysing Buffer (Gibco, A1049201) and agitated manually with gentle tapping for 2 minutes at room temperature to lyse red blood cells. The cell suspension was then diluted with 10 mL of MACS buffer (0.1% BSA (Gemini) in PBS (Gibco), filtered), followed by centrifugation at 300 × g for 5 minutes at 4°C. Supernatant was aspirated, and the pellet was resuspended in 5 mL of MACS buffer. The suspension was sequentially filtered through 70 μm and 40 μm strainers (Falcon, 352350, 352340) to remove debris. Cells were pelleted by centrifugation at 300 × g for 5 minutes at 4°C, and the supernatant was removed. The pellet was resuspended in 5 mL of MACS buffer, filtered twice through a 40 μm strainer, and counted, as well as evaluated for viability using 0.4% Trypan Blue solution (Gibco, 15250061). Cells were resuspended in the 1 mL of MACS buffer per 10^7^ total cells before proceeding to cell selection with Dclk1^+^ tagged beads.
BAL-derived (immune and epithelium) cells were harvested from 8-10-week-old male Sprague Dawley rats (n = 9; 250 g ± 25), as described above in Animal Handling and Lung Dissociation. Lungs were placed into a conical tube with 40 mL of ice-cold saline solution (PBS, 1% P/S) on ice for 60 minutes. The lungs were then maximally inflated by perfusing the trachea with sterile saline solution, gently massaging the tissue. Negative pressure was applied at the open end of the tracheal cannula by covering it with one finger and the PBS containing cells was aspirated using a 10 mL syringe. This process was repeated 5-7 times, with all the aspirated fluid being collected into a sterile 50 mL conical tube. The BAL cells were centrifuged (500 × g for 10 minutes at 4°C), supernatant aspirated, and resuspended in an equal volume of ACK lysing buffer (Gibco, A1049201) (1:1 ratio). The pellet was then agitated manually with gentle tapping for 2 minutes at room temperature before dilution with 2 mL DMEM supplemented with 10% FBS. The cells were counted and resuspended to the desired concentration for organoid seeding.
To evaluate the ability of LPM-3D medium to support the viability of distinct lung-derived compartments, representative epithelial, endothelium, immune, and mesenchymal populations were cultured under 2D conditions for 7 days. The epithelial population utilized was a heterogenous suspension of basal, ATII, and tuft cells, isolated using the previously described Dclk1-enrichment protocol. These cells were previously expanded in LPM-3D medium until non-epithelial populations were depleted. The endothelial cells utilized were passage 1, primary rat lung microvascular endothelial cells (RLMVECs; isolated from 4-6-week-old rats, VEC Technologies). Mesenchymal cells were primary neonatal rat lung fibroblasts (passage 1; isolated from 7-9-day old Sprague-Dawley pups as described previously.^136^ Immune cells were macrophages obtained by bronchoalveolar lavage (BAL) from Sprague-Dawley rats, isolated as described above in Bronchoalveolar Lavage (BAL) Cell Collection. Each population was seeded at a density of 1 × 10^5^ cells per well in standard 6-well tissue-culture-treated plates, with 3 mL of LPM-3D medium added per well. Cultures were maintained at 37°C and 5% CO2 with media changes every 24 hours. Cells were imaged by brightfield microscopy (Zeiss Axio Vert.A1) at days 1,3,5, and 7. Brightfield images were processed using the Infinity Analyze software (Teledyne Vision Solutions).
Organoid systems were generated in 0.4 μm, 24-well cell culture inserts (Falcon) using a 1 (v/v) ratio of growth factor-reduced Matrigel (Corning, Matrigel GFR, Phenol Red-Free, #356231) to cell suspension. Each transwell insert contained a total volume of 90 μL, consisting of 45 μL Matrigel and 45 μL cell suspension. All organoid suspensions (Matrigel + cells) were prepared in LPM-3D medium and maintained at 37°C. The following cell seeding densities were PD_3D (0.1 × 10^6^ cells from PD), Mixed_3D (5 × 10^4^ PD cells and 5 × 10^4^ BAL cells), and BAL_3D (0.1 × 10^6^ cells from BAL). Organoid suspensions (Matrigel + cell suspension) were prepared on dry ice, and chilled. Wide-bore pipette tips were used for sample distribution to ensure consistent handling and to avoid further physical agitation of the cells. Large bubbles in the Matrigel were removed using 27 G × ½ (0.4 mm x 13 mm; BD, 305109) sterile needles before Matrigel was allowed to solidify. Following a 30-minute incubation at 37°C, 500 μL of culture medium was added to the basolateral compartment of the Transwell system. Cell culture medium was changed every 48 hours across all conditions and media samples were saved for potential future experiments. Biological replicates were prepared for each condition (n=18, with n=6 per experiment). Organoids were maintained in culture for 10 days at 5% CO2 at 37°C. Organoids sent for single-cell RNA sequencing (scRNAseq) were dissociated on day 10 following the same methods outlined in STAR Methods section organoid dissociation for scRNAseq and other downstream applications.
Native rat lungs fixed via perfusion with 4% paraformaldehyde (PFA) at a flow rate of 60 mL/min for a minimum of 12 hours. After fixation, the lungs were transferred to a sterile petri dish for dissection. Anatomical dissection was performed to ensure that each of the five lobes and the trachea were embedded individually. Trachea and lobes were sectioned longitudinally. Tissue embedding and sectioning were performed by Yale Pathology Tissue Services (YPTS).
Organoids were fixed using one of three methods, depending on the intended downstream application. For organoids embedded and sectioned with the Transwell membrane, cell culture media was removed, and 0.5 mL of 4% PFA was added to the basolateral compartment of the Transwell insert. The plate was placed on a rocker for 16 hours without prior Matrigel digestion. After 12 hours of incubation, the Transwell membrane containing the Matrigel-embedded organoids was excised using a scalpel. The excised membrane was embedded in HistoGel (ThermoFisher, HG-4000-012) before paraffin embedding and sectioning by YPTS.
For organoids embedded in paraffin without Matrigel, samples were washed with ice-cold PBS and placed on an orbital shaker at 4°C for 60 minutes to digest as much of the ECM as possible without the use of enzymes. Organoids were collected in PBS using a P1000 wide-bore pipette coated with BSA to prevent adherence to the pipette tip. The collected organoids were centrifuged in a 15 mL conical tube coated with BSA at 200 × g for 3 minutes at 4°C. The supernatant was carefully aspirated to avoid loss of organoids, and the pellet was resuspended in 2 mL of 4% PFA. Organoids in 4% PFA were then placed on a rocker for a minimum of 4 hours at room temperature. Following fixation, PFA was aspirated, and the organoids were embedded in a 100-200 μL droplet of HistoGel (ThermoFisher, HG-4000-012) before paraffin embedding and sectioning by YPTS.
For organoids undergoing whole-mount staining, the same preparation steps as those used for organoids embedded without Matrigel were followed. However, instead of HistoGel embedding, organoids were resuspended in PBS after fixation with 4% PFA and transferred to 4-well chamber slides (ThermoFisher, Nunc Lab-Tek II) coated with 0.4% BSA in PBS, to prevent adherence. Staining and imaging was then carried out with the organoids in the chamber slides. For information on whole-mount staining, see immunohistochemistry of whole-mount organoids.
Epithelial expansion medium was prepared as previously described,^39^ with slight modifications to make the media appropriate for rat-derived cells. Briefly, Advanced DMEM/F12 basal medium (ThermoFisher, 11320033) was filtered and supplemented with the following rat growth FGF10 (50 ng/mL; R&D Systems, 7804-FG), FGF9 (50 ng/mL; R&D Systems, 273-F9), and EGF (50 ng/mL; R&D Systems, 3214-EG)(Table S1). Additionally, the medium was supplemented with small molecule inhibitors, including CHIR99021 (3 μM; GSK-3 inhibitor/Wnt activator; Cayman Chemical, 13122), BIRB796 (1 μM; p38-MAPK inhibitor; Cayman Chemical), Y27632 (10 μM; ROCK inhibitor; Cayman Chemical, 10460), and A8301 (1 μM; Activin/NODAL/TGFβ pathway inhibitor; Cayman Chemical, 9001799). Further supplements included heparin (5 μg/mL; Sigma Aldrich, 9041-08-1), insulin (10 μg/mL; Roche), and transferrin (15 μg/mL; Sigma Aldrich). The prepared medium was pre-warmed to 37°C before given to organoid systems at feeding. Antibiotics were added to prevent at concentrations of 1% penicillin-streptomycin (P/S) and 0.1% gentamicin to prevent infection.
Paraffin-embedded sections of perfused native rat lung, as well as trachea and organoid systems were immunostained to better understand protein expression localization patterns. Samples were first deparaffinized by heating at 65°C for 30 minutes, followed by sequential rinsing in a xylene and ethanol gradient (100%, 95%, and 70% ethanol). Antigen retrieval was performed by immersing the slides in antigen retrieval buffer (0.1 M citric acid, 0.05% Tween 20, pH 6) and placing them in a water bath at 75°C for 20 minutes. Slides were then cooled to room temperature for a minimum of 15 minutes.
Following antigen retrieval, slides were placed in a slide box containing PBS for 5 minutes. Tissue sections on the slides were then outlined using a hydrophobic marker (ReadyProbes™ Hydrophobic Barrier Pap Pen; ThermoFisher). Tissue sections were then rinsed 3-5 times with PBS. Two different permeabilization buffers were employed depending on the proteins being stained for. For nuclear stains (e.g. transcription factors such as Sox9), tissue sections were permeabilized using PBS supplemented with 0.2% Triton X-100 (Invitrogen, HFH10) for 15 minutes. After permeabilization, samples were blocked with blocking buffer (0.75% glycine, 5% Donkey Serum (Millipore Sigma, D9663) in PBS) for 1 hour to minimize nonspecific binding. For cell surface and cytoplasmic markers, tissue sections were permeabilized with 0.5% Tween-20 (Millipore, P1379) in PBS for 15 minutes, followed by the same blocking buffer step described previously. Primary antibodies, prepared in the same blocking buffer, were applied to the samples at their specified concentrations (Table S2) and incubated overnight at 4°C. Following incubation, slides were washed 3-5 times with PBS before the application of secondary antibodies (concentration: 500) for 1 hour at room temperature. Slides were washed again with PBS and counterstained with DAPI (ThermoFisher, 62247)(1:1000 in PBS) for 1 minute to visualize nuclei. Finally, sections were mounted using PVA-DABCO (Sigma Aldrich,10981). Details for all primary and secondary antibodies used are provided in Tables S2 and S3, respectively. All samples were imaged using either the EVOS Auto FL 2 imaging system or Stellaris 8 confocal microscope.
Organoids digested from Matrigel, as described above (see method details: organoid and tissue fixation), were transferred to a 4-well chamber slide (ThermoFisher, Nunc Lab-Tek II) coated with 0.4% BSA in PBS (to prevent adherence). Multi-well chamber slides allowed for the clear separation of condition-specific organoids. Roughly 8-12 organoids, depending on size, were transferred to each chamber using a P1000 wide-bore pipette tip coated with 0.4% BSA in PBS. The same protein-dependent permeabilization steps were performed as previously described (see method details: immunohistochemistry of paraffin embedded sections). Organoids were then blocked with 500 μL of blocking buffer (0.75% glycine, 5% Donkey Serum (Gibco, PCN5000) in PBS) for 2 hours at 37°C and placed on an orbital shaker rotating at 280 rpm. Blocking buffer was then carefully removed, ensuring not to disrupt organoids, and blocking buffer containing primary antibodies (at respective concentrations) were added to the chambers. Organoids with primary antibodies were incubated overnight at 4°C with gentle shaking on an orbital shaker (280 rpm). The next day, the primary antibody solution was removed carefully, and the samples were washed with 500 μL of PBS 3 times over, with a 10-minute incubation period between each wash. The samples were then incubated with 500 μL of secondary antibody (each at a concentration of 500) in blocking buffer for 2 hours at 37°C, with gentle rocking using the orbital shaker. Secondary antibody was then carefully aspirated, and samples were incubated with DAPI (ThermoFisher, 62247) solution (1:1000 in PBS) for 1 minute at 37°C. DAPI solution was then removed, and 250 μL of PVA-DABCO (Sigma Aldrich,10981) was added to each chamber before imaging on either the EVOS Auto FL 2 imaging system or Stellaris 8 confocal microscope.
Cytospins were generated for all three starting populations (BAL, PD, and a 1 mix of the two that subsequently became Mixed_3D). Slides were securely placed in the slide holder, ensuring proper orientation with the frosted side of the slide facing upwards. A filter card, with its absorbent side against the slide, was positioned alongside the cytofunnel, ensuring proper alignment of all components before securing the holder. Each assembled holder was then placed in the Cytospin centrifuge (Cytospin 4 Centrifuge, ThermoScientific). Each cell suspension, prepared at a concentration of 0.5 ×10^6^ cells/mL in DMEM + 10% FBS, was pipetted into each cytofunnel as a 200 μL volume. The slides were centrifuged at 1000 × g for 5 minutes. Following centrifugation, slides were removed from the centrifuge and taken out of the holders. Slides were placed cell-side up and left to air-dry for 1 hour. Spots containing cell suspensions were outlined with a hydrophobic marker (ReadyProbes™ Hydrophobic Barrier Pap Pen; ThermoFisher) before being fixed with 4% paraformaldehyde in PBS for 10 minutes at room temperature. After fixation, slides were rinsed 3 times with PBS and stored in PBS at 4°C until needed for staining (up to one week). Slides were processed as described in the previous section from the antigen retrieval step onward (see method details: immunohistochemistry of paraffin embedded sections). While BAL starting populations were stained and imaged, PD (for PD_3D) and BAL+PD (for Mixed_3D) were stained, but did not yield valuable data as the beads used for selection made it difficult to visualize cells.
Time course brightfield images of each condition (n=3) were taken at days 0, 3, 5, 7, 9, and 10 using the EVOS FL Auto 2 microscope. Stitched images of each well were generated. After image acquisition, images were analyzed using QuPath (v.0.5.1).^143^ Organoid areas in stitched images at each time point were recorded using both the ellipse and brush tools in QUPath. To facilitate accurate organoid area quantification (μm^2^), image scale was calibrated in QuPath using the embedded 1 mm scale bar present in each stitched image, which was measured by counting pixels spanning the scale bar using ImageJ and entered manually into QuPath image properties to ensure that all organoid area measurements reflected physical dimensions. Organoid count and area values were averaged across replicates per condition and timepoint. Time was converted to hours and line plots generated showing mean values ± standard deviation in R. Welch’s t-tests performed comparing conditions at each timepoint for both organoid count and area without assuming equal variance. Resulting p-values, adjusted using the Benjamini-Hochberg (BH) method to control the false discovery rate, are reported as p ≤ 0.05 (∗), ≤ 0.01 (∗∗), ≤ 0.001(∗∗∗).
Macrophage and organoid edge coordinates were extracted from immunofluorescence images using ImageJ. Ptprc^+^ macrophages and epithelial border cells of organoids were annotated using the multipoint tool. Prior to analysis, all immunofluorescence images were scaled appropriately using the embedded scale bar. For each field of view, all organoid epithelial edge points and Ptprc^+^ macrophage (x,y) coordinates were recorded. Euclidean distances were computed between each macrophage and the nearest organoid boundary. Coordinates were exported as .csv files and analyzed in R. For each macrophage, the shortest linear distance to the organoid edge was calculated using Equation 1.Equation 1d=(xmac−xedge)2+(ymac−yedge)2
The minimum value across all edge points for a given macrophage defined its nearest edge distance and all distances were retained as independent observations. Distances from multiple images (10x and 20x objectives) were pooled for downstream visualization and statistical summaries. Distance distributions are displayed in Figure S12.
Organoids were dissociated from Matrigel ECM using an adapted enzymatic digestion protocol previously described above (STAR Methods: pulmonary dissociation cell isolation and processing). The culture medium was aspirated, and 1 mL of ice-cold DPBS (Gibco, 14190144) was added to each well; the growth factor-reduced Matrigel (Corning, Matrigel GFR, Phenol Red-Free, #356231) matrix was disrupted by gently pipetting the suspension up and down approximately 20 times using a P1000 wide-bore tip coated with 0.1% BSA in PBS to prevent organoids from adhering to the walls of the pipette tip. The organoid suspension was transferred into a 15 mL conical tube pre-coated with anti-adherence solution or 0.1% BSA in PBS and centrifuged at 300 × g for 5 minutes at 4°C to pellet the organoid-cell suspension. Supernatant was aspirated, and the pellet washed with 0.5 mL of ice-cold PBS to remove as much of the residual Matrigel as possible before enzymatic dissociation. The pelleted organoids were resuspended in 1 mL of dissociation enzyme solution (DMEM HG (Gibco, 11965092), 1 mg/mL Collagenase/Dispase (Roche, 11097113001), 3 U/mL Elastase (Worthington, LS002292), and 20 U/mL DNase (Qiagen) and incubated at 37°C on an orbital shaker to provide gentle agitation and facilitate dissociation while minimizing mechanical stress. After 20-25 minutes of incubation, the suspension was gently pipetted up and down 5-10 times using a wide-bore tip before proceeding to filtration.
To terminate enzymatic activity, 10 mL of Advanced DMEM/F12 medium (ThermoFisher, 12634010) was added to the cell suspension-enzyme mixture. The cells were then filtered through a 100 μm cell strainer (Corning, 431752) and the rubber head of a 3 mL syringe plunger (BD, 309657) was used to gently push any remaining cell aggregates through the strainer. The suspension was again centrifuged at 300 × g for 5 minutes at 4°C, and the supernatant was discarded. The pellet was again resuspended in 5 mL of LPM-3D medium and filtered through a 70 μm cell strainer (Corning, 431751), using a clean syringe plunger to push remaining aggregates through. The filtered suspension was centrifuged at 300 × g for 5 minutes at 4°C one last time before being filtered through a 40 μm cell strainer (Corning, 431750) to ensure a single-cell suspension. Cells were then counted using a hemocytometer to determine concentration and viability was assessed using 0.4% Trypan Blue solution (Gibco, 15250061) before being processed for downstream applications—either scRNAseq or passaging or cryopreservation. Samples dissociated for single-cell RNA sequencing were resuspended at 1 million cells/mL in 0.1% BSA in PBS. Cells to be cryopreserved were resuspended in the desired volume of freezing medium (90% FBS (Hyclone, SH30071), 10% DMSO (Fisher Scientific, AAA132800E)) before being transferred to -80°C for no more than 24 hours and then the cryogenic dewar for longer-term storage.
Dissociated organoid and starting population cell suspensions were prepared for single-cell RNA sequencing using the Chromium Next GEM Single Cell 3’ Reagent Kits v3.1, according to the manufacturer’s instructions (10x Genomics, Pleasanton, CA). After dissociation, cells were counted, assessed for viability using 0.4% Trypan Blue solution, and then resuspended at a concentration of 1 million cells/mL in 0.1% BSA in PBS for a targeted cell recover of 2,000-10,000 cells per sample. Libraries were sequenced by the Yale Center for Genomic Analysis (YCGA) at a target depth of 50,000 reads per cell using the Illumina NovaSeq 6000 platform. Alignment was performed using Cell Ranger (v8.0.1)(10x Genomics) and the reference transcriptome Rattus norvegicus.Rnor_6.0-95. Alignment was performed with the “-include - introns” option enabled, allowing both exonic and intronic reads to be considered for gene expression quantification.
Single-cell RNA sequencing (scRNA-seq) data were processed individually for each sample using the Seurat package (v5.2.0), following best practices.^144^^145^ Raw sequencing data files were read into R using the Read10X function, and Seurat objects were created for each sample with a minimum of 3 cells and 50 features. Quality control metrics were calculated, including the percentage of mitochondrial gene content (percent.mt), total unique molecular identifier (UMI) counts (nCount_RNA), and detected features (nFeature_RNA). Sample-specific filtering thresholds for nCount_RNA and nFeature_RNA were applied as detailed in Table S4.
Initial quality control (QC) was performed to remove low-quality cells based on mitochondrial gene expression and unique feature counts. Cells with high mitochondrial percentages or low feature counts were excluded from further analysis. To ensure accurate identification and removal of low-information cells, an initial clustering step was performed at a high resolution (e.g., res = 5.0) to finely delineate clusters. This allowed us to identify and remove clusters with low gene expression or high proportions of ambient RNA contamination. Once low-information cells were filtered out, we re-clustered the remaining cells at an optimized resolution to obtain biologically meaningful clusters. QC plots for the starting BAL population can be found in Figure S5. QC plots for the starting PD population can be found in Figure S6. QC plots for the organoid sequencing data (merged) can be found in Figure S4.
Dimensionality reduction was conducted using principal component analysis (PCA), followed by Uniform Manifold Approximation and Projection (UMAP) for visualization. Clustering was performed using the Louvain algorithm implemented in Seurat. UMAP dimensionality reduction was performed using RunUMAP with the selected principal components, followed by construction of a shared nearest neighbor (SNN) graph using FindNeighbors(). Clustering was conducted using the FindClusters() function, a Louvain algorithm implemented in Seurat, with a range of resolution values to achieve optimal clustering granularity.
To identify marker genes for each cluster, differential expression analysis was performed using FindAllMarkers() with a minimum percentage threshold of 0.1 and log-fold change threshold of 0.1. See Tables S5, S6, and S7 for marker lists generated. Top markers were selected based on power metrics (ratio ∗ log-fold change), and heatmaps, as well as volcano plots were generated to visualize cluster-specific expression patterns. Heatmaps and volcano plots included in manuscript figures were generated using the ComplexHeatmap (v2.24.0)^146^ and EnhancedVolcano (v1.26.0)^147^ packages.
Cell type annotation was performed iteratively. To ensure confidence in our final cell type annotations, we consulted multiple reference datasets when assessing the most highly expressed genes of each cluster. However, prior to cell type annotation, clusters were grouped by cell class, using lineage-specific markers that have proven to be most consistent when studying rat lung biology.^148^ We used the following cell class Epcam (epithelium), Col1a1 (mesenchyme), Ptprc (immune), and Cdh5 (endothelium). Following this, cell type annotations were informed by an in-house adult rat lung atlas, created via integrated scRNAseq data from 14 different rat lung samples of both sexes.^149^ Additionally, we leveraged publically-available single-cell transcriptomic atlases, including LungMAP^150^ and PanglaoDB.^151^ These atlases provide reference datasets for both pulmonary (LungMAP) and other tissue-specific cell populations (PanglaoDB). Specifically, we utilized the LungMAP repository to cross-reference transcriptional signatures of pulmonary epithelial cell types such as ATI, ATII, basal, and secretory in adult human lung data with those in rat. While some signature genes were conserved (homologs), not all were. This is to be expected given knowledge that humans and rats evolved along different paths and some genes have evolved to have species-specific functions.^152^ Similarly, PanglaoDB, served as a useful reference repository for confirming the distinction between immune and endothelial subtypes, beyond the reliance on literature precedent. Justifications for cell type annotations (by gene expression) for our three main single-cell objects can be found in Figures S1 (BAL), S2 (PD), and S7 (Organoids: BAL_3D, Mixed_3D, PD_3D).
Specific cell type populations, such as the ‘Stressed_Progenitors’, were identified by also using gene set enrichment analysis (GSEA) of upregulated genes in that cluster of cells. This method was performed using the clusterProfiler package (v4.16.0) in R.^153^ Final annotations were validated using feature plots and differential gene expression analysis using our global (merged) Seurat object of all three organoid conditions (n=3 for each condition)(see Tables S5, S6, and S7 for specific population markers). R markdown files walking through the cleaning, clustering, and annotation process for each sample can be found on our GitHub https://github.com/RaredonLab/Edelstein2025.
To test whether differences in epithelial trajectories and organoid morphology could be explained solely by initial epithelial abundance, we calculated an Epithelial Fold Expansion (EFE) metric for each condition. For each replicate within a condition (BAL, PD, BAL_3D, Mixed_3D, PD_3D), we computed the proportion of epithelial cells relative to the total population at both day 0 and day 10 (Equation 2).⌈EFEc=p¯c,Day10p¯c,Day10whereforeachcondition(c)andreplicate(r):Equation 2|pc,t,r=#{epithelialcellsinreplicater(OrigID)attimet}#{allcellsinreplicater(OrigID)attimet},p¯c,t=1Rc,t∑r=1Rc,tpc,t,r
Here, ⌊pc,t,r denotes the epithelial fraction for replicate ⌈r of condition ⌈c at time ⌈t, and ⌊p¯c,t is the average across all replicates (⌊Rc,t) for that condition and timepoint.
To assess uncertainty, we applied a nonparametric bootstrap approach (R = 2000 resamples) to generate 95% confidence intervals for Day 0, Day 10, and EFE values. Bootstrapping was performed independently within each condition, resampling replicates with replacement to preserve biological variability.
To better characterize transcriptional programs distinguishing basal hillock cells from luminal hillock cells, we performed gene set enrichment analysis (GSEA; Broad 2022 release, v7.5.1) using the Hallmark gene sets from MSigDB.^58^
The pseudobulk BAL_3D and PD_3D/Mixed_3D comparisons were done using differential expression analysis and genes were ranked by log2 fold-change (avg_log2FC). This marker list was then cleaned, ensuring gene symbols were included as an independent column in the data-frame and any row with a missing value (NA) was removed. The resulting data-frame (marker list) was passed through the GSEA() function in the clusterProfiler package (v4.16.0),^154^^,^^155^^,^^155^ using the Hallmark gene sets for Rattus norvegicus (msigdbr) pulled from MSigDB.^151^ We defined significantly enriched pathways by adjusted p-values < 0.05.
To interpret the different components of each upregulated pathway and the implications of pathway upregulation within our systems, genes for each significant term were annotated using curated ground truth gene lists. A set of rat-specific ligands were obtained from the NICHES FANTOM5 reference database.^75^ For matrix proteins, we curated our own list by aggregating the categories “ECM Glycoproteins,” “Collagens,” and “Proteoglycans” from the Matrisome Project (MatrisomeDB).^156^ A list of transcription factors was compiled from the AnimalTFDB database.^157^ From here, we identified the number of ligands in each pathway (excluding matrix proteins), the number of matrix-associated genes, and the identity of transcription factors present in the leading-edge gene set. We then identified the number of ligands in each pathway (excluding matrix proteins), the number of matrix-associated genes, and the transcription factors present in the leading-edge gene set. In the context of these results, a positive normalized enrichment score (NES) was interpreted as being enriched in basal hillock cells while a negative NES indicated enrichment in luminal hillock cells.
Trajectory analysis and subsequent analyses of differential gene expression across pseudotime were conducted using the Monocle3, Slingshot, and TradeSeq packages. These analyses allowed for the exploration of cellular transitions within emergent populations of interest such as the hillock cells and pro-/anti-inflammatory macrophages. The combination of these methods utilized provided us with a robust way of inferring lineage associations and transcription-level differences across multiple datasets.
Monocle3 (v1.3.7) was utilized to infer cellular trajectories and identify key transcriptional changes.^158^ Previously processed Seurat objects for each emergent population were used as input. Raw gene expression counts, and associated cell metadata were extracted from the Seurat objects, and a Monocle3 cell data set (CDS) was constructed with gene annotations. Dimensionality reduction embeddings, including PCA and UMAP, were transferred from the original Seurat object to the Monocle3 CDS to maintain consistency across analyses. Parameters were set based on the Seurat subset being studied.
To minimize short, noisy bifurcation points, cells were clustered using Monocle3's clustering algorithm with a minimum branch length of 8. A Euclidean distance ratio of 2.1 and a geodesic distance ratio of 2.1 was used to optimize trajectory learning, and the graph pruning function was enabled to remove spurious connections along the inferred trajectory. Cells were ordered along the inferred trajectory to model their progression through biological states, with the defined starting cluster being the cycling proximal epithelium. Pseudotime values were extracted and added back into the original Seurat metadata.
Differential gene expression analysis was performed along the pseudotime trajectory using Moran’s I test with a statistical significance threshold set at a false discovery rate (FDR) of < 0.05. Genes were ranked based on their q-values and prioritized using Moran’s I statistic to identify key regulators of cellular transitions. Genes of interest were analyzed for their dynamic expression patterns across pseudotime. Trajectory plots were generated to visualize temporal expression trends. Monocle3 analysis scripts are available in our GitHub repository.
To independently confirm the trajectory results obtained by Monocle3, Slingshot (v2.10.0) was applied to fit principal curves to the UMAP-reduced space.^159^ Emergent populations were converted from Subset Seurat objects into SingleCellExperiment (SCE) format for compatibility with Slingshot; PCA and UMAP embeddings were carried over, and the cluster labeled “Cycling_Proximal_Epi” was used as the starting cluster for pseudotime ordering (Hillock analysis).
Slingshot inferred lineage relationships based on cluster-based lineage tracing guided by the cell type annotations. Default smoothing splines were applied to fit the trajectory curves and pseudotime values were extracted and included in Seurat metadata for comparisons with Monocle3-derived pseudotime trajectories. Differential gene expression analysis along Slingshot-derived trajectories was conducted using TradeSeq (v1.16.0), retaining genes expressed in at least 1% of cells with counts greater than 5.^160^ The number of knots for smoothing was set to 15 based on the Akaike Information Criterion (AIC) that we determined by comparing models using different numbers of knots from 3 to 15. A test of association was performed between pseudotime and genes with a significant association (FDR < 0.05), followed by start-vs-end test to evaluate gene expression changes between early vs late states in the pseudotime trajectory. Scripts for Slingshot analysis are available in our GitHub repository.
For the BAL_3D epithelium, epithelial cells were subset from the global interacted Seurat object (Condition = BAL_3D; CellClass.NodeAligned = Epithelium) and reprocessed (NormalizeData, FindVariableFeatures, ScaleData, and RunPCA). Replicates were integrated with RPCA-based integration and embeddings were generated for the integrated reduction. The same steps described above (Slingshot Analysis: Hillock Subset) were applied to the BAL_3D epithelium subset object. In this analysis, however, ‘Cycling_Proximal_Epi’ was labeled as the trajectory starting cluster and ‘ATI_Like’ was labeled as the trajectory end cluster. Pseudotime values were extracted for all points and saved in the object meta-data. For visualization, raw pseudotime values were plotted on the integrated UMAP, and a LOESS-smoothed pseudotime field (span = 0.28, degree = 2) was generated.
To summarize composition changes across the trajectory, cells were binned into pseudotime deciles and the relative proportions of epithelial subtypes (CellType.NodeAligned) were calculated per bin. These data were visualized as stacked bar plots to highlight dynamic shifts in epithelial composition within this system across pseudotime.
NICHES (v0.2.3) was used to infer ligand-receptor-mediated cell-cell communication across experimental conditions. All analyses were conducted in R (v4.2.2) using Seurat (v4.3.0) and SeuratWrappers (v.0.4.0). Global NICHES analysis was performed on an integrated Seurat object of all engineered cells (PD_3D, Mixed_3D, BAL_3D, 12 samples total; n = 3 per condition). A subset analysis was conducted on a subset of object specific populations of interest to afford a finer analysis of cell-cell connectomics.
Cells were first subset by condition and split by replicate (Orig_ID). For each replicate, cell-cell signaling was inferred using the RunNICHES function with the FANTOM5 ligand-receptor database and the Rnor_6.0 rat reference genome. We enabled only the CellToCell inference mode to focus on direct intercellular signaling. Outputs were merged across replicates and combined using JoinLayers. Cells with low signaling information (nFeature_CellToCell < 60) were removed. The resulting object was normalized, scaled, and run through PCA (100 components). UMAP embedding and unsupervised clustering were performed using the top 18 principal components (res = 0.4). To correct for batch effects, RPCA integration was performed by condition (PCs 1-16, 18-19) and by sample (PCs 1-14, 16-18). The final integrated object (eng.CTC.final) served as the basis for downstream signaling analyses. The script for this analysis can be found in our GitHub repository.
To study regenerative signaling in more detail, a focused analysis was performed on a subset of six cell populations from our condition-level Hillock_Luminal, Hillock_Basal, Rspo3+_Mes, Pdgfrb+_Pericytes, Anti_Inflamm_Mac, and Pro_Inflamm_Mac. These subsets were extracted from the global, node-aligned object. Following normalization and feature selection, PCA was performed, and PCs 5, 7, 9-12, and 14 were used for UMAP and clustering (resolution = 0.4). NICHES was then run per replicate for each condition. Cell types with fewer than two cells per replicate were excluded. Merged CellToCell results were filtered to retain cells with ≥40 features and processed using PCA (PCs 10, 13-14, 17-18, 21, 23), UMAP, and clustering. The final object (regen.subset_CTC_byCondition) was used for downstream visualization and differential signaling analysis using a composite “power” metric (avg_log2FC ∗ expression ratio). The script for this analysis can be found in our GitHub repository.
To test how mesenchymal presence may have supported the maintenance of immune cells, a subset object including only mesenchymal populations (Pdgfrb^+^_Pericytes, Rspo3^+^_Mes) and immune (Polarized_Mac) in PD_3D and Mixed_3D was created. This subset was created from the global, node-aligned object. This object was processed via our standard pipeline (scaling, feature selection, PCA, embedding, and clustering).
We visualized signaling mechanisms of interest as directed multicellular circuits (Figure 6) using the CircuitPlot() function in the NICHESMethods GitHub repository.^75^^,^^148^^,^ This function uses a custom plotting algorithm (ggCircuit) that employs trigonometric positioning, with the help of ggplot2^161^ rendering, to visualize ligand-receptor interactions (e.g., Il1a—Il1r2) as directional arrows connecting cell types defined by a metadata grouping variable (e.g., CellType) and edge aesthetics that reflect connection strength. Layout parameters such as graph.angle, offset, and unity.normalize were occasionally passed to the function, when appropriate, and all circuits were generated from the subset connectomic object (regen.subset_CTC_bySample) using annotations from CellType.regen.spec to ensure consistency in node identity and classification.
To complement NICHES analysis, CellChat (v1.6.1) was applied to our global dataset to evaluate communication to and from the hillock population.^61^^,^^62^ Whereas NICHES computes ligand-receptor signaling scores at single-cell resolution and retains sender-receiver heterogeneity, CellChat aggregates these probabilities to the level of annotated clusters.^61^^,^^62^^,^^75^^,^^62^^,^^75^ This cluster-level framework allowed us to visualize the overall number, strength, and directionality of signaling interactions between hillock cells and other populations. To conduct this analysis, normalized RNA expression values and annotated cell identities (CellType.NodeAligned) were pulled from the global object to create a CellChat object. A rat-specific ligand-receptor database was sourced.^162^ Overexpressed genes and interactions were identified and signaling probabilities were computed using the truncated mean method with 20 bootstraps. Communication networks were aggregated at the cluster level, and both pathway-level and gene-level interactions were derived. To visualize hillock-specific communication, we generated circle plots to summarize the number of statistically significant L-R pair interactions from hillock cells to all other cells and from all other cells to hillock cells. We generated chord plots to look at multiple signaling axes that were upregulated between cell types in one plot.
All statistical analyses were performed in R (v4.2.2 or later) unless otherwise states, with scRNA-seq processed in Seurat (v5.2.0) and downstream analyses using clusterProfiler (v4.16.0), Monocle3 (v1.3.7), Slingshot (v2.10.0), TradeSeq (v1.16.0), NICHES (v0.2.3), and CellChat (v.1.6.1). Sample sizes (n) represent independent biological replicates as indicated in the STAR Methods and figure legends. Animal tissue isolations used 8-10-week-old male Sprague Dawley rats (n=6; 3 for BAL, 3 for PD), organoid cultures were prepared across independent experiments (n=18 total; n=6 per experiment, Figures 1 and 2), scRNA-seq was performed on Day 10 organoids with n=3 replicates per condition (Figures 1B–1F, 1X–1Z, 2X–2Z, 3D, 4A, and 5C), and longitudinal brightfield quantification used n=3 wells per condition per timepoint (Figures 1G and 1H). Summary statistics are presented as mean ± SD unless otherwise specified (Figures 1G, 1H, and 5G) and boxplots display replicate-level distributions (Figures 5D and 5G). Organoid number and mean organoid area over time were compared using two-sided Welch’s t-tests at each timepoint with Benjamini-Hochberg (BH) correction for multiple comparisons, with significance thresholds reported in figure legends (ns, not significant, ∗, p ≤ 0.05, ∗∗, p ≤ 0.01, ∗∗∗, p ≤ 0.001)(Figures 1G and 1H). Differences in cell type proportions derived from scRNA-seq data are reported descriptively as normalized proportions across conditions, with replicate number indicated in the figure legends (Figures 1B–1F, and 2aa–2dd). Differences in Sox9^+^ cell proportions were assessed using two-proportion Z-tests on raw Sox9^+^ counts with BH correction (Figure 5D), while comparisons of RAS-like signature positivity across engineered conditions were performed using a Krusal-Wallis test (Figure 5G). Uncertainty in epithelial fold expansion (EFE) estimates was quantified using non-parametric bootstrapping (2,000 resamples) to generate 95% confidence intervals (method details; Figure S11).