Lesional gene expression profiling in cutaneous T-cell lymphoma reveals natural clusters associated with disease outcome

Cutaneous T-cell lymphoma (CTCL) refers to a heterogeneous group of non-Hodgkin lymphomas that share the common feature of infiltration of the skin with malignant T cells.²  The most common form of CTCL is historically known as mycosis fungoides (MF), and it has a variable clinical course. While many patients have a relatively indolent course, others experience rapidly progressive and often lethal disease. Like other lymphomas, prognosis is linked, albeit imprecisely, to clinical stage at presentation. Patients with stage IA disease have normal life expectancies, while those in stage IB have median survivals of 12.8 years, and median survival of stage IIB disease is only 3.2 years.³  Other variants of CTCL are less common, and include diseases with more or less aggressive clinical courses.^2,4  Sézary syndrome (SS) is the leukemic variant of CTCL with a 5-year survival of only 11%.⁵ 

In early disease, activated reactive but benign T cells greatly outnumber malignant cells in skin, and the clinical appearance is that of a patch or a thin plaque of an inflammatory skin disease.⁶  The cytokine milieu of the involved skin is dominated by type 2 T-cell cytokines, and the infiltrating cells are typically positive for cutaneous lymphocyte antigen (CLA) and CCR4.⁷  As disease progresses, the ratio of malignant cells to reactive cells increases, and thick plaques and tumors are characterized by loss of epidermotropism, increase in dermal lymphocytic infiltrate, and size and atypia of cells.⁸  The etiology of CTCL remains obscure, and we have incomplete understanding concerning the initiation of malignant disease and the role of processes such as genomic instability or antigenic stimulation in disease progression. At this time, it is impossible to predict which patients with early-stage CTCL will progress, and which will continue to manifest indolent skin-limited disease.⁹ 

Transcriptional analysis of cancer has proven to be a powerful and increasingly useful tool in biomedical research.^10-12  These patterns of gene expression have revealed unappreciated complexity in histologically homogeneous diseases like large-cell B-cell lymphoma. A goal of many studies is to link gene expression patterns revealed by transcriptional profiling to both understanding pathophysiology of disease and predicting prognosis and responsiveness to therapy.^13-15  These approaches have been applied to CTCL in a limited way to characterize leukemic cells in SS,¹⁶  to distinguish MF from benign dermatoses,¹⁷  and to predict responsiveness to interferon.¹⁸ 

In the present study, we use transcriptional analysis to gain insight into the pathophysiologic and prognostic features of CTCL. We evaluated 62 patients with various stages of CTCL for whom detailed clinical history was available. Lesional biopsies were analyzed via transcriptional profiling using oligonucleotide microarrays (Affymetrix HG-U133 chips; Santa Clara, CA). We present a robust unsupervised analysis of these samples with clustering into 3 groups, and explore the molecular pathways implied by the resultant clustering. We purposefully sought a rigorous assessment of optimal cluster number in order to better appreciate categories of disease process. We used supervised analysis to determine signatures of early- and late-stage CTCL and to evaluate responsiveness to therapy and clinical course. Taken together, these studies shed light on functional pathways in CTCL and suggest an approach to further predictive studies.

All patients were enrolled in a study protocol approved by the Dana-Farber Cancer Institute's Institutional Review Board. Informed consent was obtained in accordance with the Declaration of Helsinki. Patients were recruited from the Cutaneous Lymphoma Clinic at the Dana Farber Cancer Institute (DFCI)/Brigham and Women's Hospital (BWH). One sample was obtained through collaboration with MD Anderson. A total of 63 6-mm punch biopsies from lesional skin were collected from 62 patients between January 26, 2003, and June 1, 2005. One patient with patch/plaque disease was biopsied on 2 separate occasions from similar lesions at different body sites in order to provide methodologic verification through consistent clustering patterns. Patients enrolled were at various stages of disease, including 13 patients with limited patches and plaques covering less than 10% of body surface area (stage IA), 29 with generalized patches and plaques covering more than 10% of body surface area (stage IB), 8 with tumors (stage IIB), and 12 with generalized erythroderma (stage III). One-third of patients presented for diagnostic biopsy and were treatment naive, whereas the remainder had already started treatments with various therapies. Patients were diagnosed and staged based on clinical history and physical exam, biopsy findings, flow cytometry, and computed tomography (CT) studies. An extensive chart review was conducted to collect information on clinical parameters and outcomes. Data were stored in a study database. Mean age of patients with CTCL in this study was 60 years (range, 25-90 years). A total of 60% of patients were male.

Skin biopsies of CTCL skin lesions taken with a 6-mm punch were obtained and immediately snap-frozen in liquid nitrogen. Tissue was powdered in liquid nitrogen (Cryo-Press; Microtec Co, Chiba, Japan), and total RNA was extracted using Trizol (Invitrogen, Carlsbad, CA) according to the manufacturer's instructions. RNA was then labeled using the standard Affymetrix in vitro transcription labeling protocol, and 10 μg cRNA was hybridized to Human Genome U133A GeneChip oligonucleotide arrays, which consist of 22 283 probe sets representing approximately 13 500 named genes.

Gene expression data were processed using Robust Multi-Chip Average (RMA)¹⁹  as implemented in the R statistical package provided by Bioconductor (Seattle, WA).^20,21  RMA data processing includes background correction, quantile normalization, and expression summarization.

Genes were filtered in order to eliminate those with cross-sample variation using a variation filter²²  estimated on 14 HeLa cell lines and 5 Strategene samples (Stratagene, La Jolla, CA; Figure 1). We used this filtering to define 2 gene sets, 1 with 6997 genes (α = 0.999) and another with 4204 genes (α = 0.9999). Unsupervised analyses used the 4204-gene set, reflecting the more stringent statistical criteria required when samples and genes are allowed to cluster naturally. Supervised analyses used the 6997-gene set.

In this analysis, 2 unsupervised clustering techniques were used: hierarchical clustering (HC)²³  and self-organizing maps (SOMs).²⁴  The 2 clustering techniques were compared using consensus clustering²⁵  in order to assess the stability of the identified clusters and to evaluate sensitivity to sampling variability. Confusion matrices for gene as well as sample assignment were generated to assess agreement between clusters produced by the different clustering algorithms.

All heat maps for visualization were generated using DChip.²⁶  Fold changes were calculated based on average expression for class versus nonclass.

Statistical tests to compare predefined classes of samples based on phenotypic differences were performed using GenePattern (Broad Institute, MIT, Cambridge, MA, http://www.broad.mit.edu/cancer/software/genepattern/).^27,28  In particular, we used the Comparative Marker Selection suite, which computes statistical significance for each gene and corrects for multiple hypothesis testing (MHT).²⁹  A total of 10 000 permutations were performed for the estimation of empirical p-values, which were then corrected for MHT by the false discovery rate (FDR) procedure.³⁰  Genes with an FDR value of less than or equal to .05 were considered significant in our study.

A total of 2 different algorithms were used to evaluate contribution of gene pathways to the transcriptional differentiation of samples.

Gene Ontology (GO) is a collaborative and comprehensive gene annotation resource compiled by the Gene Ontology Consortium.³¹  GO annotations were obtained from Affymetrix and analyzed using GeneMerge.³²  GeneMerge uses defined class gene markers and calculates enrichment for GO pathways.

Gene set enrichment analysis (GSEA) was performed by generating a ranked list of genes differentiating defined classes and comparing this ranked list with specified gene sets.³³  Genes were ranked by signal to noise ratio. An enrichment score was generated for each gene set. Statistical significance was then estimated using phenotype-based permutations, with the attained P values adjusted for MHT. A total of 539 gene sets generated from several independent sources were used in our analysis. Among the resources used to compile these gene sets were (1) Biocarta (Biocarta, San Diego, CA), with 274 pathways involved in adhesion, apoptosis, cell activation, cell-cycle regulation, cell signaling, cytokines and chemokines, developmental biology, hematopoiesis, immunology, metabolism, and neurosciences; (2) GenMAPP (Gene MicroArray Pathway Profiler; www.genmapp.org), a resource for metabolic pathways, cell signaling processes, and physiologic functions, with 67 pathways; (3) Signal Transduction Knowledge Environment (STKE), a branch of the journal Science that serves as an online resource for signal transduction research and includes 29 signal transduction pathways; (4) SigmaAldrich (Sigma-Aldrich, St Louis, MO), with 12 pathways, involved in the cell cycle, mitogen-activated protein (MAP) kinase signaling, and apoptosis; (5) GO, from which were culled 23 pathways; and (6) the GEArray pathway-specific arrays made by SuperArray (SuperArray Bioscience, Frederick, MD), from which were derived 3 pathways involved in calcineurin-NFAT signaling, breast cancer and estrogen signaling, and bcl2 regulation. A total of 63 manually curated pathways involved in mitochondrial function and metabolism were also included in the analysis. An additional 59 pathways were generated from the current literature.

We used HC and SOM algorithms to differentiate biologically meaningful subsets of CTCL with similar transcription profiles. Based on consensus analysis, we determined that the most robust clustering structure included 3 distinct clusters (Figure 2). This clustering showed a high level of agreement between clusters produced by HC and SOMs, with approximately 89% of samples and more than 90% of genes assigned to the same clusters (Figure 3).

One patient was biopsied at 2 different time points and 2 different body sites to verify reproducibility of the clustering. Clustering performed with each sample separately resulted in the same cluster assignment (cluster 1) both times.

We next examined the gene sets involved in the 3-cluster structure produced by consensus clustering. Using standard differential analysis, we identified 593 genes with a P value less than or equal to .025 and fold change greater than or equal to 2 (Figure 4). We then further characterized the clusters identified by unsupervised analysis based on sample phenotype as well as gene identity.

Sample clustering appeared to correlate with clinical stage (Figure 5). Cluster 1 contained a mix of patients with stages IA, IB, IIB, and III disease. This cluster included 75% (6 of 8) of all stage IIB (tumor stage) samples. Cluster 2 contained primarily stage IA and IB disease, with only 2 patients with higher-stage disease. Cluster 3 contained only 1 stage IA patient and a relatively large proportion of stage III disease (58%, or 7 of 12 stage III patients).

Among patients with CD30⁺ disease, all 5 were in cluster 1. Interestingly, none of the 5 patients with CD8⁺ disease were found in cluster 1. Of the patients with CD8⁺ disease, 2 were distributed to cluster 2 and 3 to cluster 3. Folliculotropic disease was found to be somewhat more predominant in cluster 3, with 2 of 7 patients found in cluster 1, 1 of 7 patients in cluster 2, and 4 of 7 patients in cluster 3.

The 3 clusters were also characterized by differences in survival. Cluster 1 included 4 deaths: 3 stage IIB patients and 1 stage IB patient. Deaths occurred at 6 months after biopsy for the stage IB patient, and at 3 months, 3 months, and 4 months, respectively, for the 3 stage IIB patients. In all cases, death occurred after developing metastatic and/or leukemic involvement and failing multiple treatment regimens, including systemic therapies. Cluster 3 included 2 deaths at 5 and 9 months, respectively, both of whom were patients with SS. In comparison, the single death in cluster 2 was due to an unrelated medical condition, occurring 20 months after biopsy.

We first examined the impact of prebiopsy treatments on the clusters, and did not find a significant correlation between particular treatment regimens and clusters, suggesting that the clustering patterns were not driven by treatment signatures (Table 1). Examination of treatment response showed that clusters 1 and 3 are associated with more aggressive, less treatment-responsive disease. Recalcitrant disease across stages was more commonly found in clusters 1 and 3. In particular, 2 patients assigned to cluster 1 have since undergone bone marrow transplantation as a last resort to control their disease, and 1 patient had progressed from stage IA disease to large-cell transformation after the biopsy was taken, with subsequent persistent disease despite several treatment regimens, including denileukin diftitox (ONTAK) total skin electron beam therapy, and a suberoylanilide hydroxamic acid (SAHA) trial. By contrast, patients in cluster 2 generally exhibited extremely good response to treatment across all stages of disease. All but 3 patients are currently in remission or under good control on topical therapies, including steroids, tacrolimus, bexarotene, and/or phototherapy, administered as PUVA, narrow-band UVB, or natural sunlight. The only stage IIB patient went into remission on PUVA alone. Among stage IB patients, those who were more difficult to treat, characterized by longer time to remission, failure to respond to standard treatment regimens, and use of therapies more commonly reserved for more advanced disease, tended to fall into clusters 1 and 3. Stage IB patients whose disease was more amenable to treatment tended to fall into cluster 2. Among stage IB patients in cluster 2, 9 of 11 were maintained on intermittent topical therapies, including topical steroids, bexarotene gel, and/or phototherapy, without progression of disease. The stage IB patients not in cluster 2 included 1 patient in cluster 1 who developed metastatic disease and died despite the use of ONTAK, oral bexarotene, total skin electron beam radiation, extracorporeal photopheresis, interferon, and single- and multiple-agent chemotherapy. A total of 7 of the 18 stage IB patients in clusters 1 and 3 were placed on systemic therapies, including oral bexarotene, prednisone, extracorporeal photopheresis, interferon, ONTAK, methotrexate, chemotherapy, and SAHA.

Although differences in survival among clusters did not reach significance, differences in event-free survival, such as disease progression or large-cell transformation, did. Event-free survival curves are shown in Figure 6.

We identified dominant genes differentiating the 3 CTCL clusters by applying the t statistic and determining fold change.

Cluster 1, composed of a mix of clinical stages, showed significant up-regulation of many genes involved in immune response and T-cell activation (Table 2). Notable among these are all 3 components of the IL-2 receptor, lymphocyte-specific protein tyrosine kinase (LCK), PIK3CD, T-cell receptor α and β chains, CD3, ζ-chain–associated protein kinase (ZAP70), FYN-binding protein (FYB), T-cell receptor–interacting molecule (TCRIM), CD69, and linker for activation of T cells. SH2D1A, along with SLAMF1 (fold change, 2.99; P = .003), has been implicated in TCR activation; deficiency has been associated with the development of fatal X-linked lymphoproliferative disorder after Epstein-Barr virus (EBV) infection.³⁴  CD8⁺ T-cell markers, including granulysin, were also up-regulated in Cluster 1.

LEF1 and TCF7, both of which were up-regulated in cluster 1, are downstream members of the WNT signaling pathway. A recent study has shown that several isoforms of LEF1 and TCF7 are up-regulated in naive CD8⁺ T cells, and down-regulated by TCR activation.³⁵  Both LEF1 and TCF7 have been shown to be expressed in lymphocytic leukemias, including a subset of peripheral T-cell lymphomas and T-cell acute lymphocytic leukemia (T-ALL).^36,37  Based on activation markers, studies suggest that LEF1 and TCF7 are up-regulated in the setting of T-helper 1 (Th1) as compared with Th2 activation.³⁷ 

Also up-regulated in cluster 1 were several central memory T-cell markers and genes involved in lymphocyte homing: SELL (L-selectin), CCR7, and FUT7 (fucosyltransferase 7), an important regulator of selectin ligands, including CLA.^38,39 

Tumor necrosis factor (TNF) pathway genes, such as TNFRSF17, TNFRSF14, TNFSF4, TNFRSF7, LTB(lymphotoxin β), TNF, TRAF1, T3JAM, and TNFAIP8, are also up-regulated in this cluster. TNFRSF14, also known as herpesvirus entry mediator (HVEM), has been shown to interact with TNFSF14(also known as LIGHT) and LTA (lymphotoxin A) resulting in costimulation of T cells.^40,41 BIRC1, BIRC3, and BIRC5 were also significantly up-regulated in this cluster. BIRC3 is a member of the inhibitor of apoptosis (IAP) gene family, and inhibits apoptosis by binding to TNF receptor–associated factors TRAF1 and TRAF2.⁴²  The previous work of Tracey et al¹⁷  also emphasized TNF antiapoptotic pathway activation in MF, consistent with these results.

Interestingly, several B-cell–related genes were found to be up-regulated in this lymphocyte activation cluster. Immunoglobulin genes are among the most up-regulated by fold change, although the significance of this result (P ≤ .05) is lower than that of many other genes. POU2AF1 has been implicated in CLL.⁴³ 

Cluster 2, consisting of mostly stage IA and IB patients, showed significant up-regulation of several genes involved in keratinocyte and epidermal differentiation and proliferation. Selected genes are shown in Table 3. Among the genes up-regulated in this cluster were keratinocyte-related genes, including loricrin (LOR), late cornified envelope 2B (LCE2B), and keratin 2A (KRT2A). Members of the WNT signaling pathway were also up-regulated, including catenin β-interacting protein 1 (CTNNBIP1), transcription factor 7-like 2 (TCF7L2), and transcription factor 7-like 1 (TCF7L1).

Genes up-regulated in cluster 3 included several genes related to keratinocyte function, but also included T-cell–specific genes. Selected genes are shown in Table 4. Up-regulated genes included keratins 6 and 16, indicating keratinocyte hyperproliferation.⁴⁴ SERPINB3, CEACAM1, and WNT 5A were also significantly up-regulated. These genes have been previously associated with benign hyperplasia in the setting of psoriasis as compared with squamous cell carcinoma.⁴⁵ 

We also analyzed genes in our clusters according to established functional and molecular pathways. As our results have already shown, several pathways are easily identified ad hoc. However, we used 2 rigorous techniques to identify enriched pathways because with a large set of genes, it is often difficult to pick out patterns relevant to pathophysiology of disease, as many genes, despite relatively low fold changes, may contribute as members of an up-regulated set rather than individually. We compared GO pathways and a set of pathways curated by the Broad Institute at MIT to our dataset using the GeneMerge and GSEA tools, respectively. GeneMerge scores GO pathways according to a predefined list of cluster markers. GSEA, in comparison, uses processed expression data and class phenotype information to discover pathway enrichment.

Analysis of GO pathways gave more strength to the patterns suggested by the cluster marker gene lists. The top 15 pathways for each cluster, ranked by P value and corrected for multiple testing, are shown with the heat map in Figure 4. Cluster 1 showed up-regulation of pathways involving immune response, endogenous and exogenous antigen presentation and processing, major histocompatibility complex (MHC) receptor activity, chemokine receptor activity, TCR binding, and TNF-receptor binding. In addition, proliferative pathways are also highlighted, including mitosis, cell division, nucleus, chromosome, DNA repair, T-cell proliferation, and B-cell proliferation. In contrast, cluster 2 revealed up-regulation of several pathways involved in extracellular matrix, collagen, epidermis development, transcription, keratin filament, WNT signaling, melanin metabolism, keratinocyte differentiation, and cornified envelope. Cluster 3 showed up-regulation of many varied pathways, including angiogenesis, oxygen binding, inflammatory response, chemokine activity, chemotaxis, interleukin-1 receptor binding, and response to virus.

GSEA resulted in the identification of significantly enriched pathways in only 1 cluster—cluster 1. For all pathways containing 5 or more genes, 8 were significantly enriched in cluster 1, with FDR less than or equal to 0.05, and 3 genes with FDR less than 0.01 (Table 5). The enriched pathways involve T-cell function, TCR signaling, and T-cell activation. Also enriched was a pathway up-regulated in GATA1 mutant megakaryoblasts harvested from mice, previously implicated in leukemogenesis.⁴⁶ 

In some cases, patients had 2 biopsies available from the same lesion, 1 processed for gene expression profiling and 1 processed for histologic analysis. We examined the histology in these cases for any correlation to clustering. In particular, we were interested to determine if the relative number of T cells present in a biopsy correlated with clustering. Figure 7 shows histology and CD3 immunostaining of representative biopsies taken at the same site as for the clustering. As illustrated in the images, there appears to be variability in the number of lymphocytes with respect to cluster, suggesting that the number of T cells in a biopsy does not dominate the clustering.

Differential analysis was performed on a number of different phenotypes in order to identify genes associated with particular clinical differences. We first evaluated markers differentiating clinical stage. In a comparison of patch/plaque disease (stages IA and IB) compared with more advanced disease (stages IIB and III), we identified 480 genes with FDR less than or equal to 0.05, and 50 genes with FDR less than or equal to 0.01 (Figure 8). Genes ranked by signal-to-noise ratio contributing to differentiation of clinical stages are shown in Table 6. Genes up-regulated in earlier stage disease include several involved in epidermal development, WNT signaling pathway, and keratinocyte development. Genes up-regulated in later stage disease include several genes involved in immune response, mitosis, cell proliferation, DNA replication, apoptosis, response to virus, chemotaxis, and interferon (IFN) and TNF-alpha signaling. Although GSEA was unable to show any significant enrichment of evaluated pathways, GO analysis did show significant up-regulation in particular pathways. The top 20 up-regulated GO pathways are shown in Figure 8.

We also evaluated patients in terms of response to treatment, identifying those with disease persistence or progression despite stage-appropriate treatment as difficult to treat, and those whose disease was responsive to stage-appropriate therapies as easy to treat. Transcriptional profiling revealed a treatment response signature shown in Table 7. Changes are similar to those differentiating early- and later-stage disease: genes down-regulated in poor treatment response were involved in pathway extracellular matrix, WNT signaling pathway, epidermal development, and frizzled signaling, whereas up-regulated genes included those involved in mitosis, immune response, response to virus, apoptosis, and T-cell activation.

Comparison of patients with leukemic disease versus patients with no evidence of peripheral blood involvement identified only 3 genes with significant fold changes based on multiple testing. These genes included 2 genes up-regulated in leukemic disease, lymphocyte-activation gene 3 (LAG3; fold change 4.24) and B-factor (BF; fold change 2.36), and 1 gene down-regulated in leukemic disease, collagen XI A (COL11A1; fold change −9.84). LAG3 has been shown to be involved in lymphocyte activation, and has a close relationship to CD4; BF is a component of the alternative pathway of complement activation; and COL11A1 is an extracellular matrix component, also reported to be expressed in some T-cell lines. To further investigate these potential markers of leukemic disease, we performed reverse transcription–polymerase chain reaction (RT-PCR) to evaluate expression of these genes in peripheral blood mononuclear cell RNA samples of patients with CTCL with and without leukemic disease (as opposed to the skin lesion RNA used in our GeneChip analysis). We found no correlation between expression of these genes in blood samples and the presence of peripheral blood involvement with CTCL (data not shown).

Evaluation of SS in comparison with all other patient samples revealed significant down-regulation of only the patched homolog gene (PTCH; fold change −1.72). PTCH encodes the receptor for sonic hedgehog, and functions as a tumor suppressor.

We have previously shown that T-cell receptor excision circles (TRECs), which decrease with T-cell division, are significantly decreased in the peripheral blood of patients with CTCL, indicating increased T-cell turnover.⁴⁷  We have also previously shown that the diversity of the T-cell population appears to be decreased in CTCL, by measuring length distributions of T-cell receptor β variable regions (“spectratyping”).^47,48  We compared these parameters of T-cell immunity with clustering results in 39 of the patients of this study for whom TREC and spectratyping data were available. As shown previously, mean TRECs decreased significantly with increasing clinical stage, and mean spectratype contracture increased, consistent with loss of TCR complexity (Table 8). In this study, we compared TRECs and spectratypes by cluster and found the greatest loss of TRECs in cluster 1 and the least in cluster 2 (Table 8). Loss of T-cell diversity on spectratyping was greatest in cluster 3 and least in cluster 2.

This study shows that gene expression profiling in lesional skin biopsies provides information which may be useful in the evaluation and management of patients with CTCL. We determined the most natural clustering of 62 patients with CTCL, and examined in detail the genes and gene sets defining each group. We further correlated these clusters with other parameters including clinical stage, event-free survival, treatment resistance, TREC levels, and spectratype diversity. Our results demonstrate that grouping patients based on gene expression reveals many interesting similarities and differences with currently accepted clinical staging. This work can serve as a framework for further studies aimed at understanding the diverse spectrum of CTCL in terms of gene expression rather than clinical manifestations alone.

Using consensus analysis, we found 3 clusters to be the most statistically robust grouping of our samples, with approximately one-third of patients in each group. One cluster (cluster 2) was associated with less aggressive disease by all measures. Another cluster (cluster 3) was associated with more extensive disease, decreased event-free survival, and poor response to treatment. A third cluster (cluster 1) was associated with an even poorer prognosis and included many patients in tumor stage. A previous microarray study of lesional CTCL skin biopsies by Tracey et al¹⁷  divided subjects into 2 clusters using hierarchical clustering. These clusters were generally characterized by clinical parameters to be associated with more and less aggressive disease. They did not perform consensus clustering and did not publish any results of clustering with an explicit 3-cluster model. However, when we examined their data in a 3-cluster division for comparison with the results of our consensus clustering, we found that a 3-cluster structure for their data were highly consistent with our results, based on sample as well as gene partitioning (Figure 9).

In addition to examining individual gene differences, we examined groups of genes organized according to functional pathways. Pathways were tested using the GO and GSEA tools, which group changes in gene expression based on molecular pathways in order to uncover changes in cell functions that might otherwise be difficult to delineate because of small fold changes of individual genes. Cluster 1 demonstrated significant up-regulation of many other genes involved in lymphocyte activation, as well as up-regulation of TNF pathway genes. Gene pathways that were predominant in cluster 2, consisting of patients with relatively limited disease, tended to be associated with epidermal development and differentiation. Cluster 3, consisting of patients with more extensive disease, emphasized pathways involved in inflammation and benign epidermal hyperproliferation. Although clusters 2 and 3 appear to be clinically divergent, it is possible that they lie on a clinical spectrum, with cluster 3 tending toward more extensive skin involvement (stages IB and III, with very few stage IA patients) and with cluster 2 tending toward more limited surface area stage IA disease. Cluster 1 was found to be enriched with lymphocyte activation genes, and correlated with poorer outcome, including disease-related deaths and progression despite systemic therapies. This cluster also had lower TREC levels across all clinical stages, consistent with lymphocyte clonal expansion, as well as contracted spectratypes, indicating decreased TCR diversity. This significant correlation between our clustering and TREC/spectratype measurements supports the idea that our clustering does correlate with the biology of the disease.

Biopsies showed that the number of CD3⁺ cells in the epidermal and dermal infiltrates is highly variable across clusters, suggesting that up-regulation of T-cell activation genes is unlikely due solely to increased numbers of lymphocytes.

Markers of clinical stage in the supervised analysis were largely in line with those of the natural clusters, with lower stage disease enriched in epidermal differentiation pathways and higher-stage disease enriched in lymphocyte activation pathways, among others. More work with larger numbers of patients will be required to characterize subgroups more clearly and to isolate causative genes. However, this does appear to be a useful framework for additional study of the pathways and genes involved in clinical stage progression and response to treatment.

This study suggests that gene expression profiles of lesional skin provide information which correlates with disease progression and response to therapy. This information may prove useful in the future in combination with traditionally accepted clinical staging criteria to improve classification of CTCL and to help guide therapy, such as more aggressive treatment of those patients with gene expression profiles enriched in lymphocyte activation genes relative to epidermal differentiation and proliferation genes.

We are indebted to Marianne Tawa and Lisa Batchelder for their help with identifying and recruiting eligible patients for our study. We would also like to thank Helen Kuo for her help extracting RNA. We thank the Broad Institute for help with target preparation and hybridization.

This work was supported by a grant from the National Institutes of Health SPORE program (P50 CA093683 to T.S.K.).

National Institutes of Health

Contribution: All authors participated in designing the research; D.A. and M.D. collected samples; D.A., M.D., and J.S. collected clinical information; J.S. and D.J. examined pathology slides; J.S. and S.M. performed statistical analyses; J.S. and D.J. wrote the paper; and all authors reviewed the final version of the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: David A. Jones or Thomas S. Kupper, Department of Dermatology, Brigham and Women's Hospital, 77 Ave Louis Pasteur, HIM Room 660, Boston, MA 02115; e-mail: dajones@partners.org or tkupper@partners.org.