Spatial single-cell atlas reveals regional variations in healthy and diseased human lung

Abstract

Integration of scRNA-seq data from millions of cells revealed a high diversity of cell types in the healthy and diseased human lung. In a large and complex organ, constantly exposed to external agents, it is crucial to understand the influence of lung tissue topography or external factors on gene expression variability within cell types. Here, we apply three spatial transcriptomics approaches, to: (i) localize the majority of lung cell types, including rare epithelial cells within the tissue topography, (ii) describe consistent anatomical and regional gene expression variability within and across cell types, and (iii) reveal distinct cellular neighborhoods in specific anatomical regions and examine gene expression variations in them. We thus provide a spatially resolved tissue reference atlas in three representative regions of the healthy human lung. We further demonstrate its utility by defining previously unknown imbalances of epithelial cell type compositions in chronic obstructive pulmonary disease lungs. Our topographic atlas enables a precise description of characteristic regional cellular responses upon experimental perturbations or during disease progression.

Fig. 1. HybISS-based cell type map reveals cell type distribution and neighborhoods.

A Experimental outline including the location of sample collection from donated healthy human lungs, and the methods used for the mRNA-based cell type mapping. 1—trachea, 2—proximal bronchi, 3a—bottom part of upper left lobe, 3b—top part of upper left lobe, 3c—bottom part of lower left lobe. HybISS Hybridization-based in situ sequencing, SCRINSHOT Single-Cell Resolution IN Situ Hybridization On Tissues, RRST RNA-Rescue Spatial Transcriptomics. B UMAP of cells after leiden clustering profiled with HybISS, colored by the assigned cell type (35 cell types presented). Gen general, adv adventitial, nan not annotated. C Heatmap of relative abundance of clustered cell types between locations demonstrating their frequency across the profiled regions. Similarity in cell type composition between three distal lung regions can be assessed by hierarchical clustering dendrogram on the left side. D Representative histological images from four biological replicates of an analysed trachea (donor 4) and a distal lung (donor 1) biopsies with hematoxyllin and eosin (H&E) staining (up) coupled to the maps of cell types identified by HybISS on top of nuclei (DAPI, white) in the same sections (down). Spots represent detected transcripts, colored according to the corresponding cell type of the cell they were assigned to. Colors as in Fig. 1B. Dashed lines indicate the approximate borders of histologic compartments. SMG submucosal gland, aw airway, alv alveolar region, bv blood vessel. Scale bar 200 µm. E Cell type neighborhood enrichment graph representing cell types as nodes, and edges indicating a positive neighborhood enrichment (>2) between cell types across the profiled sections. Suggested neighborhoods are shown as bubbles. Node colors as in Fig. 1B. Source data are provided as a Source Data file.

Fig. 2. Distribution of non-epithelial cell subtypes in histological tissue compartments.

A Heatmap of average (mean) cell type proportions (%) in manually annotated histological compartments using HybISS cell type maps. Proportion calculations include non-annotated cells, not presented in the heatmap. Cell type/subtype maps on top of nuclei (DAPI, white) of HybISS (left) and SCRINSHOT (right) datasets of donor 1 trachea (B) and distal lung (C), respresentative images from four biological replicates. Dotted outlines: peribronchial and airway (AW), submucosal gland (SMG), alveolar (Alv) and arterial (Art) compartments. Cell type abbreviations: pb peribronchial, adv adventitial, alv alveolar, pv perivascular, cap capillary, gen general, art arterial, ven venous, nan not annotated. Scale bar 200 µm. D Graph of mean cell type/subtype proportions in each of the peri-epithelial compartments. Floating bars indicate Min/Max values as bounds of box, with a line on mean proportion per compartment, and indicated in percent (%), n = 3 in tracheal regions, n = 4 in proximal and distal lung regions. Values were compared in lung regions using Friedman two-sided test followed by Dunn’s multiple comparisons test. Significant differences (P < 0.05) are indicated by asterisk (*). Adjusted P values are as follows: P = 0.0370 for fibroblasts (RGCC) between alveolar parenchyma and proximal bronchi; P = 0.0370 for smooth muscle cells between alveolar parenchyma and proximal bronchi; P = 0.0370 for endothelial (CLDN5) cells between alveolar parenchyma and proximal SMG. The cell subtypes that dominated in one anatomic location are highlighted with coloured boxes. E Schematic summary of the gene expression in the parenchyma of each compartment. Source data are provided as a Source Data file.

Fig. 3. Distribution of cells and gene expression in bronchial epithelium along proximo-distal and apical-basal axes.

A Cell type frequencies in the airway epithelium according to HybISS. Left: stacked bar plot representing relative frequencies of the airway epithelial cell types across regions from different donors, using samples with representative numbers of airway cells. Right: cell type maps in the indicated regions from donor 1, respresentative images from four biological replicates. Spots represent detected transcripts, colored according to the corresponding cell type of the cell they were assigned to. Colors as in Fig. 1B. Nuclei: gray. Lu lumen, dashed line: approximate location of basal membrane. Scale bar: 50 µm. B Heatmap of the relative mean gene expression in airway epithelium of variable epithelial markers across regions (colors) in three analysed donors (numbers). The expression is normalized by gene, dividing by sum of values in each row. Superscript numbers: references of previous studies, reporting variable expression of the corresponding marker along the proximal-distal axis. Asterisk*: statistically significant expression differences of the corresponding marker between regions (P < 0.05, repeated measures ANOVA with Geisser-Greenhouse correction, followed by Tukey’s multiple comparisons test, all 161 detected genes tested). Plus + : having highest mean change. C Maps of epithelial cell types detected by SCRINSHOT in the indicated regions of the airways, respresentative images from four biological replicates. Arrows: cell clusters. Inserts in respiratory bronchiole map: (a-b) SCRINSHOT images of representative AT0 cells with either SCGB3A2 (a, orange squares) or SCGB3A1 (b, jade squares) dominating expression. (c) Zoomed area of respiratory bronchiole. Nuclei: gray. Scale bar in maps: 50 µm. Scale bar in SCRINSHOT images 10 µm. aw airway. D Area plots representing the relative apical-basal cell type distribution across regions according to HybISS (data from one representative donor 3). X axis: relative distance of cells from the basal membrane. Y axis: relative frequency of cell types. E Mean gene expression of significantly layer-variable markers along the apical-basal axis of the tracheal epithelium (n = 4). Individual data points represent biological replicates. Box bounds indicate standard deviation, line indicates mean, and whiskers indicate minimum and maximum values. Significant differences (P < 0.05) between gene dot distances are calculated using nested one-way ANOVA followed by Tukey’s multiple comparisons test, and indicated with lines and brackets. Exact P values are indicated in Source Data. Blue box* indicates the genes that are located differently from both basal and luminal values. To the right, image of HybISS-detected transcripts in tracheal epithelium, where each transcript is shown as a characteristic coloured shape. The image is representative of four biological replicates. Squared area in the image is shown magnified below. Scale bar 20 µm. Source data are provided as a Source Data file.

Fig. 4. Rare cell types and their distribution in the airways.

A Maps of rare cell types detected with SCRINSHOT from three anatomical regions (donor 4), respresentative images from three biological replicates. Scale bar 200 µm. Lu lumen. Raw SCRINSHOT signal images of cells labelled with letters (a-f) are shown in Supplementary Fig. 12B. B Heatmap of gene expression demonstrating unique and overlapping marker genes within the detected rare cell types from four donors. At least two out of four donors demonstrated each cell type in each anatomical region, except NE-GHRL population. Number of cells quantified (180): ionocytes—41, tuft-like:—27 (including tuft—24 and rare tuft-like—3), NE-GHRL—9, NE-PCSK1N—32, NE-ASCL1—40, NE-GRP—31. C Immunofluorescent staining for GHRL (green) and GRP (magenta) of two subtypes of neuroendocrine cells, as well as epithelial membrane marker CDH1 (white) on top of nuclei (DAPI, blue), distal lung of donor 4, which had largest number of neuroendocrine cells from three stained samples. Areas (a) and (b) are crops from a larger image in Supplementary Fig 12D. Scale bar: 10 µm. D Bar plot of the number of the detected cells per total airway epithelial cell number (shown in %) from four donors, error bars: standard deviation. Individual data points represent biological replicates (only donors with airway larger than 2 mm of basal membrane length). Arrow: NE-GHRLpos cells appearing only in distal lung. Source data are provided as a Source Data file.

Fig. 5. Cell type and neighborhood changes in COPD.

A Box plot of the cell type numbers in healthy and COPD (mean ± standard error, with individual values). Significant differences are highlighted with asterisk (**, adjusted P = 0.0087), according to unpaired parametric two-sided multiple t-test (20) of logit-transformed data with Holm-Sidak correction. Individual data points demonstrate samples from different donors and patients, n = 3. The direct cell type comparison without the correction also revealed significant differences between healthy and COPD samples in the following cell types: AT0 (P = 0.0004), AT1 (P = 0.0127), AT1-AT2 (P = 0.0090), basal cells (P = 0.0295), T lymphocytes (P = 0.0088), labeled by asterisks at the cell type name. B Representative (from three biological replicates) maps of alveolar epithelial cell types detected in healthy and COPD samples, nuclei: gray. Scale bar 200 µm. C UMAP plots of analyzed cells labeled according to their cellular neighborhoods (left) and their corresponding condition (right). Neighborhood annotations: Cap-Alv capillary-enriched alveoli, AT2-Alv AT2-enriched alveoli, AT0-Alv AT0-enriched alveoli, C-Epi club-enriched epithelium, G-Epi goblet-enriched epithelium, T-Epi secretory TRB-enriched epithelium, B-Epi basal cell-enriched epithelium, Imm-Par immune cell-enriched parenchyma, Str-Par stromal cell-enriched parenchyma, SM smooth muscle, SMG submucosal gland, SMG-serous serous cell-enriched submucosal gland, TRB terminal respiratory bronchiole. Arrows indicate the clusters that are predominantly composed of COPD-derived cells. D Heatmap exploring the cell type composition of each neighborhood, with cells (vertical lines) grouped by their assigned neighborhood cluster (x-axis) represented in Fig. 5C (left). Bar color represents the ratio of neighborhood enrichment with each cell type for each cell. E Maps of cell type neighborhoods in healthy and COPD samples in alveolar region with large and small airway and respiratory bronchioles. Color code as in C (left). Scale bar 200 µm. Simple arrows: AT0-enriched alveolar neighborhood, dashed arrows – immune-enriched parenchyma, double-line arrows – secretory TRB-enriched epithelium. Source data are provided as a Source Data file.