Port, M., Boltze, C., Wang, Y., Röper, B., Meineke, V. and Abend, M. A Radiation-Induced Gene Signature Distinguishes Post-Chernobyl from Sporadic Papillary Thyroid Cancers. Radiat. Res. 168, 639–649 (2007).
We investigated selected gene targets to differentiate radiation-induced papillary thyroid cancers (PTCs) from other etiologies. Total RNA was isolated from 11 post-Chernobyl PTCs and 41 sporadic PTCs characterized by a more aggressive tumor type and lacking a radiation exposure history. RNA from 10 tumor samples from both groups was pooled and hybridized separately on a whole genome microarray for screening. Then 92 selected gene targets were examined quantitatively on each tumor sample using an RTQ-PCR-based low-density array (LDA). Screening for more than fivefold differences in gene expression between the groups by microarray detected 646 up-regulated and 677 down-regulated genes. Categorization of these genes revealed a significant (P < 0.0006) over-representation of the number of up-regulated genes coding for oxidoreductases, G-proteins and growth factors, while the number of genes coding for immunoglobulin appeared to be significantly down-regulated. With the LDA, seven genes (SFRP1, MMP1, ESM1, KRTAP2-1, COL13A1, BAALC and PAGE1) made a complete differentiation between the groups possible. Gene expression patterns known to be associated with a more aggressive tumor type in older patients appeared to be more pronounced in post-Chernobyl PTC, thus underlining the known aggressiveness of radiation-induced PTC. Seven genes were found that completely distinguished post-Chernobyl (PTC) from sporadic PTC.
INTRODUCTION
Thyroid cancer is the third most common solid tumor in children and adolescents (1). Radiation exposure at a young age is known to cause thyroid cancer, predominantly of the papillary type (2, 3). The risk of thyroid tumors is increased linearly for mean doses as low as 100 mGy (4). As a result of the Chernobyl nuclear reactor accident (April 1986), an enormous increase [approximately 30-fold, especially in those <1 year old at exposure; for details see ref. (5)] in the incidence of childhood thyroid carcinoma has been observed in Belarus and Ukraine and to a lesser extent in the Russian Federation starting in 1990 (6). A comparison of these tumors with carcinomas in patients in the same age group in Italy and France revealed that the post-Chernobyl thyroid carcinomas were less influenced by gender and were nearly always papillary (7). A solid papillary histological type was found in 93% of children in the Ukraine (8), and an increased prevalence of ultrasonographic thyroid abnormalities (nodules, cysts) was found (9, 10). The tumors appeared to be more aggressive at presentation with lymph node metastasis in 60–80% of patients; other extrathyroidal spread and venous invasion were commonly observed (11, 12). The tumors were shown to be more frequently associated with thyroid autoimmunity and developed after a short latent period of 4–5 years (13). RET activation/rearrangements were found in tumors nearly 70% of the patients (6).
The exact impact of the Chernobyl accident on the development of thyroid cancer in the Ukraine and Europe has been difficult to determine. For example, between 1980 and 2000 cancer of the thyroid increased significantly in France (14). Epidemiological evidence did not support a link with the Chernobyl accident, but the Thyroid Cancer Committee recommended a national registry dedicated to thyroid cancer in young people (14). In Belgium, four children aged 2 months to 10 years at the time of the Chernobyl accident developed thyroid cancer shortly thereafter. This represented an unusual increase in thyroid cancer incidence in adolescents. The question arose as to whether the occurrence of thyroid cancer could be related to the radioactive clouds that had passed over Belgium (3). Similarly, in England, population-based data on thyroid carcinomas were obtained from the Northern Region Young Person's Malignant Disease Registry to address the question as to whether the increasing incidence might be related to the Chernobyl accident (15). Even in Belarus there was a debate about whether the reported increase in thyroid cancer was real and attributable to radiation or was an artifact due to incorrect histological diagnosis, more complete case reporting, and mass screening of children after the accident (16).
Due to these uncertainties, scientists have attempted to identify biological markers for radiation-induced thyroid cancer. Our study focused on a comparison of the same histology that had developed due to a different etiology. For this reason, RNA from tumor tissue of Belarusian patients exposed to radiation after the Chernobyl nuclear accident and RNA from a control group were hybridized on a whole genome microarray for screening of potentially up-regulated or down-regulated genes. Microarray results are semiquantitative and must be validated with another method such as real-time quantitative polymerase chain reaction (RTQ-PCR). Nearly 100 of the most promising genes were examined quantitatively for each of the tumor samples with a recently developed RTQ-PCR-based technology called low-density array (LDA) to identify the most promising genes for distinguishing between sporadic and radiation-induced PTC.
MATERIALS AND METHODS
Tumor Samples
We examined thyroid cancers from 11 patients diagnosed up to 15 years after external and internal radiation exposure due to the Chernobyl nuclear accident, with a median age of 18 years at exposure (range 15– 22 years; eight females, three males; Table 1). Patients were located within 10–50 km from the reactor during the accident. Forty-eight hours after the accident, the population living within about 20 km was evacuated. At that time the absorbed intrathyroid dose was estimated to be 0.15–1.0 Gy for 90%, 1–5 Gy for 9%, and >10 Gy for 1% of the evacuated children (17). After total thyroidectomy due to advanced carcinoma, samples of the thyroid tumor tissue were taken and stored at −80°C. Cryopreserved tumor tissues from 41 Caucasian patients (Table 1) in southeastern Germany diagnosed with thyroid carcinoma of similar histology but lacking a radiation exposure history (median age 60 years, range 15– 83 years; 23 females, 18 males) served as controls. Finally, four individuals of the control group with a median age of 24 years were selected as a closer match to the median age of 18 years of the Chernobyl group, and normalized gene expression was examined separately (Fig. 5B). The tumors were analyzed histologically by two pathologists according to the WHO classification. Only papillary thyroid histologies and no follicular subtypes were included in this study. Tumor stage was based on the TNM tumor classification system (18). The studies were approved by the local ethics committees of Otto-von-Guericke-University, Magdeburg, and all patients gave their written consent.
RNA Isolation
Cryopreserved tissues were carefully thawed in RNA-later solution (Qiagen, Hilden, Germany), homogenized (Homogenizer, Omni, Warrenton, VA), and digested (proteinase K, 20 mg/ml; Invitrogen, Karlsruhe, Germany). Total RNA was isolated (RNeasy Mini Kit, Qiagen) and the remaining DNA was digested (RNase free DNase Set, Qiagen) according to the manufacturer's instructions (Applied Biosystems, Weiterstadt, Germany). RNA was stored at −80°C until it was used for gene expression arrays and RTQ-PCR analysis.
To control the quality and purity of isolated total RNA, spectrophotometry, agarose gel electrophoresis and PCR (using β-actin primers to detect DNA contamination) were performed. Only total RNA with a ratio A260/A280 >2.0 (spectrophotometry) with 28S ribosomal bands that were present at approximately twice the amounts of the 18S RNA (agarose gel electrophoresis) and without detectable contamination of DNA (PCR) were used for gene expression arrays and RTQ-PCR.
Whole Genome Microarray Analysis (Human Genome Survey Microarray V2.0)
RNA from 10 individuals (3–4 μg RNA per tumor sample) from each group was pooled separately. This resulted in 37–40 μg RNA templates per microarray. After further RNA quality control (see above), templates were converted into cDNA using a chemiluminescence RT labeling kit. Briefly, mRNA was reverse transcribed using oligo (dT) primers and labeled with DIG-dUTP. After RNA degradation, the DIG-labeled cDNA was purified and hybridized on the Human Genome Survey Microarray V2.0. The following day, microarrays were washed and the chemiluminescence reaction was performed using the chemiluminescence kit. Chemiluminescence was detected with the AB1700 microarray reader. Biochemical reactions and data analysis were performed under the observation of experienced personnel from the Gene Expression Center in Martinsried, Germany. All materials and instruments used for microarray analyses were purchased from Applied Biosystems. Quantile normalized gene expression of 32,878 probes of the microarray for the interrogation of 29,098 genes measured in both groups were log2-transformed and plotted as a scatter plot after deletion of signals with a signal/noise ratio <2.5 (fold change of the probe's signal relative to the SD of measurements within the spot). With a signal/noise ratio ≥2.5 there is a 99.4% probability that the measurements performed are correct.
For diagnostic purposes it was assumed that genes showing a large difference in expression between the groups would be most promising. For this reason, genes with either >five- or <0.2-fold changes in expression were examined more carefully. These genes were categorized according to the biological processes and molecular functions they code for using the PANTHER database (Applied Biosystems). Differentially expressed genes detected with the microarray were collected into groups that are representative for certain functions (e.g. oxidoreductases). The significance was estimated with a mathematical tool (available at www.pantherdb.org/tools) that compared the observed number of differentially expressed genes in a certain gene category with the expected number of genes. The expected number of genes represents a fraction of the number of genes on the whole genome array that are coding for this function. This fraction corresponds to the total number of differentially expressed genes observed. Then an over-representation or under-representation in the number of differentially expressed genes of a certain category can calculated along with the associated P value.
RTQ-PCR with Low-Density Array (LDA)
The LDA allows quantification of a set of customer-specified transcripts. Aliquots of total RNA (5 μg) of each individual were reverse transcribed over 2 h using a two-step PCR protocol (High Capacity Kit). A volume of 50 μl of the resulting cDNA typically contained a 1-μg RNA equivalent. To this volume an aliquot of 50 μl 2× RT-PCR master mix was added, mixed and pipetted into a fill port of a low-density array and transferred into the 7900 RTQ-PCR instrument equipped with a specially designed thermal cycler to run the RTQ-PCR protocol with the 384-well low-density array format over 2 h.
Statistics
Gene expression microarray analyses with pooled RNA probes were performed once for screening purposes and validated using RTQ-PCR (LDA). For this purpose, 92 genes that appeared to be >fivefold up-regulated or down-regulated were selected and quantified in each tumor sample. Four gene selection algorithms were used to select the top discriminative genes from the LDA data: (1) the χ2 statistic, which determines the worth of a gene by computing the value of the χ2 statistic with respect to the class; (2) information gain, which determines the worth of a gene by measuring the information gain with respect to the class; information gain is biased in favor of genes with higher dispersion; (3) symmetrical uncertainty, which measures the worth of a gene by measuring the symmetrical uncertainty with respect to the class and compensates for information gain bias; and (4) reliefF, which is a gene-weighting algorithm that is sensitive to gene interactions. The key idea of reliefF is to rate genes according to how well their values distinguish among instances of different classes that are near each other [for details, see ref. (19)]. The mean values, geometric mean and SEM were calculated using Sigma Plot 2000 (Jandel, Erkrath, Germany).
RESULTS
Whole Genome Microarray and Validation with LDA
The Pearson coefficient of correlation in the array data was 0.921 (Fig. 1). The percentages of expressed probes/ genes (present calls) were 48.8% in the Chernobyl group and 41.4% in the controls. The number of differentially expressed probes (ratio of normalized gene expression measured in irradiated and control genes) with a ratio >2 (up-regulated) or <0.5 (down-regulated) was 3,693 and 2,566, respectively. The number of differentially expressed genes >five- or <0.2-over control was 646 and 677, respectively.
Geometric means of differential gene expression measured with the microarray (semiquantitative method) were compared with the geometric means of the same gene measured with the LDA (quantitative method). About 82% of the genes (72/88) were detected as up-regulated or down-regulated genes with both methods. However, the percentage of false positive genes (differentially expressed on the microarray while showing control values according to the LDA) was 13% (11/88). The percentage of false negative genes (control values on the microarray while being differentially expressed according to the LDA) was 3% (3/88). Contradictory measurements (same gene appears to be up-regulated with one method and down-regulated with the other method) were found in 2% of the genes examined (2/ 88).
The number of genes coding for signal transduction (P = 2 × 10−5) and its subcategory intracellular signaling cascade (P = 9 × 10−5) appeared to be over-represented in post-Chernobyl PTC (Fig. 2). Within the biological process of signal transduction, about 20% of the genes coded for G-proteins and 10% for growth factors like VEGF-A or PDGF-B or those in the EGF signaling pathway. The number of genes coding for other biological processes like fatty acid desaturation (P = 2 × 10−4) and steroid metabolism (P = 6 × 10−4) also appeared to be over-represented. All genes in these categories coded for oxidoreductases such as cytochrome b5 reductase 1 (CYB5R1) and glutathione peroxidase 2 (GPX2). The accompanying molecular functions were significantly over-represented in the categories of signaling molecule (P = 1 × 10−5), growth factor (P = 8 × 10−5) and oxidoreductase (P = 4 × 10−4). Again, the genes coding for growth factors, G-proteins and oxidoreductases predominated.
Genes with <0.2-fold change in expression were significantly over-represented in biological processes such as B-cell and antibody-mediated immunity (P = 2 × 10−74), immunity and defense (P = 1 × 10−53), signal transduction (P = 8 × 10−15), and cell communication (P = 5 × 10−14), as shown in Fig. 3. The accompanying molecular functions were significantly over-represented in the categories immunoglobulin (P = 4 × 10−83) and defense/immunity protein (P = 4 × 10−67). These changes were caused primarily by down-regulated genes coding for immunoglobulin (Fig. 3). Significantly increased numbers of down-regulated genes (although less pronounced) coding for G-proteins, cell adhesion and cytokines/chemokines were also detected for biological processes such as signal transduction and cell communication (Fig. 3).
LDA
Normalized gene expression measured with the LDA for 92 genes in every tumor sample showed an individual variation in gene expression in the range of a factor of 10 within each group examined. When comparing the geometric means of the two groups, it was shown that for certain genes (e.g. PAEP or PAGE1), an up to 1000-fold difference in the normalized gene expression between the groups was detected (Fig. 4). A ranking of the 92 genes was performed with four different mathematical algorithms (Table 2). Based on these results, 15 genes were selected that appeared most promising for differentiating between the two groups (Fig. 5A). Seven genes, COL13A1, KRTAP2-1, MMP1, SFRP1, ESM1, BAALC and PAGE1, made it possible on their own to completely differentiate between the groups without overlapping measurements (Fig. 5A). Discrimination of the groups with the seven genes was improved when examining the age-matched groups (Fig. 5B), and three additional genes (GNGT2, MST150 and RIPK4) allowed a complete differentiation between the groups (Fig. 5B).
DISCUSSION
With the microarray used, several hundred genes that were >fivefold up-regulated or down-regulated were selected. A total of 92 of these genes were chosen for diagnostic purposes and were confirmed by a recently developed high-throughput RTQ-PCR-based method called LDA. At the same time their measurements served as a control for the accuracy of the semiquantitative microarray results. The comparison of a semiquantitative and a quantitative method allowed us to compare whether up-regulation and down-regulation of the same gene could be detected similarly using the two different methods. Most of the genes (82%) were detected similarly with both methods. The overlap of similar gene expression detected with both methods is in line with experience in our laboratory using the same techniques on HSC cells exposed to a radioisotope, with an overlap of 90.3% (20).
Maenhaut et al. recently examined gene expression levels in post-Chernobyl PTC and sporadic PTC (21). In contrast to our findings, no specific gene expression signature could be identified between the tumor groups. They used a human cDNA microarray (dual-channel platform) covering 2,400 genes (Micromax, Perkin-Elmer). The single-channel microarray we used covered the whole genome (29,098 genes), including approximately 6,000 additional genes from the Celera databank, and for signal detection used chemiluminescence instead of fluorescence, which is characterized by a larger (6-log) linear dynamic range. Moreover, Maenhaut et al. isolated RNA from smaller groups (12 post-Chernobyl PTC and eight sporadic PTC) for hybridization on their microarray, and fewer genes (12 compared to 92) were examined quantitatively using RTQ-PCR and on autonomous adenoma but not on PTC. Finally, our pathologist purposely selected only pure papillary histologies and no (e.g. follicular) subtypes, which also differed from the study of Maenhaut et al. Based on molecular and cytogenetic differences between sporadic and post-Chernobyl PTC, we expected to also find differences in gene expression. In view of the fact that a large proportion [60– 70%; for details, see refs. (6, 22)] of radiation-induced cancers are characterized by RET-PTC3 translocations or AKAP9-BRAF fusion (23) that differ from sporadic papillary carcinomas (2, 24–26), it is surprising that no changes in the gene expression levels were detected when the dual-channel platform was used. In contrast, in our experiments, post-Chernobyl PTC mRNA levels were increased over sporadic PTC in HRAS (up to 4.8-fold), KRAS (3.1-fold) and NRAS (2.8-fold). This is in line with the genotypes of our post-Chernobyl PTC group (Table 1). With comparative genomic hybridization (CGH), chromosomal imbalances were detected in 30% of post-Chernobyl PTC (27). The most frequent changes in DNA copy number involved chromosomes 2, 7, 13, 16 and 21; some of these specific alterations have been reported to be associated with aggressive biological behavior (27). A recently published CGH study concluded that, in particular, DNA gains in post-Chernobyl thyroid cancers appeared to be influenced by radiation exposure (28), which should have an impact at the transcription level. Another group demonstrated the possibility of inducing RET rearrangement in vitro using X rays (29). Using array CGH, they detected a predominant pattern of subtelomeric deletions and hypothesized that chromosome 10 may be a hotspot for radiation-induced damage that could cause alterations in gene expression levels. With in situ hybridization, the presence of AXL mRNA and co-expression of AXL and GAS6 using immunohistochemistry was detected in 71–77% of post-Chernobyl PTCs (30), which again supports the argument that changes in gene expression should be expected in the post-Chernobyl group.
With the microarray used, an over-representation of up-regulated genes coding for G-proteins, growth factors and oxidoreductases was found in post-Chernobyl PTC (Fig. 2). An elevation of HRAS in the apical cell surface using immunohistochemistry was found especially in malignant thyroid tissues (31), which is in line with our findings. However, changes at the protein level are not necessarily reflected at the gene expression level, as shown after radiation exposure in cells of a murine cell line (32). With regard to post-Chernobyl PTC, the literature focuses primarily on the examination of RAS mutations that were absent in radiation-induced PTC (33) but were detected for sporadic PTC (34, 35). Increased amounts of the mRNA and protein of growth factors like VEGF-A, PDGF-B and EGFR were found in thyroid tumor cells and PTC in vitro (36) as well as in vivo (37, 38). Increased serum concentrations of VEGF-A and EGFR protein were also reported in PTC (36, 39), eventually inducing systemic effects. It is interesting to note that increased levels were especially linked with PTC characterized by lymph node metastases and more aggressive behavior (40) and patients over the age of 45 (41). Again, post-Chernobyl PTC appeared more aggressive at presentation with lymph node metastasis (60–80%); other extrathyroidal spread and venous invasion were commonly observed (11–13). Our experimental design allowed us to compare the gene expression levels in the two groups characterized by aggressive behavior. Growth factor levels appeared to be increased (VEGF-A, 5.6-fold; EGFL9, 8.9-fold; PDGFC and PDGFRB, 8.4- and 3.8-fold, respectively; IGF1R, 2.2-fold; and IGBP1, 30-fold) in post-Chernobyl PTC group compared to the control PTC group, thus supporting the hypothesis that aggressiveness, lymph node metastasis and venous invasion might be linked with growth factor mRNA levels. Likewise, it has been shown by many authors that increased mRNA or protein levels of oxidoreductases such as cyclooxygenase 2 (PTGS2), superoxide dismutase (SOD1) and, in some studies, thyroid peroxidase (TPO) are characteristic of thyroid malignancies (42–45), which was again linked with age (46) and poorer prognosis. In our experiments, mRNA levels of TPO and SOD3 measured in post-Chernobyl PTC appeared to be increased two- and 3.4-fold over the control group, respectively, again supporting the hypothesis that tumor aggressiveness could be linked to oxidoreductase mRNA levels.
Apart from these up-regulated genes in our experiments with the microarrays we used, an over-representation of down-regulated genes coding primarily for immunoglobulins was found (Fig. 3). Infiltrates of immunoglobulin-producing lymphocytes are found in papillary cancer (47–51). Some authors have speculated that the immune response appears to be important in preventing metastasis and the recurrence of thyroid cancer (49, 51). Tumor-infiltrating lymphocytes have been associated with less extensive disease at diagnosis and improved disease-free survival for adults with PTC (52). The overall down-regulation of immunoglobulins found in our post-Chernobyl PTC group could be interpreted as an expression of reduced immune defense, thus supporting tumor aggressiveness. However, these results might also reflect differences in the type and number of tumor-associated lymphocytes. It is interesting to note that another group using microarrays demonstrated a high variability and heterogeneity in gene expression when they compared sporadic PTC with normal thyroid tissue that could be attributed in particular to differing immunoglobulin gene expression (53).
Using the LDA separately on each of the tumor samples analyzed, seven genes were identified that allowed us to completely distinguish post-Chernobyl PTC from control PTC (Figs. 4 and 5 and Table 3). An association with PTC is known for MMP1 and other metalloproteinases that appear to be linked with invasion (54), PTC malignancy (55), and promoting recurrence (56). As shown for VEGF-A and COX-2 (measured in plasma), an age-related association of metalloproteinases has been found (57). In a comparison of the level of MMP1 expression in post-Chernobyl PTC with the whole control group, discrimination of the two groups was demonstrated, with the post-Chernobyl PTC showing elevated MMP1 gene expression (Fig. 5A). Discrimination could be increased by selecting four tumor samples of the control group with a median age of 24 years, so that cases and controls were of comparable age and comparable aggressiveness (Fig. 5B). Nevertheless, our data provide only a hint that age differences might be less significant in explaining our results. Sporadic papillary thyroid cancer is a histology that is seldom found in children or juveniles but appears later in life (see above). Even when performing a larger analysis, it will be challenging to control for the potential confounding variable of age. Does that really matter? After diagnosis of PTC histology, irrespective of the patient's age, pathologists can only speculate about the etiology of the tumor. It was the focus of our study to compare the same histology, which developed due to a different etiology.
Another potential confounder could be iodine deficiency, which was present in the area where Belarusians lived and absent in the area where the control group originated. Presently it appears that iodine deficiency is associated with a follicular histology type of thyroid cancer but not with the incidence of a papillary type (58–60).
Information about the remaining six gene targets with PTC could not be found in the literature (for details, see Table 3).
It is preferable to make RTQ-PCR measurements on a group of patients completely different from those used for microarraying. When we did so, four post-Chernobyl PTCs and 33 sporadic PTCs remained for the LDA measurements. However, this had no influence on the discrimination between the groups (data not shown).
It was interesting that KRTAP2-1 and PAGE1 gene expression in most of the sporadic PTCs was not detectable using LDA technology (Table 4). Additionally, progestagen-associated endometrial protein (PAEP) and ribosomal protein S4, Y-linked 1 (RPS4Y1) showed a similar behavior in this tumor entity, while granzyme H (GZMH) was undetectable in most post-Chernobyl PTCs (Table 4). The several log-scale differences in expression of these genes between the groups could have implications at the protein level, thus making the associated proteins candidates for differentiation using immunohistochemistry.
In summary, our microarray data on post-Chernobyl PTC reflected a positive stimulus for tumor growth (over-representation of up-regulated genes coding for G-proteins and growth factors) and increased aggressiveness caused by overexpressed oxidoreductases, while the immune defense (over-representation of down-regulated genes coding for immunglobulins) appeared to be weakened. A complete differentiation of post-Chernobyl PTC from controls characterized by either comparable aggressiveness (PTC of older patients, median age 60 years) or comparable age (n = 4) was accomplished by a selection of seven gene targets.
Further examinations of larger groups are needed to determine whether these findings are applicable as a diagnostic tool for identifying radiation-induced PTC.
Acknowledgments
We would like to thank R. Obermair and C. Baaske for their skillful, accurate and careful technical assistance. This study was financed by the Ministry of Defence, Germany (10Z1-S-380406).
REFERENCES
TABLE 1
Summary of the Clinical Parameters and Gene Alterations of the 11 Post-Chernobyl PTCs and 41 Sporadic PTCs Studied
TABLE 2
Selection of Most Promising Genes for Differentiating between Groups using Four Gene Selection Algorithms
TABLE 3
Characterization of the Seven Most Promising Gene Targets for Distinguishing between the Groups
TABLE 4
Genes whose Expression was Beyond the Detection Sensitivity of the LDA in >80% of the Tumor Samples of the Chernobyl Group (one gene) or the Controls (four genes)