Altered Heterogeneity of Ageing Lung Endothelium is a Hallmark of Idiopathic Pulmonary Fibrosis

Integrated Single Cell Map of Ageing Human Lungs Displays Endothelial Cell Contribution in Disease

To characterise ageing human lung endothelium in IPF and to account for its intra- and inter-lung heterogeneity within and between individuals [29], we performed comparative analysis of scRNAseq data from lung tissue samples selected only from age-matched (63± 14 years old) transplant donors and recipients with IPF. Specifically, we selected data from only those studies (cohorts) that contained samples for both conditions (transplant donors and IPF patients). From all publicly available scRNAseq datasets of human IPF lungs (4 in total, to our knowledge; [25, 27, 28, 30, 31]), only two, representing independent cohorts of subjects/individuals of different ethnic backgrounds, with controlled sample locations, fulfilled these selection criteria [25, 27] and therefore were used for generating a single dataset of ageing human lungs for integrative comparative scRNAseq analysis (Table 1; Figure 1A; for full details on dataset selection see Additional file 1: Detailed methods / Dataset Selection).

• Download figure
• Open in new tab

Figure 1. Cell type-specific composition of ageing human lung tissues from donors and IPF patients from two independent cohorts. (A)

Schematic of methodology for samples integration and comprehensive data analysis. Selected data from two studies were grouped into two cohorts (18 samples in total; see Table 1; [25, 27]) and subjected to quality control (Additional file 2: Figures S1, 2B) and merging into an integrated Seurat V4 object for subsequent analysis. The resolution was determined using the clustree R package, with 0.5 value used for this analysis and throughput the study to avoid over-clustering. Endothelial cell (EC) clusters were identified using pan-endothelial and lymphatic-specific markers. 3192 blood vessel endothelial cells (BEC) and 659 lymphatic vessel endothelial cells (LEC) were identified, and sub-clustered for their characterization and comparison in donor and IPF samples. Full details can be found in Additional file 1: Detailed methods. (B) UMAP representation of all cells and cell clusters from all 18 samples (pooled). Clusters were numbered, annotated and labelled according to their detected gene markers. Specific details about each cell population can be found in Additional file 2: Figure S2 and Additional file 3, Table S1. Two clusters (LEC and BEC) are labelled with a red ellipse. (C) Heatmap of top ten differentially expressed genes by cells in the cluster. Each column represents the average expression for a cell, hierarchically grouped by cell type. Gene expression values are scaled from 2 (yellow) to -2 (purple) across rows within each cell type. (D) Violin plots and density plots representation of expression data for pan-EC and LEC-specific markers expression in all identified clusters. Density plot expression values are scaled from between 0.04 and 1 (yellow – high expression) to 0 (purple – low expression). (E) Stacked bar chart of percentage contribution of each cluster to total lung population, split by cohort and condition - donor or fibrosis (IPF). Colour keys for cluster identity is the same as in B. Percent contributions of 3,192 BEC and 659 LEC to the total cell population in ageing human lung are shown with arrows. Specific details for each cell population can be found in the Additional file 3: Table S3. (F) UMAP representation of all cells and clusters from all 18 samples (pooled), split by condition (donor or fibrosis) and labelled by cluster (top panel), or split by cohort and labelled by condition (blue – donor; red – fibrosis; bottom panel). (G) Box and whiskers plots of EC percentage all cell types in the lung (left) or blood or lymphatic endothelial cells (BEC and LEC respectively) from total EC per sample in individual lung samples percentage from all other cell types in individual lung samples (right), split by condition (blue – donor; red – IPF/fibrosis). The lines within the boxes denotes the medians, the ‘box’ contains the 25th to 75th percentiles and the ‘whiskers’ mark the range of the data. Statistical analysis was done using Shapiro-Wilk test (p=0.162 for donor and p=0.367 for fibrosis) and unpaired t-test (p = 0.139). Means ± standard deviation (2.501±1.407 for donor and 3.959±2.388 for fibrosis). UMAP - Uniform Manifold Approximation and Projection.

The integration of the two datasets resulted in generation of a multi cohort single cell map revealing 25 distinct populations (Figure 1A-B, Additional file 2: Figure S1), based on specific gene expression patterns (Figures 1C, Additional file 2: Figure S2, Additional file 3: Table S1) that include known cell type-specific transcriptomic markers [22, 25–27]. Identities of two closely related blood and lymphatic vessel EC (BEC and LEC respectively) populations (Figure 1B) were confirmed by detailed analysis of expression of a panel of pan- EC and lymphatic EC-specific markers (Figure 1D).

Quantitative analysis of cell numbers for each of 25 populations in age-matched groups revealed IPF-associated changes in ageing human lungs in both cohorts (Figure 1E). The analysis reflected differences in location of tissue samples (Table 1), matching contribution of different cell types to the total human lung cell population and alterations in IPF lung reported in age non-matching studies [21, 22, 25–27, 30, 32] (Additional file 3: Table S2). Also, it uncovered statistically insignificant changes in proportions of ageing human lung BEC and LEC populations between cohorts and conditions (Figure 1E, Additional file 3: Table S3), which is in agreement with similar analysis in age non-matching study of IPF lungs [25]. Data integration revealed comparable spatial profiles of ageing human lung EC populations in two conditions and two cohorts (Figure 1F) and their similar contribution in donors and IPF patients (Figure 1G).

Blood Vessel Endothelium in Ageing Human Lung Exhibits Disease-Specific Heterogeneity and Alterations in Cell Numbers of Individual Subpopulations

Unsupervised clustering of ageing human lung BEC alone (n=3192) identified 12 distinct subpopulations (Figures 1A, 2A) based on their differentially expressed genes (DEG; Figures 2B, Additional file 2: Figure S3). Expression analysis of pan-EC markers and LEC-specific genes confirmed their identities as BEC (Figure 2C). The analysis of top 10 DEG for each subpopulation (110 in total) determined their selection for inclusion (or not) in 10 signatures of ageing human lung BEC subpopulations (Figures 2B, D, Additional file 2: Figure S3, S4; Additional file 3: Figure Table S4). 49 selected subpopulation-specific DEG contribute to distinct transcriptional profiles for subpopulations 1-4 and 6-11 (Figure 2D). Subpopulations zero and 5 lack obvious markers or signatures (Figures 2D, Additional file 2: Figure S5).

• Download figure
• Open in new tab

Figure 2. Cell subtype-specific EC subpopulation composition of ageing human lung blood vessel endothelial cell cluster from donors and IPF patients from two independent cohorts.

To comprehensively decipher transcriptional signatures of ageing human lung blood vessel endothelial cell (BEC) in IPF, we subtracted scRNAseq data for BEC (n=3192) from the total lung cell populations from integrated dataset (two cohorts and both conditions, see Figure 1A), and sub-clustered them. (A) UMAP representation of all cells and cell subpopulations from BEC cluster (from all 18 pooled samples) which are presented in the Figure 1. Subpopulations were labelled 0-11 according to their signatures (Additional file 2: Figure S3). Note that, although subpopulations 1, 7 and 8 visually present with two potential subgroups each on this UMAP, each of two subgroups are not significantly different from each other in the described clustering parameters (Additional file 1: Detailed methods). In subpopulation 1, this is likely to be due to the limitations of visualisation of 2D UMAP compared to 3D UMAP, which shows as an individual cluster (see also Figure 2E, Additional file 4: Video S1). In subpopulation 7 and 8, 3D UMAP reveals the shift of one subgroup towards subpopulation zero (see also Figure 2E; Additional file 4: Video S1). (B) Heatmap of top 10 differentially expressed genes (DEG) by subpopulation. Each column represents the average expression for a cell, hierarchically grouped by condition (donor or fibrosis) and cell type. Gene expression values are scaled from 2 (yellow) to -2 (purple) across rows within each cell type. Specific details about each cell subpopulation can be found in Additional file 2: Figures S3 - 5, and Additional file 3: Table S4. All cell clusters can be separated with high confidence (Additional file 3: Table S5). (C, D) The analysis of the top 10 DEG for each subpopulation (120 in total) was performed alongside tests for expression of pan-EC and LEC markers and using violin plots to test their subpopulation-specific expression for determining their selection for inclusion (or not) in 10 signatures of ageing lung BEC subpopulations. (C) Violin plots of expression of pan-endothelial and LEC markers in 12 identified BEC subpopulations. (D) Dot plots of 49 genes identified as differentially expressed by the heatmap in B, plotted alongside pan-EC and LEC markers and control genes for comparison, and used in annotation of 12 identified BEC subpopulations (see also annotation summary presented in Figure 4). Note that 39 (marked in red) of 49 identified subpopulation-specific DEG have not been previously assigned as markers of individual generic subpopulations of ageing human lung BEC. (E) Representative lineage map in a form of a three-dimensional (3D) UMAP generated by performing trajectory analysis using the R package Slingshot. Plotted lines are trajectories representing predicted relationships between 12 identified BEC subpopulations. 3D video can be found Additional file 4: Video S1. (F) Table of module scores (Additional file 1: Detailed methods) grouped by references from which reference signatures were derived. Six studies from 2013-2021 include: (i) Gillich et al., 2019 (ii) Vanlandewijck et al., 2019; (iii) Iso et al., 2013; (iv) Nukala et al., 2021; (v) Ochiya et al., 2014; (vi) Herwig et al., 2016. Abbreviations: aCAP - aerocyte capillary; gCAP - general capillary; large – large vessel; art- artery; CAP1 - capillary 1; CAP2 - capillary 2; bronc – bronchial; S vein - systemic vein; P vein - pulmonary vein; Acap - aerocyte capillary; Gcap - general capillary; PeriB – Peribronchial (Note that names given by the authors are preserved). Coral red represents a subpopulation with positive score (above 0) for signature, colour intensity is proportional to the score positivity. Additional information is available in Additional file 1: Detailed methods, Annotation of BEC subpopulation. UMAP - Uniform Manifold Approximation and Projection.

Expression analysis of the 49 subpopulation-specific DEG (Figure 2D) in the total lung revealed that the majority of them were not EC-specific, i.e. also expressed in other cell types (Additional file 2: Figure S4A). 39 of them have not been previously assigned as markers of 8 generic subpopulations of BEC identified and characterised earlier in human or mouse lung [21, 22, 25, 26, 33], and thus are newly associated with ageing human lung endothelium (Figure 2D; Additional file 2: Figure S4A, B). A pathway analysis uncovered their belonging to specific classes of proteins, including secreted molecules, cell surface receptors, transcription factors and sub-cellular distribution (Additional file 2: Figure S4B-C).

Three-dimensional pseudo-time lineage analysis revealed close relationships between identified BEC subpopulations (Figure 2E; Additional file 4: Video S1). This included a direct link between subpopulations 5 and 3, and manifestation of subpopulation zero as a “convergence point” for seven subpopulations - 1, 2, 4, 6, 9, 10, 11 and, partially, 8. The expression analysis of markers of EC sub-types proposed in other transcriptional studies in human and murine lungs [21, 26, 33–38] was applied to 12 ageing human lung BEC subpopulations identified in our study, and revealed only few distinct “marker / ageing human lung BEC subpopulation” matches (Additional file 2: Figures S6, S7). Further exploration of identities of 12 subpopulations by using these markers and also those proposed for specific EC subpopulations in HLCA and IPF cell atlas studies [21, 22, 24–26], when coupled with module scoring (Figure 2F) and cell cycle (Additional file 2: Figure S8) analyses, confirmed these findings, whilst 4 subpopulations (zero, 5, 10 and 11) remained un-matched (Figures 2F). Finally, quantitative analysis revealed a significant increase in representation of two subpopulations (1 and 2) and overall decrease in subpopulation 7 in IPF (Figure 3A-C; Additional file 3: Tables S5, S6), with specificity by the cohort (Figure 3D).

• Download figure
• Open in new tab

Figure 3. Percentage contributions of blood endothelial cell subpopulations in donors and IPF patients.

Quantitative analysis of percentage contribution of cells from each subpopulation to total blood endothelial cell (BEC) population in ageing human lung from healthy donors and IPF patients was conducted using Seurat v4. (A) UMAP of BEC split by sample condition (donor or fibrosis). (B) Table detailing the number of cells per cluster, split by condition. (C-D) Stacked bar charts of percentage contribution of each subpopulation to total endothelial cell population split by condition for (C) total integrated dataset and (D) cohort 1 (left) and cohort 2 (right). Cell numbers from each individual sample were analysed using Multiple t test. P < 0.05 (*) was considered significant. UMAP - Uniform Manifold Approximation and Projection.

Annotation of Blood Vessel Endothelial Cell Subpopulations and Identification of Dedifferentiated Cells in Ageing Human Lung

Altogether, our findings formed the foundation for the annotation of 12 ageing human lung BEC subpopulations with distinct transcriptional signatures (Figure 4A; in-depth description is provided in Additional file 1: Detailed methods, Annotation of BEC subpopulations). Integrative analysis revealed that ageing human lung BEC belong to three and 8 subpopulations of bronchial and pulmonary circulation networks respectively (with one, inflammatory, being “intermediary” as also identified in other reports [24]), and include two previously undescribed subpopulations - designated as “de-differentiated” endothelium - in pulmonary circulation (Figures 4A, 2F, Additional file 2: Figure S9). Ten subpopulations of ageing human lung BEC have distinct transcriptional signatures, whilst no obvious DEG could be identified for subpopulations zero and 5 (Figures 2B, 2D, 4A, Additional file 2: Figures S5, S11; in-depth description is provided in Additional file 1: Detailed methods, Annotation of BEC subpopulations).

• Download figure
• Open in new tab

Figure 4. Annotation of subpopulations and relationship between them in donors and IPF patients.

The profiling of markers of EC sub-types (Figure 2F, Additional file 2: Figures S6, S7, S9), when combined with the expression of pan-EC and LEC markers (Figure 2C-D) and sub-population-specific signatures identified in this study (based on expression pattern of 49 DEG; Figures 2D, Additional file 2: Figure S3), lung tissue sample collection locations (Table 1), modules scores of association with bronchial or pulmonary circulatory networks (Figure 2F, Additional file 2: Figure S9), EC differentiation status (Additional file 2: Figure S12), IPF-related alterations in cell numbers (Figure 3), high expression of IL6, SELE and ICAM1 [21, 24]; Figure 2D), cell cycle analysis data (Additional file 2: Figure S8) and findings from three-dimensional pseudo-time lineage analysis (Figure 2E, Additional file 4: Video S1), formed the foundation for the annotation of 12 identified ageing human lung BEC subpopulations with distinct transcriptional signatures (Note that in-depth description is provided in Additional file 1: Detailed methods, Annotation of BEC subpopulations). (A) Table summarizing the information taken into account for naming/annotating subpopulations, including gene signatures, sample location by cohort, blood vessel circulatory networks (bronchial and pulmonary) and cell numbers (by condition - donor or IPF/fibrosis). ^a Key colour code for subpopulations is the same as in Figures 2 and 3. ^b Sample location (by cohort) information corresponds to Table 1. DLP = distal lung parenchyma, LABS = Longitudinal apical to basal segments. ^c Cell number dot plots corresponds to data in Figure 3. Signature genes include previously reported (blue; n=10) and de novo identified in our study genes (red; n=39). (B) Dendrogram (cluster tree) and dot plot of subpopulations, based on unsupervised hierarchical clustering of 12 subpopulations of ageing human lung BEC from integrated total dataset, which was split into two separate objects by condition (donor on the left and fibrosis on the right). Stratification into two groups termed “capillary” (or “microvascular”) and “macrovascular” is based on distinct clustering of capillary sub-clusters (3, 5, 7 and 8).

Subpopulations zero, 5, 10 and 11 had similar expression levels and ridgeplot profiles (histograms) to all other BEC subpopulations for the housekeeping gene beta-2 macroglobulin (B2M; Figure 2D; Additional file 2: Figure S12A). These subpopulations also exhibited low, but similar in ridgeplot profiles to other BEC subpopulations, expression of pan-EC genes (Figure 2C-D, Additional file 2: Figure S15A, B). Subpopulations zero and 5 contained a lower average number of genes per cell compared to other BEC subpopulations (Additional file 2: Figure S12B) and did not express markers of 23 non-EC cell clusters of human lung (Additional file 3: Table S1; Additional file 2: Figure S12C) or EC progenitors (Additional file 2: Figure S12D).

Next, we used these novel identified signatures (Figure 4A) to test the hypothesis that BEC subpopulation zero (the largest one in donor group in integrated dataset; Figure 3C) is exclusively associated with ageing in the human lung. The comparison of scRNAseq datasets of four ageing and three non-ageing/young human lungs against these signatures led to annotation of BEC subpopulations in both age groups (Additional file 2: Figures S13, S14). All signatures, except for subpopulations 1 and 2 (bronchial venules IPF Endo 1 and 2, which are noticeably associated mainly with ageing IPF lung in cohorts 2 and 1 respectively; Figure 3C-D), were matching specific subpopulations in both age groups (as expected, since these were donor lung samples, Additional file 2: Figure S14E). Cell number analysis revealed that BEC subpopulation zero (identified as subpopulation A in ageing vs non-ageing/young lung analysis), which has no obvious DEG, was present within ageing human lung only (Additional file 2: Figures S13B-C, S14E), supporting our hypothesis and its annotation as an ageing-related “de-differentiated” BEC (Figure 4A).

Finally, unsupervised hierarchical clustering of 12 subpopulations of ageing human lung BEC revealed two distinct groups - “capillary/microvascular” and “macrovascular” (Figure 4B). In IPF, the relationships between subpopulations were altered, with “bronchial venule/IPF Endo 1” and “intralobular venule” subpopulations “translocating” to “capillary/ microvascular” (Figure 4B).

Cell Numbers for Bronchial Blood Vessel Endothelial Cell Subpopulations in Ageing Distal Lung Depend on Fibrosis Degree and Bronchi Number in IPF

To investigate the reason for differential cell numbers observed for two bronchial BEC subpopulations (2 and 7; bronchial venule and capillary BEC respectively) in IPF (Figure 3C, D), we focused on analysing scRNAseq data from the distal part of human ageing lungs (cohort 1; four donors and four IPF patients; Figure 5 A, B). This analysis exposed considerable variation in a percentage of subpopulation 2 cells between individual samples in IPF (5.48 - 91.6 %), suggesting inter-lung sample heterogeneity. The comparison of scRNAseq data with IPF lungs histology (based on haematoxylin and eosin staining) for each individual case, which has been made available by the authors in the supplementary materials section of the original study [27], revealed that the lowest percentage of subpopulation 2 cells was detected in sample IPF04, where there was a noticeable lack of bronchial structures. Based on these findings, we hypothesised that in IPF lung, significantly altered heterogeneity (in terms of differential cell numbers) of ageing bronchial (and also possibly pulmonary and lymphatic) endothelium reflects association with histological changes, i.e. tissue composition, and, in particular, with the increase in percentage of (i) areas occupied by peri-bronchial spaces undergoing fibrotic alterations and/or (ii) a number of bronchi, in particular, which both underpin inter- and intra-lung heterogeneity in IPF.

• Download figure
• Open in new tab

Figure 5. Spatial validation of bronchial blood endothelial subpopulations in distal lung by immunohistochemistry and quantitative high content image analysis.

(A-B) Dotplots of quantitative analysis of percentage contribution of cells from bronchial subpopulation 2 and 7 to total blood endothelial cell (BEC) population in ageing human lung from healthy donors (n=4) and IPF patients (n=4) from cohort 1 (Table 1) was conducted using Seurat v4. Donor samples are shown in blue, fibrosis samples in coral red. (C) Masson’s trichrome and (C, F) immunohistochemical staining for calcitonin receptor like receptor (CLR) was conducted on formalin fixed paraffin embedded (FFPE) tissue sections of distal lung surgical biopsy from a 58-year-old male patient with IPF (with confirmed histopathological pattern of usual interstitial pneumonia) using our anti-human CLR antibody LN1436 [43]). (C) Images of scanned full sections (scalebars represent 5000µm) and two enlarged regions of interest (ROI) 1 and 2 (see full set in Additional file 2: Figure S16; scalebars represent 100µm). Yellow dotted line represents peri-bronchial areas (between bronchi and alveolar regions) that were analysed. Black dotted line represents example measurement of width of peri-bronchial area (see Additional file 1: Detailed methods). Dot plots of quantitative analysis of (D) peri- bronchial area width and (E) number of CLR-positive nuclei from 9 analysed ROI. (D-E) Low degree fibrosis areas are shown in blue, high degree fibrosis areas in red. A cut off of 90µm for average peri-bronchial area width was used to determine high (> 90) and lower (< 90) degree fibrosis. Shapiro Wilcoxon normality test and Mann-Whitney U test were used to analyse the data. P values are indicated at the top. Dotplots presenting: (D) quantification for peri-bronchial areas width in µm and (E) a number of CLR positive nuclei per mm² (venule nuclei, left; capillary nuclei, right). (F) Two sections (S1 and S2) from the same distal lung biopsy were used for quantitative high content Image J analysis. Individual bronchi were identified within the total section under microscope and selected for downstream analysis (Additional file 1: Detailed methods). Scale bars represent 5000 µm. Indicated with red box and enlarged to the right of each scanned section are two representative examples of bronchi in individual sections (see full set in Additional file 2: Figure S17). Barcharts representing the results of quantitative image analysis of: (G) total section areas in mm², (H) number of bronchi per section and (I) percentage of CLR-positive nuclei in total section; all done for S1 and S2 sections. Colours for bronchial venules and capillaries match colours of ageing human lung BEC subpopulations 2 and 7 presented throughout paper (see Figure 4A as an example).

To test this hypothesis, we analysed a distal lung surgical diagnostic lung biopsy from a 58-year-old male patient with IPF, confirmed by a histopathological pattern of UIP. First, we identified 8 BEC subpopulations, including 2 and 7, by spatial position/localisation in ageing human lung tissue of EC by immunostaining using two pan-endothelial (for various human tissues, including lung) markers - platelet endothelial cell adhesion molecule (PECAM-1), also known as cluster of differentiation 31 (CD31; [41]) and calcitonin receptor-like receptor (CLR, encoded by CALCRL gene; [42–44]; see also Figure 1D) (Additional file 2: Figure S15). Next, image analysis allowed quantification of IPF-associated changes in cell numbers for subpopulations 2 and 7 (Figure 5C-I). The comparative analysis of low vs high degree fibrotic areas within peri-bronchial regions in distal lung parenchyma identified no differences in the number of CLR-positive nuclei in venules (vessel diameter >7µm) and a significant decrease in capillaries (vessel diameter <7µm; Figure 5C-E; Additional file 2: Figure S16), reflecting BEC from subpopulations 2 and 7 respectively, between groups.

Furthermore, comparative analysis of sections from two different areas of a distal lung parenchyma from the same IPF patient revealed the association of percentage of CLR- positive nuclei for both capillaries and venules with the number of bronchi (Figure 5F-I; Additional file 2: Figure S17).

Altogether, these findings from immunohistochemistry corroborated the scRNAseq analysis of IPF-associated changes in cell numbers of bronchial BEC subpopulations 2 and 7 in distal lung parenchyma in four IPF patients from cohort 1. Importantly, it revealed the association of changes in cell numbers of ageing human BEC subpopulations 2 (venules) and 7 (capillaries) in peri-bronchial areas with the fibrosis degree and bronchi number.

scRNAseq and Immunohistochemistry Reveal Fibrotic Alveolar Regions as Potential Niches for Dedifferentiated Blood Vessel Endothelial Cell Subpopulations in IPF

To investigate the possible location of BEC subpopulation zero in IPF, we focused on analysing scRNAseq data from the distal part of human ageing lungs (cohort 1; four donors and four IPF patients; Figure 6A). This exposed a noticeable variation in a percentage of cells between individual samples in IPF (3.79 - 73 %; Figure 6A), suggesting their inter-lung heterogeneity (similar to findings for BEC subpopulation 2; Figure 5B), and revealing a negative correlation between subpopulations zero and two (Figure 6B). Interestingly, the highest percentage of subpopulation zero cells was detected in the sample IPF04 (Figure 6C), which according to histology of individual IPF lungs (based on haematoxylin and eosin staining) for each individual case, which has been made available by the authors in the supplementary materials section of the original study ([27]), was associated with the lack of bronchial structures. Based on these findings and association of subpopulation zero with pulmonary circulation (4A; based on negative score for bronchial circulation signature, Figure 2F), we hypothesised that BEC subpopulation zero resides mainly in the alveolar regions or “niches” of the distal part of human ageing lungs in IPF, where it is likely to derive from BEC subpopulations 3, 4, 6 or 8 (Figures 2A, E, 3D and 4A).

• Download figure
• Open in new tab

Figure 6. Comparative analysis of numbers of blood endothelial subpopulations in alveolar regions of distal lung in IPF by scRNAseq, immunohistochemistry and quantitative high content imaging. (A)

Dotplot, (B) scatterplot and (C-D) stacked bar charts of quantitative analysis of percentage contribution of cells from subpopulations 0, 2, 3, 4, 6 and 8 to total blood endothelial cell (BEC) population in ageing human lung from healthy donors (n=4) and IPF patients (n=4) from cohort 1 (Table 1) was conducted using Seurat v4. (A-B) Donor samples are shown in blue, fibrosis samples - in coral red. (B) Pearson correlation analysis for percentage contribution of subpopulations zero vs subpopulation two in fibrosis samples. Line of best fit shown by grey dashed line. R - Pearson correlation coefficient. (C) Contribution of subpopulation zero and 2 to total BEC per sample coloured by subpopulation in fibrosis samples. Colours of ageing human lung BEC subpopulations match those presented throughout paper (see Figure 4A as an example). (D) Stacked bar chart showing percentage contribution of BEC subpopulation zero and combined alveolar BEC subpopulations 3, 4, 6 and 8, when analysed together as a “single alveolar pool”. Subpopulation zero is coloured by coral red, Subpopulations 3, 4, 6, 8 are coloured white. Data includes both donor and fibrosis samples. (E) Paired ridgeplots/histograms and dotplots of PECAM1 expression in donor (left) and fibrosis (right) samples from cohort 1 (Table 1) was conducted using Seurat v4. Top row, (coral ridgeplot) - subpopulation zero Middle row, (light grey ridgeplot) - subpopulations 3, 4, 6 and 8 combined, with insert reflecting data for individual subpopulations; with colours of ageing human lung BEC subpopulations matching those presented throughout paper (see Figure 4A as an example). Bottom row, (grey ridgeplot) - subpopulations zero, 3, 4, 6 and 8 combined into a “single alveolar pool”. Note that the ridgeplot and dotplot for this pool in IPF patients (right) representing the PECAM1 mRNA expression profile in alveolar region in the lung tissue, with increased number of CLR-negative BEC when compared to donors (left). (F-H) Immunohistochemical staining for CD31 was conducted on formalin fixed paraffin embedded (FFPE) tissue sections of distal lung surgical biopsy from a 58-year old male patient with IPF (with confirmed histopathological pattern of usual interstitial pneumonia) using mouse monoclonal antibody clone JC/70A. (F) Two representative regions of interest (ROI; 1 and 2) of alveolar regions of distal lung that were used for quantification are presented (see full set in Additional file 2: Figure S18). Left, low degree fibrosis, right, high degree fibrosis. Scalebar represents 100µm. Black dashed line represents example measurement of width of alveolar area analysed (for full details, see Additional file 1: Detailed methods). Dot plots of quantitative analysis of (G) alveolar area thickness in µm and (H) proportion of CD31 positive nuclei per mm2 (left) and average DAB signal intensity per CD31 positive nuclei (right). Nonfibrotic alveoli samples are shown in blue and fibrotic alveoli samples - in coral red. Nine ROIs were analysed for both nonfibrotic and fibrotic alveoli. Statistical analysis was done using Shapiro Wilcoxon test and Mann-Whitney U test. P values are indicated at the top. Note that immunohistochemistry and quantitative high content imaging profiles (H) confirm those from scRNAseq data analysis (E).

To test this hypothesis, scRNAseq data analysis was focused on distal ageing human lung samples (cohort 1), revealing that the increase in percentage of subpopulation zero cells in IPF lungs was associated with a decreased in percentage of alveolar subpopulations 3, 4, 6 or 8 BEC, when they are analysed together, as a “single alveolar pool” in ageing human lung (Figure 6D). Interestingly, subpopulation zero has expression profiles for housekeeping gene B2M (Figure 1B, Additional file 2: Figure S12A) and pan-EC markers (PECAM1/CD31 and CALCLR/CLR) (Additional file 2: Figures S15A, B) comparable to all other ageing human lung BEC subpopulations, whilst having a considerably larger proportion of cells with lower expression of the latter two genes. Furthermore, quantitative analysis of PECAM1 expression in donor and fibrosis samples confirmed similarities in ridgeplot profiles between subpopulation zero and alveolar BEC subpopulations 3, 4, 6 and 8 (in terms of presence of both PECAM1 positive and negative cells, as well as a trend for increased number of PECAM1 negative cells in fibrosis; Figure 6E, top and middle panels). Integrative analysis of alveolar BEC subpopulations (3, 4, 6 and 8) and subpopulation zero taken together, as a “single alveolar pool” in ageing human lung, revealed the trend for a reduction of PECAM1 expression (reflected by the decrease in percentage of positive cells) in fibrosis, also linked to an increased variability between individual samples (Figure 6E, bottom panel). Immunohistochemistry and image analysis of CD31 expression in distal alveoli in IPF lung revealed a statistically significant decrease in proportion of CD31-positive nuclei and CD31 signal intensity in fibrotic compared to non-fibrotic alveoli (Figure 6 F-H; Additional file 2: Figure S18). Combining analysis of PECAM1/CD31 expression by scRNAseq and immunohistochemistry revealed similar decreases in expression (compare last graphs in Figures 6E and H), thus supporting our hypothesis that BEC subpopulation zero resides mainly in the alveolar regions/niches of the distal part of human ageing lungs in IPF.

Expression of Genes Linked to Disease-Specific Signalling Pathways and Endothelial Cell Processes is Altered in Ageing Lung Blood Vessel Endothelium in IPF

Comparison of gene expression between donor and IPF lungs led to identification of 3596 significantly regulated (1600 down- and 1996 up-) DEG across 12 BEC subpopulations (Figures 7A-B, Additional file 2: Figure S19A; Additional file 3: Table S7). There were significant associations with specific signalling pathways in 7 (two bronchial and 5 pulmonary, including two “de-differentiated”) from 12 subpopulations (Figure 7C). DEG in subpopulations zero, 8 and 11 were directly linked to IPF-associated pathways, and in four others (1, 3, 5 and 11) - with fibrotic or inflammatory pathways (Figures 7C, Additional file 2: Figure S19B). Next, the utilisation of published gene set enrichment analysis (GSEA) libraries for the assessment of expression of marker genes which are associated with ten EC- relevant processes revealed positive module scores in 5-12 of the BEC subpopulations for individual processes on average (Additional file 2: Figure S20) and significant changes in several (3 bronchial and 6 pulmonary, including two “de-differentiated”) BEC subpopulations in IPF (Figure 7D-E). This reflected subpopulation-specific pattern of endothelium response to the fibrotic lung environment in all processes, ranging from two (for senescence) to 8 (for vasodilation) subpopulations amongst 9 analysed (Figures 7E).

• Download figure
• Open in new tab

Figure 7. Association of differentially expressed genes in ageing human lung blood vessel endothelial cell with specific signalling pathways and endothelial cell-relevant processes in IPF. (A)

An example of volcano plot of differentially expressed genes (DEG) (donor vs IPF/fibrosis) in blood endothelial cell (BEC) subpopulation zero. Volcano plots for other subpopulations can be found in the Additional file 2: Figure S19A. Statistical analysis was done using Shapiro Wilcoxon test and Wilcoxon rank sum test, with P value <0.05 (dotted line indicates cut-off) considered significant. (B) Doughnut plot of DEG between donor and IPF/fibrosis for each subpopulation. Outer ring shows up-regulated genes and inner ring represents down-regulated genes, with table in the centre showing numbers of DEG per subpopulation. Clusters 9 and 10 were excluded from the analysis, as they contained too few cells to analyse between conditions. (C) Bar chart detailing the results of the Ingenuity Pathway Analysis (IPA) using Core expression analysis function [77] of statistically significant altered pathways using DEG from B. Blue colour represents a negative Z score (down), and orange represents a positive Z score (up). DEG from three subpopulations (4, 6 and 7) did not significantly associate with any pathways and from two others (9 and 10) had too few cells in the IPF group for analysis. IPA was performed as described in Additional file 1: Detailed methods. (D, E) Key colour code for subpopulations is the same as in panel C and in Figures 2-4. (D) Histograms (ridgeplot) of endothelial cells differentiation (n= 5), mesenchymal transition (n= 12), inflammation (n= 567), and vasodilation (n= 36) scores in all blood endothelial subpopulations in donors (donor) and IPF (fibrosis). Ridgeplots for senescence (n= 79), apoptosis (n= 161), proliferation (n= 54), migration (n= 175), angiogenesis (n= 48), and permeability (n= 40) can be found in Additional file 2: Figure S20. Number of genes in sets are indicated in brackets. Black dotted line is added to annotate zero (cut-off point for differences). (E) Table detailing the P values for difference in module scores between donor and fibrosis. Blue colour represents down- regulation and orange - up-regulation in IPF. Statistical analysis was done using Shapiro Wilcoxon test and Mann-Whitney U test. UMAP - Uniform Manifold Approximation and Projection.

Integrated Single Cell Map and Immunohistochemistry of Ageing Human Lung Lymphatic Endothelium Reveal Its Heterogeneity and Contribution in IPF

Unsupervised clustering of subtracted ageing human lung LEC alone in integrated dataset (689 cells) identified two main populations consisting of five distinct, and to our knowledge previously uncharacterised, subpopulations in total (Figures 1A, 8A), as based on their DEG patterns (top 10 DEG for each, 50 in total, with 20 not associated with LEC formerly [22, 26] (Figures 8B-D; Additional file 2: Figure S21, 22F) and cell cycle analysis (Additional file 2: Figure S23). All signatures were specific in recognising individual LEC subpopulations (Additional file 2: Figure S21A-B). Intriguingly, ageing human lung LEC subpopulations 1 and 2 were present only in cohort 1 and subpopulations zero, three and 4 – exclusively in cohort 2 (Additional file 2: Figure S24). Quantification of cell numbers revealed a trend for a decrease (for subpopulations 1 and 2) but overall a lack of non-significant changes for all subpopulations between donor and IPF lungs in integrated (Additional file 2: Figure S24A-B; Additional file 3: Table S9) or percentage contribution when split by cohort (Additional file 2: Figure 24C-E; Additional file 3: Table S9) datasets, in contrast to cell number data on BEC (Figure 3C, D). Altogether, our findings formed the foundation for detailed transcriptional characterisation of five distinct (informed by DEG, our transcriptional signatures, sample collection site and cell cycle information) LEC subpopulations in ageing human lung (Figure 8D; in-depth description is provided in Additional file 1: Detailed methods, Annotation of LEC Subpopulations).

• Download figure
• Open in new tab

Figure 8. Cell subtype-specific composition of ageing human lung lymphatic endothelial cell cluster from donors and IPF patients from two independent cohorts.

To comprehensively decipher transcriptional signatures of ageing human lung lymphatic vessel endothelium in IPF, we subtracted scRNAseq data for ageing human lung lymphatic vessel endothelial cell (LEC) (n=659) from the total lung cell populations from integrated dataset (two cohorts and both conditions, see Figure 1A), and sub-clustered them. (A) Top, UMAP representation of all cells from total LEC cluster/population (from all 18 pooled samples; Figure 1) with pseudotime lineage analysis of 5 resulting subclusters. Bottom, UMAP of LEC clusters coloured by cohort. Clusters were labelled 0-4 according to their signatures, which are presented as (B) heatmap of top 5 differentially expressed genes by cluster. Full details about each cell sub-population can be found in Additional file 2: Figure S21-23. (C) Violin plots demonstrating expression of pan-EC markers and LEC-specific genes (PROX1, PDPN, LYVE-1 and FLT4) and CCL21 [78] was done to confirm the identities of all cell subpopulations as LEC. (D) Table summarizing the information for each subpopulation, including gene signature (including previously reported, in blue, and de novo identified in our study, in red, genes), sample location by cohort and cell numbers. ^a Key colour code for subpopulations is the same as in Figure 8A. ^b Sample location (by cohort) information corresponds to Table 1. DLP = distal lung parenchyma, LABS = Longitudinal apical to basal segments. ^c Cell number dot plots correspond to data in Additional file 2: Figure S24 A, B. (D) Dotplot of quantification for percentage contribution of subpopulation 1 and 2 to total LEC per sample. Left, subpopulation 1, right, subpopulation 2. Donor samples are shown in blue, fibrosis samples - in coral red. (F-H) Immunohistochemical staining for PDPN was conducted on formalin fixed paraffin embedded (FFPE) tissue sections of distal lung surgical biopsy from a 58-year-old male patient with IPF (with confirmed histopathological pattern of usual interstitial pneumonia) using mouse monoclonal antibody clone D2-40. (F) Image of scanned full section (scalebars represent 5000µm), with six regions of interest (ROI) 1-6 labelled with red boxes and shown enlarged in (G) (scalebars represent 100µm). Table insert to the right is the summary of information for each ROI on whether observed structure could be identified in original publications for cohorts 1 and 2 (* - [25, 27]), i.e. DLP and LABS areas of the lung respectively (see D), and which LEC subpopulations were detected in our analysis (see A). (G) |Representative images for each ROI (magnified from F; scalebars represent 100µm), including inserts (red dotted boxes) of the areas containing PDPN-positive lymphatic vessels (bottom right; scalebars represent 100µm). Coloured squares in the top left corner correspond to mapping conducted using table presented in (F). Arrows - lymphatic vessel; A - artery; V – vein; Arrowhead - small vessel; * - honeycomb cysts; B – bronchioles; P - pleura. (H) Dot plots of quantitative analysis of percentage of PDPN-positive nuclei from all identifiable nuclei per mm2 within individual ROI (n=3 per ROI), grouped by location within distal lung diagnostic biopsy of IPF (for full details, see Additional file 1: Detailed methods; images of 21 ROI used for high content image analysis are presented in Additional file 2: Figure S25). (I) Doughnut plot of DEG between donor and IPF/fibrosis for each subpopulation. Ring shows up-regulated genes, with table in the centre showing numbers of DEG per subpopulation (see also Additional file 3: Table S10). Note that no downregulated genes were found for LEC subpopulations. (J) Table detailing the P values for difference in module scores between donor and fibrosis. Blue colour represents down-regulation and orange – up-regulation in IPF. Statistical analysis was done using Shapiro Wilcoxon test and Mann-Whitney U test. UMAP - Uniform Manifold Approximation and Projection.

To investigate the reason for non-differential cell numbers observed for ageing human LEC in IPF (Additional file 2: Figure S24), compared to differential cell numbers observed for two bronchial BEC subpopulations (2 and 7), we focused scRNAseq data analysis on two LEC subpopulations (1 and 2) detected in the distal part of human ageing lungs from four donors and four IPF patients (cohort 1; Figure 8E). This exposed very small variation in a percentage of subpopulation 2 (PDPN-positive) cells between individual samples in IPF (0.71 - 5.20%; Figure 8E), contrasting findings for BEC subpopulations zero, two and 7 (Figures 5B and 6A, C). Based on this observation and on histology from 4 previous reports, jointly covering 5 different structures in ageing lung in donors or IPF patients [9, 11, 26, 46], we hypothesised that PDPN-positive LEC are distributed uniformly in all tissue regions (bronchial and alveolar) within ageing human distal lung parenchyma in IPF.

For testing this hypothesis, we immunostained human IPF lung tissue using anti-PDPN antibody D2-40 (Figure 8 F, G). The quantification analysis of immunostaining in this distal lung parenchyma sample uncovered signs of developing fibrosis and substantial intra-lung heterogeneity throughout the section areas/regions and revealed the distribution of PDPN- positive vessels (Additional file 2: Figure S25). Interestingly, the comparison to inter-lung heterogeneity, which was detected across four samples within cohort 1 (as demonstrated by histological analysis using haematoxylin and eosin staining; [27]) and cohort 2 (as demonstrated by immunohistochemistry images presented by the authors; [25]) detected the presence of all identifiable structures from both cohorts (non-fibrotic and fibrotic alveolae, bronchiole and honeycomb cyst in cohort 1; and large artery and veins, bronchiole and pleura in cohort 2) within the surgical distal lung biopsy. It also revealed uniform distribution of PDPN-positive LEC nuclei in all detected tissue regions (bronchial, alveolar and pleural) within ageing human distal lung parenchyma and honeycomb cyst regions (Figure 8 F-G; Additional file 2: Figure S25). Together with finding from scRNAseq data analysis, which showed that in ageing human lung PDPN-positive LEC subpopulation 2 was present only in cohort 1, whilst PDPN-positive LEC subpopulations zero and 4 - exclusively in cohort 2 (Figure 8C), these data suggest that our immunostaining detected all three subpopulations of PDPN-positive LEC subpopulations (zero, two and 4) in ageing human lung in IPF within our single diagnostic biopsy (Figure 8 F-G).

Next, high content image quantification of cell numbers/percentage of cells revealed uniform percentages of PDPN-positive LEC across analysed bronchial, alveolar and subpleural regions, including both fibrotic and non-fibrotic regions (ranging from 4.5% to 20.9%; Figure 8 H; Additional file 2; Figure S25). These findings support scRNAseq analysis data, which revealed lack of significant differences for these subpopulations in IPF lungs in the integrated dataset (Additional file 2: Figure S24A-B; Additional file 3: Table S9), thus supporting our hypothesis.

Finally, gene expression comparison between donor and IPF lungs led to identification of 101 significantly up-regulated DEG across 5 LEC subpopulations (Figures 8I), compared to 3596 DEG for ageing human ling BEC (Figure 7B) and significant linkage to IPF-associated pathways in subpopulations zero and one (Additional file 2: Figure S26). Finally, the utilisation of GSEA libraries identified significant alterations in ageing human lung LEC, in IPF (Figure 8J; Additional file 2: Figure S26 C-L).