Peripheral blood gene expression profiles in COPD subjects

To identify non-invasive gene expression markers for chronic obstructive pulmonary disease (COPD), we performed genome-wide expression profiling of peripheral blood samples from 12 subjects with significant airflow obstruction and an equal number of non-obstructed controls. RNA was isolated from Peripheral Blood Mononuclear Cells (PBMCs) and gene expression was assessed using Affymetrix U133 Plus 2.0 arrays. Tests for gene expression changes that discriminate between COPD cases (FEV1< 70% predicted, FEV1/FVC < 0.7) and controls (FEV1> 80% predicted, FEV1/FVC > 0.7) were performed using Significance Analysis of Microarrays (SAM) and Bayesian Analysis of Differential Gene Expression (BADGE). Using either test at high stringency (SAM median FDR = 0 or BADGE p < 0.01) we identified differential expression for 45 known genes. Correlation of gene expression with lung function measurements (FEV1 & FEV1/FVC), using both Pearson and Spearman correlation coefficients (p < 0.05), identified a set of 86 genes. A total of 16 markers showed evidence of significant correlation (p < 0.05) with quantitative traits and differential expression between cases and controls. We further compared our peripheral gene expression markers with those we previously identified from lung tissue of the same cohort. Two genes, RP9and NAPE-PLD, were identified as decreased in COPD cases compared to controls in both lung tissue and blood. These results contribute to our understanding of gene expression changes in the peripheral blood of patients with COPD and may provide insight into potential mechanisms involved in the disease.


Introduction
Chronic obstructive pulmonary disease (COPD), an inflammatory disorder that is characterized by a slowly progressive development of irreversible airflow limitation, is currently the fourth leading cause of death in the United States. Sixteen million Americans live with the disease, and there are 800 million affected individuals worldwide. Strongly associated with cigarette smoking, COPD is expected to be the third most common cause of death and fifth most common cause of disability worldwide by 2020 [1]. COPD is typically diagnosed late in life, and late in the course of disease when the patient presents with significant physiological impairment [2,3]. The need for improved early diagnosis and the identification of novel therapeutic targets for this debilitating disease has recently gained heightened interest.
Chronic obstructive bronchitis/bronchiolitis with peribronchiolar fibrosis (small airways disease), and abnormal enlargement of airspace distal to the terminal bronchioles with destruction of lung parenchyma (emphysema) are the pathological hallmarks of disease. Small airways disease and emphysema can present alone or in combination, with varying degrees of severity [4,5]. COPD is now considered primarily an inflammatory disorder involving abnormalities in both innate and adaptive immune responses. Inflammatory abnormalities in COPD include a significant increase in macrophage numbers in the lung and alveolar space, at the sites of alveolar destruction. Increased macrophage numbers may be due to increased monocyte recruitment and may result in higher secretion of inflammatory proteins leading to pathophysiological features of COPD [6].
However, systemic impairments have also been observed in patients [7].
Environmental factors contribute to varying susceptibility to COPD in the general population with the greatest environmental exposure in developed countries being tobacco smoke [8,9]. Exposure to other airborne pollutants, such as ozone, also appears to increase risk. While an increasing rate of smoking contributes to the growing incidence of COPD in developing countries as well, indoor air pollution associated with heating and cooking fuel is the major environmental risk factor, contributing to almost 3% of the global burden of disease [10]. In addition to environmental risk factors, varying genetic susceptibility to COPD exists among individuals, particularly with respect to the response to cigarette smoke [11,12]. Given the complexity of disease pathogenesis, the presence of varying levels of susceptibility in the general population and the fact that patients rarely present early in disease pathogenesis (at a time when disease-modifying therapy may be more effective) the identification of biological markers of disease susceptibility and/or progression are needed.
Numerous previous studies have sought to identify disease biomarkers in various forms, such as genetic or expression variants. DNA microarrays have been proven to be a powerful tool capable of biomarker discovery for various disease states. Multiple groups have previously applied microarray analysis to identify gene expression changes associated with COPD [13][14][15][16]. All these studies have used lung tissues obtained through invasive surgical procedures. Application of discovery approaches to samples derived from minimally-invasive procedures may provide biomarkers for diagnosis and therapeutic management of COPD. One previous study used whole blood to search for novel protein markers of COPD [16]. Here, we present a novel gene expression microarray data set generated from PBMC isolated from 24 subjects with varying levels of airflow obstruction.

Sample Collection
This study was approved by the Partners Health Care Human Research Committee. Peripheral blood, along with lung tissue, was obtained from 24 patients admitted to Brigham and Women's Hospital for suspected stage 1 lung tumors. Informed consent was provided and subjects underwent lung function testing by spirometry and completed a lung health-related questionnaire prior to surgery. Age, height, weight, sex and surgical pathology were obtained from subjects' medical charts. Predicted lung function values (FEV 1 , FVC) were calculated in SAS using the Crapo equations for Caucasians and the Hankinson equations for African-Americans. Diagnosis was confirmed by surgical pathology. A paper describing identification of COPD biomarkers identified by expression profiling from the lung tissue samples has been published previously [17].

Isolation of PBMC RNA
PBMCs were isolated by using CPT tubes (Becton Dickinson, Franklin Lakes, NJ) according to manufacturer's instructions. Approximately 8 ml whole blood was collected from each subject. Following centrifugation, cells were lysed for RNA isolation. DNAse-free total RNA was purified using the RNeasy mini kit (Qiagen, Inc, Valencia, CA) according to manufacturer's recommendations. RNA concentrations were determined by Nanodrop ND-1000 (Nanodrop Technologies, Wilmington, DE). RNA quality was assessed on an Agilent 2100 Bioanalyzer; samples with a RNA Integrity Score of > 6.0 were used in this study.

Microarray Analysis
RNA samples were used for fluorescent target generation (via in-vitro transcription), hybridized, washed, and scanned on U133 plus 2.0 GeneChips (Affymetrix, Santa Clara, CA) according to the manufacturer's instructions. Two independent versions of expression intensities were extracted from raw data files using either RMA or MAS 5.0 algorithms implemented in BioConductor. Gene annotation information was retrieved from the Affymetrix analysis portal (NetAffx http://www.affymetrix.com). Unsupervised clustering with the nonparametric bootstrap [18] was applied to check for undesirable and unanticipated structure or associations among the samples. Reliability of signal intensity measurements was determined using the Detection Call in GCOS, and analysis was restricted to probe sets reliably detected in all cases and/or all controls.
For discrete analysis (cases vs. controls), we applied two independent tests for differential expression on each version of the data set; Bayesian Analysis of Differential Gene Expression (BADGE) [19] and Significance Analysis of Microarrays (SAM) [20]. In SAM False Discovery Rate (FDR) is calculated by computing the number of significant genes for a given threshold for each permutation. The median number of significant genes from these permutations is the median False Discovery Rate. Since genes identified as significant from the permuted data are in fact false positives and as such the median number over many randomizations is a good estimate of false discovery rate. For quantitative analysis, correlation coefficients of signal intensity and lung function (FEV 1 or FEV 1 /FVC) were calculated. For each probe set, we calculated both the Pearson linear and Spearman rank correlation coefficients for both RMA and MAS5derived expression intensities using SAS.

Functional Classification
Functional classification of gene sets was performed using EASE v2.0 [21]. Affymetrix probe set IDs for the selected genes were used as the input list while probe set IDs for all filtered probe set genes served as the background set. Gene ontology categories with an EASE score of less than 0.05 were defined as significantly over-represented. Pathway analysis was performed using Ingenuity Pathway Analysis (IPA) on the set of discrete and quantitative biomarkers to identify canonical pathways that are associated with the peripheral markers determined by expression analysis. Canonical pathways with Fisher Exact test p-values less than 0.05 were identified as significantly dysregulated.

Molecular Validation
We performed quantitative real-time polymerase chain reaction (qPCR) for the genes identified as discrete and quantitative disease markers using assays from Applied Biosystems (Foster City, CA). qPCR was performed on a Agilent MX3000P (La Jolla, CA) using TaqMan chemistry, as previously described [22]. Gene expression levels were calculated according to the relative expression analysis approach using multiple endogenous controls (PPIA, GAPDH, ACTB, and HPRT1). Statistical analysis was performed on individual sample dCt values for each gene using either the parametric Student's t-test or nonparametric Mann-Whitney U-test.

Subject Demographics
The studies involved 24 subjects including 12 COPD cases with significant airflow obstruction (defined as FEV 1 < 70% predicted and FEV 1 /FVC ratio < 0.7) and 12 control subjects with normal lung function (FEV 1 > 80% predicted and FEV 1 /FVC ratio > 0.7). Cases had an average age of 63 years and average FEV1 of 47% predicted. Controls had an average age of 64 years and an average FEV1 of 99% predicted (Table 1). This population represents a subset of the population for which we have previously reported lung tissue gene expression patterns [17].

Expression Biomarker Discovery Discrete Analysis
We first extracted signal intensity data using RMA and MAS5, removed data from all probe sets that were not reliably detectable in either all cases or all control samples. We used a stringent set of conditions to identify differential gene expression in this data set, applying multiple significance testing methods (SAM & BADGE). A total of 691 probe sets were significantly different in BADGE analysis at a p-value of 0.01 or less for either RMA or MAS5 versions of data. As our data analysis approach included a combination of multiple tests and normalization approaches, we did not implement any correction on BADGE p-value. SAM analysis identified a total of 93 probe sets that were significantly different at median false discovery rate of 0 (median FDR = 0) in either RMA or MAS5 versions of the data. Ninety (97%) of the probe sets identified in SAM analysis were also identified using BADGE, and represented 47 known genes that we defined as differentially expressed in PBMC from subjects with COPD cases versus controls. Interestingly, all genes identified using these highly stringent criteria were expressed at lower levels in cases as compared to controls (Figure 1). A list of 90 probe sets identified by SAM and BADGE has been provided in Additional File 1 Table S1.

Quantitative Analysis
We also calculated Pearson's & Spearman rank correlations between expression and FEV 1 or FEV 1 /FVC for all probe sets. A total of 146 probe sets were significantly correlated with FEV 1 at p < 0.05 and 9 probe sets were   Figure  2), while at p < 0.01, the overlap was 6 probe sets, representing 2 known genes (SOX6, LMLN; both positively correlated with pulmonary function). A list of 104 probe sets significantly correlated with both FEV 1 and FEV 1 /FVC at p < 0.05 has been provided in Additional File 1 Table S2A and Table S2. There was no overlap among probe sets at p < 0.001. A total of 158 probe sets passed criteria as either discrete (90 probesets) or quantitative (104 probesets at p < 0.05) gene expression markers of COPD. Among these, 36 probe sets representing 16 known genes were significantly different in Case-Control analysis and significantly correlated with both FEV 1 and FEV 1 /FVC at p < 0.05 (Figure 3). We assessed whether differences in the distribution of tumor type between Cases and Controls contributed to the identification of these gene expression changes. The tumor types among the 24 subjects included 9 adenocarcinoma and 9 squamous cell carcinoma subjects. We applied differential expression analysis (as described for COPD cases and controls above) comparing all samples classified as adenocarcinoma versus samples classified as squamous cell carcinoma. No probe sets were identified as consistently differentially expressed between tumor types. Further, no probe sets identified as differentially expressed between tumor types in any single analysis were among the COPD biomarker gene set.

Functional Classification
In an order to identify biological systems or functions that are associated with discrete or quantitative COPD peripheral gene expression markers, we performed gene ontology assessment using EASE ( Figure 4A). We used a set of 158 probe sets that were either significantly different in cases and controls or significantly correlated with lung function and queried for over-represented ontologies using EASE. There was a consistent overrepresentation of functions relating to transcriptional activity and nucleic acid binding for all sets of COPD biomarkers. A total of 103 probe sets, or 65% of biomarkers tested for ontology (some of the probe sets lacked ontological annotation), were classified in one or more categories related to these functions. Among all 158 peripheral biomarker genes (both discrete and quantitative), 40 had an annotated molecular function, 37 had an annotated biological process and 30 had an annotated molecular function. Eighteen of 40 genes (45%; p < 0.05) were classified for the molecular function of Nucleic Acid Binding (GO: 0003676). Twelve of 37 genes (32%; p < 0.05) were classified for the biological process of DNA-dependent Transcription (GO: 0006351) and 15 of 37 genes (40%; p < 0.05) were classified for the biological process of nucleoside, nucleoside & nucleotide metabolism (GO::0008150). Among the discrete marker genes (case vs. control), 30 had an annotated molecular function and 27 had an annotated biological process. Thirteen of 30 genes (32%; p < 0.05) were classified for the molecular function of Nucleic Acid Binding (GO: 0003676). Nine of 27 genes (33%; p < 0.05) were classified for the biological process of DNA-dependent Transcription (GO: 0006351) and 11 of 27 genes (40%; p < 0.05) were classified for the biological process of nucleoside, nucleoside & nucleotide metabolism (GO::0008150). Among the quantitative marker genes (correlation), 18 genes had an annotated molecular function and 18 had an annotated biological process. Nine of 18 genes (50%; p < 0.05) were classified for the molecular function of Nucleic Acid Binding (GO: 0003676) and four of 18 genes (22%; p < 0.05) were classified for the molecular function of  Transcription factor activity (GO: 0003700). Six of 18 genes (33%; p < 0.05) were classified for the biological process of DNA-dependent Transcription (GO: 0006351) and six of 18 genes (33%; p < 0.05) were classified for the biological process of nucleoside, nucleoside & nucleotide metabolism (GO::0008150). A list of all significantly over-represented Gene Ontology (GO) classes is presented in Additional File 1 Table S3.
Pathway analysis also provided insights to canonical pathways associated with peripheral markers ( Figure  4B). Significantly over-represented pathways (p < 0.05) included those associated with cell signaling, inflammatory cell regulation, and cancer; EGF (3 biomarkers out of a total of 171 genes associated with category in IPA), Endothelin-1 (4 biomarkers out of a total of 101 genes) mTOR (4 biomarkers out of a total of 124 genes), CXCR4 (4 biomarkers out of a total of 192 genes), IL-2 (2 biomarkers out of a total of 61 genes), IL-3 (2 biomarkers out of a total of 76 genes), IL-17 (2 biomarkers out of a total of 77 genes), ILK (3 biomarkers out of a total of 191 genes), IL-8 (3 biomarkers out of a total of 193 genes), breast cancer (3 biomarkers out of a total of 106 genes) lung cancer (3 biomarkers out of a total of 101 genes), and glioblastoma (3 biomarkers out of a total of 166 genes). A list of all significant canonical pathways is presented in Additional File 1 Table S4.

Validation
We performed qPCR-based validation for a subset of genes identified as differentially expressed in COPD subjects using both discrete and quantitative analyses using all samples (n = 24). Validation analysis confirmed significant correlations between microarray-based and qPCR-based expression measures for GKAP1 (|r| = 0.25, p < 0.05) and STX17 (|r| = 0.36, p < 0.05). However, qPCR did not confirm significant differences in expression for either of these genes between cases and controls.

Discussion
Even with current advancements in medical technologies, appropriate diagnosis and management of COPD remains a major challenge. Spirometry as a measure of lung function remains the primary objective test for diagnosis of COPD, but spirometry cannot indicate whether airflow obstruction relates to emphysema, airway disease, or both processes. Additional non-or minimally-invasive approaches would be very useful for disease diagnosis and management.
In recent years, studies have attempted to identify gene expression biomarkers for COPD [13][14][15]. In those studies, genome-wide expression studies have been based on RNA derived from surgically-derived tissue samples. Although gene expression studies of lung tissues may provide useful insights into disease pathogenesis, it is not practical to consider routine COPD diagnosis from a sample that must be obtained through an invasive surgical procedure. Blood samples are less invasive, potentially provide for a larger sample size, and allow repeated sampling to monitor disease progression over time and to study therapeutic response.
Past genome-wide studies on different organ systems have shown that total RNA derived from circulating blood can distinguish between control subjects and patients with various diseases [23][24][25][26][27][28][29][30][31] including inflammatory (e.g. preeclampsia, rheumatoid arthritis, and chronic pancreatitis) and malignant (chronic lymphocytic leukemia and renal cell carcinoma) diseases [32,33]. One of the earliest demonstrations that gene expression changes in peripheral blood mononucleocytes (PBMCs) were associated with disease was demonstrated on a rat brain model, where acute neural assaults resulted in gene expression changes in PBMCs within 24 hours [34]. In the pulmonary system, Showe et al have used peripheral blood gene expression signatures to identify early-stage lung cancer in at-risk populations [35]. Karimi et al. (2006) showed that in vitro exposure of PBMC to cigarette smoke induces production of cytokines in a TLR4-dependent manner [36].
We hypothesized that peripheral blood gene expression patterns could help to improve COPD detection, diagnosis or progression. We assessed genome-wide expression patterns in RNA obtained from PBMCs isolated from a subset of 24 of the study subjects using the Affymetrix U133 Plus 2.0 microarray. Data analysis revealed novel genes that were differentially expressed in PBMCs from COPD patients. The genes we identified have not been previously implicated in COPD disease pathogenesis, and as such are likely to be true markers rather than etiological. We observed two genes, RP9 and NAPE-PLD, showing decreased expression in both lung tissue and blood of COPD subjects when compared to controls. This suggests that PBMC-derived markers may reflect processes ongoing in diseased tissues. Further, our data serves as a proof-of-principal that peripheral gene expression patterns, defined using minimally invasive samples, can be used to describe COPD.
Genome-wide linkage screens aimed to identify disease-susceptibility genes previously identified three linkage regions (chromosomes 2q33-36, 8pter-22, and 12p13-12) in the Boston Early-Onset COPD cohort [37] which includes the locus for one of the novel genes identified in our study, AT-rich domain 2 (ARID2). ARID2 is a transcriptional co-activator involved in the regulation of cardiac gene expression [38]. Among other genes displaying changes in expression between cases and controls, some have notable functions. Syntaxin 17 (STX17) expression in macrophages is regulated by Colony-stimulating factor 1 (CSF-1), a growth factor controlling the development of macrophages from myeloid progenitor cells [39]. FOXP1 is a member of winged-helix/forkhead transcription factors and is important in monocyte differentiation and macrophage function [40]. SESN1, a stress inducible sestrin regulated by p53, has been reported to be potent inhibitor of mTOR signaling and regulator of cell growth and proliferation [41].
To our knowledge only two studies have previously explored the value of genome-wide peripheral blood expression assessments in patients with COPD [16,42]; both defining serum protein levels. Hurst et al assessed paired baseline and exacerbation plasma samples from patients with COPD and identified 36 biomarkers using protein arrays [42]. They observed that although systemic biomarkers were not helpful in predicting exacerbation severity, acute-phase response at exacerbation was strongly related to monocyte activity. Pinto-Plata et al used protein array on peripheral blood from COPD patients and identified 30 biomarker clusters [16]. They identified a set of biomarkers correlated with lung function.
One major limitation of the current study is that quantitative real time-PCR (qPCR) validation indicated a potential high false discovery rate. Possible reasons for lack of validation for individual genes include expression levels below sensitivity for the assays used, poor assay specificity, alternative splice forms and inaccuracy of array data. The phenotypic heterogeneity of COPD may also be a cause of limited validation results in the current study. Regardless of the root cause of poor validation, the small size of the current study is a major limitation in the generalization of the results presented.
Another limitation of the current study is the diagnosis of lung cancer in most subjects. Recent studies have reported that genetic expression in PBMCs is altered in the context of malignancy [32,43]. Lung cancer and COPD are both typically found in smokers and the diagnosis of lung cancer can serve as an independent predictor for COPD, independent of smoking history. Even though we have previously shown any effects of the tumor on gene expression are not significant in distant, histologically normal lung tissue [17], in the case of PBMCs the presence of tumors may contribute to changes in gene expression. Even though four (PIK3C2A, JUN, FNBP1, ITPR1) of our peripheral biomarkers have been implicated in cancer pathophysiology, none of the PBMC biomarkers were differentially expressed between tumor types (among all subjects, or within cases or controls alone).
In conclusion, we used microarray technology to identify gene expression differences in PBMC obtained from COPD patients and controls. Our data contribute to the understanding of gene expression changes occurring in the blood of patients with obstructive lung disease and provide additional insight into potential mechanisms involved in the disease process. Our data suggest that PBMC may be a source of diagnostic markers. The identification and validation of markers may help to facilitate the development of non-invasive methods for diagnosis, classification of disease subtypes and/or provide a means to define response to therapeutic intervention.