Comparative Expression Profiling in Pulmonary Fibrosis Suggests a Role of Hypoxia-inducible Factor-1a in Disease Pathogenesis Argyris Tzouvelekis1*, Vaggelis Harokopos2*, Triantafillos Paparountas2*, Nikos Oikonomou2, Aristotelis Chatziioannou3, George Vilaras4, Evangelos Tsiambas4, Andreas Karameris4, Demosthenes Bouros1, and Vassilis Aidinis2 1
Department of Pneumonology, Medical School, Democritus University of Thrace, and University Hospital of Alexandroupolis, Alexandroupolis, Greece; 2Institute of Immunology, Biomedical Sciences Research Center ‘‘Alexander Fleming’’, Athens, Greece; 3Institute of Biological Research and Biotechnology, National Hellenic Research Foundation, Athens, Greece; and 4Department of Pathology, Veterans Administration Hospital (N.I.M.T.S), Athens, Greece
Rationale: Despite intense research efforts, the etiology and pathogenesis of idiopathic pulmonary fibrosis remain poorly understood. Objectives: To discover novel genes and/or cellular pathways involved in the pathogenesis of the disease. Methods: We performed expression profiling of disease progression in a well-characterized animal model of the disease. Differentially expressed genes that were identified were compared with all publicly available expression profiles both from human patients and animal models. The role of hypoxia-inducible factor (HIF)-1a in disease pathogenesis was examined with a series of immunostainings, both in the animal model as well as in tissue microarrays containing tissue samples of human patients, followed by computerized image analysis. Measurements and Main Results: Comparative expression profiling produced a prioritized gene list of high statistical significance, which consisted of the most likely disease modifiers identified so far in pulmonary fibrosis. Extending beyond target identification, a series of meta-analyses produced a number of biological hypotheses on disease pathogenesis. Among them, the role of HIF-1 signaling was further explored to reveal HIF-1a overexpression in the hyperplastic epithelium of fibrotic lungs, colocalized with its target genes p53 and Vegf. Conclusions: Comparative expression profiling was shown to be a highly efficient method in identifying deregulated genes and pathways. Moreover, tissue microarrays and computerized image analysis allowed for the high-throughput and unbiased assessment of histopathologic sections, adding substantial confidence in pathologic evaluations. More importantly, our results suggest an early primary role of HIF-1 in alveolar epithelial cell homeostasis and disease pathogenesis, provide insights on the pathophysiologic differences of different interstitial pneumonias, and indicate the importance of assessing the efficacy of pharmacologic inhibitors of HIF-1 activity in the treatment of pulmonary fibrosis. Keywords: idiopathic pulmonary fibrosis (IPF); expression profiling; tissue microarrays; hypoxia-inducible factor-1a (HIF-1a)
(Received in original form May 8, 2007; accepted in final form August 29, 2007) Supported by the Society for Respiratory Research and Treatment of Eastern Macedonia and Thrace (D.B.), European Commission Network of Excellence grant QLRT-CT-2001-01407 (V.A.), and Hellenic Ministry for Development grant GSRT-PENED-136 (V.A.). A.T. is a recipient of an annual research grant in respiratory medicine provided by GlaxoSmithKline. *These authors contributed equally to this article. Correspondence and requests for reprints should be addressed to Vassilis Aidinis, Ph.D., Institute of Immunology, B.S.R.C. Alexander Fleming, 34 Fleming Street, 16672, Athens, Greece. E-mail:
[email protected] This article has an online supplement, which is accessible from this issue’s table of contents at www.atsjournals.org Am J Respir Crit Care Med Vol 176. pp 1108–1119, 2007 Originally Published in Press as DOI: 10.1164/rccm.200705-683OC on August 30, 2007 Internet address: www.atsjournals.org
AT A GLANCE COMMENTARY Scientific Knowledge on the Subject
Despite intense research efforts, the etiology and pathogenesis of idiopathic pulmonary fibrosis remain poorly understood, which is reflected in the lack of effective treatment. What This Study Adds to the Field
Our results suggest an early primary role of HIF-1 in alveolar epithelial cell homeostasis and disease pathogenesis, provide insights on the pathophysiologic differences of different interstitial pneumonias, and indicate the importance of assessing the efficacy of pharmacologic inhibitors of HIF-1 activity in the treatment of pulmonary fibrosis.
Idiopathic interstitial pneumonias (IIPs) are a heterogeneous group of diseases comprising seven distinct clinical and pathologic entities. Idiopathic pulmonary fibrosis (IPF) and cryptogenic organizing pneumonia (COP) represent two of the most prevalent members of the disease group, with major differences in pathogenesis, clinical course, and prognosis (1). IPF is a refractory and lethal IIP characterized by fibroblast proliferation, extracellular matrix deposition, and progressive lung scarring, and comprises the histopathologic pattern of usual interstitial pneumonia (UIP). The incidence of IPF is estimated at 6.8 to 16.3 cases per 100,000 per year in the United States, and the mean survival from the time of diagnosis is 3 to 5 years regardless of treatment (2–4). Although the etiology and pathogenesis of IPF remain poorly understood, current research suggests that the mechanisms driving IPF reflect abnormal, deregulated wound healing in response to multiple sites of ongoing alveolar epithelial injury, involving increased activity and possibly exaggerated responses by a spectrum of proinflammatory and profibrogenic factors (3, 5). Expression profiling, the estimation of the expression level of thousands of genes by DNA microarrays, is a powerful tool for biologists, bioinformaticians, and statisticians in their attempt to decipher the complex organization of biological phenomena. In this context, and to identify genes and/or cellular pathways involved in the initiation and progression of IPF, we used the bleomycin (BLM)-induced animal model, the closest equivalent of the human disease. RNA lung samples were isolated at different endpoints in the development of the disease and hybridized to cDNA microarrays. After robust statistical selection of differentially expressed genes (DEGs), results were compared
Tzouvelekis, Harokopos, Paparountas, et al.: HIF-1a in IPF Pathogenesis
with all publicly available microarray datasets in IPF (6–15), both from mice and humans, thus creating a unique list of likely disease modifiers. Furthermore, gene ontology and pathway analysis revealed hypoxia signaling among the most statistically important deregulated pathways. Prompted by the meta-analysis results, we investigated the role of hypoxia-inducible factor (HIF)-1a in disease pathogenesis, in the animal model as well as in human patients, to reveal an early primary role of HIF-1a in IPF development. Some of the results of these studies have been previously reported in the form of abstracts (16, 17).
METHODS Animals All mouse strains were bred and maintained in the C57/Bl6 background for over 20 generations in the animal facilities of the Biomedical Sciences Research Center ‘‘Alexander Fleming’’ (Athens, Greece) under specific pathogen–free conditions, in compliance with the Declaration of Helsinki principles. Mice were housed at 20–228C, with 55 6 5% humidity, and a 12-hour light:dark cycle; food and water was given ad libitum. All experimentation was approved by an internal institutional review board, as well as by the Veterinary Service and Fishery Department of the local governmental prefecture. Pulmonary fibrosis was induced by a single tail vein injection of BLM hydrogen chloride (100 mg/kg body weight; 1/3 of lethal dose, 50% [1/3LD50]; Nippon Kayaku Co. Ltd., Tokyo, Japan) to 6- to 8-week-old mice as previously reported in detail (18).
Expression Profiling Total RNA from the right lobe of lung specimens was isolated by homogenization in ice-cold TRIzol reagent (Invitrogen Life Sciences, Carlsbad, CA) followed by a single passage through an RNeasy column (QIAGEN GmbH, Hilden, Germany). Isolated total RNA was reverse transcribed with Superscript Reverse Transcriptase II (Invitrogen), and the cDNA was indirectly labeled using the amino-allyl cDNA labeling method. Experimental samples were mixed with equimolar amounts of the baseline sample (which was used as a common reference sample throughout) and hybridized in quadruplicates to cDNA glass microarray slides (Riken, Yokohama, Japan) interrogating 18,816 genes. After image analysis, all microarray data were subjected to preprocessing, lowess normalization, centering, and/or averaging. To select statistically significant DEGs, and because there is no international consensus on the most appropriate method for statistical selection, we used simultaneously the two most widely used methods: a parametric and a nonparametric analysis of variance (Kruskal-Wallis), using proprietary algorithms implemented in MATLAB (version 7.1, release 14; The MathWorks, Inc., Natick, MA). Reverse transcriptase–polymerase chain reaction (RT-PCR) gene validation was performed using MMLV reverse transcriptase and an oligo-dT(15) primer (Promega, Madison, WI). Detailed information on expression profiling, including gene ontology and pathway analysis, are provided in the online supplement.
Human Subjects In total, 45 newly diagnosed patients with IIPs of two different histopathologic patterns (IPF/UIP, and COP/organizing pneumonia [COP/OP]) were recruited in our study. The diagnosis of IIPs was based on the consensus statement of the American Thoracic Society/European Respiratory Society in 2002 (1, 19). Subjects were separated according to the histopathologic pattern of the IIPs as shown in Table 1. All patients were treatment naive at study inclusion. Paraffinembedded surgical lung specimens (open lung biopsy or by videoassisted thoracoscopic surgery) from two different fibrotic regions of each individual were sampled. All patients were fully informed and signed an informed consent form in which they agreed to the anonymous usage of their lung samples for research purposes.
Tissue Microarrays, Immunohistochemistry, and Computerized Image Analysis Tissue microarrays (TMAs) were constructed from 85 tissue samples consisting of 45 lung specimens from two different histopathologic
1109 TABLE 1. DEMOGRAPHIC AND SPIROMETRIC CHARACTERISTICS OF PATIENTS WITH IPF/UIP, PATIENTS WITH COP/OP, AND CONTROL SUBJECTS Characteristics Number Sex, male/female Age, median yr Smokers/nonsmokers FVC, % pred FEV1, % pred KCO, % pred
IPF/UIP
COP/OP
Normal
25 19/6 68 (40–75) 20/5 75 6 3 84 6 4 63 6 5
20 14/6 50 (35–65) 7/13 81 6 4 86 6 3 70 6 6
40 30/10 39 (26–60) 25/15 106 6 14 104 6 17 92 6 8
Definition of abbreviations: COP/OP 5 cryptogenic organizing pneumonia/ organizing pneumonia; IPF/UIP 5 idiopathic pulmonary fibrosis/usual interstitial pneumonia. Values are expressed as mean 6 SD or median (range).
patterns of IIPs and 40 control tissues derived from the normal part of lungs removed for benign lesions. After epitope demasking, TMAs were immunostained with a number of antibodies against HIF-1a, surfactant protein A (SP-A), vascular endothelial growth factor (VEGF), p53, and DNA fragmentation factor (DFF). Signal intensities were quantified with computerized image analysis using a semiautomated system. Statistical analysis was performed using SPSS 13.0 software (SPSS, Inc., Chicago, IL). Details on these methodologies can be found in the online supplement.
RESULTS Expression Profiling
To identify genes and/or cellular pathways involved in the initiation and progression of IPF, we performed expression profiling of disease progression in the animal model of BLMinduced pulmonary inflammation and fibrosis (18, 20). In this model, and as reported previously (18), BLM administration results in progressive subpleural/peribronchial pulmonary inflammation, which subsequently diffuses into the parenchyma. Inflammation is followed by the development of mainly subpleural and peribronchial fibrotic patches, characterized by alveolar septa thickening and focal dilation of respiratory bronchioles and alveolar ducts. Concomitantly, collagen accumulation peaks 23 days post–BLM injection (Figure E4 of the online supplement). The model is very reproducible, using standardized procedures and dedicated functional readouts exhibiting minimal variation (18). RNA lung samples were isolated at 7, 15, and 23 days post–BLM administration, corresponding to the inflammatory, intermediate, and fibrotic phases of the disease (Figure E4). Similarly, as a baseline control, RNA lung samples were isolated from littermate mice 23 days after administration of saline alone. Equimolar amounts of purified RNA from five mice per endpoint were pooled, to minimize biological diversity, and fluorescently labeled using the amino-allyl indirect labeling method as described in METHODS. Identical labeled samples from the same pool were mixed with the labeled common reference sample (wt/saline) and hybridized in (technical) quadruplicates to cDNA glass microarray slides, interrogating 18,816 genes. After image acquisition and analysis, microarray data were analyzed as outlined in Figure E1, using proprietary algorithms implemented in MATLAB. Briefly, and as described in detail in METHODS, after preprocessing, lowess normalization and quality control (Figure E2), centering was applied either before or after averaging, thus producing two gene matrices. These two matrices were further analyzed with two different statistical selection methods, one parametric and one nonparametric, thus ending up with four different lists of likely DEGs. The 1,172 genes identified as differentially expressed from all methods (having therefore a very high statistical significance
1110
AMERICAN JOURNAL OF RESPIRATORY AND CRITICAL CARE MEDICINE VOL 176
and a very low false discovery rate) are shown in Table E1. The differential expression of a small number of genes (clu, Hba-a1, spp1, slc6a6, nish, mt1) was further confirmed with semiquantitative RT-PCR (at three different RNA concentrations in the linear range of the reaction) in separate pools of five experimental animals and their controls (Figure E3). Comparative Expression Profiling and Meta-analysis
To validate our list of DEGs (Table E1) in a high-throughput mode, and to compare results from different animal models as well as from human patients, we collected (through database searching and personal communications) all publicly available information from published expression profiling datasets on IPF (6–15), each one with different levels of data quality, annotation, and availability. Mouse and human Entrez-Gene IDs for all reported DEGs from the different datasets/studies were retrieved using the Ingenuity Pathways Analysis (IPA version 5.0; Ingenuity Systems, Redwood City, CA) software. Comparisons were performed separately for both human and mouse EntrezGene IDs and results were fused together to avoid exclusions due to species nonconcordance. Strikingly, and although the compared data were obtained from various models and organisms (which conceptually are governed by different pathogenetic mechanisms), using different microarray platforms (containing different genes) and statistical methods, we identified a large number of genes in common between our dataset and the published genes in pairwise comparisons (Table 2 and Table E2). Therefore, the combined gene list (Table E3) containing 296 (nonredundant) genes (common DEGs [cDEGs]) identified as differentially expressed from at least two independent studies (our own and a published one) is self-validated, has a high statistical significance, and therefore is a valuable resource of likely disease-modifying genes. Among them, 35 genes were identified from three different datasets and 6 genes from four datasets (as highlighted in Table E3), prioritizing these genes even further. To prioritize the cDEGs systematically, an in-depth metaanalysis was conducted. Initially, a very extensive literature search with automated text mining using the Biolab Experiment Assistant software (Version 2.7; Biovista, Athens, Greece) and manual (PubMed) search revealed a total of 81 genes that have been found to play a direct role in the development of the disease (Table E4). This set of genes was used as a training set for the software application Endeavour (SymBioSys, Center for Computational Systems Biology, Katholieke Universiteit Leuven, Leuven, The Netherlands), which performs computational prioritization of ‘‘test genes,’’ based on a set of ‘‘training genes’’ (21). Endeavour uses nine different data sources including both vocabulary-based (e.g., Gene Ontology [GO]) as well as other data sources (e.g., BLAST and microarray databases). The ranking of a test gene for a given data source is calculated based on its similarity with the training genes, whereas the final prioritization is calculated based on order statistics of the individual rankings (21). The statistically more significant (according to Endeavour, P , 0.05) cDEGs are shown in Figure 1. To get additional functional insights on a potential role of DEGs in the development of the disease, we then examined whether any of these genes are tumor necrosis factor (TNF) or transforming growth factor (TGF) targets, the major proinflammatory or profibrogenic factor, respectively, with a definite role in disease induction (18, 22, 23). TNF or TGF DEGs were identified from published microarray datasets (24–28) and compared with IPF DEGs (Table E1; 1,172 genes). As shown in Table 2 (and Table E2 in detail), 19 and 37 DEGs, respectively, were found to be TNF or TGF targets. Remark-
2007
ably, 14 (of 19) and 21 (of 37) TNF/TGF targets that were found as IPF DEGs were also included in the cDEG list (Table E3), highlighting the role of TNF and TGF in disease development. Five and seven of them, respectively, are also included in the statistically significant Endeavour-prioritized cDEG list (highlighted in Figure 1). Finally, in an attempt to combine expression profiling with genetic linkage studies, IPF DEGs (Table E1) were compared with possible susceptibility genes from identified quantitative trait loci for BLM-induced pulmonary fibrosis (Blmpf1 and 2; References 29–31). Eleven of 22 genes from the blmpf1/2 loci, respectively, have been identified as DEGs (highlighted in Table E1) and three of five of these were also included in the cDEG list (highlighted in Table E3). GO and Pathway Analysis
In parallel with the statistical identification of DEGs and their prioritization, and to (1) prove the validity and extend the utility of the expression data analysis even further, (2) infer deregulated biological functions from the gene expression data, and (3) define functional criteria for further gene selection, the selected genes (Table E1) were annotated in the form of the GO terms, in the categories Molecular Function and Biological Process. GO term frequencies in the selected gene list were then calculated and their statistical significance (expressed as a P value) were estimated (through their hypergeometric distribution) as reported previously, and as described in detail in METHODS. As shown in Table 3, a number of well-expected functions and processes were found to be deregulated during the pathogenesis of BLM-induced pulmonary inflammation and fibrosis, such as inflammatory response, and chemokine, cytokine, and growth factor activity. As anticipated, GO analysis indicated multiple levels of gene expression regulation during pathogenesis (RNA helicase activity, transcription corepressor activity, transcription factor binding, magnesium ion binding; RNA processing, nuclear mRNA splicing, mRNA processing). The adhesion–cytoskeleton axis was also highlighted from the analysis, as indicated (directly or indirectly) from a number of deregulated functions and processes (respectively: GTPase activity, actin binding; actin filament severing, cell matrix adhesion). Notably, oxygen transport was indicated as the most significant deregulated GO function as well as GO process, indicating hypoxia as a pathogenic insult that could lead to (or exacerbate) pulmonary fibrosis. In a similar, complementary effort, the software program IPA (Ingenuity Systems) was used for automated gene expression data integration in cellular canonical pathways, as these are (pre)defined and curated by IPA. DEGs (Table E1) were examined for their participation in IPA canonical pathways, followed by a statistical test to examine if the pathway association could be observed by chance alone. The statistically significant (P , 0.05) deregulated canonical pathways are shown in Table 3. Remarkably, integrin and hypoxia signaling were ranked first in the list of statistically significant deregulated pathways, further supporting the GO analysis results. Early HIF-1a Overexpression in BLM-induced Pulmonary Inflammation and Fibrosis
To examine the role of hypoxia in the pathogenesis of pulmonary fibrosis, as indicated by the GO/IPA analysis, we then focused on the role of the HIF-1a, the major transcription factor that mediates cellular responses to hypoxia (32). Semiquantitative RT-PCR analysis indicated that the mRNA levels of Hif-1a are found to be up-regulated upon administration of BLM and the development of pulmonary inflammation and fibrosis (Figure
Tzouvelekis, Harokopos, Paparountas, et al.: HIF-1a in IPF Pathogenesis
1111
TABLE 2. LIST OF COMPARED DATASETS Microarray Type First Author/ Reference
Organism
Bleomycin Model
Type of Array
Chip Design
Statistical Selection Sample Size
Analysis/ Clustering
Normalization
P Value
Summary
DEGs Reported
Common DEGs P < 0.1
IPF Kaminski (8)
Mouse
Intratracheal instillation of bleomycin
Oligonucleotide /6,000 genes
Mu6500 30 Lung GeneChip tissue array specimens (Affymetrix) /4 replicates
Mean hybridization intensities of all probe sets on each array were scaled to an arbitrary, fixed level.
Gene Cluster and TreeView programs/2 clusters
NA
Global analysis of gene expression in PF reveals distinct programs regulating lung inflammation and fibrosis.
468
46
Lemay (11)
Mouse
Administered bleomycin (osmotic minipumps)
Oligonucleotide /22,690 probe sets, 12,422 genes
94 Lung MOE430A tissue GeneChip specimens arrays /11 arrays (Affymetrix)
Routines from Bioconductor within the R statistical language. Robust probe level model.
NA
