Heterogeneous microenvironment analysis to explore the potential regulatory role of endothelial-mesenchymal transition in idiopathic pulmonary fibrosis.

Background: Idiopathic pulmonary fibrosis (IPF) is a chronic and progressive interstitial lung disease mainly caused by excessive proliferation of fibroblasts and activation of myofibroblasts. The cellular microenvironment is mainly composed of different types of cellular components and extracellular matrix (ECM), whose changes directly affect cellular heterogeneity, resulting in immensely complex cellular interactions. However, microenvironment study is mainly focused on the pathological process of tumors, and the microenvironment changes during IPF development remain unclear.
Methods: The current study intends to employ IPF-related single-cell sequencing and gene expression profile data to analyze the scores of different cell clusters in the IPF microenvironment, and exploit the underlying interaction between cells to illustrate the fundamental mechanism causing IPF.
Results: Our analysis revealed that the amount of endothelial cells was obviously decreased, and the amount of fibroblasts and myofibroblasts was increased during the development of IPF, suggesting a possible endothelial-mesenchymal transition (EndMT) process. Furthermore, we found that the hub genes obtained through IPF-related gene expression profile analysis may play a regulative role in the number and function of endothelial cells and fibroblasts/myofibroblasts during IPF. Conclusions: Our research represents a valuable analysis of the cellular microenvironment, and provides a novel mechanistic insight into the pathobiology of not only EndMT in IPF, but also other traumatic fibrotic disease disorders. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Idiopathic pulmonary fibrosis (IPF) is a chronic, irreversible, and fatal lung disease, which is characterized by fibroblast proliferation and excessive deposition of extracellular matrix (ECM) in lung tissue (1). With the development of the disease, lung function gradually declines, leading to respiratory failure and even death. The median survival period after initial diagnosis is 2–3 years (2,3). Clinically, lung transplantation is an effective therapeutic method to improve the quality of life and survival rate of patients with IPF (4). However, given the short supply of donor organs and some serious immune rejection responses after transplantations, only a small fraction of the patients can benefit from lung transplantation. In addition, nintedanib and pirfenidone have been approved to treat IPF by the US Food and Drug Administration (5). However, the therapeutic effect of these drugs is limited, and is accompanied by multiple adverse reactions (6,7). Therefore, there is a pressing need to better understand the underlying pathogenesis and molecular mechanism of IPF, in order to explore new treatment means (8).

The tumor microenvironment (TME) is a novel foundation inside in the bodies of tumor patients. The structure of a typical TME is composed of fibroblasts, myofibroblasts, endothelial cells, pericytes, adipocytes, immune and inflammatory cells, and ECM, which can directly affect the growth, migration, and differentiation of cancer cells (9,10). Normally, the stroma maintains a stable environment within tissue and acts as a barrier to tumor formation. However, when cells become cancerous, the surrounding matrix changes to support cancer development. Likewise, the development of fibrotic disease is a result of the interaction of multiple cell subsets with distinct genetic and phenotypic characteristics (11). Changes in cellular diversity contribute to pathological pulmonary fibrosis (12). In the TME, among fibroblasts, pericytes, and bone marrow-derived mesenchymal cells, endothelial cells adopt cancer-associated fibroblast (CAF) phenotypes through the endothelial-mesenchymal transition (EndMT) process (13). It has been reported that mesenchymal cells, immune cells, and endothelial cells synergistically regulate the production, degradation, deposition, and remodeling of the ECM, leading to the development of liver fibrosis. Interactions between fibrotic components create a unique ‘fibrotic niche’ microenvironment (14). A previous study has shown that M1 bone marrow-derived macrophages (BMDMs) recruit and modify the activation of endogenous macrophages and natural killer (NK) cells by regulating the liver microenvironment, possibly leading to hepatic stellate cell (HSC) apoptosis, thereby hindering liver fibrosis (15). Furthermore, it has also been reported that the transformation of bone marrow-derived monocytes to fibroblasts is a key step in the pathogenesis of renal fibrosis and is regulated by the inflammatory microenvironment. The axis of cluster of differentiation eight positive (CD8+) T cell and interferon (IFN)-γ-CD4+ T cell as an important microenvironment for the monocyte-to-fibroblast transition, which negatively regulates renal fibrosis (16). However, the direct link between microenvironments and IPF has not yet been confirmed. Thus, there is a pressing need to explore the connections between different cell types in the IPF microenvironment to better understand the underlying pathogenesis and molecular mechanism of IPF.

Transcriptomics is increasingly used to explore disease mechanisms. Through transcriptomics and next-generation sequencing, we can obtain an unprecedented level of detail in understanding cell phenotype and function by investigating genes expressed in specific physiological and pathological states (17). The increased availability of well characterized human tissues and the emergence of high-throughput transcriptome analysis technology have promoted a new era of IPF research. Transcriptome research has revealed many new molecules and pathways highly related to the pathogenesis of IPF, including the role of matrix metalloproteinase (MMP), developmental pathway, micro ribonucleic acid (RNA) and the importance of alveolar epithelial and myofibroblast regulatory networks in IPF. At the same time, using single cell and tissue microenvironment to improve the cellular and spatial resolution of transcriptomics is very important to study the occurrence of IPF (18). In this study, we explored the composition of the IPF microenvironment and the potential interactions of its subcellular populations by analyzing IPF single-cell sequencing data and transcriptome gene expression data.

In the past, there have been many studies on IPF, but the research on pulmonary microenvironment is still rare. However, the existing bioinformatics articles simply explore the pathogenesis of IPF through the transcriptome characteristics of large tissue homogenate, failing to capture the complexity of microenvironment. Our research starts from the heterogeneous microenvironment of IPF, and uses a unique perspective to understand how cells interact in the remodeled heterogeneous microenvironment of IPF through the changes of cell number and phenotype. We aim to discuss the pathogenic mechanism of IPF from a novel perspective, and to lay a theoretical foundation and provide data support for drug development and the discovery of pathological triggers. We present the following article in accordance with the MDAR reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-22-1438/rc).

Methods

IPF related data set acquisition and pre-processing

The Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) is an international public database that provides high-throughput microarray and next-generation sequence-functional genomic datasets submitted by research communities (19). We obtained three IPF-related gene expression profile datasets from the GEO database: GSE150910, GSE92592, and GSE110147 (20-22), and the corresponding clinical information of each dataset is listed in . We downloaded the original data files and experimental information files of the above datasets. The GSE150910 dataset contains 103 IPF samples, 103 control samples, and the GSE92592 dataset contains 20 IPF samples and 19 control samples. The reference genomes for transcriptome sequencing of the GSE150910 and GSE92592 datasets are human gencode.v27 and human hg19, respectively. We used R4.0.3 software (http://cran.r-project.org) combined with the GenomicFeatures package (23) to annotate the gene names of the two different versions. The GSE110147 dataset contains 22 IPF samples and 11 control samples. We used R4.0.3 software to annotate the probe expression profile into the expression profile corresponding to the gene name, and the median expression of multiple probes corresponding to a single gene name was taken.

Table 1

Clinical characteristics of IPF patients and controls in the GSE150910, GSE110147, and GSE92592 datasets

Characteristics	GSE150910	GSE110147	GSE92592
Age (years)
<60	88	NA	NA
≥60	117	NA	NA
NA	1	NA	NA
Gender
Male	102	NA	NA
Female	104	NA	NA
Sample size
Total	206	33	39
IPF	103	22	20
Control	103	11	19
Smoke
Yes	105	NA	NA
No	83	NA	NA
NA	18	NA	NA

IPF, idiopathic pulmonary fibrosis; NA, not available.

Single-cell RNA sequencing data acquisition and pre-processing

We obtained the single-cell RNA sequencing dataset, GSE135893 (12), from the GEO database, and used the Seurat package in the R4.0.3 software to read the single-cell sequencing sample data for quality control (quality control threshold: nFeature_RNA >1,000, and mitochondrial-related genes account for less than 25%). The total expression level of each cell in the corrected data was further normalized to 10,000, and logarithmic normalization was subsequently performed. The normalized data were subjected to the identification of genes with highly variable expression between cells to obtain 2,000 genes. We used the ScaleData function to perform linear regression scaling based on the 2,000 highly variable genes, and then applied the principal component analysis (PCA) method for dimensionality reduction analysis. We further combined the harmony package to correct for batch effects between different samples in the PCA dimensionality reduction analysis results. The FindNeighbors function was used to construct a nearest-neighbor graph of the dimensionality reduction analysis results after eliminating the batch effect, which was then clustered, and the results were visualized.

IPF microenvironment analysis

Bulk gene expression

The xCell package (24) in R4.0.3 software was used to calculate the microenvironmental cell fractions of both the IPF and control samples for the three gene expression profiles, respectively. The combined the ggpubr package was applied to visualize the fractions with boxplots. P<0.05 was considered to indicate statistically significant differences.

Single-cell gene expression

The Seurat package (25) in R4.0.3 software was used to combine protein tyrosine phosphatase receptor type C (PTPRC), epithelial cell adhesion molecule (EPCAM), platelet and endothelial cell adhesion molecule 1 (PECAM1), actin alpha 2, smooth muscle (ACTA2), insulin-like growth factor binding protein 6 (IGFBP6), and other cell markers, and annotate the above single-cell clustering results as immune cells, epithelial cells, endothelial cells, myofibroblasts, fibroblasts, and other stromal cells, respectively. All of the above proteins were confirmed to be specific biomarkers of each cell type (26-30), and the proportions of these cells in the IPF and control groups were visualized using the ggplot2 package (31).

IPF differentially expressed genes (DEGs) and enrichment analysis

To identify DEGs between the IPF and control groups, we performed DEGs analysis on the GSE150910 dataset, with the largest sample size in the gene expression profile analyzed using the limma package (32) in R4.0.3 software [DEGs screening threshold: P<0.05 after correction by Benjamini-Hochberg (BH) method, |logFC| >1]. We then employed the gseGO function of the clusterProfiler package to perform gene set enrichment analysis (GSEA) (33) on the obtained DEGs combined with the BPs (Biological Processes) background in a logFC descending order (result screening threshold: P<0.05). Finally, BPs related to endothelial cells and fibroblasts in the enrichment results were visualized.

Correlation between IPF endothelial cells and multi-scale embedded gene co-expression network analysis (MEGENA)

In this study, we used R4.0.3 software combined with the cor.test function to analyze the Pearson’s correlation between the endothelial cell score and all GSE150910 dataset gene expression values. The screening threshold was as follows: P<0.05 and the absolute value of the correlation >0.3. We then used the MEGENA package (34) to perform multi-scale chimeric gene co-expression network analysis on the gene expression profiles corresponding to the screened genes, and identified the hub genes of the co-expression module as well as its sub-modules with the largest number of genes. Finally, we investigated their expression in the GSE150910 dataset.

Transcription factors (TFs) analysis

The upstream transcription factors of the hub genes were searched in ChIPBase v2.0 (https://rna.sysu.edu.cn/chipbase/), and cytoscape 3.9.0 (www.cytoscape.org/) was used to draw the transcription factor regulatory network corresponding to the hub genes.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Statistical analysis

In this study, the limma package was used for differential gene analysis, the difference significance between groups was calculated by Bayes corrected t-test, and then the difference significance P value was corrected by BH method. Correlation analysis was performed using Pearson. R4.0.3 software was used for statistical analysis, and P<0.05 was considered statistically significant.

Results

Advances in single-cell and transcriptome sequencing have facilitated exploration of biological systems with unprecedented resolution. Firstly, the composition of cell types from lung tissues of healthy control and IPF patients in the GSE150910 (n=206), GSE110147 (n=33), and GSE92592 (n=39) datasets was analyzed to identify and construct developmental ‘trajectories’ for lung fibrosis. Compared with the control group, there were a total of 20 different cell types with altered cell scores in the IPF group, including lymphoids (B-cells, CD8+ T-cells, CD8+ central memory T-cells (Tcm), class-switched memory B-cells, memory B-cells, naive B-cells, plasma cells, type2 T-helper (Th2) cells), stem cells (platelets), myeloids (activated dendritic cells (aDC), dendritic cells (DC), immature dendritic cells (iDC), macrophages M1, monocytes, neutrophils), and stromal cells (chondrocytes, endothelial cells, fibroblasts, lymphatic (ly) endothelial cells, microvascular (mv) endothelial cells), as shown in . Given that the relationship between immune-inflammatory cell (lymphoids and myeloids) changes and pulmonary fibrosis in transcriptome sequencing has been reported previously (35,36), we explored the underlying relationship between the changes to other cell types. As shown in , there were four types of cell scores increased in all three datasets, including fibroblasts, B-cells, chondrocytes, and plasma cells. Interestingly, the scores of endothelial cells, ly endothelial cells, and mv endothelial cells were significantly reduced. In addition, the amount of epithelial cells exhibited different change tendencies in these three datasets.

Figure 1

Cell score analysis of different cell types associated with EndMT in the lung tissue microenvironment of IPF patients and healthy controls. The microenvironment scores of the 20 major cell types screened in the training [GSE150910 (NIPF=103, NCtrl=103)] (A) and validation [GSE110147 (NIPF=22, NCtrl=11)] (B) and GSE92592 [(NIPF=20, NCtrl=19)] (C) sets were analyzed, and the ggpubr package was used for boxplot visualization of the microenvironment scores. aDC, activated dendritic cells; CD8+ Tcm, CD8+ central memory T-cells; DC, dendritic cells; iDC, immature dendritic cells; ly endothelial cells, lymphatic endothelial cells; mv endothelial cells, microvascular endothelial cells; Th2 cells, type 2 T-helper cells; IPF, idiopathic pulmonary fibrosis.

EndMT is a process in which endothelial cells lose their specific markers and obtain mesenchymal cell phenotype. This process is accompanied by the reduction of endothelial cells, leading to the proliferation and activation of fibroblasts (37). The fibroblast scores increased in all three datasets. We analyzed the cell fraction changes in the microenvironment through transcriptome sequencing data, which indicated that EndMT might occur during IPF.

Analysis of IPF microenvironment in single-cell sequencing data

In order to verify the hypothesis of IPF-associated EndMT, we further selected a single-cell sequencing dataset (GSE135893) to investigate the cellular heterogeneity and identify different cell populations. We performed annotation analysis on 33,042 cells from lung tissues of healthy control patients and 60,217 cells from lung tissues of IPF patients (total =93,259 cells). The cells were divided into five types, including immune cells, epithelial cells, endothelial cells, myofibroblasts, and fibroblasts (), which were labeled with biomarkers including CPTPRC, EPCAM, PECAM1, ACTA2, and IGFBP6. Furthermore, we used the ggplot2 package to visualize the proportions of the aforementioned cells in the control and IPF groups, respectively. As shown in , endothelial cells accounted for 9.9% in the control group, but decreased to 7.8% in the IPF group. Meanwhile, the fraction of fibroblasts obviously increased from 0.48% to 1.6% in the IPF group. In addition, the down-regulated proportion of endothelial cells was highly correlated to the increased amplitude of fibroblasts in IPF, and the proportion of myofibroblasts in the IPF group was also increased from 0.6% to 1.1%. By calculating the proportion of cells using single-cell sequencing data, we reported that EndMT occurred during pulmonary fibrosis.

Figure 2

Analysis of the IPF microenvironment using sing-cell sequencing data. (A) The cell markers (PTPRC, EPCAM, PECAM1, ACTA2, and IGFBP6) were used to annotate the single-cell clustering results as corresponding different cell types; the size of the dots corresponds to the percentage of cells expressing the marker in the cluster, and the dot color depth corresponds to the average expression level of the marker in the cluster. (B) Standard cell markers were used to label cell populations (cell types were classified as immune cells, epithelial cells, endothelial cells, myofibroblasts, fibroblasts, and other stromal cells) and color-coded by their associated major cell types were visualized with tSNE plot. (C) A histogram was used to show the proportion of the six major cell types in the IPF and control groups, color-coded by relevant taxa. (D) Pie charts were used to display the proportions of the six major cell types in the IPF and control groups, color-coded by relevant taxa. IPF, Idiopathic pulmonary fibrosis; tSNE, t-distributed stochastic neighbor embedding.

IPF DEGs and enrichment analysis

In order to explore the correlation between gene regulatory functions and EndMT in the IPF microenvironment, we first used the dataset with the largest sample size, GSE150910 (NIPF=103, NCtrl=103), to identify the DEGs. As shown in , the results showed a total of 327 up-regulated genes and 187 down-regulated genes in the IPF group. We further used these DEGs for GSEA enrichment analysis to reveal changes in BPs associated with endothelial cells and fibroblasts. The control group mainly exhibited correlations with endothelial cell migration and differentiation, including blood vessel endothelial cell migration, endothelial cell differentiation, and negative regulation of endothelial cell proliferation (). Meanwhile, BPs related to fibroblast function, including the cellular response to fibroblast growth factor stimulus and collagen fibril organization deposition, were mainly enriched in the IPF group ().

Figure 3

DEGs analysis and GSEA in the IPF (n=103) and control (n=103) groups in the GSE150910 dataset. (A,B) IPF DEGs analysis was performed on the GSE150910 dataset with the largest sample size in the gene expression profiles, in order to generate a hierarchical clustering heat map and volcano plot of the significantly DEGs (IPF DEGs analysis using limma package, BH-corrected P<0.05, |logFC| >1). (C,D) Using the gseGO function of the clusterProfiler package, GSEA was performed on the background of DEGs binding BPs, and the BPs related to endothelial cells and fibroblasts were visualized. (E,F) Functional enrichment network map and heat map displaying the genes involved in endothelial and fibroblast-related BPs and enriched therein. DEGs, differentially expressed genes; GSEA, gene set enrichment analysis; IPF, idiopathic pulmonary fibrosis; BPs, biological processes.

Based on the enrichment results, a functional enrichment network map and heat map were constructed to investigate genes that may regulate the BPs of endothelial cells and fibroblasts, as shown in . We found 15 genes that were enriched in endothelial cell-related BPs, including ACVRL1, CLDN5, TMEM100, SLC40A1, ARHGEF26, CLEC14A, KLF4, EFNB2, JCAD, PPARG, RGCC, CAV1, COL4A1, HEY1, and ID1, indicating that they may play an important role in regulating endothelial cell migration and differentiation. Also, there were 13 genes enriched in fibroblast-related BPs, including COL1A1, CXCL13, POSTN, SULF1, COL14A1, COMP, COL3A1, GREM1, SFRP1, SULF2, FGFBP1, TNC, and SCGB1A1, which may be involved in the regulation of fibroblast activation and collagen deposition. The above results suggested that in the occurrence and development of IPF, genes enriched in endothelial cell BPs may promote endothelial cell function and phenotype changes by regulating endothelial cell migration and differentiation, leading to fibroblast proliferation, further transformation into myofibroblasts (mesenchymal cells), and collagen deposition.

Analysis of MEGENA co-expression network and expression of hub genes

To further dissect the complex cooperative regulatory network between genes regulating endothelial cell function and phenotype in the IPF microenvironment, we used the cor.test function to analyze the Pearson correlation between endothelial cell fraction and all gene expression values in the GSE150910 dataset. A co-expression network map with multiple modules was obtained using MEGENA (). The network diagram consisted of 21 modules of different sizes, of which the largest was the No. 4 module, and CFAP73 and LRRC74B were the key nodes of this module. Module 4 contains six sub-modules, including c1_11, c1_12, c1_13, c1_14, c1_15, and c1_16.

Figure 4

IPF endothelial cell correlation and MEGENA analysis. (A) The cor.test function was used to analyze the Pearson correlation between endothelial cell score with all GSE150910 dataset genes expression values (screening threshold: P<0.05 and absolute value of correlation >0.3). The MEGENA package was also used for co-expression network analysis. Each black circle represents a module, and the nodes with labels are the hubs of the modules. (B) Enlarged view of co-expression module No. 4 and its sub-modules with the highest number of genes. (C) Expression of hub genes screened in module 4 and its sub-modules in the IPF and control groups. (D) Correlation between the hub genes expression levels and the endothelial cell score. IPF, idiopathic pulmonary fibrosis.

We further constructed an enlarged map of the co-expression network of module 4 (containing 98 genes) and its sub-modules (), and screened a total of nine hub genes, including CFAP73, LRRC74B, C12orf74, PTPRT, TUBA4B, STMND1, DNAH2, RSPH10B2, and GPR87. To explore the role of hub genes in IPF and their relationship with endothelial cells, we investigated the expression of these nine hub genes in IPF in the GSE150910 dataset, and found that the expressions of these nine hub genes was significantly up-regulated in the IPF group, which were negatively correlated with endothelial cell scores (). From this, we inferred that hub genes may be involved in regulating the endothelial cell function and phenotype changes in the IPF microenvironment, and their expressions may play an important role in EndMT.

TFs analysis of hub genes

We further searched for the TFs of the above 9 hub genes in ChIPBase v2.0, and drew the TFs regulatory network corresponding to the hub genes, as shown in . We can observe that FOXA1, ZBTB7A, CTCF, MYC, SP1, and CEBPB are the most common TFs that regulate the above hub genes. Among them, FOXA1 and ZBTB7A are involved in the transcriptional regulation of most hub genes, including CFAP73, LRRC74B, C12orf74, and STMND1, GPR87, TUBA4B, DNAH2.

Figure 5

The TFs regulatory network corresponding to the hub genes. TFs, transcription factors.

Discussion

Fibrotic tissue is very intricately and well organized into distinct tissues, including epithelial, endothelial, stromal, and immune cells. Numerous fibrosis diseases, such as pulmonary fibrosis, arise from abnormal activation in fibroblasts, which is mainly generated by inherent normal fibroblasts in lung tissues and epithelial or endothelial cell transformation. Owing to the limitation of traditional transcriptome analysis of bulk lung tissues, the hierarchy and heterogeneity of pulmonary cells within these lineages remain unclear.

EndMT is a cellular transformation process in which endothelial cells lose their endothelial characteristics and acquire a mesenchymal cell-like phenotype (38). EndMT was originally considered to be an important biological process in embryonic development and some cardiovascular diseases (39). Recently, it has been found that EndMT is closely related to the development of fibrotic diseases. When IPF occurs, endothelial cells acquire a mesenchymal phenotype and present typical markers of myofibroblast differentiation, such as α- smooth muscle actin (α-SMA), vimentin and collagens, while reducing the expression of vascular endothelial cadherin (VE-cadherin) (40). It has been found that in the bleomycin-induced pulmonary fibrosis model, under the action of transforming growth factor-β (TGF-β) combined with Ras signal activation, endothelial cells undergo EndMT process to generate a large number of fibroblasts (41). It has also been reported that EndMT caused by heat shock protein B (small) member 1 (HSPB1) deficiency contributes to the development of pulmonary fibrosis, suggesting that HSPB1-targeted therapy may be suitable for the treatment of a range of fibrotic diseases (42). Another study showed that vildagliptin improved pulmonary fibrosis by inhibiting EndMT (43). Therefore, EndMT is considered to be one of the key factors in the occurrence of IPF. However, there is no currently relevant research that explores the relationship between EndMT and IPF from the perspective of cell heterogeneity in the fibrotic microenvironment.

Considering that epithelial-mesenchymal transition (EMT) has been widely studied in IPF (44,45), we also focused on the changes in epithelial cell scores. Interestingly, the trends in the number of epithelial cells differed between the three datasets. It has been reported that epithelial cells proliferate significantly in the early stage of IPF (46). Another study has demonstrated that in the late stage of fibrosis, the number of epithelial cells decreases, while that of fibroblasts increases (47). Therefore, we speculated that the occurrence and development of EMT might be related to the different stages of IPF. In contrast, EndMT could persist throughout IPF.

In this study, we analyzed the IPF microenvironment through single-cell sequencing data and transcriptome gene expression data. Firstly, we found that the endothelial cells scores were significantly reduced, whereas fibroblasts scores were increased in all three datasets. It was inferred that EndMT might occur in the IPF microenvironment. In order to verify our speculation, we further performed IPF microenvironment analysis using single-cell sequencing data combined with PCA, and found that the proportion of endothelial cells was significantly down-regulated, while the proportions of fibroblasts and myofibroblasts were significantly increased. Moreover, the proportion of endothelial cells decreased was similar to that of fibroblasts increase. Therefore, we speculated that EndMT might occur in the occurrence and development of IPF.

When the amount of cells changes, the function of the cells also changes considerably (48). We performed GSEA to further explore the changes in the function and phenotype of endothelial cells during EndMT. Our results indicated that BPs related to endothelial cell function, including blood vessel endothelial cell migration, endothelial cell differentiation, and negative regulation of endothelial cell proliferation, were mainly enriched in the control group. BPs related to fibroblast function, including cellular response to fibroblast growth factor stimulus and collagen fibril organization, were mainly enriched in the IPF group. This indicates that the function and phenotype of endothelial cells is also changed in the EndMT of IPF.

Gene co-expression network analysis can effectively identify functional co-expressed gene modules associated with complex diseases (34). We performed MEGENA analysis to explore key modules related to endothelial cell function and identify hub genes. We found that there were nine hub genes in the most critical modules related to endothelial cell function regulation, including CFAP73, LRRC74B, C12orf74, PTPRT, TUBA4B, STMND1, DNAH2, RSPH10B2, and GPR87. We verified the expression of these genes in the dataset and found that their expression was significantly increased in the IPF group and negatively correlated with the endothelial cell scores, indicating that these genes are likely involved in the occurrence of EndMT in IPF. The previous view that PTPRP is closely related to cell adhesion, PTPRT interacts with endothelial cadherin (E-cadherin), and this interaction leads to dephosphorylation of E-cadherin, which affects the stability of the junction complex (49). During the EndMT process, endothelial cells lost intercellular adhesion, changed the cytoskeleton organization and acquired stronger migration and invasion ability (50), indicating that PTPRP may participate in the EndMT process by affecting cell adhesion and promote the development of IPF. The G protein-coupled receptor GPR87 is overexpressed in a variety of cancers, and studies have reported that GPR87 plays a key role in the proliferation of lung cancer cells, and indicated its potential as a new target for lung cancer therapy (51-53). GPR87 overexpression promotes EMT, significantly reduces the expression of the epithelial marker E-cadherin and increases the expression of the mesenchymal cell marker N-cadherin, and enhances the invasive and metastatic phenotype of cancer cells, which is associated with poor patient prognosis (54). TUBA4B can act as a tumor suppressor in various cancers (55,56), and low expression of TUBA4B is a predictor of poor prognosis and regulated cell proliferation in non-small cell lung cancer (57). The key genes obtained from the above screening may become new potential targets for the pathogenesis of IPF, and their correlation with IPF should be explored in future studies.

We further analyzed the upstream TFs of these hub genes, and found that FOXA1 and ZBTB7A may be involved in the transcriptional regulation of most hub genes. Studies have shown that long non-coding RNA 19 attenuates lipopolysaccharide-induced pulmonary fibrosis by regulating miR-423-5p/FOXA1 (58), and demonstrated increased expression of FOXA1 in isolated lung fibrotic tissues (59). ZBTB7A has been reported to play an important role in tumor metastasis (60), and ZBTB7A promotes tumor invasion by downregulating the expression of E-cadherin (61). High expression of TGF-β can combine with E-cadherin and β-catenin to inhibit ZBTB7A (62), and Smad4, a key component of TGFβ signaling, was found to interact with ZBTB7A (63). It indicated that TFs may participate in the signal transduction of IPF by targeting hub genes.

In this study, we discussed the heterogeneity of cells in terms of the IPF microenvironment for the first time, and demonstrated that EndMT may occur throughout the IPF process.

The article’s supplementary files as