Overview of the gene expression profiles

We identified 4,149 probes whose expression varied from their median values in these samples by more than 2-fold in at least two of the 75 arrays and analyzed them by two-dimensional hierarchical clustering [18]. The resulting sample dendrogram shows that the samples separate into two main branches (Figure 1A) that in part stratify patients by their clinical diagnosis. The branch lengths in the tree are inversely proportional to the correlation between samples or groups of samples. The diversity in gene expression among the patients with scleroderma is greater than previously shown (Figure 1B) [16], [17] as distinct subsets of scleroderma are evident in the gene expression patterns. Some of these delineate existing classifications, such as the distinction between limited and diffuse, while others reflect new groups. One subset of dSSc patients cluster on the left branch (red) and has gene expression profiles that are distinct from both healthy controls and patients with lSSc (Figure 1B–1C), while a second subset of dSSc skin clusters in the middle of the dendrogram tree (black), and a third set clusters with healthy controls. We found lSSc samples formed a group in the middle portion of the dendrogram and could be associated with a distinct, but heterogeneous gene expression signature that also showed high expression in a subset of dSSc patients (Figure 1H). LSSc samples are partially intermixed with normal controls on the right boundary and with dSSc on the left boundary of the tree, illustrating that their gene expression phenotype is highly variable (Figure 1A). Samples taken from individuals with morphea also grouped together with a gene expression signatures that overlapped with those of dSSc and lSSc (Figure 1). Although nodes can be flipped, we have left the nodes of the dendrogram as originally organized by the clustering software, which places nodes with the most similar samples next to one another. Although, the assignment of samples into particular clusters (Table 3) would not change even if nodes were flipped.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Gene expression signatures in scleroderma.

4,149 probes that changed at least 2-fold from their median value on at least two microarrays were selected from 75 microarray hybridizations representing 61 biopsies. Probes and microarrays were ordered by 2-dimensional average linkage hierarchical clustering. This clustering shows that the dSSc, lSSc, morphea samples form distinct groups largely stratified by their clinical diagnosis. A. The unsupervised hierarchical clustering dendrogram shows the relationship among the samples using this list of 4,149 probes. Samples names have been color-coded by their clinical diagnosis: dSSc in red, lSSc in orange, morphea and EF in black, and healthy controls (Nor) in green. Forearm (FA) and Back (B) are indicated for each sample. Solid arrows indicate the 14 of 22 forearm-back pairs that cluster next to one another; dashed arrows indicate the additional 3 forearm-back pairs that cluster with only a single sample between them. Technical replicates are indicated by the labels (a), (b) or (c). 9 out of 14 technical replicates cluster immediately beside one another. B. Overview of the gene expression profiles for the 4,149 probes. Each probe has been centered on its median expression value across all samples analyzed. Measurements that are above the median are colored red and those below the median are colored green. The intensity of the color is directly proportional to the fold change. Groups of genes on the right hand side indicated with colored bars are shown in greater detail in panels C–H. C. Immunoglobulin genes expressed highly in a subset of patients with dSSc and in patients with morphea, D. proliferation signature, E. collagen and extracelluar matrix components, F. genes typically associated with the presence of T-lymphocyes and macrophages, G. Genes showing low expression in dSSc, H. Heterogeneous expression cluster that is high in lSSc and a subset of dSSc. In each case only a subset of the genes in each cluster are shown. The precise location of each gene in the cluster can be viewed in Supplemental Figure S1.

https://doi.org/10.1371/journal.pone.0002696.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Cluster assignments using the scleroderma intrinsic genes.

https://doi.org/10.1371/journal.pone.0002696.t003

Multiple distinct gene expression programs are evident in each subgroup. Some of these recapitulate the major themes in our prior microarray study of dSSc skin [16] while others reflect gene expression programs not previously observed. A subset of these biological themes and selected genes are discussed below. The entire figure with all gene names is available in the supplementary material (Supplementary Figure S1; Supplementary Data File S1).

Immunoglobulins typically associated with B lymphocytes and plasma cells are expressed in a subset of the dSSc skin biopsies (Figure 1C). Previously we found dense clusters of infiltrating B cells in dSSc by immunohistochemistry (IHC), indicating that these genes may be from infiltrating CD20+ B cells rather than from a small number of infiltrating plasma cells [16].

Previous studies have identified infiltrating T cells in the skin of dSSc patients [19], [20], [21], [22], [23], although an association between T cell gene expression and dSSc was not observed in our prior study [16]. A new result from this study is genes typically associated with T cells are more highly expressed in a subset of the patients (Figure 1F). These genes include the PTPRC (CD45; Leukocyte Common Antigen Precursor), which is required for T-cell activation through the antigen receptor [24], [25], [26], as well as CD2 [27], [28] and CDW52 [29] that are expressed on the surface of T lymphocytes. Also found were CD8A, Granzyme K, Granzyme H, and Granzyme B that are typically expressed in cytotoxic T lymphocytes [30], [31], [32], [33], [34], and CCR7, which is expressed in B and T lymphocytes [35]. Genes induced by interferon (IFIT2, GBP1), genes involved in antigen presentation (HLA-DRB1, HLA-DPA1 and HLA-DMB) and CD74, the receptor for Macrophage Inhibitory factor (MIF), are also present [36], [37], [38], [39], [40]. Genes typically associated with the monocyte/macrophage lineage, B cells and dendritic cells (DCs) were also found in this cluster including Leukocyte immunoglobulin-like receptor B2 and B3 (LILRB2 and LILRB3; Figure 1F) [41], [42]. Finally, genes specific to the monocyte/macrophage lineage such as CD163 are expressed in this group of genes (data not shown) [43].

Genes typically associated with the process of fibrosis were co-expressed with markers of T lymphocytes and macrophages. These genes showed increased expression in the central group of samples that included patients with dSSc, lSSc and morphea (Figure 1E). Included in this set of genes were the collagens (COL5A2, COL8A1, COL10A1, COL12A1), and collagen triple helix repeat containing 1 (CTHRC1), which is typically expressed in vascular calcifications of diseased arteries and has been shown to inhibit TGF−β signaling [44], [45] (Figure 1E). Also found in this cluster was fibrillin-1 (FBN1). The phenotype of the TSK1 mouse, a model of scleroderma, results from a partial in-frame duplication of the FBN1 gene and defects in FBN1 are the cause of Marfan's syndrome (OMIM: 154700).

A surprising result in this study is the differential expression of a ‘proliferation signature’ (Figure 1D). The proliferation signature is defined as genes that are expressed only when cells are dividing [46]. We have previously shown that the proliferation signature, originally identified in breast cancer [47], [48], is composed almost completely of cell cycle-regulated genes [49]. Genes showing increased expression in this cluster include the cell cycle-regulated genes CKS1B, CDKS2, CDC2, MCM8, and E2F7 [49]. The existence of a proliferation signature is consistent with reports demonstrating that a subset of cells in dSSc skin biopsies show high levels of tritiated thymidine uptake indicative of cells undergoing DNA replication [50], [51]; a more recent study has shown increased expression of the cell cycle-regulated gene PCNA in a perivascular distribution [52]. IHC of dSSc skin biopsies with the proliferation marker KI67 also shows proliferating cells primarily in the epidermis (see below).

Another cluster of genes is expressed at low levels in the dSSc skin biopsies but at higher levels in all other biopsies, however it is not clearly associated with a single biological function or process. Included in this cluster are the genes WIF1, Tetranectin, IGFBP6, and IGFBP5 found in our original study [16] with similar patterns of expression (Figure 1G).

Since the skin of lSSc patients does not show any clinical or histologic manifestations at the biopsy site, it was possible that the skin of those patients would not show significant differences in gene expression when compared to normal controls. In fact, lSSc skin showed a distinct, disease-specific gene expression profile. This novel finding demonstrates that microarrays are sensitive enough to identify the limited subset of SSc even when discernable skin fibrosis was not present. There is a signature of genes that is expressed at high levels in a subset of lSSc patients, and variably expressed in dSSc and normal controls (Figure 1H). Included in this signature is the urotensin 2 receptor (UTS2R). The ligand for this receptor, urotensin 2, is considered to be one of the most potent vasoconstrictors yet identified [53], [54], [55]. This finding raises the intriguing possibility that this vasoactive peptide is involved in the vascular pathogenesis of lSSc.

We previously demonstrated that skin biopsies from patients with early dSSc show nearly identical patterns of gene expression at a clinically affected forearm site and a clinically unaffected back site, and the gene expression profiles are distinct from those found in healthy controls [16]. This finding is confirmed in this larger cohort of patients analyzed on a different microarray platform. 14 of 22 forearm-back pairs cluster immediately next to one another indicating that these samples are more similar to each other than to any other sample (Figure 1A). An additional 3 forearm-back pairs grouped together with only a single sample between them (Figure 1A). In total, 17 of 22 (77%) forearm-back pairs show nearly identical patterns of gene expression. This result holds true for patients with lSSc even though neither the forearm or back biopsy sites in lSSc patients are defined as clinically affected [16]. We found 9 out of 14 technical replicates clustered next to one another. The five technical replicates that did not cluster together are likely misclassified as a result of noise in the genes selected by fold change.

Classification of scleroderma using the intrinsic genes from skin biopsies

A list of genes selected by their fold change alone is not ideal for classifying samples because they emphasize differences between samples rather than the intrinsic differences between patients [47], [56]. To select genes that captured the intrinsic differences between patients, we exploited the observation that the forearm-back pairs from each SSc patient show nearly identical patterns of gene expression to select the ‘intrinsic’ genes in SSc. We selected 995 genes with the most consistent expression between each forearm-back pair and technical replicates, but with the highest variance across all samples analyzed [47], [56] (Supplementary Data File S2). Each of the 995 intrinsic genes was centered on its median value across all experiments, and the data clustered hierarchically in both the gene and experiment dimension using average linkage hierarchical clustering. The dendrogram summarizes the relationship among the samples and shows their clear separation into distinct groups (Figure 2A). As a direct result of our gene selection, all forearm-back pairs cluster together and all technical replicate hybridizations cluster together when using the intrinsic genes. Sample identifiers have been colored according to the patient diagnosis: dSSc is red, lSSc is orange, morphea and EF are black, and normal controls are green (Figure 2A). The dendrogram has been colored to reflect the signatures of gene expression that are an inherent feature of the biopsies (Figure 2A–C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Cluster analysis using the scleroderma intrinsic gene set.

The 995 most ‘intrinsic’ genes selected from 75 microarray hybridizations analyzing 34 individuals. Two major branches of the dendrogram tree are evident which divide a subset of the dSSc samples from all other samples. Within these major groups are smaller branches with identifiable biological themes, which have been colored accordingly: blue for diffuse 1, red for diffuse 2, purple for inflammatory, orange for limited and green for normal-like. Statistically significant clusters (p<0.001) identified by SigClust are indicated by an asterisk (*) at the lowest significant branch. A. Experimental sample hierarchical clustering dendrogram. Black bars indicate forearm-back pairs which cluster together based on this analysis. B. Scaled down overview of the intrinsic gene expression signatures. C. Limited SSc gene expression cluster. D. Proliferation cluster. E. Immunoglobulin gene expression cluster. F. T-lymphocyte and IFNγ gene expression cluster. The full figure with all gene names can be viewed in Supplemental Figure S2.

https://doi.org/10.1371/journal.pone.0002696.g002

The gene expression signatures further subdivide samples within existing clinical groups. We find a consistent set of genes that are highly expressed in a subset of the dSSc samples, which occupy the left branch of the dendrogram tree (Figure 2D). These groups have been labeled diffuse 1 (Figure 2A; blue branches) and diffuse 2 (Figure 2A; red branches) as they consistently cluster as two separate groups (c.f. Figures 1 and 2) and have distinct signatures of gene expression. The most consistent biological program expressed across the diffuse 1 and diffuse 2 scleroderma samples is that of proliferation (Figure 2D). We refer to this group broadly as Diffuse-Proliferation. A second group contains dSSc, lSSc and morphea samples on a single branch of the dendrogram tree (Figure 2A, purple branches). The genes most highly expressed in this group are those typically associated with the presence of inflammatory lymphocyte infiltrates as described above and this group has thus been labeled the Inflammatory group (Figure 2F). A third group contains primarily lSSc samples (Limited, orange branches, Figure 2A), which has low expression of the proliferation and T cell signatures but has high expression of a distinct signature found heterogeneously across the samples (Figure 2C). A branch of samples comprised primarily of healthy controls (green branches, Figure 2A) also contains samples from one patient with a diagnosis of dSSc and a patient with lSSc. This group has been labeled Normal-like, since the gene expression signatures in these samples more closely resemble and cluster with normal skin.

Significance and reproducibility of intrinsic clustering

To examine the robustness of these groups, we performed two separate analyses: Statistical Significance of Clustering (SigClust)[57] and consensus clustering [58]. SigClust analysis was performed with the 995 intrinsic genes. At a p-value<0.001 we find five statistically significant clusters. The four major groups of diffuse-proliferation, inflammatory, limited and normal-like groups are each found to be statistically significant (Figure 2A); samples of patient dSSc8 form a statistically significant group of their own in the SigClust analysis (Table 3). Thus, the major groups identified in the hierarchical clustering using the 995 intrinsic genes are statistically significant and cannot be reasonably divided into smaller clusters with the current set of data. The two branches within the diffuse-proliferation group do not reach statistical significance in this analysis even though there are identifiable differences in their gene expression profile.

To perform a second validation of the intrinsic groups, we used consensus clustering [58], which performs a K-means clustering analysis on randomly selected subsets of the data by resampling without replacement over 1,000 iterations using different values of K. Figure 3A shows the consensus over 1,000 iterations for K = 4, 5 and 6. To determine the number of clusters present in the data, we examined the area under the Consensus Distribution Function (CDF; Figure 3B). The point at which the area under the CDF ceases to show significant changes indicates the probable number of clusters (Figure 3C). The largest change occurs between three and four clusters with a slight change between four and five clusters (Figure 3C).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Robustness of sample classification.

The robustness of the sample classifications was analyzed by consensus clustering, which uses multiple iterations of K-means clustering with random restart. 500 subsets of the data were sampled without replacement. The results of consensus clustering and Principal Component Analysis (PCA) applied to the 75 arrays and 995 intrinsic genes are shown. A. Consensus matrices are shown for K = 4, 5 and 6. Cluster numbers are shown and cluster assignments are summarized in Table 3. B. Empirical consensus distribution function (CDF) plots corresponding to K = 2,3,4…10. The ideal number of clusters can be identified when the area under the curve shows minimal increases with increasing K. C. Proportion increase Δ(K) in the area under the CDF. D. PCA was performed using TIGR MeV software; principal components 1 and 2 are plotted in 2-dimensional space. Samples (points in space) have been colored according Figure 2. Normal-like are green, limited orange, diffuse-proliferation in red and inflammatory in black. Circles indicate groups of samples distinguished by the top two principal components. E. Principal components 1 and 3 were plotted in two-dimensional space and show distinction between two groups within the diffuse-proliferation, normal-like and limited scleroderma.

https://doi.org/10.1371/journal.pone.0002696.g003

Based on this analysis and the SigClust analysis, we propose that there are approximately four to five statistically significant clusters in the data. The statistically significant cluster assignments from both SigClust and consensus clustering are summarized in Table 3. These are (1) Diffuse-proliferation comprised completely of patients with dcSSc, (2) Inflammatory, which includes a subset of dSSc, lSSc and morphea, (3) Limited, characterized by the inclusion of lSSc patients and (4) Normal-like, which includes five of six healthy controls along with two dSSc and one lSSc patients. Notably, three samples are not consistently classified into the primary clusters. These are: dSSc2 which is assigned to the either the diffuse-proliferation, normal-like or into a single cluster by itself, dSSc13 which is assigned to either diffuse-proliferation or the limited groups, and the patient EF which clusters either on the peripheral edge of the diffuse-proliferation cluster or is assigned to a cluster by itself.

To determine how sensitive the clustering was to the selection of the intrinsic genes, we analyzed the clustering results using a larger list of 2071 intrinsic genes and compared that clustering to that obtained with 995 intrinsic genes (Supplemental Figure S3). Although we find slight differences in the ordering of the samples, the major subsets of diffuse-proliferation, inflammatory, and limited are again identified. The Normal-like group is split onto two different branches using this larger set of genes. Samples that show inconsistent clustering are from patient dSSc2, dSSc8, dSSc13, and the single array for patient EF. The samples for each of these patients were also inconsistently classified in the SigClust and consensus clustering analysis using the 995 intrinsic gene set.

Principal Component Analysis (PCA) was used to confirm the sample grouping found by hierarchical clustering. PCA is an analytic technique used to reduce high dimensional data into more easily interpretable principal components by determining the direction of maximum variation in the data [59]. The 995 intrinsic genes were analyzed by PCA using the MultiExperiment Viewer (MeV) software [60]. The first and second principal components that capture the most variability in the data (Figure 3D), and the first and third principle components (Figure 3E) have been plotted in 2-dimensional space. The 2D projection shows that the samples group in a manner similar to that found by hierarchical clustering analysis: normal controls and limited samples group together and the two different groups of diffuse scleroderma group together (Figure 3D). Notably, the first and second principal components (Figure 3D) separate the diffuse-proliferation, the inflammatory and the normal-like/Limited groups. When the first and third principal components are analyzed (Figure 3E) we find that the distinction between dSSc group 1 and dSSc group 2 is clearly delineated, as is the distinction between normal-like and limited. The PCA analysis provides further evidence, in addition to our hierarchical clustering analysis, that the gene expression groups are stable features of the data.

Biological processes differentially expressed in the intrinsic groups

In order to systematically investigate the biological processes found in the gene expression profiles of SSc, we created a module map using Genomica software [61], [62] (Figure 4A). A module map shows arrays that have co-expressed genes that map to specific gene sets. In this case, each gene set represents a specific biological process derived from Gene Ontology (GO) Biological process annotations [63], or from previously published microarray datasets [49], [64].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Scleroderma Module Map.

A. Module map of the Gene Ontology (GO) Biological Processes differentially expressed among the scleroderma samples is shown. Each column represents a single microarray and each row represents a single GO Biological process. Patient samples are organized as described in Figure 2. Only modules that were significantly enriched (minimum 2-fold change, p<0.05) on at least 4 micoarrays are shown. The average expression of the gene hits from each enriched gene set is displayed here. Only gene sets that show significant differences after multiple hypothesis testing were included. Select GO biological processes are shown. The entire figure with all biological processes can be viewed in Supplementary Figure S4. B. Module map using gene list created from an experimental identification of all cell cycle-regulated genes [49] and genes found to be expressed in specific lymphocyte subsets [64].

https://doi.org/10.1371/journal.pone.0002696.g004

Modules with significantly enriched genes (p<0.05, hypergeometric distribution) and corrected for multiple hypothesis testing with an FDR of 0.1% are shown (Figure 4A). Expressed among the group diffuse-proliferation are the biological processes of cytokinesis, cell cycle checkpoint, regulation of mitosis, cell cycle, DNA repair, S phase, and DNA replication, consistent with the presence of dividing cells. Decreased in this group are genes associated with fatty acid biosynthesis, lipid biosynthesis, oxidoreductase activity and decreased electron transport activity. The decrease in genes associated with fatty acid and lipid biosynthesis is notable given the loss of subcutaneous fat observed in dSSc patients [4].

Expressed in the inflammatory group are biological processes indicative of an increased immune response, including the GO biological processes of immune response, response to pathogen, humoral defense, lymphocyte proliferation, chemokine binding, chemokine receptor activity, and response to virus (Figure 4A). Biological processes of icosanoid and prostanoid metabolism, which represents synthesis of prostaglandin lipid second messengers, have been associated with immune responses [65], found to be highly expressed in rheumatoid arthritis [66], [67], [68] and associated with severity in collagen-induced arthritis in mice [69], [70]. Also expressed in this group are processes associated with fibrosis including trypsin activity, collagen and extracellular matrix. The full figure with all differentially expressed biological processes is available as Supplemental Figure S4.

In order to better define the proliferation signature observed, we created gene sets representing the genes periodically expressed in the human cell division cycle as defined by Whitfield et al. [49] (Figure 4B). Gene sets were created that included the genes with peak expression at each of the five different cell cycle phases, G1/S, S, G2, G2/M and M/G1 [49]. The enrichment of each of these five gene sets was statistically significant (p<0.05 using the hypergeometric distribution) and more highly expressed in the diffuse-proliferation group (Figure 4B).

To better characterize the lymphocyte infiltrates we generated gene sets representing lymphocyte subsets using results reported by Palmer and coworkers [64]. Using isolated populations of lymphocytes and DNA microarray hybridization, the genes specifically expressed in different lymphocyte subsets were identified. Subsets included T cells (total lymphocyte and CD8+), B cells, and granulocytes. We found four of these gene sets, B cells, T cells, CD8+ T cells and granulocytes, to have a statistically significant over-representation in the inflammatory group (Figure 4B). This suggests the gene expression signature expressed in this group is determined by the presence of infiltrating lymphocytes and specifically implies the infiltrating cells include T cells, B cells and granulocytes. Although we do not have a gene expression signature representative of macrophages or dendritic cells in this analysis, the macrophage marker CD163 is highly expressed in this group, suggesting innate immune responses may play an important role in disease pathogenesis.

Immunohistochemistry (IHC)

In order to verify that the gene expression reflected increased numbers of infiltrating lymphocytes or proliferating cells, we performed IHC for T cells (anti-CD3), B cells (anti-CD20) and cycling cells (anti-KI67). Summarized in table 4 is a full enumeration of marker positive cells counted from representative fields of all biopsies analyzed by IHC, with the observer blinded to disease state. All IHC images for all three markers are available as supplemental figure S5. Analysis of biopsies from each of the major intrinsic groups confirmed the results found in the gene expression signatures. The presence of infiltrating T cells was confirmed in the inflammatory group (Table 4). The largest numbers of T cells were found in perivascular and perifollicular distributions, as well as in the dermis, of two dSSc patients (dSSc5, dSSc6) assigned to the inflammatory group (Table 4). IHC was also performed on skin biopsies from two patients with morphea (Morph1, Morph 3) and each showed large numbers of infiltrating T cells. Only a small number of T cells were observed in two healthy controls analyzed (Nor2 and Nor3). A slight increase in T cells was observed in a perivascular distribution in the four patients assigned to diffuse-proliferation (dSSc1, dSSc2, dSSc11, dSSc12; Table 4), which had a lower expression of the T cell signature.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Immunohistochemical staining for KI67 and CD3 in the intrinsic subsets.

https://doi.org/10.1371/journal.pone.0002696.t004

Few CD20+ B cells were observed in the SSc skin biopsies. The immunoglobulin gene expression signature was observed in eight diffuse patients (dSSc1, dSSc3, dSSc6, dSSc7, dSSc8, dSSc10, dSSc11, dSSc12) and one limited patient (1SSc7; Figure 2F). Of the six patients analyzed by IHC (dSSc1, dSSc2, dSSc5, dSSc6, dSSc11, dSSc12), two samples (dSSc1and dSSc12) showed small numbers of CD20+ B cells.

The presence of the proliferation signature is correlated with an increase in the mitotic index or number of dividing cells in microarray studies of cancer [46], [47], [48], [49], [71]. To confirm the presence of proliferating cells in the dSSc skin biopsies, we performed IHC staining for KI67, a standard marker of cycling cells. Analysis of skin from healthy controls (Nor2, Nor3), morphea (Morph1, Morph3), and diffuse patients in the inflammatory group (dSSc5, dSSc6), showed no proliferating cells in the dermis, and a small numbers of proliferating cells surrounding dermal appendages and in the epidermal layer (Table 4). In contrast, analysis of the skin from four patients in the diffuse-proliferation subgroup (dSSc1, dSSc2, dSSc11 and dSSc12) showed higher numbers of proliferating cells primarily in the epidermis (Table 4). Therefore, we conclude that the proliferation signature is likely the result of an increased number of proliferating cells in the epidermal compartment of the SSc skin biopsies. The identity of these cells is very likely to be keratinocytes.

Intrinsic gene expression maps to identifiable clinical covariates

To map the intrinsic groups to specific clinical covariates, Pearson correlations were calculated between the gene expression of each of the 995 intrinsic genes and different clinical covariates. Shown are the results for three different covariates: the modified Rodnan skin score (MRSS; 0–51 scale), a self-reported Raynaud's severity score (0–10 scale), and the extent of skin involvement (dSSc, lSSc and unaffected). Each group was analyzed for correlation to each of the clinical parameters listed in Table 1; only the significant associations are shown. Figure 5A shows the gene expression patterns of the 995 intrinsic genes with each row representing a microarray and each column representing a gene. Pearson correlation coefficients were calculated between each of the clinical parameters and the expression of each gene. The moving average (10-gene window) of the resultant correlation coefficients is plotted for MRSS (Figure 5B), Raynaud's severity (Figure 5C) and degree of skin involvement (Figure 5D). Areas of high positive correlation between a clinical parameter and the expression of a group of genes indicate that increased expression of those genes is associated with an increase in that clinical covariate; a negative correlation indicates a relationship between a decrease in expression of the genes and an increase in a clinical covariate.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Correlation between gene expression and clinical covariates.

A. Shown is the color-coded heatmap of the 75 arrays and 995 intrinsic genes. The graph on the right of the heat map shows disease duration for each sample. Disease duration was set to zero for normal controls and morphea samples. B. Pearson correlations were calculated between skin score and the expression values for each gene in the list. The moving average of the Pearson correlation (10-gene window) was plotted. Regions of high negative and high positive correlations to the three different clinical parameters are indicated (regions I–III shaded grey). C. Moving average of the Pearson correlation coefficients (10-gene window) between the self-reported Raynaud's severity score and the expression of each gene, D. Moving average of the Pearson Correlations (10-gene window) between extent of skin involvement and a diagnosis vector (see Methods) for dSSc(red), lSSc (orange) and healthy controls (green). E. Box plot of disease duration for dSSc patients. The patients included in the diffuse-proliferation group had disease duration of 8.4±6.4 years. The dSSc patients that fell into the inflammatory or normal-like groups have disease duration of 3.2±3.9 yrs (p<0.12, t-test). F. Genes that ideally discriminate ‘Diffuse 1’ and ‘Diffuse 2’ groups were selected using Significance Analysis of Microarrays (SAM). 329 genes were selected with an FDR<1%. Pearson correlation coefficients were calculated between each clinical parameter and the expression for each gene and plotted as a 10-gene moving window.

https://doi.org/10.1371/journal.pone.0002696.g005

Areas of high positive or high negative correlation are highlighted in three different panels (Figure 5B–D, Regions I–III). Each of the three clinical covariates shows high positive correlations to a subset of gene expression signatures. Most notably, the MRSS skin score shows a high positive correlation to the ‘proliferation signature’ (Figure 5B, region II) with correlations ranging from 0.5 and 0.6. This signature is highly expressed in diffuse-proliferation samples but has low expression in the inflammatory group. The Raynaud's severity score has a high positive correlation to genes expressed at higher levels in the limited group and heterogeneously expressed in patients with dSSc (Figure 5C, region I). Not surprisingly, the genes highly correlated with MRSS also show a high positive correlation with diffuse skin involvement (Figure 5D, region II). While this signature associates with diffuse skin involvement, it is important to note that a subset of dSSc skin biopsies do not express this signature and have low skin scores. Similarly, the genes that have a high positive correlation with Raynaud's severity and a high positive correlation with the limited group (Figure 5C), which typically has more severe vascular involvement, are uncorrelated with the diagnosis of dSSc and are expressed at low levels in healthy control samples (Figure 5D, region III). Moving averages of the Pearson correlation between the intrinsic genes and other clinical covariates (digital ulcers, ILD, or GI involvement) were also calculated but did not reveal significant regions of positive or negative correlation to the gene expression profiles (data not shown).

One initial hypothesis was that there would be an obvious trend in the gene expression data reflecting the progressive nature of SSc in some patients. To examine this more carefully, disease duration in years since first onset of non-Raynaud's symptoms is plotted along the X-axis of the heat map (Figure 5A, right panel). The mean disease duration for the diffuse-proliferation group is 8.4±6.4 yrs, whereas mean disease duration for the inflammatory group, which includes dSSc and lSSc, is 6.5±6.1 yrs. Using a Student's t-test with a two-tailed distribution we find that this difference is not statistically significant. To test the hypothesis that a subset of the patients was grouping by disease duration, we analyzed the disease duration between the dSSc patients in the diffuse-proliferation group and the dSSc patients that were classified as either inflammatory or normal-like (Table 3). The diffuse-proliferation group has a mean disease duration of 8.4±6.4 years, and the dSSc patients in the inflammatory and normal-like groups have a mean disease duration of 3.2±3.9 yrs (Figure 5E, p = 0.12, t-test). The difference in the means between these two groups is clear, but outliers in each reduce the significance of the result. Dropping the two outliers results in p = 0.0042 (unequal variance two-sample t-test, two-sided)). Therefore, we conclude that there is a significant association between disease duration and the intrinsic groups for dSSc samples.

Since no obvious clinical covariate was identified that differentiated the dSSc group 1 from dSSc group 2, we selected the genes that most differentiated the two groups. Genes were selected that differentiated group 1 from group 2 using a non-parametric t-test implemented in Significance Analysis of Microarrays (SAM) [72]. 329 genes were selected that were differentially expressed between these two groups with an FDR of 0.19% (Figure 5F; Supplementary Data File S3). We analyzed these 329 genes for correlation to clinical covariates. Three clinical covariates were found associated with these two groups. The genes highly expressed in the dSSc group 2 (9 patients) are highly correlated with the presence of digital ulcers (DU) and the presence of interstitial lung disease (ILD) at the time the skin biopsies were taken. In contrast, dSSc group 1 (2 patients, both male) did not have DU or ILD at the time of biopsy. Although this grouping could result simply from stratification by sex, it also may reflect a true difference in disease presentation. Only 18 of the 329 genes map to either the X or Y chromosomes and thus are expected to be differentially expressed, suggesting the remainder may represent biology underlying these groups.

A subset of genes is associated with increased modified Rodnan skin score

To identify genes associated with MRSS we selected the subset of genes most highly correlated with each covariate from the intrinsic list using Pearson correlations. 177 genes were selected from the 995 intrinsic genes that had Pearson correlations with MRSS >0.5 or <−0.5. We then used this list of 177 genes to organize the skin biopsies by average linkage hierarchical clustering (Figure 6; Supplementary Data File S4). We find that both forearm and back skin biopsies from 14 patients with dSSc (mean MRSS of 26.34±9.42) clustered onto a single branch of the dendrogram. All other samples, including the forearm-back pairs of 4 patients with dSSc (mean MRSS 18.11±6.45) clustered onto a separate branch of the dendrogram. Using a two-tailed Student's t-test we find that the difference in skin score between the two groups of dSSc is statistically significant (p = 0.0197).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Genes correlated with MRSS.

We selected the genes from the 995 intrinsic list that had a correlation greater than 0.5 or less than −0.5 to the MRSS. This list of 177 genes was then used to organize the skin biopsies. Forearm-back pairs from 14 patients with dSSc (mean MRSS of 26.34±9.42) clustered onto one branch of the dendrogram tree. The forearm-back pairs of 4 patients with dSSc (Mean MRSS 18.11±6.45) clustered onto a different branch of the dendrogram tree. The difference in skin score between these two groups is statistically significant (p<0.0197).

https://doi.org/10.1371/journal.pone.0002696.g006

From this analysis, 62 genes were expressed high levels and 115 genes were expressed at low levels in the patients with the highest skin score. Genes highly expressed include the cell cycle genes CENPE, CDC7 and CDT1, the mitogen Fibroblasts Growth Factor 5 (FGF5), the immediate early gene Tumor Necrosis Factor Receptor Superfamily member 12A (TNFRSF12A) and TRAF interacting protein (TRIP). Since skin score is considered to be an effective measure for disease outcome, this 177-gene group may contain genes that could be further developed into surrogate markers for skin score.

Quantitative Real Time PCR

In order to validate the gene expression in the major groups found in this study, we performed quantitative real time PCR (qRT-PCR) on three genes selected from the intrinsic subsets (Figure 7). These included TNFRSF12A, which is highly expressed in the dSSc patients and shows high expression in patients with increased MRSS (see Figure 6), WIF1, which shows low expression in SSc its decreased expression is associated with increased MRSS, and CD8A, which is highly expressed in CD8+ T cells and is highly expressed in the inflammatory subset of patients. A representative sampling of patients from the intrinsic subsets was analyzed for expression of these three genes. Each was analyzed in triplicate and standardized to the expression of GAPDH.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Quantitative Real Time PCR analysis of representative biopsies.

The mRNA levels of three genes, TNFRSF12A (A), CD8A (B) and WIF1 (C) were analyzed by Taqman quantitative real time PCR. Each was analyzed in two representative forearm skin biopsies from each of the major subsets of proliferation, inflammatory, limited and normal controls. In the case of TNFRSF12A, patient dSSc11 was replaced by patient dSSc10, which cluster next to one another in the intrinsic subsets and show similar clinical characteristics (Table 1). Each qRT-PCR assay was performed in triplicate for each sample. The level of each gene was then normalized against triplicate measurements of GAPDH to control for total mRNA levels (see materials and methods). The relative expression values are displayed as the fold change for each gene relative to the median value of the eight samples analyzed.

https://doi.org/10.1371/journal.pone.0002696.g007

Each gene is shown with the fold change relative to the median value for the eight samples analyzed. TNFRSF12A shows highest expression in the patients with dSSc and the lowest in patients with limited SSc and normal controls. The three patients with highest expression are dSSc and include the proliferation group (Figure 7A). CD8A shows highest expression in the inflammatory subgroup as predicted by our gene expression subsets (Figure 7B). WIF1 shows highest expression in the healthy controls with approximately 4–8 fold relative decrease in patients with SSc (Figure 7C). The most dramatic decrease is in patients with dSSc with smaller fold changes in patients with lSSc.

cRNA synthesis, microarray hybridization and data processing

200 ng of total RNA from each biopsy was converted to Cy3-CTP (Perkin Elmer) labeled cRNA, and Universal Human Reference (UHR) RNA (Stratagene) was converted to Cy5-CTP (Perkin Elmer) labeled cRNA using a low input linear amplification kit (Agilent Technologies). Labeled cRNA targets were then purified using RNeasy columns (Qiagen). Cy3-labeled cRNA from each skin biopsy was competitively hybridized against Cy5-CTP labeled cRNA from Universal Human Reference (UHR) RNA pool, to 44,000 element DNA oligonucleotide microarrays (Agilent Technologies) representing more than 33,000 known and novel human genes in a common reference design [81]. Hybridizations were performed for 17 hours at 65°C with rotation.

After hybridization, arrays were washed following Agilent 60-mer oligo microarray processing protocols (6× SSC, 0.005% Triton X-102 for 10 min. at room temperature; 0.1× SSC, 0, 005% Triton X-102 for 5 min at 4°C, rinse in 0.1× SSC). Microarray hybridizations were performed for each RNA sample resulting in 61 hybridizations. Fourteen replicate hybridizations were added, resulting in a total of 75 microarray hybridizations.

Microarrays were scanned using a dual laser GenePix 4000B scanner (Axon Instruments). The pixel intensities of the acquired images were then quantified using GenePix Pro 5.0 software. Arrays were visually inspected for defects or technical artifacts, and poor quality spots were manually flagged and excluded from further analysis. Only spots with fluorescent signal at least twofold greater than local background in both Cy3- and Cy5- channels were included in the analysis. Probes missing more than 20% of their data points were excluded, resulting in 28,495 probes that passed the filtering criteria. The data were displayed as log₂ of the LOWESS-normalized Cy5/Cy3 ratio. Since a common reference experimental design was used, each probe was centered on its median value across all arrays.

Selection of intrinsic genes

An intrinsic gene identifier algorithm was used to select a set of intrinsic scleroderma genes. Detailed methods on the selection of intrinsic genes are available in [47] and http://genome-www.stanford.edu/breast_cancer/molecularportraits/. A gene was considered ‘intrinsic’ if it showed the most consistent expression between forearm-back pairs and technical replicates for the same patient, but had the highest variance in expression across all samples analyzed. The intrinsic gene identifier computes a weight for each gene, which is inversely related to how intrinsic the gene's expression is across the samples analyzed. A lower weight equals a higher ‘intrinsic’ character. A total of 34 experimental groups were defined, each representing the 34 different subjects in our study. Replicate hybridizations for a given patient were assigned to the same experimental group.

In order to estimate False Discovery Rate (FDR) at a given intrinsic weight, the analysis was repeated on data randomized in rows (i.e. across each gene). The FDR at a given weight was estimated by determining the number of genes that received the same weight or lower in the randomized data. 995 genes were selected that had an intrinsic weight <0.3; in randomized data 39±7 genes (calculated from 10 independent randomizations) had a weight of 0.3 or less, resulting in an FDR of approximately 4%. We found that a cutoff of 0.3 balanced the number of genes selected with an acceptable FDR, while retaining reproducible hierarchical clustering of technical replicate samples. Although it is possible to select a more or less restrictive list of genes with FDRs of 5% (weight <0.35; 2071 genes), 3.4% (weight <0.25; 425 genes) or 2.4% (weight <0.20; 171 genes), these smaller lists of genes resulted in less reproducible hierarchical clustering suggesting overfitting (Supplemental Figure S3).

Hierarchical clustering

Average linkage hierarchical clustering was performed in both the gene and experiment dimensions using either Cluster 3.0 software (http://bonsai.ims.u-tokyo.ac.jp/mdehoon/software/cluster/software.htm) or X-Cluster (Gavin Sherlock; http://genetics.stanford.edu/sherlock/cluster.html) using Pearson correlation (uncentered) as a distance metric [18]. Clustered trees and gene expression heat maps were viewed using Java TreeView Software (http://jtreeview.sourceforge.net/) [82].

Robustness and statistical significance of clustering

The statistical significance of clustering was assessed using Statistical Significance of Clustering (SigClust) [57] and Consensus Cluster [58]. SigClust tests the null hypothesis that the samples form a single cluster. A statistically significant p-value indicates the data came from a non-Gaussian distribution and that there is more than one cluster. Two different p-values were used to identify significant clusters, p<0.01 and p<0.001. The statistical significance of the clusters was first assessed at the root node of the tree derived from hierarchical clustering with the 995 intrinsic genes. If the cluster was statistically significant, the next node further down the tree was tested. The process ended when a cluster had a p-value greater than the established cutoff.

In addition, we analyzed the 995 intrinsic genes using Consensus Cluster [58]. Consensus Cluster is available through GenePattern (v.1.3.1.114; [83]). Assessment of sample clustering was performed by consensus clustering with K clusters (K = 2, 3, 4…10) using 1000 iterations with random restart. Samples that clustered together most often in each of the K clusters received a correlation value. The resulting consensus matrix was visualized as a color-coded heat map with varying shades of red, the brighter of which corresponded to higher correlation among samples (Figure 3A). Summary statistics are shown, including the empirical consensus distribution function (CDF) vs. the consensus index value (Figure 3B). Also shown are the proportion change (ΔK) under the CDF for each K = 2, 3…10 (Figure 3C). Consensus Cluster assignments for each sample are summarized in Table 3.

Principal Component Analysis

Principal Component Analysis was performed using Multiexperiment Viewer (MeV) software version 4.0.01 (http://www.tm4.org/mev.html; [60]). Data was loaded into MeV as a tab delimited text file of log2-transformed Cy3/Cy5 ratios. For PCA analysis [59], missing data were first estimated using K-nearest neighbors (KNN) imputation with N = 4.

Module Maps

Module maps were created using the Genomica software package [61], [62]. Gene sets containing all human Gene Ontology (GO) Terms were obtained from the Genomica web site (http://genomica.weizmann.ac.il/; Human_go_process.gxa, created Nov. 20, 2006). Additional custom gene sets representing the human cell division cycle [49] and lymphocyte subsets [64] were created specifically for this study. The human cell division cycle gene set was created from the genes found to periodically expressed in human HeLa cells [49]. Genes found to show peak expression at the five different cell cycle phases G1/S, S, G2, G2/M and M/G1 were each put into their own independent gene list. Gene sets representing different lymphocyte populations, T cells (total population, CD4+, CD8+), B cells, and granulocytes, were derived for this study from the genes expressed in isolated lymphocyte subsets by Palmer and coworkers [64].

All 75 microarray experiments and 28,495 DNA probes were included in the module map analysis. The 28,495 probes were collapsed to 14,448 unique LocusLink Ids (LLIDs) [84]. Only gene sets with at least three genes but fewer than 1000 genes were analyzed. A gene set was considered enriched on a given array if at least 3 genes from that set were considered to be significantly up-regulated or down-regulated (minimum 2-fold change, p<0.05, hypergeometric distribution) on at least four microarrays. Each gene set was corrected for multiple hypothesis testing using an FDR correction of 0.1%.

Correlation to clinical parameters

Pearson correlations were calculated between each clinical parameter and the gene expression data in Microsoft Excel. Pearson correlations between the diagnosis of dSSc, lSSc and healthy controls and the gene expression data were calculated by creating a ‘diagnosis vector’. The diagnosis vector was created by assigning a value 1.0 to all dSSc samples and 0.0 to all remaining samples for the dSSc vector; lSSc and healthy controls were treated similarly creating a vector for each. Pearson correlations were calculated between the gene expression vector and the diagnosis vector for dSSc, lSSc and healthy controls. Correlations between the gene expression and clinical data were plotted as a moving average of a 10-gene window.

Immunohistochemistry (IHC)

IHC was performed on paraffin embedded sections at the University of California, San Francisco in the Immunohistochemistry and Molecular Pathology core facility. All immunostaining was completed via a semi-automated protocol utilizing an automated immunostainer (DAKO Corp, Carpenteria, CA). Slides were heated, deparaffinized and then hydrated. Protease digestion was completed followed by antigen retrieval via pressure cooker as per standard protocols. After an endogenous peroxidase block with 3% H₂0₂, slides were loaded on to the automated immunostainer. A primary antibody cycle of 30 min was followed by a secondary antibody cycle using the ENVISION+ system. Color development was completed using DAB followed by counterstaining with Gills #2 Hematoxylin. Specific conditions for the antibodies utilized were as follows: anti-CD20 (DAKO) was used at 1∶600 for 30 minutes in citrate buffer (pH 6.0); anti-CD3 (DAKO) at 1∶400 for 30 minutes in Tris buffer (pH 9.0), and anti-Ki67 (MiB1; DAKO) was used at 1∶1000 for 30 minutes in Tris buffer (pH 9.0). Marker positive cells were enumerated by tissue compartment in equal sized images of n skin biopsies, with the observer blinded to disease state and array results of the specimens (Table 4).

Data Access

The full dataset, figures in both red/green and blue/yellow format, as well as searchable versions of figures 1 and 2 are available at website maintained by the authors: http://whitfieldlab.dartmouth.edu/SScSubsets/. Raw data is available from the UNC Microarray Database (https://genome.unc.edu/) and has been deposited to NCBI's Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/; accession Number: GSE9285).