This article is mentioned in:

Introduction

Idiopathic pulmonary fibrosis (IPF) is a chronic progressive fibrosing pneumonia of the lung with no known etiology or pathogenesis. The disease is characterized by a damaged alveolar network that is gradually replaced with fibrous scars, resulting in lung deformity and organ failure. It is a fatal disease that kills patients within 2–5 years of diagnosis. The 5-year survival rate is 20% (1). The most common symptoms of IPF include dyspnea, persistent cough and fatigue. However, these symptoms are not specific to IPF and can be found in other lung conditions, such as chronic obstructive pulmonary disease, pulmonary tuberculosis and bronchiectasis, making diagnosis difficult due to lack of definitive diagnostic indicators. Thus, these challenges necessitate the identification of associated genes to develop diagnostic models and targeted therapies.

Disrupted RNA methylation and its associated downstream signaling pathways contribute to development and progression of various diseases, such as gastric (2) and breast cancer (3). Several types of RNA methylation have been discovered, including N6-methyladenosine (m6A), 7-methylguanosine and 5-methylcytosine (m5C), the latter being one of the most common modifications of RNA (4,5). Using regulatory genes involved in RNA methylation, prognostic signatures can be created. Studies (6,7) have been conducted to establish prognostic signatures using RNA methylation regulatory genes for cancer. Yi et al (6) identified a prognostic signature for head and neck squamous cell carcinoma based on m6A regulatory genes. Furthermore, Huang et al (7) developed a m5C-related signature to predict prognosis in cutaneous melanoma. RNA m5C methylation modification serves a role in the regulation of several cancers, including lung, gastric, liver, bladder, prostate and breast cancer, influencing their development, occurrence and invasive behavior (8–13). Altered m5C methylation levels can contribute to tumor progression and a poor prognosis in cancers. Understanding of the role of RNA methylation in IPF remains limited. However, research (14–16) in this area is expanding and yielding promising results about the potential consequences of altered RNA methylation in diverse pathologies. Further research is required to confirm and validate previous findings about this disease. The present study aimed to identify and validate the function of m5C-associated genes in IPF diagnosis and typing.

Recent research by Zhou et al (14) using the GSE150910 dataset investigated the role of m5C regulatory factor in IPF diagnosis; m5C regulatory factor mediates RNA methylation modification patterns and immune microenvironment infiltration characteristics, implying that this genetic factor may be used as an immunotherapy agent. The present study used data from patients with IPF from the Gene Expression Omnibus (GEO) to investigate the relationship between m5C-associated genes and the occurrence of IPF. The present study created a bleomycin (BLM)-induced IPF mouse model to validate these findings. The present results may serve as a foundation for future research into development and progression of IPF, as well as the search for new and effective therapeutic targets for PF treatments such as immunotherapy.

Materials and methods

Data collection and processing

A total of three datasets, GSE150910, GSE173355 and GSE124685 were chosen from the GEO database (hncbi.nlm.nih.gov/geo/). The GSE150910 dataset included 288 lung tissue samples, with clinical data for 103 IPF and 103 normal control (CTRL) samples. GSE173355 dataset contained 37 lung tissue samples, including 23 IFP and 14 CTRL samples. Both datasets had read count data from whole genome detection downloaded using the detection platform [GPL24676 Illumina NovaSeq 6000 (Homo sapiens)]. Finally, the GSE124685 dataset included 84 lung tissue samples, with clinical data for 49 IPF and 35 CTRL samples. Fragments per kilobase of exon model per million mapped fragments data from genome-wide detection was downloaded via the detection platform [GPL17303 Ion Torrent Proton (H. sapiens)]. The training datasets were GSE150910 and GSE173355, while the validation dataset was GSE124685. Because GSE150910 and GSE173355 contained different batches of gene expression data, sva package (version 3.38.0) in R3.6.1 (17) was used to eliminate batch effects and artifacts from the merged data. A total of 243 samples were collected, including 126 IPF and 117 CTRL samples.

Identification of differentially expressed genes (DEGs)

The limma package (version 3.34.7) in R3.6.1 (18) was used to investigate significant differences in m5C gene expression between patients with IPF and CTRL group. A false discovery rate (FDR) <0.05 was used to indicate significance. Spearman coefficients between pairs of DEGs were calculated using the R cor function and results were visualized with the heatmap package (version 1.0.8) (19,20).

Protein-protein interaction (PPI) and enrichment of functional pathways

PPI network of DEGs was established using STRING (version 11.0, string-db.org/) (21) while keeping linkage pairs with interaction scores ≥0.4, and the interaction network was visualized using Cytoscape (version 3.9.0) Display (22). The nodes in the network were annotated using DAVID (version 6.8), based on Gene Ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (23,24), with FDR <0.05 indicating significant enrichment.

Construction of diagnostic model based on m5C-related genes

One-way logistic regression analysis in R package rms (version 6.3–0) (25) was utilized to screen DEG levels determined in the combined dataset samples in step with the lars package (version 1.2) (26) using the least absolute shrinkage and selection operator (LASSO) algorithm for optimization of the m5C-related genes.

Construction of the diagnostic model

In the combined training set, the support vector machine method in R3.6.1 e1071 (version 1.6–8) (27) was used to build the disease diagnostic classifier based on m5C-associated genes (Core: Sigmoid Kernel; Cross: 10-fold cross-validation). The receiver operating characteristic (ROC) curve in R3.6.1 pROC (version 1.12.1) (28) was used to predict model accuracy in training and validation cohorts.

KEGG signaling pathway analysis

KEGG and gene data were obtained from the Gene Set Enrichment Analysis (GSEA) MSigDB database (gsea-msigdb.org/gsea/msigdb/index.jsp). The combined sample genome-wide expression level data were determined for each KEGG signaling pathway by gene expression levels, again using GSVA (version 1.36.3) in R3.6.1 language (29). The distribution of quantified KEGG signaling pathway values in different subtype groupings was compared with FDR <0.05 set as the threshold for significant differences.

Subtype grouping-related KEGG enrichment analysis

GSEA database was used to identify KEGG signaling pathways that were significantly associated with subtype grouping based on genome-wide expression levels in combined samples.

Screening DEGs associated with subtype groupings

In the combined sample expression profiles, limma (version 3.34.7) package in R3.6.1 (18) was used to identify DEGs between different subtype subgroups, with |log2FC| >1 and FDR <0.05 as threshold criteria. KEGG and GO analyses were used to investigate the biological functions of genes using DAVID (version 6.8) (23,24). FDR<0.05 was used as the threshold to select DEGs. T test in R3.6.1 language was used to compare the expression levels of immune checkpoint genes and Human leukocyte Antigen (HLA) family genes in different subtypes.

BLM-induced PF in mice

A total of 12 male C57BL/6 mice (age, 6–8 weeks; weight 16–20 g) were obtained from Beijing HFK Bioscience Co., Ltd. All mice were raised in a pathogen-free environment at 22–25°C, relative humidity of 50–60% and a 12/12 light-dark cycle. The mice had unrestricted access to food and water. The animal experiments were authorized by the Ethics Review Committee of Fujian Medical University (approval no. IACUC FJMU 2023-0327) and all procedures followed ARRIVE guidelines (30). BLM was acquired from Hisun Pfizen Pharmaceuticals Co., Ltd. The mice were anesthetized by intraperitoneal injection of 1% pentobarbital sodium (50 mg/kg) to prepare for BLM injection via the intratracheal route. The skin of the throat was cut lengthwise along the trachea and the muscles and fascia were bluntly separated using bending forceps to expose the trachea. A total of 12 mice were randomly divided into control group (50 µl 0.9% saline, administered intratracheally once) and BLM group (3.5 mg/kg BLM dissolved in saline to 50 µl, intratracheally injected once; both n=6). Overall health and welfare were regularly monitored, including their weight, respiratory rate, activity levels and food intake. Humane endpoints were weight loss >20%, respiration rate >100 breaths/min or abnormal behaviors. However, none of the mice died or reached the humane endpoints. After 21 days, mice were euthanized via cervical dislocation following anesthesia with 1% pentobarbital (50 mg/kg) by intraperitoneal injection and the lung tissues were harvested and collected for further analysis.

Histopathological analysis

The right lungs were fixed with 10% formalin at room temperature for 24 h, embedded in paraffin wax and cut into lung sections with a thickness of 3 µm. The sections were then stained with hematoxylin for 4 min and eosin for 20 sec and Masson trichrome (hematoxylin for 5–10 min, Ponceau S acid fuchsin stain for 5–10 min, and aniline blue for 5 min) at room temperature. Image acquisition was performed using a light microscope (CX-31, Olympus Corporation) at ×100 magnification.

Reverse transcription-quantitative (RT-q)PCR

Total RNA was extracted from the lung tissue using the SteadyPure RNA extraction kit (Accurate Biology, China), following the manufacturer's instructions. qPCR assay was then carried out using the EVO M-MLV RT kit and SYBR Green (both Accurate Biology) according to the manufacturer's instructions. Pre-denaturation was set at 95°C for 30 sec. Additionally, denaturation, annealing and extension were completed in 40 cycles at 95°C (5 sec) and 60°C (30 sec) respectively. The relative mRNA change was normalized using the 2−ΔΔCq method (31). Primer pairs are summarized in Table SI.

Statistical analysis

The statistical analysis was performed with R software (version 3.6.1). CIBERSORT (cibersort.stanford.edu/index.php) (32) was used to calculate relative abundance of immune cell subtypes in the training cohort and the proportions of immune cell subtypes in the IPF and CRTL groups were compared using the Kruskal-Wallis test. Spearman coefficients between the expression of m5C genes used in the model and immune cell types with significantly different distributions was calculated using the cor function in R. ConsensusClusterPlus (version 1.54.0) (33) in R was used to analyze samples for subtypes concerning diagnostically relevant m5C gene levels. The m5C scores of each sample were calculated using Gene Set Variation Analysis (GSVA) (version 1.36.3) (29). Kruskal-Wallis test was used to compare distribution of immune cells and m5C scores across subtypes in the IPF and CTRL groups. R package ESTIMATE was used to determine the immune, stromal and ESTIMATE scores of IPF samples and the differences in distribution between the immune cell subsets and ESTIMATE scores were analyzed using the Kruskal-Wallis test in R. The association between clinical traits and disease subtypes was investigated using Fisher's exact test. Additionally, an unpaired intergroup t test was used to compare gene expression levels across subtype groups. P<0.05 was considered to indicate a statistically significant difference.

Results

m5C regulators expressed differentially between IPF and CTRL samples

Following data processing and quality control, 126 IPF and 117 CTRL samples were examined (Fig. 1). A total of 23 m5C genes were significantly differentially expressed, of which methyl-CpG binding domain protein1 (MBD1), Single-strand selective monofunctional uracil-DNA glycosylase1, Aly/REF export factor (ALYREF), Uracil-DNA glycosylase (UNG), Zinc finger and BTB domain-containing protein38 (ZBTB38), NSUN6, ZBTB33, Methyl CpG binding protein2 (MECP2), UHRF1 and TDG were significantly upregulated in IPF compared with CTRL (Fig. 2A). Highest positive correlations were between TET3 and ZBTB4, TET2 and MECP2, and MECP2 and ZBTB4, while the strongest negative correlations were between NSUN5 and TDG and TET3 and ZBTB4 (Fig. 2B).

Figure 1. Study design. IPF, idiopathic pulmonary fibrosis; NCBI GEO, national center for biotechnology information gene expression omnibus; CTRL, control; PPI, protein-protein interaction; DEG, differentially expressed gene; q, quantitative; GVSA, gene variation and structural analysis; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 2. Differentially expressed m5C genes. (A) Expression levels of m5C genes differ significantly between IPF and CTRL samples. (B) Correlation plot of the expression levels of significantly differentially expressed m5C genes. (C) Interaction network of differentially expressed m5C proteins. Thicker lines represent higher interaction scores. (D) Differentially expressed protein enrichment. *P<0.05, ***P<0.001. IPF, idiopathic pulmonary fibrosis; m5C, 5-methylcytosine; CTRL, control; FDR, false discovery rate; FC, fold change.

PPI network and functional pathway enrichment analysis

A PPI network (Fig. 2C) was created to investigate interactions between proteins encoded by DEGs. This produced a total of 132 linkage pairs. To analyze systematic characterization and biological functions of m5C proteins in the network, GO and KEGG analysis was performed using DAVID with FDR <0.05 set as the threshold. In total, 17 biological processes and four KEGG signaling pathways were identified. IPF-associated DEGs were involved in GO biological processes, including ‘depyrimidination’, ‘RNA methylation’, ‘base-excision repair’, ‘negative regulation of transcription from RNA polymerase II promoter’ and ‘negative regulation of transcription, DNA-templated’. KEGG pathway analysis revealed that IPF-associated DEGs were enriched in ‘base excision repair’, ‘microRNAs in cancer’, ‘metabolic pathways’ and ‘cysteine and methionine metabolism’ (Fig. 2D).

Construction of a diagnostic model centered on m5C-related genes

To assess the role of m5C-associated genes in IPF diagnosis, 23 m5C-related DEGs were screened using univariate Cox regression. A total of nine genes, including NTHL1, NSUN5, DNMT3B, MBD3, NSUN6, ZBTB33, UHRF1, TDG and TET2, were identified (Fig. S1). The final optimized m5C-related genes were then screened using the lars package and the LASSO algorithm in R (Fig. S2). Based on feature selection of LASSO algorithm, the ZBTB33 gene was excluded. In the combined training set, NSUN6, UHRF1, TDG and TET2 were highly expressed in patients with IPF, whereas NTHL1, NSUN5, DNMT3B and MBD3 were downregulated; similar expression patterns were also found in the independent validation dataset GSE124685. Area under the curve was >0.8 in both the training and validation datasets, indicating the potential effectiveness for IPF diagnosis (Fig. S3).

Diagnostic m5C gene and immune correlation analysis

Studies (34,35) have demonstrated that there is a fibro-inflammatory response during formation of PF and immune cells may accumulate around PF lesions. CIBERSORT was used to calculate immune infiltration in the training dataset, comparing distribution between IPF and CTRL groups. The proportion of immune cell infiltration varied among the 13 types (naive, memory, B cell plasma, T cell CD4 naive, T cell CD4 memory resting, T cell follicular helper, T cell regulatory Tregs, NK cell resting, Monocyte, Macrophage M0, Mast cell activated, Eosinophil, and Neutrophil). Resting memory CD4 T cells were strongly and positively associated with NSUN6, TDG, and TET2 levels and negatively correlated with NSUN5 (Fig. 3).

Figure 3. Diagnostic m5C gene and immune correlation analysis. (A) Distribution of 22 immune cell types. (B) Heatmap showing immune cell correlations for diagnostic m5C genes with significantly different distributions. *P<0.05, **P<0.01, ***P<0.001 vs. CTRL. m5C, 5-methylcytosine; IPF, idiopathic pulmonary fibrosis; NK, natural killer; CTRL, control.

Sample subtype analysis based on diagnostic m5C genes

All samples were classified into two subtypes based on the levels of the eight diagnostic m5C genes. Subtype 1 samples had significantly higher m5C score than subtype 2 (Fig. 4A and B). There was a significant difference between the subtypes based on sex but age and smoking had no significant effect (Fig. 4C; Table I). Immune cell infiltration varied by subtype, with 15 types of immune cell differently distributed (Fig. S4A). ESTIMATE scores also differed significantly between subtypes (Fig. S4B), with subtype 2 showing higher resting memory CD4 T cell infiltration than subtype 1. Immune checkpoint and HLA family gene analysis revealed 16 HLA family genes with significantly different expression: Beta-2-Microglobulin, Human Leukocyte Antigen-A (HLA-A), Human Leukocyte Antigen-C (HLA-C), HLA-DMA, HLA-DOA, HLA-DOB, HLA-DBP1, HLA-DQA1, HLA-DQA2, HLA-DQB1, HLA-DRA, HLA-E, HLA-F, HLA-G, Transporter2 ATP-Binding Cassette Sub-Family B (TAP2) and TAP Binding Protein (TAPBP) (Fig. S5A). A total of seven immune checkpoint DEGs, CD27, CD274, CD40, CD70, CD86, Cytotoxic T-lymphocyte-associated protein4 (CTLA4) and Hepatitis A virus cellular receptor2 (HAVCR2), were also discovered (Fig. S5B).

Figure 4. Subtype analysis. (A) Clustering of sample subtypes. (B) Comparison of m5C score. (C) Expression levels of the eight DEGs according to clinical factors. (D) KEGG signaling pathways with significantly different distributions. (E) Heat map of DEGs. Subtype classification of DEGs significantly associated with (F) biological processes and (G) KEGG signaling pathways. FDR, false discovery rate; m5C, 5-methylcytosine; DEG, differentially expressed gene; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table I. Clinical information of samples.

Table I.

Clinical information of samples.

	Subtype
Characteristic	1 (n=106)	2 (n=20)	P-value
Age, years
≤60	49	6	0.2212
>60	56	14
Sex
Female	46	4	0.0479
Male	60	16
Smoking status
Yes	57	14	0.3051
No	40	5

Genome-wide expression-based GSVA quantification of KEGG analysis

A total of 163 KEGG signaling pathways with significantly different distribution were screened by comparing distribution of quantified KEGG signaling pathway values across subtypes. The logFC values of differences were ordered from smallest to largest and the top 10 KEGG were identified (Fig. 4D). KEGG signaling pathways of subtypes 1 and 2 differed significantly: Subtype 1 was primarily involved in ‘glycan biosynthesis’, ‘dorsoventral axis formation’ and ‘ECM receptor interaction’, whereas subtype 2 was primarily enriched in metabolic aspects such as ‘glycine, serine and threonine metabolism’, ‘retinol metabolism’ and ‘alpha-linolenic acid metabolism’.

Subtype grouping-related KEGG enrichment analysis

Using GSEA, the KEGG signaling pathways that were significantly associated with subtype groupings were screened and 12 significantly associated KEGG signaling pathways were identified (Fig. S6). In general, the higher the absolute value of NES, the lower the P-value, implying that the higher the enrichment of the functional gene set, the greater the confidence in the analysis results.

Screening significantly differentially expressed genes associated with subtype grouping and enrichment analysis of GO biological processes and KEGG signaling pathway

A total of 2,870 DEGs were identified (Fig. 4E). DAVID was then used to identify genes that were significantly differentially expressed for GO biological processes and KEGG signaling pathways. In total, 55 biological processes and 12 KEGG signaling pathways were screened (Fig. 4F and G), and the top 10 in each category were visualized based on FDR value (smallest to largest). ‘Ribosome’ and ‘oxidative phosphorylation’ showed the most significant difference.

Normal and fibrotic mouse lung tissue

Histopathological staining revealed BLM-induced PF in mice. Mouse lung tissue showed a high concentration of fibroblasts, severe and extensive interstitial fibrosis, alveolar structure destruction and irregularly shaped and small alveolar cavities. In Masson-stained sections, a large blue-stained collagen deposit could be seen, and the alveolar structure was destroyed, indicating a more typical form of lung interstitial fibrosis (Fig. 5A).

Figure 5. Lung fibrosis tissue in mice. (A) Normal and fibrotic lung tissue. Magnification, ×100. (B) Expression of eight diagnostic 5-methylcytosine -related genes in mouse lung tissue. Control vs. BLM. **P<0.01, ***P<0.001, ****P<0.0001. HE, hematoxylin and eosin; NSUN, nucleolar RNA methyltransferase; NTHL, nth-like; DNMT3B, DNA methyltransferase3 beta; MBD, methyl-CpG binding domain; UHRF, ubiquitin-like with plant homeodomain and ring finger; TDG, thymine DNA glycosylase; TET2, ten-eleven translocation2; BLM, bleomycin

Expression of m5C methylation-related genes in lung fibrosis tissue samples

NSUN6, UHRF1, TDG and TET2 were significantly upregulated, while the expression levels of NSUN5, NTHL1, DNMT3B and MBD3 were downregulated in fibrotic compared with normal lung tissue (Fig. 5B).

Discussion

IPF is a progressive fatal lung disease whose symptoms frequently overlap with other conditions, making diagnosis difficult. Identifying factors involved in its development and progression can facilitate development of reliable markers capable of accurately predicting the prognosis of patients with IPF (4). RNA modifications such as m6A, m5C and m1A are increasingly recognized to serve critical roles in diverse biological processes, and their dysregulation has been linked to development and progression of several diseases, particularly cancer, by influencing gene expression (5,13,36,37). There has been a great deal of research on m6A regulators and their role in the etiology of IPF (15,16), but to the best of our knowledge, there have been no studies on m5C. The present study aimed to develop a reliable and effective diagnostic model for IPF using m5C-associated regulatory genes. The differential expression and interactions of 23 IPF-associated m5C regulatory genes were investigated using GEO database, resulting in the development of a diagnostic model. The model identified two subtypes of patients with IPF. Furthermore, CIBERSORT was used to compare differences in immune cell distribution and immune checkpoint and HLA family gene expression.

The present study identified eight differentially expressed IPF-associated m5C regulatory genes. These included NTHL1, DNMT3B, MBD3, UHRF1, TDG, NSUN5 and 6 and TET2. The model efficacy was assessed in training and validation datasets, and each of the genes was further validated using a PF mouse model. The results were consistent with the information obtained from GEO: NSUN6, UHRF1 and TDG were significantly upregulated, while the expression levels of NSUN5, NTHL1, DNMT3B and MBD3 were downregulated in fibrotic compared with normal lung tissue. TET2 demethylates both DNA and RNA (38). It preserves the stemness of trophoblast stem cells by converting 5-methylcytosine to 5-hydroxymethylcytosine and its inhibition slows trophoblast stem cell proliferation and promotes epithelial-mesenchymal transition (EMT) (39). EMT process involves transformation of epithelial cells into mesenchymal-like cells. In ~33.3% of cases of PF, myofibroblasts undergo EMT (40). Furthermore, TET3 has been shown to be upregulated in IPF fibroblasts, which contradicts the present findings (41). Macrophage polarization in the alveoli has also been found in PF caused by low levels of DNMT3B, IL-4 and transforming growth factor β1 (TGF-β1) (42). This suggests DNMT3B may play a protective role against fibrosis. This is consistent with the present findings, which showed that DNMT3B was downregulated in the IPF samples compared with CTRL. Furthermore, development of fibrosis has been linked to DNMT1, which regulates DNA methylation in fibroblasts and alveolar epithelial cells (43). NTHL1 encodes a DNA repair enzyme that primarily removes oxidatively damaged bases from DNA molecules. Inflammation in the lungs can cause oxidative DNA damage (44), which NTHL1 can help to mitigate. However, aberrations in the NTHL1 repair system cause pathological oxidative damage in the lung, eventually leading to PF pathogenesis (45,46). One study found UHRF1 is a mediator of KRAS-driven oncogenesis in lung adenocarcinoma (47). The protein encoded by TDG is a key DNA repair enzyme that is primarily responsible for the removal of modified pyrimidine bases, specifically demethylated thymine. It is key for maintaining genomic stability and regulating the DNA damage response. A study discovered a link between DNA damage and development of PF, implying that TDG could play an important role in repairing damage (48). NSUN5 and NSUN6 are members of the NSUN family, encoding RNA methyltransferase proteins (49). Previous studies (50,51) on NSUN5 have concentrated on its role in neurological disorder. NSUN6 is associated with development of pancreatic (52) and colorectal cancer (53). However, no link has been found between NSUN5, NSUN6 and PF to date. Compared with mice without macrophage MBD2 deficiency, mice with MBD2 deficiency in macrophages are protected from BLM-induced pulmonary fibrosis and exhibit significantly lower levels of TGF-β1 and M2 macrophages (54). MECP2 is associated with pathological fibrosis in the heart by promoting fibroblast proliferation and fibrosis via Dual Specificity Phosphatase5 downregulation (55). Here, MECP2 was highly expressed in IPF samples compared with CTRL, implying that the gene may also regulate lung fibrosis, as confirmed by Cheng et al (56). The present findings indicated that ZBTB4 may regulate development of PF via m5C methylation, as it was significantly differentially expressed in the IPF and CTRL groups. This demonstrated that the model based on these m5C-related genes, combined with disease subtype classification, may be useful for diagnosing IPF.

IPF is characterized by increased infiltration of inflammatory cells; whether this is a primary cause of IPF or an epiphenomenon remains unknown (34). Furthermore, immune dysregulation is hypothesized to contribute to IPF development (35). In the present study, the proportion of inflammatory cells such as T and B cells, macrophages, monocytes, natural killer (NK) cells, neutrophils and eosinophils was significantly different between IPF and CTRL groups. Macrophages are abundant in healthy lungs (57). Alveolar macrophages (AMs) maintain homeostasis by removing cell debris, including apoptotic cells, regulating wound healing and initiating anti-pathogen immune responses. Monocytes are recruited and stimulated to differentiate into macrophages in response to lung injury (58). However, M2 macrophages predominate in the lung of patients with IPF/Usual Interstitial Pneumonia and animal studies indicate that M2 macrophages may be a useful target for treating and preventing PF (59,60). These cells produce a high-affinity IL-13 receptor, IL-13Rα2, which interacts with IL-13 and increases production of TGF-β1, promoting fibrosis (61). Although patients with IPF had significantly higher levels of M0 macrophages, the present study found no difference in number of M2 macrophages between the IPF and CTRL groups. Little is understood about how NK cells and M0 macrophages contribute to PF. The present findings may provide novel insight into how immunological infiltration occurs in PF. Resting memory CD4+ T cells are hypothesized to serve as reservoirs for viruses such as HIV and, if activated, they can promote HIV infection (62). Resting memory CD4+ T cell levels were significantly lower in patients with IPF and they showed strong positive correlations with TET2 and TDG, as well as a negative correlation with NSUN5. TET2 is involved in development of IPF (63). Thus, combination of resting memory CD4+ T cells and TET2 may offer a novel approach to IPF diagnosis.

Based on expression of m5C-associated genes, samples were classified into two subtypes. Subtype 1 had significantly higher m5C scores compared with subtype 2 samples. Sex was the only clinical characteristic that distinguished subtypes. The distribution of immune cells revealed that subtype 1 contained significantly more M2 macrophages than subtype 2. By contrast, resting memory CD4+ T cells, which have not previously been linked to IPF tissues, were highly expressed in subtype 2. However, currently, IPF is not clinically classified into distinct subtypes. The present study found differences in the expression of m5C-related genes across subtypes. Future studies should enroll patients and measure the expression levels of these m5C genes and assess correlation with the clinical characteristics associated with each subtype, resulting in a more complete understanding of potential stratification in IPF.

Subtype 1 had significantly lower levels of CD27 and CD70 expression than subtype 2. CD27, a member of the TNF receptor family, is found almost exclusively on naive CD4 T cells (64). CD27- and CD28-expressing T cells are associated with various inflammatory lung conditions, including IPF, and CD27 levels are associated with lung function parameters in patients with IPF (65). A key aspect of the pathophysiology of fibrosis is the secretion of extracellular matrix proteins, which is reduced when CD70 binds to fibroblasts (58). CD70 is a type II transmembrane member of the TNF family found in both fibroblasts and lung tissue. CD70 and CD27 work together as a ligand-receptor system to regulate T cell co-stimulation and interaction with fibroblasts (65). Thus, CD27-CD70 interaction may be a promising target for fibrosis treatment. CD40 regulates a range of processes, including innate and adaptive immune responses. Several inflammatory pulmonary conditions, including acute lung injury, bronchial asthma interstitial pneumonia and acute respiratory distress syndrome, are associated with CD40-CD40L interactions (66). CD40 also decreases inflammation in the early stages of IPF, making it a potential target for slowing progression of PF. Furthermore, CD274, CD86, CTLA4 and HAVCR2 expression were significantly higher in subtype 2 compared with subtype 1. Programmed cell death 1 ligand 1 (PD-L1), also called CD274 or B7-H1, is a member of the B7 family of immune regulatory molecules. PD-L1 is highly expressed in numerous types of tissues, including the heart and lungs. Wang et al (67) found that PD-L1 knockout in septic mice decreases plasma levels of TNF-α and IL-6 and alveolar edema. This suggests that PD-L1 serves a protective role against lung inflammation. Immune checkpoint inhibitor-associated pneumonia (CIP) frequently develops in patients receiving PD-L1 inhibitors in clinical practice (68–70). CIP frequently presents with interstitial lung changes, which can progress to PF in severe cases (70). The present findings support the protective effects of PD-L1, which may be applicable in IPF. The immunoglobulin superfamily includes B lymphocyte-activating antigen B7-2, also known as CD86, which is expressed as a type I membrane protein in dendritic and Langerhans cells and binds to T cell surfaces, activating CD28 and CTLA-4 (71). CD86 negatively regulates T cell activation and decreases immune responses when bound to CTLA4 (72,73). HAVCR2, also known as TIM3, is an immune checkpoint receptor protein that inhibits antitumor immune responses. Its overexpression by AMs exacerbates PF (68). HLA-G expression on mast cells helps to counteract fibrosis (74). Certain morphological features, such as male sex, lower baseline forced vital capacity and diffusing capacity of the lung for carbon monoxide (DLCO) and older age, have been identified in retrospective studies as risk factors for progressive fibrosing interstitial lung disease progression and mortality (75–77). Patients with IPF are more likely to be male and the present study found a significant difference in the sex distribution of subtypes, with subtype 2 having a higher proportion of male patients than subtype 1. Subtype 2 is hypothesized to have a worse prognosis than subtype 1 because it expresses lower levels of CD27 and CD70 and has higher levels of CD274, CD86, CTLA4 and HAVCR2 compared with subtype 1. This requires additional clinical and experimental verification. The association between m5C methylation factors and immune checkpoints also warrants further investigation.

The present study had limitations. First, due to individual differences and other confounding factors, findings based on existing databases may not be accurate. The present results should be verified with a well-designed multicenter clinical study based on patient clinical data. The small sample sizes in GSE150910 and GSE173355 datasets are a second limitation. Finally, because IPF is rare, insufficient specimens have been collected to confirm the predicted DEGs and the disease diagnosis model; this should be validated in future studies.

To the best of our knowledge, the present study is the first to develop an effective prognostic signature for IPF and the findings highlighted the significance of m5C-related genes in diagnosing and typing IPF, as well as expression of immune cells and immune checkpoint-associated genes in patients with IPF and different disease subtypes.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by the National Natural Science Foundation of China (grant no. 82171566), Joint Funds for the Innovation of Science and Technology, Fujian Province (grant no. 2020Y9086) and Foundation of Fujian Provincial Department of Finance (grant no. 2021XH005).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

LC designed and supervised the study. LT analyzed data. LT and WS wrote the manuscript. WS, JW and YL performed bioinformatics analysis. LT and LC confirm the authenticity of all the raw data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

The animal experiments were approved by the Ethics Review Committee of Fujian Medical University (approval no. IACUC FJMU 2023-0327) and all steps were carried out in accordance with ARRIVE guidelines.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References