Increased IL-1 beta gene expression in peripheral blood leukocytes is associated with increased pain and predicts risk for progression of symptomatic knee osteoarthritis

Patients

Three independent cohorts of patients were examined: two from NYUHJD (designated “Learning Cohort” and “NYUHJD Validation Cohort”) and the “The Duke Validation Cohort.” Studies were approved by the Institutional Review Boards (IRB) of the NYU School of Medicine and Duke University, and all participants provided written informed consent prior to initiation of the studies.

NYUHJD Learning Cohort

Patients were diagnosed with knee OA by the referring physicians, according to the 1986 American College of Rheumatology (ACR) classification, and met clinical plus either radiographic or laboratory criteria for the diagnosis of idiopathic OA of the knee. Demographics, comorbidities, medications, number and location of joints affected (by both patient and physician report) were collected. We excluded patients with any other form of arthritis; body mass index (BMI) >33; any disorder requiring the use of systemic corticosteroids; history of bilateral knee replacements; major co-morbidities including diabetes mellitus, non-cutaneous cancer within 5 years of screening, chronic hepatic or renal disease, chronic infectious disease, congestive heart failure; and hyaluronan and/or corticosteroid injection to the affected knee within 3 months of screening. Non-OA controls had no clinical signs or symptoms of any arthritis. We assessed pain, functional status and quality of life with the Western Ontario and McMaster Osteoarthritis Index (WOMAC) (14), Short Form Health Survey (SF-36) (15), and Multidimensional Health Assessment Questionnaire (MDHAQ) (16); within the MDHAQ are modules for the patient global assessment of their disease and for visual analog scale (VAS) pain assessment.

The NYUHJD Validation Cohort

As part of an NIH-funded study, a separate cohort of 86 patients were followed for 24 months who met ACR clinical symptomatic criteria (17) and radiographic criteria [Kellgren-Lawrence (KL) grade ≥1] (18) for symptomatic knee OA. Exclusion criteria were as noted above. All patients underwent bilateral standardized weight bearing fixed-flexion posteroanterior (PA) knee radiographs using the SynaFlexer™ X-ray positioning frame (Synarc). We also screened 41 age-matched healthy controls, and among these recruited 20 subjects who had KL radiographic score <1 and no knee pain, of whom 12 were available for the analyses reported here. All patients were examined by an NYUHJD investigator every 6 months in this 24-month study. Radiographic assessments at baseline and 24 months include bilateral (signal and non-signal knee) KL determination, quantitative measurement of medial and lateral joint space width (JSW), medial and lateral JSN, osteophytes, medial tibial/lateral femoral subchondral sclerosis, and medial tibial attrition using the OARSI atlas, by two musculoskeletal radiologists blinded to patient information. Kappas for inter-rater agreement were 0.85 and 0.77 for KL scores of the right and left knees, respectively, and >0.85 for most other radiographic outcome measures. Concordance correlation coefficients were >0.90 for JSW measurements. JSW was measured at the narrowest portion of the joint space via electronic calipers and a medical monitor; the average value between the two readers was used in the statistical analysis. We collected genomic DNA and RNA from peripheral blood mononuclear cells (PBMCs) from each OA patient and non-OA control.

Duke Validation Cohort

A third independent population was selected from among 159 participants (118 female, 41 male) enrolled in the NIH-sponsored Strategies to Predict Osteoarthritis Progression (POP) study (N=56) (19). Excluding patients with BMI ≥33, this study provided data on 36 patients. The POP study patients were recruited from the Durham, NC area under a protocol approved by the Duke University IRB and in accordance with the policies of the Duke University Medical Center. Participants were recruited who met the ACR criteria for symptomatic OA of at least one knee (18). Exclusion criteria were similar to above, as reported (19). For gene expression studies by TaqMan QPCR analysis from the Duke patients, PAXgene samples were utilized.

Sample collection

Blood samples in the NYUHJD cohorts were collected in pyrogen-free heparinized tubes for isolation of PBMCs and plasma, as well as serum collection tubes. PBMCs were isolated as described (20). PAXgene tubes were collected for RNA in all three cohorts. Blood was processed within 30 minutes of collection.

Total RNA isolation and gene expression analysis

Total RNA from PBMCs was isolated using Qiagen RNeasy column. Total RNA from PAXgene tubes was isolated using a PreAnalytiX blood RNA isolation kit (Qiagen Biosciences). Isolated total RNA was stored in aliquots at −70°C. It should be noted that in contrast to the PBMCs isolated from the NYUHJD patients by Ficoll-Hypaque density centrifugation, which were mononuclear leukocytes, neutrophils are the predominant PBL cell type that contribute to mRNA isolated from PAXgene tubes.

Real Time PCR (QPCR)

Total RNA (1 μg) was primed using oligo (dT) 18 primers and cDNA synthesized using the Clontech cDNA synthesis kit (Clontech). Predesigned TaqMan primer sets were purchased from Applied Biosystems. Real-time PCR reactions were performed as described (21).

Labeling and hybridization of microarray

Five micrograms of total RNA were used for the first strand synthesis using a T7-(dT)24 oligomer. Complementary RNA was synthesized using ENZO kit (Affymetrix) and purified using the Qiagen RNeasy kit, fragmented at 95°C for 35 min for target preparation, and hybridized against Human U133A microarray (Affymetrix).

Normalization of microarray data

Raw expression data from array scans were pre-processed using the Affy package from Bioconductor (22), normalized using the RMA method of Irizarry et al (23), and analyzed further using the R language for statistical computing (24). Genes with less than 70% reliable expression measurement (“present” call in Affymetrix) were filtered out, leaving 8700 genes for analysis.

Statistical methods

Significance Analysis of Microarrays (SAM) implemented in the R package samr was used to find genes differentially expressed between non-OA controls and OA subjects (25). SAM uses the False Discovery Rate (FDR) approach to control for multiple comparisons (26).

In order to identify and validate OA subclasses, the training/test approach was used. In the training set, complete-linkage hierarchical clustering was used to identify OA subclasses based on the expression of a pre-selected set of cytokine genes. Cluster and TreeView packages were used for cluster analysis (27). Nearest centroid classification was used to classify test set samples into subclasses of OA. Briefly, nearest centroid classification method computes a centroid for each OA subclass in the training set, based on the expression of the pre-selected set of cytokine genes, and assigns each test set sample to the OA subclass whose centroid it is closest to in squared distance.

The “nearest shrunken centroid” method implemented in the R package pamr (Prediction Analysis of Microarrays) was used to construct a gene expression classifier of OA and that of OA subclasses (28). The nearest shrunken centroid classification makes one important modification to standard nearest centroid classification. It shrinks each of the class centroids toward the overall centroid for all classes by an amount called the threshold. This shrinkage consists of moving the centroid towards zero by threshold, setting it equal to zero if it hits zero. Shrinkage makes the classifier more accurate by reducing the effect of noisy genes and performs automatic gene selection. The value of the threshold is determined by cross-validation in pamr. The predictive accuracy of the gene expression classifier of OA subclasses was estimated using the test set.

The method of pre-validation (PV) was used to assess the significance of the microarray biomarker of OA while controlling for age, BMI and other demographic characteristics in a logistic regression model (29). PV outputs a prediction for each patient based on the model that is estimated without using that patient's data, thus reducing bias associated with reuse of the data.

Box-Cox transformations were used to normalize variables prior to using statistical tests and regression models. All statistical analyses were performed using the R language for statistical computing (24) and Bioconductor packages (22). Linear (logistic) regression methods were used to compare the groups (OA^IL-1, OA^nl). Spearman's correlation coefficient r was used for different correlation analysis.

Study cohorts

Three independent cohorts of patients were examined: two from NYUHJD (“Learning Cohort” and “NYUHJD Validation Cohort”) and the “Duke Validation Cohort.” The strategy for identifying differentially expressed PBMC genes in the NYUHJD Learning Cohort, and validating expression of those genes in the NYUHJD Validation and Duke Validation Cohorts, is outlined in Figure 1.

NYUHJD Learning Cohort

We recruited 45 OA patients and 25 non-OA healthy controls into the NYUHJD Learning Cohort: OA patients were significantly older (p <0.001) and had higher BMI than non-OA healthy controls (p=0.017). Approximately 25% of the subjects were males and 80% were non-Hispanic among OA patients as well as non-OA healthy controls. To address differences in age and BMI we selected a sub-cohort (designated the Matched Sub-cohort) of 19 OA patients and 17 healthy controls, frequency-matched on age, gender, ethnicity and BMI, for further analysis.

PBMC gene expression profile as a diagnostic combinatorial biomarker: Pre-validated Microarray Predictor

Using the nearest shrunken centroid classification (28), we constructed a microarray combinatorial biomarker diagnostic of OA based on the NYUHJD Learning Cohort microarray data. We performed this analysis based on the entire NYUHJD Learning Cohort and repeated it based on the matched sub-cohort. We fitted a logistic regression model of OA with the microarray biomarker alone and in combination with BMI, age and ethnicity as independent variables. To reduce bias associated with re-use of the data, we used the method of pre-validation when assessing the significance of the microarray predictor (29). The pre-validated microarray predictor was significant when alone (p=0.0002) and while controlling for the subject characteristics (p=0.0017). Significant predictors included the differential gene expression profile microarray [odds ratio (OR) 8.51, 95% confidence interval (CI) 2.03-35.7], age (OR 1.10, 95% CI 1.02-1.19) and BMI (OR 1.23, 95% CI 1.02-1.49).

Significant genes

In the NYUHJD Learning Cohort, we identified 332 upregulated genes and 193 downregulated genes in OA using a FDR of 5% and with at least a 1.25-fold change. Using the matched sub-cohort, we identified 828 upregulated genes and 448 downregulated genes in OA by the same criteria. Seventy-three percent of the upregulated genes and 62% of down-regulated genes identified using data on all subjects were also found to be significant based on the matched sub-cohort data. Among the most significant genes were junB, CDC42, Granulin, MMP-9 and cyclin D2. Using the method of cross-validation, we estimated that this gene expression pattern has sensitivity of 89% and specificity of 76% for an OA diagnosis [see Supplemental Tables 1 and 2 (ST1 and ST2)].

PBMC inflammatory gene expression identifies OA subclasses

Given the role that inflammation in joint tissues has been shown to play, we determined whether PBMCs were “activated,” reflecting exposure to inflammatory stimuli as these cells traversed diseased synovium and bone. Therefore, to identify subclasses of OA, we randomly divided the NYUHJD Learning Cohort OA samples into a separate training set (n=21) and test set (n=20). The samples had been pre-stratified to ensure that the two sets were frequency matched on sex, ethnicity and BMI. In the training set, we used complete-linkage hierarchical clustering to identify subclasses of OA based on the expression of 21 pre-selected cytokine genes (see Table 1). Two subclasses of OA were identified in the training set: cytokine overexpressors (OA^IL-1) and cytokine underexpressors (OA^nl). OA^IL-1 class had elevated levels of IL-1β (up 6.56-fold), IL-8 (up 2-fold), COX-2 (up 2.75-fold), and chemokines GRO2 (up 5.37-fold), macrophage inflammatory protein-1α (MIP-1α) (up 6-fold) and MIP-1β (up 2.75-fold). Based on the gene expression of the 21 cytokines, we computed the centroids C1 and C2 of OA classes in the training set and labeled each OA test sample as OA^IL-1 or OA^nl according to whether it was most correlated with C1 or C2 (nearest centroid classification). Therefore, we defined OA classes using cluster analysis based solely on the expression of the 21 cytokine genes. Presented in Figure 2A is the hierarchical clustering diagram of all OA samples based on the expression of the 21 pre-selected cytokine genes. The centroids of the two classes of OA and on-OA controls are depicted in Figure 2B.

Cross-validated gene expression (excluding 21 cytokine genes) predictor of OA classes

Using the method of nearest shrunken centroids, we constructed a 10-fold cross-validated predictor of these cluster-defined subclasses of OA based on the training set and excluding the 21 cluster-defining inflammatory genes, in order to find whether other genes or probe sets can stratify subsets of OA patients. The resulting microarray predictor of OA subclasses consisted of eight genes, including cytokine genes such as IL-1β (a different probe set) and TNF inducible proteins [see Supplemental Table 3 (ST3)]. We evaluated the performance of the resulting 8-gene predictor of OA subclass using the test set. The training, cross-validated and test set error of this gene expression predictor of OA subclass was estimated to be over 96%.

Analysis of differentially expressed cytokines by real time RT-PCR

TaqMan real-time QPCR was performed to confirm the findings based on microarrays. The pre-designed primers were purchased from Applied Biosystems. In this study, we analyzed expression of three different house-keeping genes, actin (act), glycerolaldehyde-3 phosphate dehydrogenase (GAPDH) and ribosomal protein large PO (RPLPO). Among the three, only GAPDH showed no significant difference between the non-OA control and OA group (data not shown). The “inflammatory” genes IL-1β, TNFα, IL-6, IL-8, COX-1 and COX-2 were analyzed. Relative expression levels of each gene were calculated using the delta-delta comparative C_T method (21). As observed by microarray data (Figure 3), significantly elevated levels of IL-1β, TNFα, IL-6, IL-8 and COX-2 expression were observed in the OA^IL-1 as compared to the OA^nl subclass (Figure 3A). Relative expression levels of inflammatory genes quantitated by QPCR in OA PBMCs (NYUHJD Learning Cohort) were correlated with pain assessment scores in the linear regression model while controlling for age, BMI, gender and ethnicity. Among the inflammatory genes, only IL-1β and IL-6 levels were significantly associated with pain (r=0.38, p=0.02 and r=0.44, p=0.04, respectively).

Elevated inflammatory gene expression in OA PBMCs confirmed in two independent cohorts

To further confirm the elevated levels of both genes as diagnostic biomarkers for a specific OA subset, we used two different additional independent cohorts of OA patients 1) The Duke Validation Cohort consisted of 36 OA patients and no controls; the Duke OA patients had higher mean BMI than the patients in the other cohorts (Table 2). 2) The NYUHJD Validation Cohort consisted of 86 OA patients and 12 non-OA controls (Table 2). As noted, we defined OA^IL-1 and OA^nl subclasses in the validation cohorts based on the samples IL-1 expression using the 90% quantile method.

As observed in NYUHJD Learning Cohort, elevated expression of IL-1β and TNFα was confirmed in both the NYUHJD and Duke Validation Cohorts (Figure 3B-C). For the NYUHJD Validation Cohort, OA subjects were defined as cytokine overexpressors (OA^IL-1 class) if they had elevated expression (at least 2-fold) of IL-1β. Some characteristics of the OA^IL-1 and OA^nl subclasses of the NYUHJD Learning Cohort were also confirmed in the NYUHJD Validation Cohort. Specifically, pain assessment scores were significantly higher in the OA^IL-1 vs OA^nl subclass by all three measures (WOMAC pain, 52.16 ± 23.67 vs. 37.05 ± 19.81, p= 0.0034; WOMAC total, 161.78 ± 64.01 vs. 112.18 ± 55.84, p=0.0004; VAS pain, 61.80 ± 25.80 vs. 42.77 ± 28.17, p=0.0069) in the NYUHJD Validation Cohort. Additionally, IL-1 levels were significantly associated with pain assessment scores by all three measures in the linear regression model while controlling for age, BMI, gender and ethnicity (adjusted p=0.0009, 0.01 and 0.02 for pain VAS, WOMAC total and WOMAC pain, respectively) in the NYUHJD Validation Cohort.

The OA^IL-1 subset is at increased risk for radiographic progression

Eighty-six patients in the NYUHJD Validation Cohort completed a 24-month study of radiographic progression. Comparing standardized fixed-flexion radiographs read independently by two radiologists (LR, JB), we determined change in JSW (both in millimeters and as >30% change from baseline) and KL score between visit 0 and 24 months. For KL scores, a change in one grade was considered to represent progression. Baseline IL-1 mRNA expression was used to segregate OA patients into two groups, OA^IL-1 and OA^nl, compared for progression based on JSW at baseline and >30% JSN at 24 months. OA^IL-1 patients exhibited greater JSN at 24 months than the OA^nl group (0.76 vs 0.15 mm, p <0.02) (Table 3). Since the mean JSW at baseline differed between OA^IL-1 (3.8±1.7 mm) and OA^nl (2.1±2.1 mm) groups, we separately analyzed those OA patients who had ≥3 mm baseline JSW. In this sub-analysis (Table 3), JSN at baseline remained higher in the OA^IL-1 group compared to OA^nl group, and OA^IL-1 patients exhibited greater JSN at 24 months than the OA^nl group (0.89 vs 0.38 mm, p=0.0507). Logistic regression analysis adjusted for BMI, age and gender also showed increased progression in the OA^IL-1 group, as assessed by decrease in JSW (i.e., JSN) >30%, observed in 19 of 86 patients [OR=3.2 (1.16-8.66), p <0.03). Thus, using percent change in JSW (JSN >30%), we could identify a top quartile of “fast progressors” who were identified at baseline by increased PBMC expression of IL-1β. Additionally, using the same logistic regression model, the change of KL score in the OA^IL-1 group was significantly greater (p<0.004) compared to the OA^nl group over 24 months. After controlling for BMI, age and gender, OA^IL-1 also exhibited clinical features that differed from OA^nl patients, with higher WOMAC pain, decreased physical function and higher VAS pain scores (p=0.0069).

PBMC overexpression of TNFα also predicts OA progression

Analysis of the 86 “completer” patients also revealed that patients who are high expressors of IL-1β (OA^IL-1) exhibit increased PBMC gene expression of TNFα (OA^TNFα) (significant positive correlation r=0.78, p <0.001). We therefore also stratified patients based on PBMC TNFα expression (QPCR) into OA^TNFα (overexpressors) and OA^nl(TNF) (whose TNFα expression was comparable to controls). Similar to the OA^IL-1 group, the OA^TNFα group demonstrated increased risk of progression in the signal knee compared to OA^nl(TNF) group, as assessed by JSN greater than 30% (OR=8.9; 95% CI 2.57–30.81; p<0.0005).