MassARRAY genotyping of mutations in HA, HMIN, and HWIDE tumors

Comparative genomic hybridization analysis

We used the Agilent 1M array comparative genomic hybridization (CGH) platform for global copy number analysis (www.agilent.com/genomics) according to the manufacturer's instructions. The data were normalized, and then the Cy5 to Cy3 ratio for each probe was expressed as log₂(Cy5/Cy3). Normalized data were then analyzed using RAE, a computational segmentation framework that adapts parameters to individual tumors, and then identifies statistically significant regions of interest across tumors (19). Candidate copy number alterations are based on a false discovery rate (FDR) Q-value. Regions with FDR ≤10% were plotted. Analysis and scores were calculated as described (19). These regions were then compared with known copy number alterations in PTC and FTC from published array CGH data. Clustering analysis of the copy number data was done using GenePattern.

Gene expression profiling

All gene expression analyses were performed with the HG-U133A Plus version 2.0 (Affymetrix, Santa Clara, California) array as described by the manufacturer. The microarray data were quantile-normalized (20), and the gene expression values were estimated using the robust multi-array average method (21). Moderated t statistics (22) were used to test whether genes were differentially expressed between the groups of interest. P value < .001 was considered statistically significant. Genes with P values < .001 were used for hierarchical clustering analysis. The expression values for each gene were centered and scaled before analysis. Complete-link algorithm was used for clustering, and Pearson correlation coefficient was used as the similarity distance.

Comparisons of expression profiles for HA vs HMIN, HA vs HWIDE, and HMIN vs HWIDE were performed (Supplemental Tables 3–5) and a hierarchical clustering analysis for genes with P values < .001 and a fold change differential of ×1.8 performed. For the HA vs HWIDE comparison, select genes of significance were internally validated using real-time PCR with the iCycler System (Bio-Rad Instruments, Hercules, California) using SYBR green detection. Primers for each gene were designed using Primer-Blast (www.ncbi.nlm.nih.gov/BLAST/) and are shown in Supplemental Table 6 with optimum annealing temperature. Quantitative RT-PCR analysis was performed using the ΔΔCt method (23). As a reference gene, the housekeeping gene GADPH was amplified using the primers forward 5′-AAGGTGAAGGTCGGAGTCAA-3′ and reverse 5′-AATGAAGGGGTCATTGATGG-3′ and an annealing temperature of 60°C.

Concept module mapping with Oncomine, Gene Set Enrichment Analysis (GSEA), and Ingenuity Pathway Analysis

Concept module mapping was performed as follows. The differential gene expression signature identified from our analysis for HA vs HWIDE, HA vs HMIN, and HMIN vs HWIDE (P value < .001 and differential fold change of ×1.8) was imported into Oncomine (http://www.oncomine.org) to search for associations with molecular concept signatures derived from independent cancer profiling studies. We report statistically significant overlaps of our gene expression signature with the top-ranking gene expression signatures of clinical outcome using percentile cutoffs (10%). Q-value is calculated as previously described (24). GSEA was performed with GSEA software version 2.0.7 (www.broadinstitute.org/cancer/software/gsea). We assessed the significance of the gene sets with the following parameters: number of permutations = 1000, and permutation type = phenotype with an FDR Q-value cutoff of 25%. The most differentially expressed genes from statistically significant gene sets were identified with the leading edge subset that consists of genes with the most contribution to the enrichment score of a particular gene set. For Ingenuity Pathway Analysis, the same differential gene expression signature for HA vs HWIDE was imported into IPA (http://ingenuity.com/) and then correlated with genes associated with vascular invasion to search for interacting genes.

Comparison of transcriptomes of HWIDE, PTC, and FTC tumors

Using GEO (gene expression omnibus: http://www.ncbi.nlm.nih.gov/geo/), expression array profiles for PTC (25) and FTC (26) were obtained and then compared with the transcriptomes for the HWIDE profile. Datasets were identified using the same Affymetrix platform (HG-U133A and 2.0) and are shown in Supplemental Table 7.

Using published gene expression profiles for normal thyroid tissue, the differential gene expression profiles for PTC, FTC, and HWIDE were compared with normal thyroid. To reduce any potential batch effect, the expression data were first normalized, and then the differential gene expressions identified with Significance Analysis of Microarray (SAM) (27) using MultiExperiment Viewer (MeV) with 100 permutations. The complete lists of genes upregulated and downregulated for PTC, FTC, and HWIDE compared with normal are shown in Supplemental Tables 8, 9, and 10. Differentially expressed genes were used to identify the key pathways.

Global copy number diversity in Hurthle tumors

Regions of chromosomal gain and loss are shown in Supplemental Table 11 (a more detailed list is shown in Supplemental Table 12 and Supplemental Gene Tables 1–6). This table also lists genes and microRNAs for each significant chromosomal region. Genomic regions showing amplification were chromosomes 4p16, 5p15-5q35,6p25, 7p15-7q36, 8p21-23, 10p13, 12p13-q24, 16q23, and 17p13-q25. Chromosomal regions of loss included 4q24, 6p23, 7p15, 9q33, 13q14, and 16q23. In particular, there were very large regions of amplification on chromosome 5 (174 gene region and 650 gene region), chromosome 7 (128 gene region and 776 gene region), chromosome 12 (995 genes), and chromosome 17 (471 gene region, 322 gene region, 59 gene region, 252 gene region). The RAE plot showing amplified and deleted genes of interest is shown in Supplemental Figure 1. A summary of chromosomal regions of gain and loss in HCC compared with PTC and FTC is shown in Figure 2. Regions of gain are shown in red, and regions of loss are shown in green.

Open in a separate window

Figure 2.

Summary of chromosomal regions of gain and loss in HCC compared with PTC and FTC. Regions of gain are shown in red, and regions of loss are shown in green.

A number of chromosomes demonstrated increases in number, including chromosomes 5, 7, 10, 12, 17, and 20, suggesting the presence of either chromosomal instability or selection for increased copies of these chromosomes. Clustering analysis of the copy number data revealed 3 main groups of Hurthle tumors (Figure 3). These were tumors where large regions of amplification were predominant (subtype 3), tumors where deletions were predominant (subtype 1), and tumors with a mix of focal chromosomal alterations (subtype 2). Four widely invasive tumors were subtype 3, 5 were subtype 2, and 1 was subtype 1. The 3 tumors that were clinically more aggressive with distant metastases were of subtype 2. We analyzed all 22 chromosomes to determine whether individual chromosomal gains/losses could differentiate HA, HMIN, and HWIDE. The differential pattern for all chromosomes closely resembled the clustering pattern observed for all 22 autosomes. We were unable to identify any pattern of chromosomal gains/losses that could specifically differentiate HA from HMIN and from HWIDE.

Open in a separate window

Figure 3.

Diversity of the Hurthle cell tumor copy number landscape. Clustering of copy number alterations for each tumor. Amplifications (red) and deletions (blue) are indicated across the 22 autosomes. Three main groups of Hurthle tumors are identified: tumors where large regions of amplification are predominant (subtype 3 shown on the right), tumors where deletions are predominant (subtype 1 shown on the left), and tumors with more focal chromosomal alterations (subtype 2).

Transcriptomes of Hurthle cell tumors reveal differentially active pathways

Principal component analysis using the Partek Genomics Suite (http://www.partek.com/software) was carried out to determine the similarity of HA, HMIN, and HWIDE tumors (Figure 4A). The HMIN tumors as a group were very similar to the HA tumors. In contrast, the HWIDE tumors clustered away from the HA and HMIN groups, although there were 3 tumors similar to the HMIN group. This indicated that the HWIDE tumors were quite different in their transcriptional profiles from each other and from both HA and HMIN.

Open in a separate window

Figure 4.

A, Principal component analysis of expression data derived from Hurthle cell tumors. The subclasses HA, HMIN, and HWIDE are noted in the color legend. B, Heatmap showing differential gene expression between HA, HMIN (Min Inv Ca), and HWIDE (Widely Inv Ca). HMIN shows a similar profile to HA, whereas HWIDE shows a more divergent expression profile. Histogram from clustering analysis is shown at the top. Relative expression is noted in the legend.

Unsupervised hierarchical clustering (Figure 4B) also showed that the 3 groups of tumors clustered separately from each other with some differences between the HA and HMIN tumors but a more marked difference between HWIDE and the other 2 groups. Supplemental Table 13 shows the number of differentially expressed genes between HA, HMIN, and HWIDE tumors. The complete list of differentially expressed genes with a P value < .001 for HA vs HMIN, HA vs HWIDE, and HMIN vs HWIDE are shown in Supplemental Tables 3–5. Because the greatest number of gene expression differences were in the HA vs HWIDE comparison (henceforth called the HWIDE signature), we focused further analysis on this dataset. This comparison showed 403 gene expression differences; 324 genes were underexpressed in the HWIDE tumors and 79 genes overexpressed compared with the HA tumors. Select genes were validated by RT-PCR. Genes selected for validation were PFKFB2, CTNNB1, PRO1073, SULT131, MALAT1, ITCH, RHOU, TXN, ZNF567, SHB, AYTL1, and GATA6. Primers and PCR conditions used for RT-PCR validation are shown in Supplemental Table 6 and Supplemental Figures 2 and 3.

Genetic programs involving vascular invasion and β-catenin in Hurthle cell tumor subsets

Concept module mapping of the HWIDE signature using Oncomine showed that there were several strongly significant gene sets that were enriched (Supplemental Table 14). The top program was a CTNNB1 (β-catenin)-driven signature resulting from expression of β-catenin (31) (Q-value = 1.5 × 10⁻²³) (Figure 5 and Supplemental Table 15). This implicates β-catenin in HCC oncogenesis.

Open in a separate window

Figure 5.

Correlation of HWIDE signature with select gene sets in Oncomine database. The significant gene overlap is represented by −log of the P value. CTNNB1 transfection, outlined in red, has the most significant P value with a −log P value of 27.7. The complete concept datasets are shown in Supplemental Table 13.

To determine the molecular context of β-catenin pathway enrichment, we made use of Ingenuity Pathway Analysis. Analysis of genes involved in vascular invasion with the HWIDE signature showed that β-catenin is intimately involved in processes regulating vascular invasion. Supplemental Figure 4, A and B, shows the results from this analysis. Another concept highlighted from Oncomine analysis was the temsirolimus-sensitive signature (31) (genes relating to temsirolimus sensitivity), which had 60 overlapping genes with a significant Q-value of 1.2 × 10⁻⁸. Temsirolimus is a mammalian target of rapamycin (mTOR) inhibitor and recent phase I studies at our institute have shown patients with distant metastases from widely invasive HCC responded to this drug (unpublished data). The gene set was also analyzed using GSEA and the top datasets with significant overlapping genes are shown in Supplemental Table 16. There is also significant overlap with genes upregulated in HepaRG cells (liver cancer) expressing active forms of mTOR (23-gene overlap; P = 5.08 × 10⁻⁶). This again was evidence that mTOR may be important in HCC pathogenesis. Significant overlap with genes downregulated in poorly differentiated thyroid cancer (31-gene overlap; P = 1.6 × 10⁻⁷) and genes downregulated in anaplastic thyroid cancer (18-gene overlap; P = 3.4 × 10⁻⁴) were also found. Concept analysis of the HMIN gene signature (HMIN vs HA) using Oncomine (Supplemental Table 17a) and GSEA (Supplemental Table 17b) showed different concepts compared with the HWIDE signature. This provides some evidence that the progression of HCC from HA to HMIN and HWIDE may not be a stepwise progression model. In contrast, concept analysis of HWIDE vs HMIN showed similar concepts to the HWIDE vs HA differential expression set (Supplemental Table 18, a and b). Additional concepts of interest identified included sorafenib, torcetrapib, pazopanib, and foratinib sensitivity signatures.

We would like to thank Agnes Viale and the Memorial Sloan Kettering Genomics Core for excellent technical support.

This work was supported by Dr Peter Scardino (Department of Surgery, MSKCC) and by Mr Richard Pechter of the Pechter Foundation. Special thanks goes to Mr Pechter and Dr Scardino for providing encouragement and inspiration for this work. Funding was also provided by the Louis Gerstner Foundation (T.A.C.).

Disclosure Summary: All authors have no conflicts of interest.