The lung tissues from six PF-PH and five healthy subjects were obtained from patients undergoing an open lung biopsy or lung transplant in Shanghai Pulmonary Hospital in the period from 2020 to 2022. The diagnosis of PF-PH was established from the European Society of Cardiology and the European Respiratory Society guidelines (1). All of the experimental procedures using human tissues were approved and supervised by the Ethics Committee of Shanghai Pulmonary Hospital (numbers: K22-137Y). Written informed consent was obtained from all participants.

We used “pulmonary arterial hypertension” or “idiopathic pulmonary fibrosis,” the two typical diseases for pulmonary vascular remodeling and PF, as keywords to search in the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and obtain the expression profiles. Database inclusion criteria were as follows. First, the sequencing methods of the selected database are expression profiling by array, and the dataset must contain both control and disease groups. Second, the tissue sources used for sequencing were all human lung tissue. Third, the number of samples should be sufficient to ensure the accuracy of each analysis, and the database must provide either processed or raw data that can be used for reanalysis. Finally, we selected GSE53845, GSE113439, GSE110147, and GSE15197 from the GEO database. Subsequently, we used R 4.2.0 software to convert probes in GSE53845 and GSE113439 into gene symbols. We provide all involved R programs in Supplementary Table 1 .

PF- and PH-related modules were gained by weighted gene co-expression network analysis (WGCNA) (17), and WGCNA was carried out using SangerBox (sangerbox.com). The gene cluster map and the correlation heat map between modules and phenotypes were obtained.

We selected the modules highly correlated with PF and PH and screened the overlapping shared genes in the modules that are positively associated with PF and PH. ClueGO (18) explored the potential role of these shared genes in PF and PH and made a biological analysis of these shared genes. A p-value <0.05 was considered significant. Protein–protein interaction (PPI) network analysis used STRING (string-db.org), and the “MCODE” algorithm in Cytoscape software was used to realize visualization.

To further identify potential functional genes for PF and PH, differential expression analysis was performed using GEO2R between patients and healthy controls in GSE110147 or GSE15197. The screening threshold of differentially expressed genes (DEGs) was set at a log2 |fold change (FC)| >1 and a p-value <0.05. Overlapping DEGs in PF and PH databases were acquired by Venn. A set of immunity-related genes was achieved from the ImmPort database (19) and crossed with identified common genes by Venn mapping.

The normalized GSE53845 and GSE113439 expression data were analyzed by CIBERSORT (https://cibersort.stanford.edu/) to evaluate enrichment for immune invasions and to obtain an immunocyte infiltration matrix. Furthermore, we displayed the percentage of each immunocyte in the samples in a histogram, constructed a heatmap of 22 kinds of immunocytes, and compared the levels of 22 immunocytes between patients and controls using SangerBox. Spearman’s rank correlation analysis was used to analyze the relationship between immunocytes and feature genes. A p-value <0.05 was considered statistically significant.

Lung tissue samples were obtained from patients with PF-PH who underwent a lung transplant and healthy donors (controls) and stored at −80°C until use. Total RNA was extracted from the lung tissue using TRIzol reagent. RNA purity was determined using the NanoDrop 2000 spectrophotometer at 260/280 nm (ratio = 1.9–2.1). Then, reverse transcription using the PrimeScript™ RT reagent Kit (TaKaRa, RR0371, Japan) was performed according to the manufacturer’s instructions. RT-PCR was then performed with the SYBR Green-based RT-PCR (TOYOBO, QPK201, Japan). For RT-PCR analysis, the expression level of genes was represented as a fold change using the 2^−△△Ct method, and a p-value <0.05 was considered statistically significant. Primer sequences are listed in Supplementary Ta ble S2 .

Human pulmonary microvascular endothelial cells (PMECs), human pulmonary artery smooth muscle cells (PASMCs), and human fibroblasts were purchased from ScienceCell (Shanghai, China) and cultured in endothelial cell medium, smooth muscle cell medium, and Dulbecco’s modified Eagle’s medium, respectively, and incubated with 5% CO₂ at 37°C. Hypoxia-induced PMECs were acquired in anoxic incubators with 5% O₂ at 37°C. PASMCs and fibroblasts were induced to proliferate with 10 ng/mL platelet-derived growth factor-BB (PDGF-BB) (PeproTech, USA) and 5 ng/mL transforming growth factor-β1 (TGF-β1) (PeproTech, USA), respectively.

Cells were transfected with 20 μM siRNA IGF1 (Genomeditch, Shanghai, China) using Lipo2000 Transfection Reagent (Invitrogen, USA). The control groups were treated with equal concentrations of negative control sequences to eliminate nonspecific effects.

Cell Counting Kit-8 (CCK-8, Dojindo, Japan) was used to measure cell proliferation. Cells were inoculated into a 96-well plate at a density of 3 × 10³ cells per well. After 48 h of siRNA IGF1 treatment under the stimulation of proliferation induction, cells were treated with 10 μL CCK-8 solution at 37°C for 60 min and then tested. The absorbance at 450 nm was then detected, with a relative proliferation level obtained by normalization from five independent experiments.

In this study, 2 × 10³ treated cells were planted in each well of a 96-well plate. After 48 h of treatment, cell proliferation was measured by EdU Kit (BeyoClick™, China). Hoechst (blue) labeled the nucleus and EdU (green) labeled the proliferating cells, and these were observed and counted under the fluorescence microscope.

The Annexin V apoptosis kit (Dojindo, Japan) was used to detect the apoptosis of cells. In this study, 2 × 10⁵ cells were planted into each well of a 6-well plate with or without treatment. After 48 h, cells were digested with trypsin, and cells were washed with phosphate buffered saline (PBS) three times. Then, 500 μL of 1× Annexin V binding solution was added to the washed cells to prepare a cell suspension, and 3 μL of Annexin V-FITC and PI were added, respectively. The samples were mixed thoroughly and were left in the dark for 15 min. Finally, the apoptosis was analyzed by flow cytometry.

The cells were analyzed for fluorescence by Calcein-AM/PI staining (Calcein-AM/PI Double Stain Kit, YESEN, China). Cells were seeded in 24-well plates (2 × 10⁵ cells) for 48 h with or without treatment. Dyeing solution was added, and the cells were incubated at 37°C for 30 min. Finally, the cells were analyzed using a fluorescence microscope.

In order to further verify the ability of protein IGF1 to target drugs, molecular docking was carried out. Firstly, a database was selected in The NCGC Pharmaceutical Collection, which contained 7,929 FDA-approved drugs. The crystal structure of 2OJ9 was downloaded from the PDB database as the receptor protein structure for molecular docking. Secondly, the Discovery Studio LibDock program was used for molecular docking screening, and the top 20 compounds of LibDock score were selected for further study. Subsequently, for the top 20 compounds, the ligands and protein needed for molecular docking were prepared by AutoDock Vina software (http://vina.scripps.edu/), and for the target protein, the crystal structure obtained from the PDB database (https://www.rcsb.org/) needed pretreatment. Finally, the target structure was molecularly docked with the active ingredient structure, and the affinity value of the target structure is the binding ability of the target structure and the active ingredient structure using vina in pyrx software (https://pyrx.sourceforge.io/).

SangerBox (sangerbox.com) was used for all statistical analyses. Differences between the two groups were compared using parametric tests for normally distributed variables. All of the statistical differences were p-value <0.05.

Based on the screening criteria mentioned in the method, we selected four GEO datasets numbered GSE53845, GSE110147, GSE113439, and GSE15197. Furthermore, GSE53845 (20) and GSE113439 (21) served as discovery queues for WGCNA. A concise description of these samples was summarized in Supplementary T able S3 . Additionally, GSE110147 and GSE15197 served as validation queues for DEG analysis.

We identified eight relevant modules for GSE53845 through WGCNA, with each color representing a different module. A heat map of the module–trait relationship was then drawn based on the correlation coefficients to assess the association of each module with disease ( Figures 1A, C ). Among them, the “red” module (r = 0.81, p = 4.9e-12) had a high correlation with PF and was selected as the module with the highest positive correlation with PF, including 1,260 genes. Similarly, in GSE113439, module “pink” (r = 0.89, p = 1.8e-9) was the module with the highest positive correlation with PH, including 5,013 genes ( Figures 1B, D ).

Among PF and PH, 313 genes in the most relevant positive correlation module overlapped and were defined as gene set 1 ( Figure 1E ). To explore the potential function of gene set 1, Gene Ontology (GO) enrichment was analyzed using ClueGO. The first significantly enriched GO terms for Biological Process (BP) are “lymphocyte activation,” which accounted for 53.12% of the total GO terms ( Figures 2A, B ), confirming that this pathway may be essential in both PF and PH. Then, we constructed a protein-level PPI network for dataset 1 ( Figures 3A-D ). Four clusters were extracted using MCODE analysis, and cluster 1 contains 60 nodes and 759 edges (score = 25.73). We performed functional enrichment analysis on the genes of the four clusters respectively and selected the first 10 immune-related pathways ( Figure 3E ). GO enrichment analysis showed that the genes in cluster 1 were mainly related to the lymphocyte activation pathway; this was consistent with the enrichment of dataset 1. Therefore, this gene cluster belongs to the gene part shared by PF in PH.

We performed DEG analyses of GSE110147 and GSE15197 to validate our results. Among the 2,579 DEGs in GSE110147, 1,213 were significantly downregulated and 1,366 were significantly upregulated. In GSE15197, 3,497 DEGs were identified, including 1,328 downregulated and 2,169 upregulated genes. Further cross-analysis of GSE110147 and GSE15197 revealed that there were 32 co-downregulated genes and 165 co-upregulated genes in PF and PH, defined as gene set 2 and expressed by Venn diagram ( Figures 4A, B ). There were eight overlapped genes, ACTR2, COL5A2, COL6A3, CYSLTR1, IGF1, RSPO3, SCARNA17 and SEL1L, in GS1 and GS2 ( Figure 4C ). Notably, COL5A2, IGF1, and COL6A3 are contained in cluster 1, and RSPO3 is contained in cluster 3. At the same time, the expression level of the overlapped genes in the four datasets was visualized ( Figures 4E-H ).

In the enrichment assay, immune-related biological functions were significantly enriched in PF and PH. Therefore, we further performed CIBERSORT to predict immune cell infiltration between patients with the disease and control groups. The percentage of each of the 22 immune cells in each sample was shown in the bar and heat graphs of GSE53845 ( Figures 5A, B ) and GSE113439 ( Figures 5D, E ). The box diagram of the differences in immunocyte infiltration when samples were classified into disease and control groups showed significant differences of T cells CD4 memory resting, T cells CD4 memory activated, T cells regulatory (Tregs), T cells gamma delta, natural killer (NK) cells resting, Macrophages M1, Macrophages M2, and Neutrophils in PF patients compared with those in controls ( Figure 5C ). Compared with the control group, there were significant differences in T cells CD8, T cells CD4 memory activated, T cells follicular helper, NK cell activated, Mast cells resting, Eosinophils, and Neutrophils in patients with PH ( Figure 5F ). As can be seen from the results, T cells CD4 memory activated was a common significantly different immune feature to PF and PH. In addition, IGF1 and CYSLTR1 were the shared genes between 1,792 immune-related genes retrieved from the ImmPort database and eight common genes ( Figure 4D ). IGF1 and CYSLTR1 were also related to several immune cells ( Figure 6 ).

To further demonstrate the expression of overlapping genes in PF-PH patients, we performed RT-PCR in six lung tissue samples from PF-PH patients and five healthy donors. The results showed that IGF1 and COL5A2 expression was distinctly increased in PF-PH samples compared with healthy samples ( Figure 7 ), which is consistent with previous findings in DEGs of PF and PH.

In view of the upregulation of IGF1 in lung tissue of patients with PF-PH, we used siRNA IGF1 to reduce the expression of IGF in PMECs, PASMCs, and fibroblasts to further study the effect of IGF1 on cell function. Hypoxia, PDGF-BB, and TGF-β1 were used to induce the proliferation of PMECs, PASMCs, and fibroblasts, respectively. As shown in Figures 8A–C , the cell proliferation was increased obviously after stimulation and significantly inhibited after adding siRNA IGF1. In addition, we used EdU Kit to verify the cell proliferation, and the results were consistent ( Figures 8D–F ). Then, using flow cytometry and AM/PI Kit, we identified that siRNA IGF1 could reduce the apoptosis rate and death rate promoted by hypoxia, PDGF-BB, and TGF-β1 in PMECs ( Figures 9A, D ), PASMCs ( Figures 9B, E ), and fibroblasts ( Figures 9C, F ), respectively.

The LibDock docking results showed that 2,720 compounds in 7,929 were successfully docked to protein IGF1. The higher the LibDock score, the closer the interaction between molecules and proteins. We selected the top 20 compounds of LibDock score ( Table 1 ) to enter the next research. Subsequently, the lower the binding ability, the more stable the binding between the ligand and the receptor (the lower the binding energy, the better the binding) in vina docking results. Therefore, among the first 20 compounds, we screened out the compounds with binding energy less than -7.0 kcal/mol by vina (22), namely, deslanoside ( Figure 10A ), flavitan ( Figure 10C ), ivermectin ( Figure 10E ), and posaconazole ( Figure 10G ), which have strong binding activity with IGF1. The results also show the interaction between these drugs and IGF1. Deslanoside forms hydrogen bonds (HBs) with ARG1062, GLU985, and ARG973 of IGF1 and forms alkyl bonds with ARG1062 and ARG1054 ( Figure 10B ). Six HBs were formed between flavitan and IGF1, namely, SER1059, SER1063, ARG1054, GLU1115, GLU985, and GLU974, and also formed alkyl with ILE1130 and pi-alkyl with ARG1062 ( Figure 10D ). Ivermectin forms HBs with GLU985 and ASP1056 of IGF1 and alkyl with ARG1054 ( Figure 10F ). Posaconazole forms HBs with ILE1136, ARG1062, and SER1059 of IGF1, alkyl with ARG1054, pi-pi stacked with TYR1131, and halogen with GLU974 and GLU985 ( Figure 10H ).