Our study flow chart is shown in Figure 1 . All expression profiling data generated with arrays used in this study were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). Our inclusion criteria were as follows: more than 30 lung samples from Homo sapiens and concurrent inclusion of both IPF patients and healthy controls. According to the inclusion criteria, five GEO datasets including GSE32537 (healthy controls=50, IPF=119), GSE17978 (healthy controls=20, IPF=38), GSE53845 (healthy controls=8, IPF=40), GSE110147 (healthy controls=11, IPF=27), and GSE10667 (healthy controls=15, IPF=31) were identified. After background calibration and normalization of the data, the differentially expressed genes (DEGs) were screened by the Limma package according to the following criteria: |log2-fold change| > 1 and adjusted P value <0.05.

Gene Ontology (GO) is an international standardized gene function classification system to comprehensively describe the properties of genes and gene products in an organism, including biological process (BP), molecular function (MF) and cellular component (CC) terms (14). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a well-known database for pathway-associated investigation of genes (15). To explore functions in IPF, we performed GO and KEGG enrichment analyses of the DEGs identified between IPF and healthy control samples. Disease Ontology (DO) is a database established in 2003 that includes common and rare diseases (16). In this study, we also performed DO enrichment analysis to explore the diseases favored by the DEGs.

Among the included datasets, GSE32537 had the largest number of samples, so we selected it as our training set. First, LASSO regression and SVM-RFE machine learning algorithms were performed to identify candidate hub genes for a diagnostic model (17). Then, we constructed the diagnostic model using hub genes and calculated the risk score of the model using multiple logistic regression analysis. The model was presented as a nomogram. ROC curve analysis was subsequently performed to evaluate the diagnostic value of the nomogram. To further validate the reliability of the diagnostic model, we integrated GSE17978 and GSE53845 into one dataset, creating one dataset with more than 100 samples, as our main validation cohort. To eliminate batch effects in the merged dataset, we applied the Bioconductor “SVA” R package. We also used the GSE110147 and GSE10667 datasets as external validation cohorts. Calibration curve (CC) analysis, decision curve analysis (DCA) and clinical impact curve (CIC) analysis were also performed to evaluate the clinical diagnostic value of the model.

In this study, we used the Cell Type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm to identify the proportions of 22 immune cell types in lung samples from IPF patients and healthy controls (18). Correlation analyses between infiltrating immune cells and hub genes were also performed to explore which kinds of infiltrating immune cells were the main participants in the development of IPF. Differences were considered significant when the CIBERSORT output P value < 0.05. We also performed correlation analysis of the different types of infiltrating immune cells in IPF (19).

Our animal experiments were approved by the Ethics Committee of Sun Yat-sen University. All C57BL/6 mice were fed in a colony room in the Sun Yat-sen University Animal Center with a 12:12 hour light/dark cycle. We injected mice with bleomycin (Teva Pharmaceutical; 3 U/kg) or an equal volume of PBS intratracheally when they were eight weeks old. On day 21 after injection, all mice were harvested, and lung samples were collected for further analysis.

Mouse lung samples were embedded in paraffin and cut into sections. Then, we performed hematoxylin and eosin (H&E), Masson trichrome and Sirius Red staining to analyze pulmonary fibrosis. As previously described (20), the method for immunofluorescence (IF) staining of mouse lung samples was as follows: lung samples were incubated with appropriate primary antibodies ( Supplementary Table S1 ) overnight at 4°C after dewaxing and antigen retrieval. Then, we stained the sections with a secondary antibody ( Supplementary Table S1 ) that recognized the primary antibody for 40 min at room temperature. Stained sections were imaged with a Zeiss 800 laser scanning confocal microscope. The immunohistochemistry (IHC) method involves several steps. First, antigen retrieval and blocking were performed to enhance the visibility of target antigens. Next, primary antibodies, followed by secondary antibodies ( Supplementary Table S1 ), were sequentially incubated with lung tissues. Subsequently, diaminobenzidine (DAB) was applied as a chromogen to stain the target antigens. The presence of a brown signal is considered indicative of positive staining. As suggested in other study (21), the analysis of occupied area was conducted by selecting brown regions using a consistent threshold value within the macro function of ImageJ software (NIH). The findings were expressed as the percentage of the area occupied by positive staining relative to the total area in each specimen.

RNA was extracted from lung samples using TRIzol reagent (Molecular Research Center, Inc.) and reverse transcribed into cDNA with the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, K1622). Then, we performed quantitative real-time PCR (qPCR) using LightCycler480 SYBR Green I Master Mix (Roche, 4887352001-1). 18S served as the internal control. Differences were considered statistically significant with a P value < 0.05 using an independent-sample t test or the Mann-Whitney U test. All primer sequences used in this study are listed below: mouse Aspn forward 5’-TCCTCTGACAAGGTTGGACT-3’, and reverse 5’-AGAGAGTTGTCGTCATCATCGT-3’; mouse Sfrp2 forward 5’-CGTGGGCTCTTCCTCTTCG-3’, and reverse 5’-ATGTTCTGGTACTCGATGCCG-3’; mouse Slco4a1 forward 5’-CGATCTGCACAGCTACCAGAG-3’, and reverse 5’-GCTGACGAAGGTAAGGCATAG-3’; And mouse 18s forward 5’-GTGACGTTGACATCCGTAAAGA-3’, and reverse 5’-GCCGGACTCATCGTACTCC-3’.

Mouse lung samples were prepared using RIPA buffer (Beyotime, P0013B). The total protein concentration was assessed using the Pierce BCA Protein Assay (Thermo Fisher, 23227). After electrophoresis, proteins were transferred to PVDF membranes (Millipore). Then, the target proteins were immunoblotted with specific antibodies. The signal intensity of protein bands was visualized with a chemiluminescent substrate (Millipore). All protein bands from three independent blots were quantified using ImageJ software. The corresponding primary and secondary antibodies utilized are listed below: rabbit anti-Asporin (Invitrogen, PA5-28124), rabbit anti-SFRP2 (Affifinity Biosciences, DF4451), rabbit anti-SLCO4A1 (XY-Bioscience, XY12713), rabbit anti-GAPDH (Cell Signaling Technology, D16H11), and anti-rabbit IgG HRP-linked Ab (Cell Signaling Technology, 7074).

All statistical analyses performed in this study were conducted with R software, version 4.1.3 (http://www.r-project.org). An independent-sample t test was applied to validate the differences between two groups if the continuous variables were normally distributed; otherwise, we used the Mann-Whitney U test. A P value < 0.05 was considered statistically significant.

First, we downloaded RNA sequences of lung samples from the GEO database (GSE32537, including 119 IPF subjects and 50 healthy controls). After bioinformatic analysis with Limma, a total of 412 DEGs were identified between the IPF and healthy control samples, of which 283 were upregulated and 129 were downregulated ( Figure 2A ; Supplementary Figure S1A ).

Then, we investigated the potential biological functions of the DEGs through GO enrichment analysis. We found that the DEGs were mainly enriched in the extracellular matrix (ECM), extracellular structure and external encapsulating structure organization for BP terms. With regard to CC terms, these genes were mainly involved in collagen-containing ECM and external side of plasma membrane. In the MF category, they were strongly related to ECM structural constituent ( Figure 2B ). The GO enrichment analysis revealed that IPF was tightly linked to the pathological ECM deposition process. Moreover, DO enrichment analysis showed that these DEGs were mainly enriched in lung diseases. Together, these results indicated that these DEGs were reliable for subsequent IPF research ( Figure 2C ). Through KEGG enrichment analysis, we also found that IPF was strongly related to cytokine-cytokine receptor interaction and focal adhesion signaling pathways ( Figure 2D ).

We identified the candidate hub genes of a diagnostic model by LASSO regression combined with SVM-RFE machine learning algorithms. LASSO regression identified seven potential genes (ASPN, SFRP2, SLCO4A1, IL1R2, MMP7, FCN3, and CP) ( Figure 3A ), and the SVM-RFE algorithms screened eight potential genes (ASPN, SFRP2, FCN3, COL14A1, MGAM, SLCO4A1, CD24, and SPATA18) ( Figure 3B ). By taking the intersecting genes, we identified four candidate hub genes (ASPN, SFRP2, FCN3, and SLCO4A1) ( Figures 3C, D ). Each gene showed reliable classification between IPF and healthy control samples by ROC curve analysis ( Supplementary Figures S1B–E ).

To further verify the differential expression of candidate hub genes between IPF and control samples, we then established a bleomycin-induced pulmonary fibrosis mouse model and detected the expression of SLCO4A1, ASPN and SFRP2 between fibrotic lung samples and control samples from model mice ( Figure 4A ). Since FCN3 is a pseudogene in mice, we were unable to validate it in the mouse model (22). As shown by H&E, Masson trichrome, and Sirius Red staining, there was significant collagen deposition and fibrosis in the lung samples from bleomycin-induced model mice compared with those from PBS-treated mice ( Supplementary Figure S2A ). By qPCR and western blotting, we found that the mRNA and protein levels of ASPN and SFRP2 were upregulated in lung samples from bleomycin-induced model mice, while SLCO4A1 was downregulated ( Figures 4B–E ). Similarly, IF staining also produced the same findings ( Figure 4F ). Overall, our hub genes showed significant differential expression between fibrotic lungs and control lungs, and their change trends were consistent with those in the meta-GEO cohort including GSE32537, GSE17978, GSE53845, GSE110147 and GSE10667 ( Supplementary Figures S2B–D ).

To explore the relationships between candidate hub genes and infiltrating immune cells, the CIBERSORT algorithm was used. We found that all four candidate hub genes showed significant correlations with immune cells ( Supplementary Figure S3A–D ), including neutrophils, regulatory T cells, and resting mast cells. Among the different types of infiltrating immune cells, the expression of all four hub genes had the strongest link with a decline in the neutrophil level. The expression of ASPN (P<0.001, R=-0.63) and SFRP2 (P<0.001, R=-0.62) had a negative relationship with neutrophils ( Figures 5A, B ), while SLCO4A1 (P<0.001, R=0.75) and FCN3 (P<0.001, R=0.58) had a positive relationship with neutrophils ( Figure 5C ; Supplementary Figure S3E ). Ly6G has been widely acknowledged as a key marker for identifying neutrophils in murine models (23). In order to further substantiate the relationship between hub genes and neutrophils, we conducted an analysis of the correlation between the expression of hub genes and Ly6G using IHC analysis in a bleomycin-induced mouse model ( Figure 5D ). Our findings revealed a significant negative relationship between the expression of Aspn and Sfrp2 and the expression of Ly6G ( Figures 5E, F ), consistent with the aforementioned results. Conversely, the expression of Slco4a1 demonstrated a positive relationship ( Figure 5G ). Collectively, these results reinforce our previous findings and suggest a potential decrease in neutrophil levels within the immune microenvironment of IPF.

Considering that we could not validate the differential expression of FCN3 in animal models, we ultimately identified three hub genes (ASPN, SFRP2 and SLCO4A1) to construct a diagnostic model for IPF. The model that incorporated the above three hub genes was developed using multivariate logistic regression analysis and is presented as a nomogram ( Figure 6A ). Furthermore, we obtained a risk score based on the expression levels of the three hub genes, where the risk score = -5.6339 + (the expression level of ASPN × 0.8534) + (the expression level of SFRP2 × 0.7400) + (the expression level of SLCO4A1 × -0.8030). As shown in a CC, the nomogram for the classification of IPF and healthy control samples showed good agreement between prediction and reality ( Figure 6B ). DCA demonstrated that patients could benefit from the diagnostic model developed with the three hub genes at threshold probabilities from 0 to 1 ( Figure 6C ). Similarly, a CIC showed that the predicted number of high-risk patients was very close to the actual number of high-risk patients, further affirming that our model had clinical diagnostic value for IPF ( Supplementary Figure S3F ). Through ROC curve analysis, we found that the AUC of the diagnostic model was 1.000 for the training set ( Figure 6D ). To further validate the reliability of the diagnostic model, we merged two external datasets into a validation cohort (GSE17978 and GSE53835) that included 78 IPF lung samples and 28 healthy control lung samples. The AUC was 0.962 ( Figure 6E ). In regard to the other external cohorts, the AUCs of GSE110147 and GSE10667 were 0.985 and 0.925, respectively ( Supplementary Figure S3G, H ). Overall, it was suggested that our diagnostic model had a strong ability to distinguish IPF patients from healthy controls.

As shown above, all three hub genes had good diagnostic value for IPF and showed significant correlations with immune cells, especially neutrophils. To further explore the immune microenvironment in IPF, we performed immune cell infiltration analysis between IPF and healthy control samples using the meta-GEO cohort generated by merging five datasets, namely, GSE32537, GSE17978, GSE53845, GSE110147 and GSE10667. The proportions of different immune cells are shown in Supplementary Figure S4A , and we found that there was a significant difference in immune cells between IPF and healthy control samples ( Supplementary Figure S4B ), including differences in naïve B cells (P<0.001), memory B cells (P<0.001), plasma cells (P<0.001), naïve CD4 T cells (P<0.001), resting memory CD4 T cells (P<0.001), activated CD4 memory T cells (P<0.001), follicular helper T cells (P=0.024), gamma delta T cells (P=0.003), resting NK cells (P<0.001), monocytes (P<0.001), M0 macrophages (P=0.039), resting dendritic cells (P<0.001), activated dendritic cells (P=0.012), resting mast cells (P<0.001), eosinophils (P<0.001) and neutrophils (P<0.001). As shown above, the proportions of most infiltrating immune cells involved in activating an adaptive immune response were increased, and the frequencies of a majority of innate immune cells were reduced in IPF. By analyzing the correlations of different infiltrating immune cells, it was shown that resting memory CD4 T cells had the strongest negative correlation with CD8 T cells (R=-0.57) and neutrophils showed a negative correlation with resting mast cells in IPF lung samples ( Supplementary Figure S4C ). All of the above results showed that the changes in diverse infiltrating immune cells might have roles in the pathogenesis of IPF.