Abstract

Idiopathic pulmonary fibrosis (IPF), the most frequent form of irreversible interstitial pneumonia with unknown etiology, is characterized by massive remodeling of lung architecture and followed by progressive loss of lung function. However, the key regulatory genes and the specific signaling pathways involved in the onset and progression of IPF still remain unclear. The present study is aimed at investigating the key role of long noncoding RNAs (lncRNAs) and transcription factors (TFs) involved in the pathogenesis of IPF through the integrated analysis of three gene expression profiles from the GEO dataset (GSE2052, GSE44723, and GSE24206). A total of 8483 differentially expressed genes (DEGs) including 988 upregulated and 7495 downregulated genes were filtered. Subsequently, following the intersection of these DEGs, 29 overlapping genes were identified and further analyzed using a bioinformatics approach. Furthermore, the protein-protein interaction (PPI) network was used to obtain 18 modules of related genes. The hub genes were identified through hypergeometric testing, which were closely associated with ubiquitin-mediated proteolysis, the spliceosome, and the cell cycle. The significant difference was observed in the expression of these key genes, such as lncRNA MALAT1, E2F1, and YBX1, in the peripheral blood of IPF patients when compared with those normal control subjects by real-time polymerase chain reaction (RT-PCR) analysis. This study indicated that lncRNA MALAT1, E2F1, and YBX1 may be key regulators for the pathogenesis of IPF.

1. Introduction

Idiopathic pulmonary fibrosis is a chronic and progressive lung tissue damage of unknown etiology, which is characterized by the abnormal proliferation of activated fibroblasts/myofibroblasts and excessive deposition of collagen in the extracellular matrix (ECM) from adjacent alveoli to the lung parenchyma. IPF has a poor prognosis and high mortality rate with the postdiagnosis median survival rate of only 20% to 30% and the median survival of approximately 3 to 5 years [1, 2]. Due to the complexity and heterogeneity of IPF, its incidence and mortality rate, which has a positive relationship with advanced age, have shown a steadily increasing trend worldwide [3]. Although the pharmacotherapy of IPF has made certain progress over the past 5 years, the therapeutic efficacy is unsatisfactory because of the variable and unpredictable course of IPF and large individual differences [4].

Increasing studies related to transcriptome, including both protein-coding mRNAs and noncoding RNAs (ncRNAs), have provided novel insights into the molecular mechanism of IPF pathogenesis. Among them, ncRNAs implicated in multiple fibrotic diseases have been divided into short and long ncRNAs (lncRNAs) based on its length of nucleotide sequences. Multiple studies have shown that lncRNAs (≥200 nucleotides) contribute to the pathogenesis and progression of IPF and gain more attention [5, 6]. However, varying proportions of transcripts that can be detected and the accuracy of measurements of changes in low-abundance transcripts reduce detection accuracy of lncRNAs in transcriptome-related lung fibrosis research. In addition, transcripts detected and measured are very different in different microarray platforms. These factors imply that some lncRNAs may be overlooked and false-positive or false-negative results may be generated [7, 8]. Based on publicly available microarray expression datasets in the Gene Expression Omnibus (GEO) database, an in-depth bioinformatics analysis of lncRNAs may provide a comprehensive understanding of not only transcriptional regulation but also posttranscriptional regulation. Hence, using bioinformatics methods to analyze the comprehensive gene network, this study was performed to identify the biological processes and pathways of differentially expressed genes (DEGs) that are involved in the pathogenic mechanism of IPF. These results may be useful in elucidating the critical regulatory mechanism of IPF from a systematic perspective and providing the relevant effective interventions to attenuate or reverse the process of lung fibrosis.

2. Materials and Methods

2.1. Microarray Data Information

NCBI-GEO (http://www.ncbi.nlm.nih.gov/geo/) is a free database repository comprising microarray/gene profile, next-generation sequencing, hybridization array, and chip data. All data were derived from GEO datasets GSE2052, GSE44723, and GSE24206. The microarray data of GSE2052 were based on GPL1739 Platforms (Amersham Biosciences CodeLink Uniset Human I Bioarray, University of Pittsburgh, PA, USA) and included 15 IPF and 11 control lung tissues (submission date: 09 December 2004) [9, 10]. The GSE44723 data were based on GPL570 Platforms (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA) and included 10 pulmonary fibrosis and 4 normal lung tissues (submission date: 10 April 2013) [11]. The GSE24206 data were based on GPL570 Platforms (Affymetrix Human Genome U133 Plus 2.0 Array, Affymetrix, Santa Clara, CA, USA) and included 17 IPF and 6 normal lung tissues (submission date: 01 November 2011) [12]. The total RNA of the samples was extracted to analyze the genomic profile of the RNA. All data came from expression profiling with microarrays conducted for Homo sapiens.

2.2. Identification of Differential Gene Expression in IPF

The original data from these datasets, including SOFT-formatted family files and Series Matrix Files, were downloaded for analysis. DEGs were identified with the R package limma (http://bioconductor.org/packages/release/bioc/html/limma.html). Unsupervised hierarchical clustering was performed to center the normalized and log2-scaled expression values on the median by using Cluster 3.0 (Fig. S1–S3). After pretreatment of the genes that came from more than one probe set, the DEGs identified with cutoff criteria of and by the classic -test were considered statistically significant.

2.3. Gene Ontology and KEGG Enrichment Analysis of DEGs

Functional and pathway enrichment analyses of candidate DEGs were performed with the online bioinformatics Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.ncifcrf.gov) (version 6.7), which can integrate biological data and comprehensively annotate the biological functional information of genes. Gene Ontology (GO) analysis can provide annotation of DEGs regarding biological processes (BPs), molecular functions (MFs), and cellular components (CCs) and allows further analysis of the bioprocesses of these genes. The Kyoto Encyclopedia of Genes and Genomes (KEGG) provides high-level functions and biological system information derived from large-scale molecular datasets generated with high-throughput experimental technologies. DAVID was applied to analyze the function of DEGs, and a value of less than 0.05 was considered statistically significant.

2.4. Construction of Protein-Protein Interaction (PPI) Network and Module Analysis

Functional analysis of interactions between the candidate DEG-encoded proteins can provide a new perspective on the pathogenesis and development of IPF. The protein-protein interaction network (PPI) of DEGs was constructed with the Search Tool for the Retrieval of Interacting Genes (STRING) online database (http://string-db.org) (version 11.0) considering combined scores of interaction greater than 0.4 to indicate statistical significance, and the network was visualized in the form of modules by using ClusterONE Cytoscape plug-in (version 1.0) [13]. Cytoscape (version 3.6.1) is powerful bioinformatics software that is utilized to visualize molecular interaction networks. Then, GO and KEGG enrichment analyses of genes in the module were conducted by using DAVID.

2.5. Selection of the Key Genes

Using the Molecular Complex Detection (MCODE) (version 1.4.2) plug-in of Cytoscape, the hub genes were selected by means of clustering the dense connection domain based on the topology of a given network [14]. The GO and pathway enrichment analyses of hub genes were performed with the ClueGO (2.5.1) plug-in of Cytoscape [15]. Subsequently, the biological pathway relationship network of these hub genes was constructed with the Biological Networks Gene Ontology tool (BiNGO) (version 3.0.3) plug-in of Cytoscape [16]. Using the hypergeometric test of the empirical Bayes approach, key genes were obtained through calculation. A value of less than 0.05 was considered statistically significant.

2.6. Subjects and Blood Samples

Considering the particular and complex nature of IPF, not all patients can undertake the invasion operation including bronchial and surgical lung biopsies to obtain the lung tissue samples. Moreover, obtaining healthy control samples would be not only extremely difficult but also restricted by ethical concerns. Due to the feasibility and convenience of obtaining blood samples, we validated the expression levels of candidate genes in the peripheral blood samples of all subjects. IPF patients () were diagnosed at the Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University. Healthy physical examinees () were selected as the control group. The cohort of 40 subjects provided written informed consent in compliance with the code of ethics of the World Medical Association. The collection and usage of the blood samples were approved by the Medical Research Ethics Committee of Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University (the Scientific Research Project 2018XE0109-1).

2.7. Validation of Key Genes by RT-PCR Analysis

Purification of RNA from blood samples from 20 IPF patients and 20 normal control subjects was performed using the TRIzol™ LS Reagent (Invitrogen, USA). The RNA was reverse-transcribed using the PrimeScript™ RT reagent kit with ɡDNA Eraser (TAKARA, Japan) according to the manufacturer’s recommendations. The cDNA from each sample was used as a template with GAPDH as an internal reference. The specific primer sequences used to amplify the 4 key candidate genes are listed in Table S1. Real-time PCR (RT-PCR) was performed using the StepOnePlus™ Real-Time PCR System (Thermo Fisher Scientific, USA). The results are represented as the means of 3 repetitions and were quantified via the 2-ΔΔct method. The mRNA levels of key genes between the IPF and normal lung tissues were compared using a paired -test () using GraphPad 6.0 (GraphPad Software, La Jolla, CA, USA). The data were presented as the (Table S2). Counting data was assessed using a test. Multiple-group comparison was assessed using one-way analysis of variance (ANOVA) followed by the Bonferroni multiple comparison test. Comparison of two groups was assessed by a two-tailed -test.

3. Results

3.1. Identification of Differentially Expressed Genes in IPF

After normalization and standardization of the raw data from these three GSE2052, GSE44723, and GSE24206 datasets (Figures 1(a)1(c)), we identified a total of 8483 aberrantly expressed genes, including 988 upregulated and 7495 downregulated genes, in IPF tissues compared to normal lung tissues (Table 1). There were 29 overlapping genes between the GSE2052, GSE24206, and GSE44723 datasets according to the Venn diagram, including 29 overlapping genes between GSE2052 and GSE44723, 268 overlapping genes between GSE44723 and GSE24206, and 389 overlapping genes between GSE2052 and GSE24206 (Figure 1(d)).

3.2. GO and KEGG Enrichment Analyses of Differentially Expressed Genes

The biological processes associated with the DEGs were determined by using the DAVID online bioinformatics database. As shown in Table 2, the top 6 GO results revealed that the significantly enriched BPs of IPF DEGs were mainly concentrated in the cell adhesion, biological adhesion, regulation of cell proliferation, and so on. The top 6 significantly enriched MFs were mainly concentrated in the calcium ion binding, cytokine binding, chemokine activity, and so on. The top 6 significantly enriched CCs were mainly concentrated in the extracellular region part, extracellular space, extracellular matrix, and so on. The top 6 significantly enriched KEGG pathways were mainly concentrated in the ECM-receptor interaction, cytokine-cytokine receptor interaction, and so on.

3.3. Construction and Enrichment Analysis of Modules

The IPF DEGs were used to construct the PPI network by using the STRING online database, and a total of 18 modules were obtained with the ClusterONE Cytoscape plug-in (Figure 2). To obtain functional and pathway enrichment information, the genes involved in these 18 modules were further analyzed by using DAVID as shown in Figures 3(a)3(d). The top 10 modules of significantly enriched BPs were mainly concentrated in the protein polyubiquitination, ciliary basal body-plasma membrane docking, Golgi vesicle transport, and so on. The top 10 modules of significantly enriched CCs were mainly concentrated in the ubiquitin ligase complex, microtubule organizing center part, microtubule-associated complex, and so on. The top 10 modules of significantly enriched MFs were mainly concentrated in the ubiquitin-protein transferase activity, structural constituents of the cytoskeleton, microtubule motor activity, and so on (Table 3). The KEGG pathways of these 18 modules were concentrated in the ubiquitin-mediated proteolysis, spliceosome, purine metabolism, Fanconi anemia pathway, and so on (Table 4).

3.4. Selection and Analysis of Key Genes

The biological network of differentially expressed IPF genes was constructed by using the BiNGO plug-in of Cytoscape, and the results revealed that most of the DEGs were significantly enriched in mitochondrial translation, cellular macromolecule metabolic process, cellular process, and so on (Figure 4(a)). ClueGO, another plug-in of Cytoscape, can annotate and visualize the pathway networks of DEGs integrating GO terms as well as KEGG pathways. The results from ClueGO revealed that most of the DEGs were significantly enriched in the glutathione metabolism, Fanconi anemia pathway, etc. (Figure 4(b)).

Subsequently, the key genes were obtained through calculation of the hypergeometric test. The 30 miRNAs and 4 lncRNAs enriched in 13 modules and the 44 transcription factors (TFs) enriched in 10 modules are presented in (Figures 4(c) and 4(d)). According to the enrichment scores, the corresponding relevant noncoding RNAs (ncRNAs) were closely associated with ubiquitin-mediated proteolysis module m1, spliceosome module m2, cell cycle modules m14 and m18, and endocytosis module m12, which included long noncoding RNAs (lncRNAs) MALAT1 (, ), FENDRR (, ), RNU1-1 (, ), and TUG1 (, ). The transcription factors (TFs) identified based on the enrichment scores were closely associated with GPR signaling pathway module m3, ECM-receptor interaction module m4, glutathione metabolism module m5, neuroactive ligand-receptor interaction module m9, endocytosis module m12, cell adhesion module m13, nucleotide excision repair module m17, homologous recombination module m16, and cell cycle modules m14 and m18, which included E2F1 (, ), TP53 (, ), YBX1 (, ), E2F4 (, ), SP1 (, ), BRCA1 (, ), CREB1 (, ), and CIITA (, ). As shown in Table 5, these key genes, such as MALAT1, RNU1-1, FENDRR, TUG1, E2F1, TP53, SP1, YBX1, BRCA1, E2F4, CREB1, and CIITA, play significant functional roles in their associated modules, suggesting that these genes may play roles in cell cycle regulation, methylation, acetyltransferase activity, and the splicing cycle. According to the integrated analysis results, these key genes of lncRNAs and TFs might play pathogenic roles in the occurrence and progression of IPF.

3.5. Validation of the lncRNAs and TFs in IPF with qRT-PCR

Demographic and clinical features of the IPF patients and the healthy control group are listed in Table 6. Patients with IPF smoked fewer cigarettes than in the control group. Moreover, the clinical features of the IPF group were decreased in pulmonary function. To validate the results obtained through integrated analysis of the three datasets related to IPF, the relative expression of the key genes was analyzed by RT-PCR (Fig. S4). We found that 3 of the 4 candidate genes have statistically significant differences between the IPF and normal groups (MALAT1, E2F1, and YBX1 with , FENDRR with ).

4. Discussion

Chronic and progressive airway remodeling is a major characteristic of IPF with unknown etiology. Although accumulating evidence reveals that activated fibroblasts have important effects on the pathogenesis and progression of IPF, the underlying molecular mechanisms involved in the regulation of IPF remain unclear. Previous findings of gene regulation on IPF have mainly focused on protein-coding genes which can delay but do not inhabit the development of fibrosis. Recently, with the development of high-throughput sequencing technology, epigenetic researches provide new insights into the underlying molecular and etiological mechanisms of IPF. Epigenetics, such as functional ncRNAs, refers to heritable changes in DNA and chromatin that influence gene expression other than changes in DNA sequence and has gradually become the research hotspot. Multiple studies have indicated that lncRNAs can influence the pathological process involving the structural remodeling of pulmonary architecture and eventually lead to respiratory failure. As multifunctional adaptor molecules, lncRNAs play multifunctional roles in the regulation of gene expression by regulating mRNA decay, splicing, and gene looping by binding to DNA, proteins, and certain other RNAs [17, 18]. In this study, we integrated three publicly available microarray datasets (GSE2052, GSE44723, and GSE24206) and found that differential expression of 8483 genes comprised 988 upregulated and 7495 downregulated genes. Consistent with the results of previous studies on the molecular mechanism of IPF, we found that DEGs were mainly concentrated on the extracellular matrix and these biological functions were mainly related to cell adhesion, proliferation, cytoskeleton development, and cytokine interaction. After a series of bioinformatics analysis, the regulatory network consisting of key lncRNAs and transcription factors (TFs), which may contribute to the pathogenesis of IPF, was ultimately obtained. We found that the biological functions of these key genes, which were related to epithelial-mesenchymal transition (EMT), mainly focused on mitochondrial translation, RNA processing, and ubiquitin-mediated proteolysis. We performed a comprehensive literature search and judged by integrating degrees, closeness, and betweenness centrality of the regulatory network ultimately identifying 2 lncRNAs and 2 TFs (MALAT1, FENDRR and E2F1, YBX1, respectively). Subsequently, we further validate the expression levels of these key genes related to the regulation of pulmonary fibrosis in blood samples between the IPF and control groups using real-time polymerase chain reaction (RT-PCR). As a result, differential expression of these genes including downregulated YBX1 and upregulated MALAT1 and E2F1 reached statistical significance except for FENDRR between two groups.

The research on epithelial-to-mesenchymal transition (EMT) related to a fibrotic process has received an increased attention in recent years. Considering the possible false-negative or false-positive results between the individual studies and the research sample size, we integrated and analyzed the potential lncRNAs and TFs related to the pathogenesis and progression of IPF from the public microarray data. Metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) located on chromosome 11q13.1, also known as nuclear-enriched abundant transcript 2 (NEAT2), is involved in various biological functions including molecular scaffolds for ribonucleoprotein complexes, transcriptional regulator for genes, and regulation of cell cycle [19]. Substantial research has confirmed that MALAT1 has many important physiological and pathological function in a wide range of diseases such as various solid cancers, septic lung injury, myocardial or renal ischemia-reperfusion injury, cardiac fibrosis, liver fibrosis, and silica-induced pulmonary fibrosis [2023]. Furthermore, MALAT1 also play an important role in EMT related to the pulmonary fibrosis [24]. Although MALAT1 was reported to be mainly localized in the nucleus, it could transfer from the nucleus to the cytoplasm during the G2/M-phase cell cycle [25]. E2F1 located on chromosome 20q11.2, which was screened out and validated in this study, belongs to TFs of the nuclear factor of the E2F family and participates in the cell cycle G1/S phase regulation mediating both cell proliferation and apoptosis [26]. Many studies have confirmed that E2F1 could activate the expression of stromal markers related to EMT such as vimentin and fibronectin and facilitate the pathogenesis processes such as fibrosis and tumor progression [27]. YBX1 located on chromosome 1p34.2, which is another screened and validated transcription factor in this study, belongs to a member of cold-shock protein family and acts as an important regulator related to cell proliferation and cell cycle [28]. Some studies have reported a correlation between abnormal expression of YBX1 and EMT markers such as vimentin and N-cadherin [29]. The above results suggest that genes screened and validated in this study might act as key regulators of the pathogenesis and progression of IPF.

However, the limitations of this study are as follows. First, although key genes related to the pathogenesis of IPF have been screened and validated through integrating three datasets and performing a series of bioinformatics approaches, expression levels of these key genes need further validated experiments such as western blotting (WB) and immunohistochemistry analysis (IHC). Second, differential gene analysis is one of the crucial data analysis strategies for expression profiling of IPF in GEO datasets. However, the three datasets combined and analyzed in this study from the GEO database microarray and platform were not unified. Meanwhile, the sample sizes of these three datasets were relatively small and imbalanced. The potential selection bias and information bias were inevitable. Therefore, the accuracy and reliability of candidate genes could be improved greatly by integrating more various types of datasets. Third, verification of the expression levels of candidate genes in clinical samples is far from enough. Further functional verification of these candidate genes was necessary to perform by loss-of-function and gain-of-function experiments in vivo and in vitro. Lastly, the verification and discussion of the underlying molecular mechanisms of these candidate genes involved in the pathogenesis and progression of IPF will be necessary to confirm through a chromatin immunoprecipitation assay (CHIP) or dual-luciferase reporter gene assay and so on. Regardless of the limitations mentioned above, this study provided preliminary evidence for the candidate genes related to the pathogenesis of IPF. As a time-saving and cost-saving method for analysis of biomedical data, we took an extensive bioinformatics data-mining approach from different microarray platforms to obtain candidate lncRNAs and TFs in IPF. This study provided a framework and broad application prospects for exploring pathological molecular networks related to IPF.

5. Conclusion

This study provides reliable and comprehensive perspectives on the pathogenesis and progression of IPF; potential lncRNAs and TFs related to the pathogenesis of IPF were obtained through bioinformatics analysis. Ultimately, the 3 key genes that were found to show abnormal expression in IPF compared to normal lung tissues may be considered as biomarkers for the diagnosis and treatment of IPF, which should be verified in subsequent studies.

Data Availability

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Ethical Approval

The collection and usage of the blood samples were approved by the Medical Research Ethics Committee of Traditional Chinese Medicine Hospital Affiliated to Xinjiang Medical University (the Scientific Research Project 2018XE0109-1).

Conflicts of Interest

The authors issue the statement with no conflicts of interest in this work.

Authors’ Contributions

Conceptualization was handled by F.W. and F.L. Investigation, formal analysis, resources, and writing (original draft preparation) were taken care by F.W. and P.L. Writing (review and editing) was done by all authors. Fan Wang and Pei Li contributed equally to this work.

Acknowledgments

This research received a specific grant from the Natural Science Foundation of Xinjiang Uygur Autonomous Region (No. 2019D01A06) and Xinjiang Uygur Autonomous Region Graduate Research and Innovation project (No. XJ2019G175).

Supplementary Materials

Fig. S1: hierarchical clustering using differentially expressed genes across all samples from GSE2052. Fig. S2: hierarchical clustering using differentially expressed genes across all samples from GSE44723. Fig. S3: hierarchical clustering using differentially expressed genes across all samples from GSE24206. Fig. S4: the relative expression of the key genes was analyzed by RT-PCR. Table S1: sequences of primers for candidate genes. Table S2: the RT-PCR results of four candidate genes. (Supplementary Materials)