Identification of potential biomarkers for idiopathic pulmonary fibrosis and validation of TDO2 as a potential therapeutic target

Abstract

BACKGROUND

Idiopathic pulmonary fibrosis (IPF) is a progressive interstitial lung disease with a high mortality rate. On this basis, exploring potential therapeutic targets to meet the unmet needs of IPF patients is important.

AIM

To explore novel hub genes for IPF therapy.

METHODS

Here, we used public datasets to identify differentially expressed genes between IPF patients and healthy donors. Potential targets were considered based on multiple bioinformatics analyses, especially the correlation between hub genes and carbon monoxide diffusing capacity of carbon monoxide, forced vital capacity, and patient survival rate. The mRNA levels of the hub genes were determined through quantitative real-time polymerase chain reaction.

RESULTS

We found that TDO2 was upregulated in IPF patients and predicted poor prognosis. Surprisingly, single-cell RNA sequencing data analysis revealed significant enrichment of TDO2 in alveolar fibroblasts, indicating that TDO2 may participate in the regulation of proliferation and survival. Therefore, we verified the upregulated expression of TDO2 in an experimental mouse model of transforming growth factor-β (TGF-β)-induced pulmonary fibrosis. Furthermore, the results showed that a TDO2 inhibitor effectively suppressed TGF-β-induced fibroblast activation. These findings suggest that TDO2 may be a potential target for IPF treatment. Based on transcription factors-microRNA prediction and scRNA-seq analysis, elevated TDO2 promoted the IPF proliferation of fibroblasts and may be involved in the P53 pathway and aggravate ageing and persistent pulmonary fibrosis.

CONCLUSION

We provided new target genes prediction and proposed blocking TGF-β production as a potential treatment for IPF.

Key Words: Idiopathic pulmonary fibrosis; Lung function; Overall survival; Transforming growth factor-β; TDO2 inhibitor

Core Tip: This study is unique in several aspects: (1) We identified six hub genes for idiopathic pulmonary fibrosis (IPF) and determined through quantitative real-time polymerase chain reaction; (2) Multi-omics analysis proved that TDO2 was upregulated in IPF patients and promoted the IPF proliferation of fibroblasts; (3) TDO2 may be involved in P53 pathway and aggravate aging and persistent pulmonary fibrosis; and (4) TDO2 inhibitor effectively suppressed transforming growth factor-β-induced fibroblast activation.

INTRODUCTION

Idiopathic pulmonary fibrosis (IPF) is a chronic progressive fibrotic lung disease characterized by exertional dyspnoea and diminished lung function. Its progression seriously affects the quality of life of patients, leading to respiratory failure and death in severe cases. In addition, IPF is a fatal disease characterized by abnormal lung epithelial cells and excessive accumulation of lung interstitial matrix[1]. Clinical symptoms are usually progressive and manifest as shortness of breath, dry cough, and dyspnoea. Smoking exacerbates the loss of lung function in patients with IPF[2]. In addition to smoking, other factors, including viral and bacterial infections, genetic factors, age and sex, have different effects on IPF patients[3-7]. IPF mainly occurs in male patients and may develop earlier in men than in women for various biological reasons[8]. In terms of age, IPF predominantly affects elderly individuals, and its prevalence increases with age[9]. These microinjuries induce abnormal epithelial-fibroblast communication[10,11]. Moreover, slight injuries induce matrix-producing myofibroblasts and the accumulation of extensive extracellular matrix (ECM) with remodelling in the lung interstitium[12].

The median annual survival after the diagnosis of IPF is 2 to 3 years[13]. Although two antifibrotic drugs (pirfenidone and nintedanib) are currently available, they can only slow the progression of the disease instead of curing IPF[14]. There is no effective treatment for the illness, and a major need for new therapies has not been satisfied[15]. Therefore, identifying and intervening in key genes of IPF and exploring new therapeutic approaches are essential. The diagnosis of IPF usually requires respiratory physicians, radiologists, and pathologists to review various clinical features, imaging results, and biopsy results of patients in a group discussion to arrive at the final diagnosis. Standard imaging assessment of IPF with high-resolution computed tomography provides diagnostic and predictive information[16]. The degree of fibrosis and cellularity correlated with forced vital capacity (FVC) and diffusing capacity of carbon monoxide (DLCO), as well as predicted mortality, can be observed. The most common indicators of lung function associated with prognosis are total vital capacity, FVC, and DLCO[17,18]. The decreases in FVC and DLCO in IPF patients reflect disease progression and predict mortality[19,20].

Previous studies on IPF still need to be extended, and further identification of additional gene-targeted therapies for IPF remains challenging. In this study, several novel genes were investigated to explore the treatment of IPF. We adopted quantitative real-time polymerase chain reaction (qRT-PCR) to identify common hub genes in multiple sets of IPF data and discussed the correlation between these genes and lung function and the overall survival rate. Afterwards, a hub gene, TDO2, was used as an example to demonstrate the expression level. This gene was increased in IPF and validated by western blotting. After inhibitor treatment, transforming growth factor-β (TGF-β)-induced fibroblast activation was effectively inhibited.

MATERIALS AND METHODS

Data collection and preprocessing

All raw data in this study were retrieved from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). GSE53845 (platform: Agilent microarray GPL6480) consists of lung tissues from 8 healthy subjects and 40 IPF patients[21]. In GSE47460[22], only samples with IPF from the interstitial lung disease population were included. The Agilent GPL14550 platform consists of 91 normal and 122 IPF samples, and the Agilent GPL6480 platform consists of 17 normal and 38 IPF samples. There were 110 men and 50 women with IPF, of whom 58 had never smoked, 96 had smoked, 2 were current smokers and 4 had an unclear smoking status. The age range of individuals with these IPF data was 37 to 82 years. GSE24206 (platform: Affymetrix GPL570) consists of 6 healthy and 17 IPF samples[23]. GSE110147 (platform: Affymetrix GPL6244) consists of 11 normal and 22 IPF samples[24]. The GEOquery package[25] was used to download the series matrix files of the databases above in R (v4.0.2). Soft formatted family files were downloaded to correctly map the probe ID to the gene symbol. GSE136831 contains lung tissue from 32 IPF and 28 control patients.

Differentially expressed gene screening

We screened differentially expressed genes (DEGs) between IPF patients and controls using the R package of the Microarray Data Linear Model (limma, version 3.50.3)[26]. The significant DEGs were identified according to the thresholds of adjusted P < 0.05 and fold change (FC) > 1.5. The common DEGs in the datasets were visualized by ggVennDiagram[27].

Functional and pathway enrichment analyses

Gene Ontology (GO) enrichment analyses of hub genes were performed by the clusterProfiler package with a background set of all Entrez IDs mapped to a GO pathway[28]. There are three functional categories, specifically biological processes (BPs), cell components (CCs), and molecular functions (MFs). Because GO is organized in a parent-child structure, a parent term can have a large percentage of overlap with its child terms. Concerning this problem, a simplified approach in the clusterProfiler implements was used. The organism reference was set as “org.Hs.eg.db”. The terms meeting P < 0.05, pAdjustMethod = "BH" and P adjust < 0.05 were selected. Enrichment analyses of the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome were performed using Database for Annotation, Visualization and Integrated Discovery (DAVID)[29]. After Benjamini-Hochberg adjustment, terms meeting P < 0.01 were selected.

Protein-protein interaction network establishment and hub gene exploration

The protein-protein interaction (PPI) network of identified DEGs was built using the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org, v11.5) to determine the functional interactions between the products of the key genes. The condition for constructing the PPI network was an interaction combined score > 0.4. The interaction information network was downloaded and visualized through Cytoscape software (v3.9.1). The Molecular Complex Detection (MCODE) plugin was applied to select the significant subnetworks from the PPI network (degree cut-off ≥ 2, node score cut-off ≥ 0.2, K-core ≥ 2, and max depth = 100, score ≥ 3)[30]. The hub genes were identified through a plugin from Cytoscape named cytoHubba[31]. The cytoHubba plugin quantifies the importance of nodes in biological interaction networks through 11 node ranking methods. Five of these methods were adopted in this study: Degree, maximum neighbourhood component (MNC), density of maximum neighbourhood component (DMNC), maximal clique centrality (MCC), and edge percolated component (EPC). The top 30 Hubba nodes were ranked by these methods, and common hub genes were taken as the true final hub genes.

Overall survival analyses

Disease status, survival days, age, and survival status were extracted from GSE70866[32]. Clinical data from 176 IPF-diagnosed cases were used to analyse survival changes in specific genes in IPF. According to the mRNA expression from one of the hub genes, IPF patients were divided into high-expression (upper tertile) and low-expression (lower tertile) groups. The survival analysis of hub genes was performed in R using the functions “survfit” and “coxph” of the “survival” package. Survival curves were generated using the “ggsurvplot” function of the “survminer” package. The final results were visualized by ggplot2[33].

Clinical correlation

Several factors, including age, sex, smoking history, carbon monoxide diffusion capacity (DLCO), and FVC, have been identified as predictors of poor survival in IPF patients[34]. Varying degrees of decline in FVC and DLCO predict higher mortality in patients with IPF[35]. In this study, clinical metadata were analysed using clinical datasets of GEO accession No. GSE47460[22]. Linear models were created for selected gene expression and clinical variables such as DLCO and FVC percentage using the lm function in the R stats package. Scatterplots were made using the ggscatter function in the ggpubr package. The Pearson correlation coefficient and significance were calculated in R using the stat_cor function in the ggpubr package.

Prediction of target microRNAs

We used three online microRNA (miRNA) databases (miRWalk, DIANATools, and miRDB) to predict miRNAs of hub genes with the default parameters. The miRNA was selected if it targeted a gene in all three databases. Subsequently, Cytoscape was used to construct the mRNA-miRNA coexpression network.

Transcription factors-miRNA prediction

Both miRNAs and transcription factors (TFs) are important in gene transcription and expression. Thus, clarifying the regulatory relationship between miRNAs and mRNAs and discovering key TFs can provide more insight into the pathological mechanisms of IPF. The "Enrichment analysis" module in TransmiR v2.0[36] was used to identify significant TFs that may regulate the list of miRNAs targeting TDO2.

Mouse model of Bleomycin-induced pulmonary fibrosis

SPF C57/6J mice were purchased from Beijing Weitonglihua Experimental Animal Co., Ltd. This experiment was approved by the Experimental Animal Welfare Ethics Review Committee of Henan University of Traditional Chinese Medicine (Review No. DWLL202110014). Mice were fed for one week in an SPF barrier environment in the Animal Experimental Center of Henan University of Chinese Medicine. Mice were randomly divided into two groups. For the phosphate-buffered saline (PBS) group, mice were intratracheally administered 50 μL of saline at Day 0. All mice were subjected to an anaesthesia experiment with a small animal anaesthesia machine (Ruiwode Life Technology Co., Ltd., Shenzhen, China). The light source was fixed at the pharyngeal skin, and the surgical instrument was fully exposed to the pharynx. At this time, a bright spot could be observed to open and close continuously with the spontaneous respiration of mice, and a tracheal intubation was inserted into the bright spot. For the Bleo group, mice were given intratracheal instillation of 50 μL of bleomycin (5 mg/kg). On Day 21, mice were euthanized by intraperitoneal injection of excessive pentobarbital sodium. Lung tissue was homogenized in TRIzol (Invitrogen) and stored at -80 °C for subsequent RNA isolation.

Cell culture and treatment

The A549 cell line and IMR-90 cell line were purchased from ATCC. A549 cells were cultured in F-12K medium (Boster Biological Technology) containing 10% foetal bovine serum (FBS), 100 U/mL penicillin, and 100 μg/mL streptomycin. IMR-90 cells were cultured in DMEM (Biological Industries) containing 10% FBS, 100 U/mL penicillin, and 100 μg/mL streptomycin. The MLE-12 cell line was purchased from the Shanghai Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences (Shanghai, China). MLE-12 cells were cultured in DMEM/F-12 medium (Corning) containing 10% FBS, 100 U/mL penicillin, and 100 μg/mL streptomycin. A549, MLE-12, and IMR-90 cells were incubated in a humidified incubator with 95% air and 5% CO₂ at 37 °C. A549 cells, MLE-12 cells, and IMR-90 cells were seeded in six-well plates and cultured overnight. Cells were treated with TGF-β (10 ng/mL) for 24 h[37] and then exposed to 680C91 (20 μm) for 24 h.

Quantitative real-time polymerase chain reaction assay

Primers were designed using PrimerBank (https://pga.mgh.harvard.edu/primerbank/index.html) and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). Total RNA was extracted using RNAiso Plus (TaKaRa) according to the manufacturer’s instructions. Reverse transcription (RT) was performed using HiScript^® II Q RT SuperMix for qPCR (Vazyme, Nanjing, China). The reactions were performed using a QuantStudio 6 real-time fluorescence quantitative PCR System (Life Technologies, Singapore). The procedure was detailed as follows: Lung tissue was homogenized on ice with TRIzol. After complete homogenization, the homogenized solution was allowed to settle at room temperature for 5 min. The homogenate was then mixed with 200 μl/mL chloroform, allowed to stand at room temperature for 3 min and then centrifuged at 12000 × g for 15 min at 4 °C to separate the phases. The aqueous RNA phase was removed and dispensed into a new 2 mL microcentrifuge tube. The upper aqueous phase (RNA solution) was transferred to a new EP tube, mixed with isopropyl alcohol equal to the supernatant, and placed at 12000 × g for 10 min at 4 °C, and the supernatant was discarded. Then, 75% ethanol was added to 1 mL/mL TRIzol and gently shocked at 4 °C and 12000 × g for 10 min, and the supernatant was discarded. The pellets were allowed to air dry for approximately 5 min and then dissolved in DEPC water to measure the RNA concentration. The TOYOBO ReverTra Ace qPCR RT Kit was used for RT. RNA was incubated for 10 min at 65 °C and immediately placed on ice. The RT reaction contained 1 μg of RNA, 2 μL of 5 × RT Buffer, 0.5 μL of RT Enzyme Mix, 0.5 μL of Primer Mix, and up to 10 μL of DEPC water and was subsequently incubated at 37 °C for 15 min and 98 °C for 5 min. The values of the target genes were normalized using the 2^−ΔΔCT method with the values of the housekeeping gene 18S. Technical repeats were performed three times. Notably, all mRNA levels are expressed as the mean ± SD.

Western blot assay

After treatment with TGF-β and 680C91, cells were lysed with RIPA buffer in ice. Protein samples of equal concentrations were separated using a 10% SDS-PAGE gel and electrotransferred onto PVDF membranes. Membranes with proteins were blocked with 5% skim milk, followed by incubation with primary and secondary antibodies. TDO2 (Cat No.:15880-1-AP), α-smooth muscle actin (α-SMA) (Cat No.:14395-1-AP), COL-I (Cat No.:14695-1-AP), and GAPDH (Cat No.:10494-1-AP) were detected using the Bio-Rad Imaging System. These antibodies were purchased from Proteintech Group, Inc.

Statistical analysis

Independent-samples t tests were performed using SPSS 22.0 (IBM Corporation, Armonk, NY, United States), and P < 0.05 was considered significant. Pearson’s correlation coefficient r was calculated using the R function stat_cor.

RESULTS

Identification and enrichment analysis of DEGs

Four sets of raw GEO data (GSE53845, GSE47460, GSE24206, and GSE110147) from five platforms were used to analyze the DEGs between the normal and IPF groups. Data sources are summarized in Supplementary Table 1. After background correction and data standardization, DEGs were determined between IPF and normal lung tissues. DEG analysis was performed by the R package “limma”. Based on the threshold of adjusted P < 0.05 and FC > 1.5, 2157 DEGs were identified in GSE53845, 2861 DEGs were screened out from GSE47460 under the GPL6480 platform, 2113 DEGs were obtained in GSE47460 under the GPL14550 platform, 2525 DEGs were detected in GSE24206, and 7146 DEGs were selected in GSE110147. Venn analysis showed that the DEGs of GSE110147 had less than 50% consistency with other datasets (Supplementary Figure 1). Therefore, GSE110147 was discarded in subsequent analyses. A total of 254 common genes were screened by intersecting the remaining four DEG datasets (Figure 1A).

Figure 1 Identification and enrichment analysis of differentially expressed genes in idiopathic pulmonary fibrosis data compared with control data. A: Analysis of common differentially expressed genes was performed between idiopathic pulmonary fibrosis and normal lung tissues with thresholds of P < 0.05 and fold change > 1.5 based on GSE53845, GSE47460 (GPL14550 and GPL6480), and GSE24206; B-D: Three categories of Gene Ontology analysis: Biological processes, cell components, and molecular functions; E: Kyoto Encyclopedia of Genes and Genomes and Reactome pathway analysis.

The biological roles and functional pathways of the 254 DEGs were revealed by GO enrichment analysis. The three functional categories were BPs (Figure 1B), CCs (Figure 1C), and MFs (Figure 1D). The results showed that DEGs were significantly involved in “extracellular matrix organization”, “response to transforming growth factor beta”, “regulation of insulin-like growth factor receptor signalling pathway”, “regulation of neuron maturation”, “collagen-containing extracellular matrix”, “extracellular matrix structural constituent”, “serine-type endopeptidase activity”, “calcium-activated potassium channel activity”, and “Wnt-protein binding”. The results of the KEGG and Reactome pathway analyses (Figure 1E) identified by DAVID showed that the DEGs were significantly correlated with “extracellular matrix organization”, “collagen degradation”, “regulation of insulin-like growth factor transport and uptake by insulin-like growth factor proteins”, and “molecules associated with elastic fibres”.

PPI network analysis and identified hub genes

For further elucidation of the interaction of common DEGs, interaction networks were constructed using STRING (STRING v11.5, https://string-db.org/). The TSV file containing information on PPI network interactions was downloaded and visualized by Cytoscape (V3.9.1). The network included 69 nodes and 59 edges, with an average number of neighbours of 2.0. The MCODE plugin in Cytoscape was used to detect the most significant hub gene modules and cluster scores under the default parameters. Three modules were identified with the MCODE plugin (Supplementary Figure 2) according to the following filter criteria: Cluster 1 (score: 4.000, 4 nodes and 6 edges), Cluster 2 (score: 4.000, 4 nodes and 6 edges), and Cluster 3 (score: 3.000, 3 nodes, and 3 edges). Afterwards, five algorithms (MCC, DMNC, MNC, Degree, and EPC) were used to identify hub genes in the cytoHubba plugin (Figure 2A-E). The top 30 genes in each algorithm were considered hub genes (Supplementary Tables 2-6). Finally, 15 genes were detected in all five algorithms: RET, TGFB3, HSD17B6, SULF1, ASPN, SDR16C5, ALDH1A3, COL3A1, COL1A2, KCNMA1, COL5A2, KCNMB4, AOX1, KCNN3, and KCNN4 (Figure 2F).

Figure 2 Hub genes identified by maximum neighbourhood component, Degree, density of maximum neighbourhood component, maximal clique centrality, edge percolated component, venn diagram of hub genes. A: Maximum neighbourhood component; B: Degree; C: Density of maximum neighbourhood component; D: Maximal clique centrality; E: Edge percolated component; F: Venn diagram. The 15 hub genes were RET, TGFB3, HSD17B6, SULF1, ASPN, SDR16C5, ALDH1A3, COL3A1, COL1A2, KCNMA1, COL5A2, KCNMB4, AOX1, KCNN3, and KCNN4.

Correlation between gene expression and lung functions

The above analysis showed significant enrichment of genes in epithelial cell proliferation. We expect a variety of candidate targets that can be used to treat IPF. Combined with the results of GO molecular functional enrichment analysis, these findings showed that tryptophan 2,3-dioxygenase activity and amino acid binding were significantly enriched (P < 0.05), and TDO2 was involved in these pathways. Figure 1D does not show this pathway because it was not ranked high enough. Moreover, in conjunction with the query results from DGIdb (https://dgidb.org/search_categories), which indicated that TDO2 has a potential drug effect, these data were added to the downstream analysis. Moreover, TDO2 has been confirmed to promote tumour cell proliferation and differentiation in the progression of oesophageal squamous cell carcinoma[38]. No relevant studies have been reported on the role of TDO2 in IPF. Therefore, 15 genes identified by PPI and TDO2 were included in the subsequent analysis.

The potential role of these genes in IPF development was examined by tracking expression changes in public datasets. We used linear models to assess changes in gene expression and lung function. The results showed that changes in 16 genes were directly associated with lung function. Twelve (TDO2, RET, TGFB3, SULF1, ASPN, ALDH1A3, COL3A1, COL1A2, KCNMA1, COL5A2, KCNN3, KCNN4) of the 16 genes were upregulated in IPF patients in all datasets. Additionally, these 12 genes significantly decreased in terms of lung function (decrease in DLCO and FVC, as shown in Figure 3A-L and Supplementary Figure 3A-L respectively). The changes in AOX1 expression were conflicting. We observed a downregulation in IPF in the GSE53845 and GSE24206 data (logFC of -0.94 and logFC of -1.27, respectively) and an upregulation in IPF in the GSE47460 data (logFC of 0.9). Therefore, although AOX1 is negatively correlated with DLCO and FVC (Figure 3M and Supplementary Figure 3M), we have excluded this gene from the subsequent analysis. Three (HSD17B6, SDR16C5, KCNMB4) of 16 genes were downregulated in IPF patients in all datasets, revealing a decrease in lung function (decreased DLCO and FVC, as shown in Figure 3N-P and Supplementary Figure 3N-P respectively).

Figure 3 Correlations between lung gene expression and forced vital capacity and diffusing capacity. A-M: TDO2, RET, TGFB3, SULF1, ASPN, ALDH1A3, COL3A1, COL1A2, KCNMA1, COL5A2, KCNN3, KCNN4 and AOX1 expression was negatively correlated with diffusing capacity of carbon monoxide (DLCO); N-P: HSD17B6, SDR16C5, and KCNMB4 were significantly positively correlated with DLCO. The Pearson correlation metric was computed by using the ‘stat_cor’ function in R. DLCO: Diffusing capacity of carbon monoxide.

Survival analysis

To explore the prognostic value of hub genes in IPF, we performed a survival analysis by Cox models. The data were from GSE70866[32]. Sample grouping information and the average expression and difference statistics of hub genes in each group are shown in Table 1. Survival analysis showed that TDO2, HSD17B6, SULF1, SDR16C5, COL3A1, KCNMA1, KCNMB4, KCNN4, and ALDH1A3 were independent prognostic factors for poor survival in IPF (Figure 4). Patients with high expression of these genes are likely to suffer unfavourable outcomes. There was no significant difference between the expression of other hub genes and survival in IPF patients.

Figure 4 Survival analysis of overlapping hub genes. Overall survival curves show that higher expression of TDO2, HSD17B6, SULF1, SDR16C5, COL3A1, KCNMA1, KCNMB4, KCNN4 and ALADH1A3 is associated with a worse survival rate. A: TDO2; B: HSD17B6; C: SULF1; D: SDR16C5; E: COL3A1; F: KCNMA1; G: KCNMB4; H: KCNN4; I: ALADH1A3.

Table 1 Primer sequences.

Gene name	Forward	Reverse
TDO2	AACATGCTCAAGGTGATAGCTC	GAACCGAGAACTGCTGTACCA
SULF1	TGTGTTCCACCGTTCGGTC	CACATCCTGGTCGTCAGTGAG
KCNN4	GCTCAACCAAGTCCGCTTC	GTGATCGGAATCAGCCACAGT
COL3A1	CTGTAACATGGAAACTGGGGAAA	CCATAGCTGAACTGAAAACCACC
HSD17B6	GGAGCGTGTTGGAGACAGAG	GAGGTTCACTTGAAAGATAGGCA
KCNMA1	TCACGGAACTCGCTAAGCC	AATGTGCGTCCCACTGTTTTT
KCNMB4	ACCAACCCCAAGTGCTCCTAT	GAATGGCTGGGAACCGATCTC
ALDH1A3	ATCAACAACGACTGGCACGAA	CACATCGGGCTTATCTCCTTC
SDR16C5	TTGAGTGTTTTGGAGGCCCTA	AACTGCAATGCTAAGAGCCTT
18S	GTAACCCGTTGAACCCCATT	CCATCCAATCGGTAGTAGCG

Validation of hub gene expression by qRT-PCR

The mRNA expression levels of the hub genes associated with low survival rates were detected through qRT-PCR. The sequences of the primers are shown in Table 1. The results showed that the mRNA expression levels of six genes (TDO2, COL3A1, KCNN4, SULF1, KCNMA1, ALDH1A3) were significantly increased in the bleomycin-induced pulmonary fibrosis group compared with the control group (P < 0.01, Figure 5). The direction of gene dysregulation was consistent with the test data analysis. These results demonstrate the validity of these genes for further analysis.

Figure 5 mRNA expression levels of TDO2, COL3A1, KCNN4, SULF1, HSD17B6, KCNMA1, KCNMB4, ALDH1A3, and SDR16C5 in the lung tissue vs PBS (n = 5 for PBS group, n = 5 for Bleo group). The statistical test used was the t test. Data are expressed as the mean with SD. ^aP < 0.01, ^bP < 0.05.

Prediction of target miRNAs and identifying genes of interest

Previous studies have shown that miRNAs are widely involved in disease regulation by targeting key genes[39]. The typical function of miRNAs is to bind to their 3′ untranslated regions to regulate various target genes. TFs can upregulate or downregulate miRNAs, corresponding to positive and negative feedback phenomena. Identification of the TF-miRNA-mRNA coregulatory network can help determine which miRNAs play a crucial role in the pathogenesis of hub genes involved in IPF and how they interact with TFs to regulate IPF-related genes. We used three online miRNA databases (miRWalk, DIANATools, and miRDB) to predict miRNAs targeting the above nine hub genes. The targeted miRNA was selected when it targeted a gene in all three databases. Finally, 367 miRNAs were selected, and 401 mRNA-miRNA pairs were obtained (Supplementary Table 7). The interaction network of mRNAs and miRNAs, comprising 376 nodes and 401 edges, was constructed by Cytoscape (Supplementary Figure 4). Finally, TDO2 was selected to further explore the regulatory network.

TFs and miRNAs regulate each other to form feed-forward loops or feedback loops, where a TF regulates a miRNA or a miRNA inhibits a TF[40,41]. In this study, TransmiR was adopted to predict a list of significant TFs that may regulate miRNAs, and it provided 3730 highly reliable TF-miRNA regulations[36]. TDO2 was regulated by seven miRNAs (hsa-let-7g, hsa-mir-6878, hsa-mir-4270, hsa-mir-4441, hsa-mir-7974, hsa-let-7a-2, hsa-mir-6754). Based on the TFs in TransmiR prediction with a threshold of P < 0.05, the number of targeted miRNAs for regulation was 7. Nine molecules corresponding to this condition were TCF12, SPI1, HIF1A, RUNX1, ETS-related gene (ERG), MYC, EP300, RELA, and STAT1. On this basis, the TF-miRNA-mRNA network of TDO2 was constructed by Cytoscape (Figure 6).

Figure 6 Transcription factors-miRNA-mRNA regulatory network. Transcription factors, miRNA, and mRNA are represented by hexagons, circles, and diamonds, respectively.

TDO2 in IPF and alveolar fibroblasts

TDO2 was upregulated in IPF in all trained datasets, as shown in Figure 7A. TDO2 expression was upregulated in IPF, while lung function was decreased (decreased DLCO and FVC, Figure 7B and C and Supplementary Figure 5A). In addition, single-cell RNA-seq data containing 32 IPF and 28 control lung tissues (GSE136831)[42] were used to validate TDO2 expression between IPF and normal lung tissues of patients. The parameters and cell type identification in UMAP were consistent with those in a study of CD38-mediated pulmonary fibrosis[42]. In addition, TDO2 expression was significantly enriched in alveolar fibroblasts (Figure 7D-F). This result is consistent with the proliferation of fibroblasts and myofibroblast accumulation in IPF. The above results were also supported by GSE135893 and GSE122960 (Supplementary Figure 5B-E). A UMAP plot with cells labeled by disease identity in GSE136831 has been added to Supplementary Figure 5F.

Figure 7 TDO2 increased in Idiopathic pulmonary fibrosis and was validated in scRNA-seq data. A: Upregulation of TDO2 in Idiopathic pulmonary fibrosis of trained data (e.g., in GSE47460). The statistical test used was the t test; B and C: Correlations between lung TDO2 expression and diffusing capacity of carbon monoxide (B) and forced vital capacity (C) in GSE47460. The Pearson correlation metric was computed by using the ‘stat_cor’ function in R; D: UMAP of 8,942 randomly selected resident lung parenchymal cells from GSE136831; E: UMAP with cells labelled by normalized TDO2 expression; F: Differential gene expression analysis using the log-fold change expression vs the difference in the percentage of cells expressing the gene comparing alveolar fibroblast cluster vs other cluster cells (Δ Percentage Difference). Labelled genes have a log2-fold change > 1, Δ Percentage Difference > 20% and adjusted P value from Wilcoxon rank sum test < 0.05. DLCO: Diffusing capacity of carbon monoxide; IPF: Idiopathic pulmonary fibrosis; FVC: Forced vital capacity.

Reduction of TGF-β-induced fibrosis by TDO2 deficiency

To identify the effect of altered TDO2 gene expression on fibrosis, we treated cells with TDO2 inhibitors. Cells treated with TGF-β were the pulmonary fibrosis model group. After 24 h of TGF-β treatment, cells were treated with the TDO2 inhibitor 680c91 to explore fibrotic changes. The same experiments were performed in three cell lines, including fibroblasts (IMR-90) and epithelial cells (A549 and MLE-12). The attenuated fibrosis was supported by the decrease in α-SMA and COL-I expression (Figure 8). These in vitro data reveal that decreased TDO2 expression inhibited TGF-β-induced fibroblast activation, indicating that TDO2 is a potential therapeutic target for fibrosis.

Figure 8 Western blot analysis of α-smooth muscle actin and COL-1 in A549, MLE-12, and IMR-90 cells. Western blot quantification results in A549, MLE-12 and IMR-90 cells. A and D: A549; B and E: MLE-12; C and F: IMR-90. In line with the homogeneity of variance, the least significance difference method was used when assuming homogeneity of variance. ^aP < 0.01 vs phosphate-buffered saline; ^bP < 0.05 vs transforming growth factor-β. PBS: Phosphate-buffered saline; α-SMA: α-smooth muscle actin; TGF-β: Transforming growth factor-β.

DISCUSSION

The progression of IPF seriously affects the lives of patients, while effective treatments are limited. To address the shortage of IPF-targeted therapeutic candidate genes, this study collected four sets of GEO public data to find common DEGs between IPF and normal lung tissues. Information on the common DEGs has been added to Supplementary Table 8. Since screening large public data ensures more accurate and generally applicable hub genes, data containing extensive IPF and control samples were selected for screening. For more robust data results, further screening was performed with clinical lung function indicators and survival analyses. Based on the PPI network and five cytoHubba algorithms, 15 hub genes were identified for the analysis. The relationship of these genes with lung function (such as DLCO and FVC) and overall survival was examined. Nine genes (TDO2, HSD17B6, SULF1, SDR16C5, COL3A1, KCNMA1, KCNMB4, KCNN4, and ALDH1A3) significantly decreased lung function in IPF and caused poor survival. The selected genes were validated by scRNA-seq data, qRT-PCR and western blotting in three cell lines, indicating these results are consistent with the actual situation. Compared with the validation in a single cell line, the results in multiple cell lines provide better support for the therapeutic function of TDO2 inhibitors.

GO analysis showed that the common DEGs were primarily enriched in extracellular matrix organization (BP), collagen-containing extracellular matrix (CC), and extracellular structure constituent (MF). These biological pathways are involved in the degradation, regeneration, and remodelling of fibrosis[43]. IPF is the most common idiopathic interstitial pneumonia. However, nothing, other than a lung transplant, has thus far increased survival[44]. The antifibrotic drugs on the market are pirfenidone and nintedanib, which are used to slow the progression of the disease. Thus, further exploration of other therapeutic drugs is warranted. Among the screened genes, SULF1 promotes histone H4 acetylation, enhances the effects of HDAC inhibitors, and inhibits tumorigenesis in hepatocellular carcinoma[45]. HDAC inhibitors can downregulate COL3A1 expression in primary IPF lung fibroblasts[46]. Therefore, targeting histone deacetylases in IPF may also be a future therapeutic approach[15].

MiRNAs are endogenous noncoding RNAs with regulatory functions and a length of approximately 22 nt[47]. These molecules play important roles in regulating the expression of genes related to the growth and development of organisms and the occurrence of diseases. Studies have shown that each miRNA may have dozens or hundreds of target genes[48]. With the continuous discovery of specific miRNAs, the study of target genes in some disease pathways may lead to the discovery of new therapeutic methods. In this study, a regulatory pathway that may have the potential to regulate the expression of TDO2 was identified through the TF-miRNA-mRNA network. In this way, the underlying molecular mechanisms of IPF were revealed, providing potential targets for developing new therapeutic agents to treat IPF patients. The results showed nine TFs targeting seven miRNAs to coregulate TDO2. Among them, HIF1A induced overexpression of the p53/hypoxia pathway in the tissues of IPF patients[49]; EP300 was highly enriched at enhancers of several genes in multiple profibrotic pathways[50]; immunostaining of RELA (v-rel reticuloendotheliosis viral oncogene homologue A) was increased in type 2 alveolar epithelial cells of mice treated with bleomycin[51]. The function of the TF ERG is dysregulated in ageing. A previous study showed that ERG dysfunction exacerbated pulmonary vascular ageing and persistent fibrosis[52]. We speculated that the disorder of ERG caused positive feedback, such as mechanisms involving miR-4270 and miR-4441 to regulate the expression of TDO2, so specific inhibition of ERG could be assessed in follow-up research, and combined treatment may be the future direction of pulmonary fibrosis treatment.

This study also has some limitations. We used five common sets from IPF data, but the intersection of common DEGs between one group of datasets (GSE110147) and the other groups was less than 50%. This group was removed in subsequent analysis, and the DEGs in the remaining groups were screened using thresholds of adjusted P < 0.05 and FC > 1.5. The selected common DEGs had identical direction of dysregulation, indicating that the validity of the thresholds and the data are reasonable. With the progress of scientific research, more IPF-related data can be generated, and additional datasets or novel computing methods can be adopted to find relevant genes affecting IPF. In addition to human lung tissue data, mouse data can be analysed and verified. With the diversification of sequencing technologies, multiomics data analysis can be adopted to validate genes. For example, single-cell RNA sequencing technology can examine the expression of specific genes in different cell populations and the existence of tissue preferences. Some genes were differentially expressed in at least one cell type, while some genes related to ECM were largely expressed in fibroblast subpopulations[53]. By investigating the expression of TDO2 in each cell cluster in scRNA-seq data (GEO series accession No. GSE135893 and GSE122960[54]), we found higher expression in fibroblasts than in other cell clusters. The specific mechanism of TDO2 in fibroblasts should be further investigated.

Moreover, when using STRING to construct PPI networks, we used different interaction score thresholds to obtain different network structures. Important candidate genes in common DEGs might be missed, and these genes would not appear in the network due to insufficient research. We reviewed these genes again and found that TDO2 has potential drug effects. Previous research has shown that TDO2 knockdown inhibits colorectal cancer progression[55]. In vitro experiments have shown that knocking down or blocking TDO2 expression in myofibroblast subsets can effectively reverse T-cell immunosuppression in oral squamous cell carcinoma[56]. No relevant studies have been reported on the role of TDO2 in IPF. Therefore, we included TDO2 for follow-up analysis. Fortunately, the results of the network analysis in this study validated independent data from public databases and qRT-PCR data from our laboratory. However, the disadvantage is that we performed the experiment once, and there is a lack of repeated experiments to prove the repeatability of these findings. Additionally, it is important to explore whether the other group receiving more or less bleomycin has any effect on the reproducibility of the results. Representative images of rat body weight, survival rate, lung wet-to-dry weight ratio, lung hydroxyproline content, mouse hematoxylin and eosin staining and Masson’s trichrome staining of lung sections are more conducive to increasing the persuasiveness of our results, which will be further studied in the future. In summary, we demonstrated that the TDO2 gene is significantly upregulated after bleomycin treatment and proposed blocking TGF-β production as a potential treatment for IPF. A more comprehensive study of TDO2 in animal models would be worthwhile.

CONCLUSION

This study identified novel hub genes to explore the treatment of IPF, examined the upregulated expression of TDO2, and investigated the TF-miRNA-gene regulatory network of TDO2. After inhibitor treatment, TGF-β-induced fibroblast activation was effectively inhibited. TDO2 appears to be an effective treatment for IPF. Therefore, the molecular mechanisms of these TFs and TDO2 deserve further exploration. The findings provide new target gene prediction, and we propose blocking TGF-β production as a potential treatment for IPF.

ARTICLE HIGHLIGHTS

Research background

Idiopathic pulmonary fibrosis (IPF) is a progressive interstitial lung disease with high mortality rate. Therefore, exploring potential therapeutic targets to meet the unmet needs of IPF patients is of great significance.

Research motivation

To explore potential targets for IPF treatment, as well as potential pathways and therapeutic methods.

Research objectives

Explore novel hub genes for IPF therapy. Whether TDO2 can be a therapeutic target for IPF.

Research methods

We used public datasets (GSE53845, GSE47460, GSE24206, and GSE110147) to identify differentially expressed genes between IPF patients and healthy donors. Potential targets were considered based on multiple bioinformatics conditions, especially the correlation between hub genes and carbon monoxide diffusing capacity, forced vital capacity, and patient survival rate. The mRNA levels of the hub genes were determined through quantitative real-time polymerase chain reaction. Transforming growth factor-β (TGF-β) induced pulmonary fibrosis mouse model and the expression of TDO2 was observed before and after the addition of inhibitor.

Research results

This study identifies novel hub genes to explore for IPF treatment. TDO2 was upregulated in an experimental mouse model of TGF-β-induced pulmonary fibrosis and a TDO2 inhibitor effectively suppressed TGF-β-induced fibroblast activation.

Research conclusions

TDO2 could be a potential target for treatment of IPF.

Research perspectives

More complete work with animal models is needed, as the TDO2 gene is apparently upregulated by bleomycin treatment.