Single-cell RNA sequencing reveals profibrotic roles of distinct epithelial and mesenchymal lineages in pulmonary fibrosis

Single-cell RNA sequencing provides new insights into pathologic epithelial and mesenchymal remodeling in the human lung.


INTRODUCTION
Pulmonary fibrosis (PF) is a form of chronic lung disease characterized by progressive accumulation of extracellular matrix (ECM) in the peripheral lung, accompanied by destruction of functional alveolar gas exchange units (1). PF is a clinical phenotype that reflects the end stage of chronic interstitial lung diseases (ILDs) of both known and unknown causes. The most severe form of PF, idiopathic PF (IPF), is a relentlessly progressive disorder of unknown cause that leads to respiratory failure and death or lung transplantation within 3 to 5 years of diagnosis (2). Spatially heterogeneous expansion of interstitial collagen and ECM (i.e., "temporal heterogeneity") is a defining feature of the histopathology of IPF (3). Genetic evidence suggests that the lung epithelium plays a central role in disease susceptibility and initiation (4)(5)(6), while epithelial-mesenchymal interactions have been proposed to play a central role in driving pathologic epithelial remodeling and ECM expansion in PF (7). While there have been advances in identifying factors that regulate fibrosis in experimental models (8), to date, there has been limited progress toward developing an integrated understanding of the central mechanisms driving pathologic epithelial remodeling and ECM expansion in the human lung. Bulk tissue-based genomic studies have provided insights into the molecular programs that characterize PF lungs (9)(10)(11)(12)(13), but the inherent spatial heterogeneity of disease pathology and largescale changes in cell types present in diseased tissue have limited the ability of these studies to identify key mechanistic mediators. Singlecell biology approaches including single-cell RNA sequencing (scRNAseq) offer the potential to overcome these limitations inherent to bulk tissue-based analyses (14,15), and throughout the past decade, advances in microfluidic technologies now enable simultaneously transcriptomic profiling of thousands of individual cells (16)(17)(18). Initial proof-of-concept studies have demonstrated feasibility of these approaches in the airway epithelium (19,20) and peripheral lung in PF (21)(22)(23)(24), with a focus on macrophages phenotypes (21,24) and type II alveolar epithelial cells (22,23) in a small number of samples.
To develop a more comprehensive understanding of the diversity of cell types, states, and molecular programs in PF lungs, we performed scRNA-seq of single-cell suspensions generated from 20 PF and 10 nonfibrotic control lungs. Unbiased analyses revealed a spectrum of pathologic epithelial cell states and programs that are found in the fibrotic distal lung. We identified multiple distinct epithelial and mesenchymal lineages contributing to ECM expansion in PF lungs. In addition to myofibroblasts, we identified a unique HAS1 hi ECMproducing population, which is markedly enriched in lungs from patients with IPF and localizes to peripheral and subpleural regions. We also found a previously undescribed KRT5 − /KRT17 + epithelial cell population expressing collagen and other ECM components that is conserved across a subset of histopathologic patterns of PF. Together, these high-resolution transcriptomic data and the identification of multiple previously undescribed pathologic cell types provide remarkable insights into the cellular architecture of the human lung and the fundamental mechanisms driving disease pathology in PF.

RESULTS
To determine the cellular populations and mediators shared across different forms of PF, we generated single-cell suspensions from peripheral lung tissue of explanted lungs from patients with IPF (n = 12), chronic hypersensitivity pneumonitis (n = 3), nonspecific interstitial pneumonia (NSIP; n = 2), sarcoidosis (n = 2), unclassifiable ILD (n = 1), and nonfibrotic controls (declined donors; n = 10; tobacco users, 8 of 10) (table S1) and performed scRNA-seq using the 10x Genomics Chromium platform (see Materials and Methods and Fig. 1A). The samples were collected and processed at two different sites (table S2, detailed metadata from each sample); however, both sites collected cases and controls. In an effort to maximize our ability to identify rare cell populations, we jointly analyzed data from all samples. We defined inclusion criteria for cells based on observations from the entire dataset, removed low-quality cells accordingly, then performed dimensionality reduction, and unsupervised clustering of the 114,396 recovered cells using the Seurat (25,26) package in R (see Materials and Methods and fig. S1). Notably, we did not observe over batch effects driven by processing site or sequencing batch in our dimensionality reduction and visualization ( fig. S2). This may   (E) Cell type of origin and disease state informed expression of selected biomarkers and putative mediators of PF. (F) Heatmap depicting relative expression (normalized and scaled z-scored) of known ECM components previously shown to be increased in PF lungs, each cluster is downsampled to 100 cells for visualization. All of the genes except COL6A1, COL6A2, LAMA3, LAMA5, and LAMB2 were differentially expressed in at least one cell type. NK cells, natural killer cells; pDCs, plasmacytoid dendritic cells; cDCs, classical dendritic cells; cHP, chronic hypersensitivity pneumonitis. be due in part to the normalization and variance stabilization approach used in Seurat V3 (26). Using canonical lineage-defining markers to annotate clusters, we defined 31 cell types/states in the lung (see  Materials and Methods, Fig. 1B, table S3, and figs. S2 and S3). Most of cell types were identified both in diseased (PF) and control lungs as well as across diagnoses ( Fig. 1C and table S4). Comparing cells originating from PF lungs to controls, we identified hundreds of differentially expressed genes in each cell type (see Materials and Methods, Fig. 1D, and tables S5 and S6). As expected, this analysis identified a number of previously undescribed expression changes associated with PF, and because of the analysis being performed at a cell type resolution, we cataloged patterns of shared and cell typespecific differentially expressed genes ( fig. S4). We also used the heterogeneity in diagnoses to identify gene expression changes between patients diagnosed with IPF and those with a non-IPF diagnosis (table S5). We next examined a small number of known biomarkers for PF (27,28) and genes linked to PF disease mechanisms (Fig. 1E) (6,29,30). All of these genes except for SMAD3 were significantly differentially expressed in at least one cell type (table S5). There was considerable heterogeneity among the cell type-specific expression patterns of these genes. For example, MMP7, one of the most highly validated biomarkers for IPF, has expression changes across many cell types, whereas CDKN2A is more restricted. Together, these data suggest that multiple distinct epithelial and mesenchymal cell types are involved in pathologic tissue remodeling in PF. Turning our analysis to genes encoding for ECM components, we identified multiple cell types expressing such genes that have previously been reported to be increased in IPF lungs (31). All but five of these genes were differentially expressed in at least one cell type, and as expected, these ECM genes were largely expressed in fibroblast subpopulations. Further, while subsets of these genes (basement membrane components from Fig. 1F and fig. S5) were expressed in smooth muscle cells, endothelial cells, and type 1 alveolar epithelial (AT1) cells and AT2 cells ( Fig. 1F and fig. S5), this analysis also identified a previously undescribed population of KRT5 − /KRT17 + epithelial cells that expressed COL1A1, FN1, and other pathologic ECM components (Fig. 1F).
Genetic studies have suggested a central role of epithelial cells in mediating IPF risk (4,6,32,33); thus, we next performed more focused analyses of epithelial cells. Overlaying canonical lineage markers onto unbiased clustering of 37,325 EPCAM + cells, we identified 12 distinct cell types/states (Fig. 2, A and B, and fig. S6). These included transcriptional programs consistent with classifically defined lung epithelial cells (basal cells, ciliated cells, and AT1 and AT2 cells) and also demonstrated four secretory cell subtypes, putative transitional states between secretory/ciliated and AT2/AT1 cells, a previously undescribed KRT5 − /KRT17 + population, and a cluster of proliferating cells. Among these four secretory cell types, two expressed classical secretory lineage marker SCGB1A1 and airways mucins (MUC5AC and/or MUC5B), while the other two were defined by the expression of SCGB1A1 and SCGB3A2 or SCGB3A2 only ( Fig. 2A and fig. S7). Quantification of cell types from transcriptomic data demonstrated significantly increased proportions of basal cells, SCGB3A2 + secretory cells, and KRT5 − /KRT17 + cells and significantly decreased proportions of AT1 cells and transitional AT2 cells (Fig. 2C, fig. S6, and table S4).
To understand the spatial distribution of these secretory lineages within PF and control lungs, we performed RNA in situ hybridization (ISH). In areas of cystic and fibrotic remodeling in PF lungs, we observed a notable pattern characterized by discrete regions of re-modeled epithelium with a near-exclusive population by a single secretory phenotype (Fig. 2, D and E). In control lungs, low-level MUC5B expression was observed in SCGB1A1 + cells in large airways; a subset of SCGB1A1 + cells coexpressed SCGB3A2. SCGB3A2 + cells lacking SCGB1A1 or MUC5B expression were found in a subset of airways in control lungs (Fig. 2F). Quantification of secretory cell subsets in matched formalin-fixed, paraffin-embedded tissue by RNA-ISH and automated image analysis revealed a significant increase in MUC5B + cells (Fig. 2H); SCGB1A1 + and SCGB3A2 + cells were observed in a patchy distribution and were found in highly variable density across different analyzed regions. In contrast, AT2 cells (quantified by SFTPC expression) (Fig. 2J) and AT1 cells (quantified by AGER expression) (Fig. 2K) were significantly less frequently found in PF lungs. These patterns are consistent with the epithelial proportions quantified using the scRNA-seq data (Fig. 2C) The proportions of secretory cell subtypes differed significantly between PF and control lungs, with a relative increase in MUC5B + cells and a decrease in SCGB1A1 + /SCGB3A2 + secretory cells (Fig. 2L). A smaller proportion of all Secretoglobin Family 3A Member 2 (SCGB3A2) cells coexpressed SCGB1A1 or MUC5B in PF lungs (Fig. 2M). Analysis of gene expression programs discriminating between the SCGB3A2 + /SCGB1A1 + and SCGB3A2 + populations demonstrated that the SCGB3A2 + cells expressed a subset of alveolar programs including NKX2-1, HOPX, CAV1, surfactant genes, and major histocompatibility complex-II ( fig. S7).
The role of AT2 cells as a progenitor population for AT1 cells during lung injury and repair has been well established (34,35). In both control and PF lungs, we identified a population of cells that express some features of both AT2 and AT1 cells (Fig. 2, A and B), resembling an indeterminate population reported in previous work (22). A pseudo-time analysis-used to identify trajectories between cell types/states-performed using the Slingshot (36,37) package in R suggested that this transitional AT2 cell population represents a state during the differentiation trajectory from AT2 to AT1 with increasing AT1 markers along the trajectory (see Materials and Methods and Fig. 3, A and B). Unexpectedly, among both control and PF lungs, a proportion of transitional AT2 cells express SCGB3A2 ( Fig. 2A and fig. S8), and trajectory analyses demonstrated SCGB3A2 + cells up-regulating AT1 and AT2 programs (Fig. 3, A and B), leading us to hypothesize that SCGB3A2 + cells are capable of generating AT1 cells by differentiating to transitional AT2 cells. These results are consistent with a recent report suggesting that both club cells and AT2 cells may generate AT1 cells via a common progenitor state in mice following experimental lung injury (38). To test the robustness of these trajectories and account for the fact that trajectory analyses are sensitive to choices made during dimensionality reduction, we analyzed the concordance of three different approaches [Slingshot based on Uniform Manifold Approximation and Projection (UMAP), Slingshot based on principal components analysis (PCA), and linear ordering based on PCA; fig. S9]. We found the ordering of cells between methods to be highly concordant ( fig. S9). Last, to confirm these states exist in vivo, we performed RNA-ISH using tissue sections from sequenced lungs to localize SCGB3A2, SFTPC (AT2 marker), and AGER (AT1 marker) and identified a putatively transitional state coexpressing SFTPC and AGER in both control (Fig. 3C) and fibrotic lungs (Fig. 3D). Consistent with transcriptomic data, a subset of these SFTPC + /AGER + dual-positive cells expressed low levels of SCGB3A2 in PF samples and were rarely observed in control lungs ( Fig. 3, C to F). Quantification of colocalization of SCGB3A2, SFTPC, and AGER demonstrated a larger proportion of SFTPC + cells as a proportion of all SCGB3A2 + cells and increased SCGB3A2 + cells as a proportion of SFTPC + cells (Fig. 3, E and F). Together, these data demonstrate that there is a spectrum of distal lung epithelial cells expressing a continuum of secretory, AT2, and AT1 programs and raises the possibility that multiple distal lineages may converge upon a common transitional phenotype that can give rise to AT1 cells during injury repair.
We also identified a previously undescribed epithelial cell population that expressed COL1A1 and other pathologic ECM components and was found nearly exclusively in PF lungs ( Fig. 2C and fig. S6). Furthermore, these cells were identified across all histopathologic  patterns of PF, although they were infrequently recovered from lungs with NSIP pathology (fig. S10 and table S2). Using RNA-ISH, we found that KRT5 − /KRT17 + cells intimately overlay foci of high collagen expression in the distal PF lung and coexpressed COL1A1 (Fig. 4, B to F); KRT17 + basal cells were observed in airways of control lungs but lacked coexpression of COL1A1 (Fig. 4G). Notably, KRT17 + /COL1A1 + cells could be identified in a lung biopsy specimen from an asymptomatic individual with a family history of IPF ( Fig. 4F) (27), was the predominant cell type-expressing integrin V6 (Fig. 4H), previously implicated in transforming growth factor  activation (40), and highly expressed CDKN2A (encoding for p16), suggesting cell cycle arrest and/or a senescent phenotype. KRT5 − /KRT17 + cells expressed canonical basal cell transcription factor TP63, but these cells lacked expression of the basal cell marker KRT5 (Fig. 4, I and J). Enrichment analyses of genes up-regulated in KRT5 − /KRT17 + compared to basal cells indicated increased expression of cell adhesion and motility programs, increased p53 and integrin signaling pathways, and increased expression of signaling targets downstream of the mitogen-activated protein kinase, glycogen synthase kinase 3, and protein kinase C- pathways (Fig. 4K). Given that these cells shared characteristics with airway basal cells and the distal/alveolar epithelium, the origin of these cells was unclear. To this end, we again performed Slingshot (see Materials and Methods) using only cells from PF lungs. The pseudo-time analysis suggested that the KRT5 − /KRT17 + cells may be derived from transitional AT2 cells, which themselves can arise from either AT2 or

of 15
SCGB3A2 + cells (Fig. 5, A and B). Note that these methods do not have sufficient resolution to assign individual cells to a specific trajectory but rather depict global patterns of changing gene expression that occur along an inferred differentiation path. Similar to our previous analysis, we tested the robustness of this trajectory by comparing three different approaches. Again, we found our cell order to be largely consistent between trajectory analysis techniques ( fig. S12), and the trajectory was recapitulated when using only cells from one processing site or a single individual ( fig. S13). We also compared the results to an alternative, recently developed method, TooManyCells, that is intentionally designed to be orthogonal to traditional pseudo-time analyses and does not rely on dimensionality reduction; this method also recapitulated our results ( fig. S13) S9). Both SOX4 and SOX9 are also associated with the trajectory (Fig. 5F) and are known to be key regulators of epithelial cells in the lung (42)(43)(44). Furthermore, we observed an increase in SOX4 expression in transitional AT2 cells, high expression of SOX4 in SCGB3A2 + cells, and high expression of SOX9 in KRT5 − /KRT17 + cells (Fig. 5F). We also observed a significant up-regulation of NR1D1 in AT2 and transitional AT2 cells from control lungs, while it was near undetectable in PF samples ( Fig. 5F and table S5). We next wanted to compare how these key transcription factors and some of the target genes differed along the KRT5 − /KRT17 + trajectory in PF samples compared to the AT1 trajectory in controls. To this end, we calculated a pseudo-time trajectory from AT2 cells to AT1 cells using only cells from control lungs ( Fig. 5F and fig. S16). The observed expression differences suggest that up-regulation of SOX4 and SOX9 combined with down-regulation of NR1D1 in AT2, SCGB3A2 + , and transitional AT2 cells may drive aberrant transcriptional programming, leading to the generation of ECM-producing KRT5 − /KRT17 + cells instead of the AT1 cells seen in a homeostatic lung. We also directly examined the relationship between expression of SOX4 and SOX9 and putative target genes in cells, where both genes were expressed in AT2, transitional AT2, or KRT5 − /KRT17 + cells (see Materials and Methods, Fig. 5G, and fig. S17). We identified significant associations between SOX4 and SOX9 and a number of genes with known function in PF, including the known senescence marker Cyclin Dependent Kinase Inhibitor 1A (CDKN1A) (Fig. 5G and fig. S17) (45). Last, we used RNA-ISH to demonstrate colocalization of KRT17 and SOX9 in the tissue from patients with PF (Fig. 5H). Together, these data provide direct evidence of an epithelial role in collagen/ECM production in PF lungs. This stands in contrast to genetic lineage-tracing studies performed in mouse models of PF (34,35). Furthermore, the cytokeratin profile of this KRT5 − /KRT17 + niche appears distinct from that characterizing "lineage-negative epithelial progenitor" cells previously described in mouse models (46,47). These discordant findings may reflect both biological differences in the distal lung epithelium between humans and mice and mechanistic differences between human PF and the bleomycin mouse model, as well as the inherent limitations of lineage-tracing studies using genetic Cre-lox systems and single marker-based approaches to mesenchymal cell annotation, which may differ between mice and human lungs (48). Focusing next on mesenchymal and other stromal cells (see Materials and Methods), we identified four discrete populations of fibroblasts in addition to smooth muscle cells and mesothelial cells (Fig. 6A). While the overall proportion of fibroblasts was higher in PF lungs compared to controls as anticipated (table S4), there was also a specific enrichment in fibrotic lungs of ACTA2 + myofibroblasts, a PLIN2 + lipofibroblast-like group, and a previously undescribed HAS1 hi fibroblast population ( fig. S18 and table S4). Each fibroblast subtype had a unique gene expression signature (Fig. 6, B and C), but there was also a collection of dysregulated genes shared across the fibroblast subpopulations (Fig. 1F). Further investigation revealed that the HAS1 hi population was composed entirely of cells from IPF lungs. Analysis of genes that were up-regulated in HAS1 hi cells indicated enrichment for pathways associated with cellular stress, interleukin-4 (IL4)/IL13 signaling (fig. S19), and programs previously implicated in epithelialmesenchymal transition in other systems (TWIST and SNAI1) ( Fig. 6C and table S5). To determine the localization of these different fibroblast populations, we performed multiplexed immunofluorescence/ ISH using specific markers enriched in each fibroblast population (Fig. 6, D to H). As anticipated, there was notable expansion of LUM + cells (a pan-fibroblast marker) in PF lungs compared to control (Fig. 6, D to F). -Smooth muscle actin (SMA) + /LUM + myofibroblasts were found coexpressing COL1A1 in subepithelial regions around large airways (Fig. 6, D and E) in PF and rarely around airways in control lungs (Fig. 6F). HAS1 + cells were found near exclusively in subpleural regions, where they colocalized with COL1A1 in IPF lungs (Fig. 6, G and H). In PF lungs, HAS1 + cells were found at increased depth from the pleural surface compared to control lungs, suggesting a possible invasive phenotype (Fig. 6I). PLIN2 + fibroblasts were found diffusely within interstitial regions and around alveoli (Fig. 6, J and K). Together, these data suggest that a global dysregulation of fibroblasts occurs in PF along with the expansion of specific pathogenic subtypes in spatially distinct niches driving pathologic ECM deposition and expansion (Fig. 6L).
To determine potential mediators of pathologic mesenchymalepithelial communication, we performed an interactome-based analysis identifying putative ligand-receptor (LR) binding pairs in epithelial and mesenchymal cells (see Materials and Methods) and constructed interaction networks for mesenchymal-driven (Fig. 6L) and epithelialdriven (Fig. 6M) signaling (fig. S20 and table S10). Both networks were significantly enriched for genes differentially expressed between PF and control lungs (see Materials and Methods). This analysis implicated matrix-driven signaling through integrin receptors as the central mechanisms through which fibroblast lineages interact with epithelial cells in PF lungs. In contrast, a more complex network involving multiple growth factors, cytokines, and chemokines signaling through integrins, Wingless (Wnt) co-receptors, and other pathways was observed for epithelial-driven signaling.
In addition, to assess the robustness of the previously unreported cell types identified in this study, we compared our results with two previously published studies that performed scRNA-seq on patients with IPF (21,24). We recapitulated the cell types identified by these two studies (table S11) and were able to recover the same cell types identified in this study in the data from the other groups (figs. S21 and S22). Last, the major findings of this study are consistent with those reported in a companion study from an independent group (49) using an entirely separate dataset and processing pipeline.

DISCUSSION
Although previous studies in experimental models have identified a number of pathways and mediators of lung fibrosis, the central mechanisms underlying pathologic epithelial remodeling and ECM expansion in the human lung have remained poorly understood. In this study, we provide a single-cell landscape of the PF lung, with a focus on epithelial and fibroblast cell populations, and identify several previously unrecognized cell types/states.
These granular data reveal a greater degree of transcriptomic heterogeneity within the lung epithelium than that recognized in the classical model of the lung epithelium, with considerable implications for understanding of lung disease mechanisms. First, these data reveal that the expression of the airway mucin MUC5B is not restricted to a SPDEF hi , "Goblet"-like cell (MUC5AC + cells), but rather is predominantly expressed in a distinct SCGB1A1 + secretory population, consistent with observations in an independent recent report focusing on the normal human lung (48). Further, we identify distinct genomic programs in SCGB1A1 + /SCGB3A2 + versus SCGB3A2 + secretory cells, the latter of which is characterized by relative up-regulation of alveolar programs. Histologic quantification identified that a proportion of SCGB3A2 + secretory cells coexpressed SFTPC (5.4% on control versus 13.7% in PF lungs). The convergence of alveolar and secretory programs is reminiscent of "bronchiolo-alveolar stem cells", which have been described in the murine lung (50)(51)(52)(53), but whose existence in the human lung has been uncertain. Rare identification of SCGB3A2 + / SFTPC + /AGER + cells raises the possibility that an analogous multipotent progenitor exists in the human lung; experimental validation of this hypothesis will be required.
In both control and PF lungs, we identified a transitional AT2 cell state characterized by simultaneous expression of lower levels of canonical AT2 programs compared to AT2 cells in addition to a subset of AT1 markers, including CAV1 and AGER. SFTPC + /AGER + cells were detected in both control and PF lungs by RNA-ISH, indicating that this state exists in vivo. This appears to share some degree of homology to transitional states recently described in the mouse lung (38,54). Further studies will be required to determine whether this reflects an intermediate state along the AT2 → AT1 differentiation axis, a stable and discrete epithelial niche, or a facultative state adopted by any AT2 cell in the proper context.
A putative role of lung epithelial cells producing pathologic ECM has been controversial, with ex vivo experiments, immunohistochemical studies, and some lineage-tracing studies in various experimental models supporting this hypothesis (55)(56)(57)(58)(59)(60), while other well-designed and detailed studies have not found evidence of epithelial-derived fibroblasts/myofibroblasts (34). The data we present provide some harmonization of prior seemingly discordant findings. Specifically, KRT5 − /KRT17 + cells retain canonical epithelial lineage-defining genes (EPCAM, NKX2-1, and cytokeratins) while coexpressing a small subset of "typical" mesenchymal genes, including COL1A1, FN1, and CDH2 (N-cadherin), while lacking generally recognized stromal surface markers and transcription factors. Notably, these data demonstrate that certain markers that have been used as lineage-defining in prior studies (HOPX for AT1 cells, VIM for mesenchymal cells, and S100A4 for mesenchymal/immune cells) lack a high degree of cell type specificity in the diseased human lung, potentially confounding interpretation of histology-based colocalization studies. Further studies will be required to determine the origin and function of KRT5 − /KRT17 + cells. Our computational trajectory analyses lead us to speculate that this phenotype may reflect aberrant differentiation of a transitional state on the AT1 differentiation axis, a hypothesis supported by recent lineage-tracing studies in mice (54); however, further studies are required to provide more direct evidence of this using human cells. Our trajectory analyses suggest SOX4 as a regulator of this trajectory. SOX4 has also been implicated in adoption of mesenchymal phenotypes in the context of lung cancer (61)(62)(63), raising the possibility that this cell state may also play a role in other lung diseases.
These data also provide new insights into the role of distinct fibroblast cell types/states in PF pathogenesis. Specifically, we identify multiple fibroblast phenotypes expressing pathologic ECM. These data suggest that while activated myofibroblasts produce high levels of collagen and other ECM components, other transcriptionally distinct fibroblast groups may also play a significant role. These distinct fibroblast subtypes localize to distinct regions of the lung; myofibroblasts in subepithelial regions around airways and areas of cystic remodeling, while HAS1 hi fibroblasts appear restricted to the immediate subpleural region; PLIN2 + and other LUM + fibroblasts are found diffusely in parenchymal regions. The stability of these distinct phenotypes over time in vivo cannot be assessed by this study; however, these data underscore the concept that a diversity of fibroblast phenotypes is found in PF lungs. Understanding the functional roles of these heterogenous populations may yield novel therapeutic targets for PF.
In addition to these notable findings, there are several limitations that warrant discussion. First, most control lungs (8 of 10) were current or former smokers and had undergone some period of mechanical ventilation. While this may be a particularly relevant control group given that tobacco use was similarly common among PF lung, it should be recognized that nonfibrotic control lungs declined for organ donation are unlikely to be "healthy" and changes in cell type-specific gene expression with regard to inflammatory pathways, for example, should be interpreted with caution. Second, PF samples were collected from lungs at the time of organ transplantation, reflecting advanced disease; whether similar changes in cell types present and gene expression programs would be found at early disease stages is not known. Third, these samples were collected exclusively from peripheral regions of the lung and cannot address a potential role of the more proximal airway epithelium, inflammatory cells, and resident stromal populations.
Together, our results provide substantial insight into the complexity, heterogeneity, and plasticity of the peripheral lung in human disease, building upon molecular atlasing efforts in the diseased (21,23,24) and healthy lung (48). This high-resolution overview identifies the genes, pathways, and programs that characterize pathologic lung remodeling in PF. Future studies investigating the origin, behavior, and function of these cell types, subtypes, states, and pathologic gene expression programs should provide additional insights into the foundational mechanisms regulating homeostasis and disease in the human lung.

MATERIALS AND METHODS
The code for genomic analyses in this paper is available at https:// github.com/tgen/banovichlab/.

Subjects and samples PF tissue samples were obtained from lungs removed at the time of lung transplantation at two lung transplant centers [Vanderbilt University
Medical Center (VUMC) and the National Thoracic Institute (NTI)]. Nonfibrotic control tissue samples were obtained from lungs declined for organ donation. For both PF and control lungs, tissue sections were taken from multiple peripheral (within ~2 cm of the pleural surface) regions in a given lung. For PF lungs, representatively diseased areas were selected on the basis of preoperative chest CT, while for control lungs, the most normal-appearing region was identified by gross inspection and selected for biopsy. For PF lungs, diagnoses were determined according to American Thoracic Society/European Respiratory Society consensus criteria (3). All studies were approved by the local Institutional Review Boards (Vanderbilt IRB nos. 060165 and 171657 and Western IRB no. 20181836). Surgical lung tissue sections were obtained from the ILD registry/repository at VUMC from patients with IPF who had undergone surgical lung biopsy for diagnosis of their disease (n = 3). Transbronchial lung biopsy specimens were obtained from siblings or children of patients with familial IPF with no known personal history of lung disease (n = 3). Details of this cohort have previously been reported (2). The selected slides were from individuals with mild interstitial lung abnormalities noted on chest CT ( fig. S11).

Tissue processing
Biopsies of multiple regions from each lung sample were digested in an enzymatic cocktail [collagenase I/dispase II (1 g/ml) tissue or Miltenyi Multi-Tissue Dissociation Kit] using a gentleMACS Octo Dissociator (Miltenyi Inc.). Adjacent tissue was fixed in 10% formalin for 24 to 72 hours and used for tissue localization studies. Tissue lysates were serially filtered through sterile gauze, 100-and 40-m sterile filters (Fischer). Single-cell suspensions then underwent cell sorting using serial columns (Miltenyi Microbeads, CD235a and CD45) at VUMC or fluorescence-activated cell sorting at the Translational Genomics Research Institute (TGen). CD45 − and C45 + populations were mixed 2:1 in samples processed at VUMC and were used as input for generation of scRNA-seq libraries. At TGen, calcein acetoxymethyl was used to stain live cells, and 10,000 to 15,000 total live cells were sorted directly into the 10x reaction buffer and transferred immediately to the 10x 5′ chip A (10x Genomics).

scRNA-seq library preparation and next-generation sequencing
scRNA-seq libraries were generated using the 10x Chromium platform 3′ v2 or 5′ library preparation kits (10x Genomics) following the manufacturer's recommendations and targeting 5000 to 10,000 cells per sample. From 12 of the samples, multiple 10x Genomics libraries were generated from distinct biopsies to account for regional heterogeneity (table S2). Next-generation sequencing was performed on an Illumina NovaSeq 6000 or HiSeq 4000. Reads with read quality less than 30 were filtered out, and CellRanger Count v3.0.2 (10x Genomics) was used to align reads onto GRCh38 reference genome. We have included QC (quality control) metrics from sequencing and demographic information of each library in table S2.

Image analysis and quantification
For cell quantification in tissue, RNA-ISH slides from control (n = 4) and IPF (n = 5) lungs were split into quadrants, and one image including an airway was taken in each quadrant. Images were captured at ×20 magnification using the Keyence BZ-X710 all-in-one fluorescence microscope and BZ-X Viewer software. Each quadrant image was composed of 5 × 5 overlapping images, which were stitched using Keyence BZ-X Analyzer software and then merged using HALO image analysis software v3.0.311.293 (Indica Labs), resulting in an area of at least 20 mm 2 analyzed for each slide. Cell number and target expression were automatically quantified using HALO v3.0.311.293 software (Indica Labs) and their FISH-IF v1.2.2 algorithm. Areas of high background were manually annotated and subsequently excluded by analysis algorithm. Data are presented as total cells per area analyzed or percentage of total cells.

Quantification of maximum subpleural depth of HAS1 Hi cells
Randomly selected tissue sections of diseased (n = 6) and control (n = 2) lungs underwent RNA-ISH for hyaluronan synthase 1 (HAS1) and COL1A1 as described above. The section was manually scanned at ×10 magnification for the expression of HAS1, and regions where there was HAS1 expression were imaged ×100 and ×400 magnification, with 400× images focused on capturing cells that were furthest from the pleura. All 100× images were then analyzed in ImageJ by cross-referencing to 400× images to identify the HAS1 + cells maximally distant from the pleura. The shortest straight-line distance to the pleura was then measured and recorded. This was repeated on each 10× image to obtain an average maximum depth from the pleura.

Dimensionality reduction, clustering, and visualization
Seurat v3 was used to perform dimensionality reduction, clustering, and visualization for the scRNA-seq data (3,4). Individual sample output files from CellRanger Count were read into Seurat v3 to generate a unique molecular identifier count matrix that was used to create a Seurat object containing a count matrix and analysis. All Seurat objects were combined into a merged dataset, and a percentage of mitochondrial genes were calculated for each sample in the merged object. Cells containing less than 1000 identified genes or more than 25% of reads arising from mitochondrial genes were removed ( fig. S1). SCTransform with default parameters was used to normalize and scale the data, and dimensionality reduction was performed using PCA on the top 3000 most variable genes. We also performed SCTransform with the batch_vars parameters to correct for batch effects between different processing sites (TGen versus Vanderbilt) and different sequencing batches. However, applying batch correction did not change our clusters markedly ( fig. S2). To determine the optimal number of principal components (PCs) for UMAP visualizations and to avoid overfitting, we attempted to identify an optimal number of PCs that kept the relative distance between points on UMAP-1 and UMAP-2 stable, i.e., the UMAP plot remained stable between across a range of PCs used to generate the UMAP plot. To this end, we used Mantel randtest to calculate correlation between relative location of points on the UMAP plot between two adjacent numbers of included PCs for both UMAP-1 and UMAP-2. Correlation coefficient values were plotted, and we manually selected PC ranges where both UMAP-1 and UMAP-2 plateaued. The smallest PC number within the plateau was chosen for input into the final UMAP (figs. S1 and S3).

Cell type annotation and doublet removal
Markers specific for major cell types PTPRC + (immune cells), EPCAM + (epithelial cells), PECAM1 + /PTPRC − (endothelial cells), and PTPRC − /EPCAM − /PECAM1 − (mesenchymal cells) were used to split Seurat clusters into four subgroups ( fig. S3). Each subgroup object underwent the same dimensionality reduction, clustering, and visualization approach as described above. Each subgroup object was then further split into clusters and manually annotated with known cell type markers (table S3). Doublet cells were identified manually as expressing markers for different cell types, and the final object was created by merging all four annotated, doubletremoved subgroups.

Differential expression analysis
To identify differentially expressed genes between cell types, we used a negative binomial model as implemented in the Seurat FindMarkers function, comparing each individual cell type to all other cells within the major cell type cluster (immune, epithelial, endothelial, and mesenchymal). Genes were considered differentially expressed if the adjusted P value was lower than 0.01 (table S5). To identify genes that were differentially expressed between PF and control lungs, we took each cell type independently, for all cell types with a minimum of 50 cells in both PF and control lungs, and used the negative binomial model to test for differences in expression with a log fold change (FC) cutoff of 0.25. Genes were considered differentially expressed if the adjusted P value was less than 0.01 ( Fig. 1D and table S5). To test for robustness of the differentially expressed analysis and assess for batch effects, we applied latent.vars function embedded in Seurat FindMarkers to assign processing site, flow cell, or processing site and flow cell as latent variables. However, we did not find significant differences between the different processing sites or flow cells being used as latent variables, compared to the original analysis where no latent variables were specified (table S6).

Cell trajectory analysis
The R package Slingshot (37) was used to perform a pseudo-timebased cell trajectory analysis. The slingshot wrapper function was performed with the UMAP dimensionality reduction and cluster labels as in Seurat objects to identify the trajectory. The slingshot trajectories are calculated such that each trajectory has a single start and end point. For clarity and brevity, multiple trajectories were plotted on the same graph (Figs. 3A and 5, A and B). Robustness of the trajectories was assessed by comparing the UMAP-derived pseudo-time to a Slingshot pseudo-time calculated on the first two PCs and linear ordering of cells based on the first or second PC only (figs. S9 and S12).

Identifying genes associated with trajectory analysis
To identify genes along the course of the trajectory, a general additive model (GAM) was used to regress each gene on the pseudo-time variable. The top 400 significant genes, after correcting for multiple testing [false discovery rate (FDR), <0.01], were chosen for heatmap based on loess GAM P values, and heatmap was plotted using the function plotHeatmap in the R package clusterExperiment (64). Loess plot for individual genes of interest was generated using ggplot2 geom_smooth function (smoothing method: "loess") on the trajectory variable timeline resulted from the slingshot wrapper function above.

TFBS analysis
For motif analysis, findMotifs.pl (package HOMER v4.10) (65) was used on promoter sequences [700-base pair (bp) upstream and 100-bp downstream of transcription start-site (TSS)] of genes in the trajectory heatmap (Fig. 5, D and E, and table S9) using HOMER pre_built human promoter database and default parameters. To identify SOX and NR1D1 motif location, HOMER annotatePeaks. pl was performed on promoter regions of the targeted genes using the motif files generated from findMotifs.pl function.

Correlation between SOX genes and putative targets
We assessed the correlation of putative targets with SOX4 and SOX9 in transitional AT2 and KRT5 − /KRT17 + cells. Within each individual, the mean of each gene (after removing zeros) was calculated in each cell type. The strength of the correlation was assessed across all individuals with detectable expression of both genes using a linear regression, and those with a nominal P value of less than 0.05 were plotted color-coded by PF and control ( Fig. 5G and fig. S17).

Pathway enrichment analysis
Pathway enrichment was performed using PANTHER Pathways through the WebGestalt 2019 (66). Differentially expressed genes (FDR, <0.05) with an absolute increase in proportion of cells expressing a given gene of >0.1 and a logFC of >0.4 (HAS1 hi versus other mesenchymal cells) were selected as input for enrichment analysis. For analysis of KRT5 − /KRT17 + cells versus basal cells, all genes up-regulated in KRT5 − /KRT17 + cells (logFC, >0.5; FDR, <0.1; n = 227 genes) were selected as input for enrichment analysis. Gene Ontology biological process enrichment, PANTHER Pathways, and Kinase target enrichment analyses were performed using WebGestalt 2019.

Interactome analysis
The scaled gene expression matrix for diseased cells from a Seurat subgroup object containing epithelial cells and mesenchymal cells (excluding HAS1 hi cells, which localized to subpleural regions and were not in spatial proximity with nonmesothelial epithelial populations) was filtered to contain the top 20% highly expressed genes for each cell type using the iTALK function rawParse (67). These top genes were used as input for the iTALK FindLR function to find LR pairs between all cell types. Hundreds of LR pairs were found between epithelial cells and mesenchymal cells (table S10). Next, we tested whether highly coexpressed LR pairs (top 20%) were enriched for differentially expressed genes between PF and control samples. First, we generated an empirical null distribution for each cell type interaction [excluding Perilipin 2 (PLIN2 + ) fibroblasts and KRT5 − / KRT17 + epithelial cells, which did not contain sufficient numbers of control cells for differential expression testing] by randomly sampling a matched number (table S10) of LR pairs from all LR pairs identified for that interaction from the iTALK database 1000 times and identifying the number of significantly differentially expressed genes within each iteration. P values were calculated as the proportion of permutations that exceeded the number of differentially expressed genes within the top 20% most highly coexpressed pairs. For highly coexpressed ligands from mesenchymal cells and corresponding receptors in epithelial cells, significant enrichment (P < 0.05) was found for LR pairs between fibroblasts and each epithelial cell type, as well as between myofibroblasts and each epithelial cell type, excluding ciliated cells (P = 0.083). For highly coexpressed ligands from epithelial cells and corresponding receptors in mesenchymal cells, significant enrichment was found for LR pairs between fibroblasts and each epithelial cell type excluding ciliated (P = 0.071), differentiating ciliated (P = 0.53), and proliferating epithelial (P = 0.168) cells, as well as between myofibroblasts and each epithelial cell type, excluding ciliated (P = 0.068) and differentiating ciliated (P = 0.179) cells. For Fig. 6 (M and N), the top five most highly coexpressed pairs (ranked by the product of mean ligand expression and mean receptor expression) for each epithelial cell-mesenchymal cell interaction were visualized using Cytoscape (68). Differentially expressed genes within the visualized LR pairs are denoted by red text.

Data comparison with current publications
To assess the robustness of the novel cell types identified in this study, we obtained data from two previously published studies that performed scRNA-seq on patients with IPF. These data were obtained from Gene Expression Omnibus (GEO) accessions GSE122960 (24) and GSE128033 (21). The sample data from each study were then processed in Seurat following the same QC procedure as the data presented in this study. After QC and normalization, individual samples [17 samples for Reyfman et al. (24) and 18 samples for Morse et al. (21)] were combined together into their own respective Seurat objects per study. Seurat's Transfer Anchors function was then independently applied to each combined unlabeled Seurat object using the dataset from this manuscript as a reference for cell type prediction and annotation. Cell type proportion (table S11) was calculated per cell type per study, excluding HAS1 high cell type since TransferAnchor was unable to detect this cell type in Morse et al. (21) and Reyfman et al. (24) dataset. For the novel cell types KRT5 − /KRT17 + , SCG3A2 + , transitional AT2, and PLIN2 + , we generated heatmaps of the 20 most discriminating markers (from the differential expression analysis described above), along with key marker genes, across the cells in our study and the transfer anchors annotated cells from Reyfman et al. (24) and Morse et al. (21). These results show the high degree of similarity between the cells across studies (figs. S21 and S22).

Data and code accessibility
Raw and processed 10x Genomics data can be found on GEO using accession number GSE135893. The code used to analyze the data can be found at https://github.com/tgen/banovichlab/.

SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/6/28/eaba1972/DC1 View/request a protocol for this paper from Bio-protocol.