Iodine-131 dose-dependent gene expression: alterations in both normal and tumour thyroid tissues of post-Chernobyl thyroid cancers.

BACKGROUND: A strong, consistent association between childhood irradiation and subsequent thyroid cancer provides an excellent model for studying radiation carcinogenesis.
METHODS: We evaluated gene expression in 63 paired RNA specimens from frozen normal and tumour thyroid tissues with individual iodine-131 (I-131) doses (0.008-8.6 Gy, no unirradiated controls) received from Chernobyl fallout during childhood (Ukrainian-American cohort). Approximately half of these randomly selected samples (32 tumour/normal tissue RNA specimens) were hybridised on 64 whole-genome microarrays (Agilent, 4 × 44 K). Associations between I-131 dose and gene expression were assessed separately in normal and tumour tissues using Kruskal-Wallis and linear trend tests. Of 155 genes significantly associated with I-131 after Bonferroni correction and with ≥2-fold increase per dose category, we selected 95 genes. On the remaining 31 RNA samples these genes were used for validation purposes using qRT-PCR.
RESULTS: Expression of eight genes (ABCC3, C1orf9, C6orf62, FGFR1OP2, HEY2, NDOR1, STAT3, and UCP3) in normal tissue and six genes (ANKRD46, CD47, HNRNPH1, NDOR1, SCEL, and SERPINA1) in tumour tissue was significantly associated with I-131. PANTHER/DAVID pathway analyses demonstrated significant over-representation of genes coding for nucleic acid binding in normal and tumour tissues, and for p53, EGF, and FGF signalling pathways in tumour tissue.
CONCLUSION: The multistep process of radiation carcinogenesis begins in histologically normal thyroid tissue and may involve dose-dependent gene expression changes.

One of the most important health consequences of the 1986 Chernobyl nuclear power plant accident was a dramatic increase in thyroid cancer incidence among those who were children or adolescents at the time (Ron, 2007; Cardis and Hatch, 2011). While numerous epidemiological studies have established that this increase is primarily related to iodine-131 (I-131) thyroid dose received from the accident (Davis ; Cardis ; Tronko ; Brenner ; Zablotska ), the mechanisms of radiation-related thyroid carcinogenesis remain poorly understood. Most early post-Chernobyl molecular studies focused on evaluation of mutations in mitogen-activated protein kinase pathway including RET/PTC rearrangements and BRAF mutations (Nikiforov ; Thomas ). These studies suggested that although RET/papillary thyroid cancer (PTC) rearrangements are common in radiation-related thyroid cancer, they are also present in spontaneous thyroid cancer diagnosed at a young age and a large proportion of radiation-related cancers harbour no known mutations. Therefore, other more specific alterations induced by ionising radiation might exist ().

In recent years, additional attempts have been made to identify molecular changes in thyroid tissue specific to radiation exposure. These studies took advantage of biological materials available through the Chernobyl Tissue Bank (CTB) and high throughput technologies including genome-wide gene expression and DNA copy number variation (, 2007; Port ; Boltze ; Stein ; Hess ). Although each transcriptome study has reported a set of genes that discriminated post-Chernobyl thyroid cancers from spontaneous thyroid cancers, there is little consistency across studies at a gene level. This is likely due to small sample sizes, poor control of confounding factors, lack of methodological validation in independent samples, and different analytic approaches. The two studies that used exposed and unexposed cases matched on age and ethnicity, and validation of promising targets by qRT−PCR, provided intriguing findings. The study by Hess ), identified a gain of chromosome band 7q11 in exposed PTC cases compared with non-exposed PTC cases with young age of onset, whereas the study by Dom identified a set of genes permitting differentiation between normal thyroid tissue of I-131-exposed and non-exposed cases. These findings require further validation in independent populations as well as substantiation of the dose−response relationship based on individual dose estimates, because an assumption that all exposed cases received the same dose could be misleading and result in false-positive or false-negative associations.

We recently identified 11 genes with evidence of differential dose−expression relationship using specimens from well-characterised PTC cases who underwent thyroid surgery in the Ukrainian-American (UkrAm) cohort following standardized thyroid screening of ∼13 000 Ukrainian residents (<18 years at the time of the accident) with individual radioactivity measurements taken shortly after the accident (Abend ). We hypothesised that if dose-related gene expression patterns in tumour tissue truly reflect an important event in radiation carcinogenesis, they should differ from patterns observed in normal tissue. Although such different patterns were identified, we did observe significant dose-related changes in gene expression not only in tumour but also in contralateral normal thyroid tissue. This motivated us to validate additional genes with the evidence of dose-related expression in either normal or tumour thyroid tissue.

In the current study, we first conducted an initial screen in half of the cases to identify promising gene candidates that exhibited I-131-related expression in normal or tumour thyroid tissue based on whole-genome mRNA microarrays (phase I). We then validated the top candidates for normal or tumour tissue in the remaining cases using qRT−PCR (phase II). Our findings provide evidence that long-lasting dose-related changes in gene expression may be present in histologically normal thyroid issue and potentially point to early events in a multistep process of radiation carcinogenesis.

Materials and methods

Patients and tissue samples

Detailed description of the study cases has been provided previously (Abend ). We included 71 PTC cases, diagnosed in the UkrAm cohort between 1998 and 2008 at the Laboratory of Morphology of Endocrine System of the Institute of Endocrinology and Metabolism (IEM, Kiev, Ukraine), that had tumour and/or normal tissue RNA specimens in the CTB (http://www.chernobyltissuebank.com/). All tissue specimens were taken intraoperatively after patients signed informed consent forms approved by the institutional review boards (IRBs) of the IEM and the Coordinating Center of the CTB project (Imperial College Research Ethics Committee, London, UK). Annual review of the entire project was also provided by the IRB of the National Cancer Institute, USA.

Detailed operating procedures for the collection, documentation, and processing of frozen tumour and normal thyroid tissue specimens are available from the CTB website (http://www.chernobyltissuebank.com/) and were developed jointly with the Laboratory of Morphology of Endocrine System of the IEM and the Wales Cancer Bank.

Dosimetry

Dosimetric methods have been described elsewhere (Likhtarev et al, 2003, 2005, 2006). Briefly, individual I-131 thyroid doses and their uncertainties were estimated from the combination of thyroid radioactivity measurements, data on dietary and lifestyle habits, and environmental transfer models using a Monte-Carlo procedure with 1000 realisations per individual (Likhtarev ). For analysis, we used the arithmetic mean of each individual's 1000 realisations as the best estimate of I-131 dose corrected for thyroid masses typical of the Ukrainian population (Brenner ).

RNA extraction and quality control

Full details of RNA extraction can be obtained from the CTB website (http://www.chernobyltissuebank.com/). In brief, frozen thyroid tissue is homogenised using a tissue lyser (Qiagen, Hilden, Germany). RNA is extracted using Qiagen column-based systems. Standard 20 μl aliquots containing 5 μg of total RNA are frozen at −80 °C. Quality and quantity of isolated total RNA is measured spectrophotometrically (NanoDrop, PeqLab Biotechnology, Erlangen, Germany) while RNA integrity is assessed by the 2100 Agilent Bioanalyser (Life Science Group, Penzberg, Germany). For analysis, we used only RNA specimens with a ratio of A260/A280⩾2.0 (Nanodrop) and RNA integrity number (RIN)⩾7.5 for whole-genome microarray or RIN ⩾5.5 for qRT−PCR analyses (IMGM Laboratories, Martinsried, Germany).

Although the CTB provided 137 RNA specimens for 71 individuals with PTC, we were able to use only 126 paired (tumour/normal) RNA specimens corresponding to 63 individuals (Figure 1). Eleven RNA specimens from eight individuals were excluded because of either missing complementary tissue specimens (n=5) or failing our quality criteria (n=6).

Figure 1

Flow diagram showing our study design, included samples, gene expression experiments, and statistical/bioinformatics analyses.

Phase I: whole-genome microarray experiments

Genome-wide expression profiling was carried out using the Agilent oligo microarray (4 × 44 K format; Agilent Technologies, Waldbronn, Germany) combined with a one-color-based hybridisation protocol on 32 normal and 32 tumour RNA specimens from 32 randomly selected individuals (about half the sample set, Figure 1) as described in detail elsewhere (Abend ).

Phase I: statistical analysis of microarray data

We analysed the gene expression in normal and tumour thyroid tissues separately using quintile normalised log2-transformed probe signals of the normal and tumour tissues as an outcome. We used the non-parametric Kruskal−Wallis test (P kruskal) to compare gene expression across three dose categories (⩽0.30, 0.31–1.0, >1.0 Gy) with cut-off points approximately corresponding to tertiles of dose distribution among cases, and linear regression models that included continuous dose in a linear term (P linear). Only those gene transcripts that had a call ‘present' in at least 50% of RNA specimens from tumour or normal tissue were included in the analysis of gene expression (∼15 000).

Phase I: bioinformatic analysis of microarray data

All genes associated with dose at P kruskal<0.05 or P linear<0.05 in normal (n=3099) or tumour tissue (n=2233), respectively, underwent separate gene set enrichment analyses using PANTHER pathway software (http://www.pantherdb.org/) that group genes with similar biological function based on their annotation. To evaluate reproducibility of gene set enrichment analyses and to compare our findings with Dom ), we repeated analyses using DAVID database (Huang , 2009b).

Selection of gene candidates for independent validation

To adjust for multiple comparisons we employed Bonferroni correction. In total, 832 genes in tumour or normal tissue withstood such correction (P kruskal or P linear<10−6). To reduce the number of the promising candidates further, we used two approaches: (1) selecting genes with the lowest corrected P-values (P kruskal or P linear<10−25) and (2) selecting genes with corrected P-value <10−6 and absolute dose−response slope >1 (corresponding to two-fold increase in expression per dose category). These steps resulted in the selection of 155 genes. We further narrowed the candidates to 91 gene based on available inventoried TaqMan assays for qRT−PCR. We added three genes located at the chromosomal region 7q11 (AUTS2, MLXIPL, and LATL that had P kruskal/P linear values of 0.014/0.013, 0.016/0.012, and 0.042/0.019, respectively), because there was evidence in the literature for their involvement in radiation carcinogenesis (Hess ). We also added one housekeeping gene (GAPDH) to check its suitability for normalisation purposes. Thus, 95 gene candidates were selected for validation by qRT−PCR in phase II.

Phase II: quantitative RT−PCR experiments

To validate our phase I findings, we evaluated gene expression by qRT−PCR (TaqMan primer probe assays) on the remaining 28 normal and 30 tumour tissue RNA specimens using a low-density array (LDA). Because of RNA consumption for validation of differentially expressed genes in our previous study (Abend ) not all remaining 31 normal and 31 tumour RNA specimens were available for the current study. A 0.75-μg RNA aliquot of each RNA sample was reverse transcribed using a two-step PCR protocol (High Capacity Kit). 50 μl cDNA (equivalent to about 0.25 μg RNA) was mixed with 50 μl 2 × RT−PCR master mix and pipetted into two of eight fill ports of the LDA. Cards were centrifuged twice (1200 rpm, 1 min, Multifuge3S-R, Heraeus, Germany), sealed, and transferred into the 7900 qRT−PCR instrument. The qRT−PCR was run for 2 h following the qRT−PCR protocol for 384-well LDA format. All technical procedures for qRT−PCR were performed in accordance with standard operating procedures implemented in our laboratory in 2008 when the Bundeswehr Institute of Radiobiology became accredited according to DIN EN ISO 9001/2008. All chemicals for qRT−PCR using TaqMan chemistry were provided by Life Technologies (Darmstadt, Germany).

We used a previously established upper limit of the linear-dynamic range of our qRT−PCR, CT ⩽30 (Abend ). CT values were normalised relative to the median gene expression of the examined gene. This approach to normalisation was more robust compared with the use of recommended housekeeping gene expression (18S rRNA or GAPDH, Life Technologies homepage).

Finally, we compared mean differential gene expression (tumour relative to normal tissue) of whole-genome microarray data from phase I individuals with mean differential gene expression of qRT−PCR data from phase II individuals for the same gene to check for the reliability of the methods employed. Mean differential gene expression in phase I and phase II data was highly correlated (r2=0.86) and had an overall agreement of about 93% (data not shown), supporting the reliability of the methods.

Phase II: statistical analysis of qRT−PCR data

To confirm phase I findings, the phase II analyses used only individuals not included in phase I. Normalised CT values (corresponding to normalised gene expression values) of all genes were normally distributed. We analysed normalised gene expression values y in tumour and normal tissue jointly in linear mixed models with individual negative log2-transformed CT values as the outcome variable. The models included separate I-131 dose terms for tumour and normal tissues, and were adjusted for age at thyroid surgery (three categories), sex, and oblast or state of residence (Chernigov, Zhytomyr, and Kiev),

for subjects i (i=1, 2,…, 30) on the jth sample (j=1, 2 for tumour and normal tissues), where μ is the overall mean expression level and ɛ is the normally distributed error term. In model (1), the dose effect in tumour specimens is quantified by dosetumour, and the dose effect in normal tissue specimens is given by dosenormal. To evaluate the dose effect within each tissue type, we used I-131 dose in two ways: (1) in three independent dose groups, which leads to a 2 degree of freedom (d.f.) Wald test P-value for the null hypothesis H0 of no dose effect, and (2) in three ordered dose groups, corresponding to a 1 d.f. test for H0. Fold changes in expression associated with I-131 dose were computed as two to the power of the slopes, that is, 2^(-dosenormal) and expressed on a linear scale in our tables and figures. Parameters for the mixed models were estimated using the restricted maximum likelihood method incorporated in PROC MIXED (SAS, 2005; SAS 9.1.3).

Results

Characteristics of PTC cases

Of 63 cases included in our study, 56% were females and 54% were residents of the Chernigov oblast. Age at the time of the accident ranged from 0 to <18 years (mean 7.9 years) and cancers were diagnosed 12.5–21.6 years after the accident (mean 16.5 years). Overall mean I-131 thyroid dose was 1.25 Gy, ranging between 0.008 and 8.6 Gy, while means for the three dose categories were 0.11, 0.57, and 2.62 Gy, respectively. The most common histological subtype of PTC was mixed (48%) and the remainder consisted of follicular (25%), classic papillary (19%), and solid (8%) subtypes. The mean of the largest tumour diameter was 16.0 mm, with a range from 6.0 to 45.0 mm.

Whole-genome microarray

Of 19 596 gene mRNAs (41 079 transcripts) spotted on the whole-genome microarray, on average 73.4% (range: 63.3–91.0%) were distinguishable from background (expressed). The total number of gene transcripts significantly associated with I-131 dose either in normal or in tumour tissue specimens (Bonferroni corrected P kruskal or P linear<10−6) was 832; of these 95 gene candidates were selected for validation by qRT−PCR as described in Materials and Methods.

qRT−PCR

Of 95 genes assayed, the qRT−PCR data were available for 74 genes in normal tissue and 79 genes in tumour tissue because either no gene-specific amplification plots developed or plots were detected in less than half of the samples. For eight and six genes, the I-131 dose-related expression in normal or tumour tissue, respectively, was significant based on a categorical or ordinal trend test (Table 1, Figure 2). Expression of NDOR1 gene was significantly associated with dose both in normal and tumour thyroid tissues. The strongest association with I-131 dose, more than a two-fold increase or decrease in gene expression per dose category, was observed for ABCC3 and UBA3 genes in normal tissue and for SCEL and SERPINA1 genes in tumour tissue. Details of gene ontology and function for 13 genes significantly associated with I-131 dose in normal or tumour tissue are available in Supplementary Table 1.

Table 1

Summary statistics are shown for genes with significant dose−expression relationship based on qRT−PCR measurements

			Fold change per dose categorya
Gene	Cytoband	N	1 vs 0	95% CI		2 vs 0	95% CI		P 2 d.f.b	P 1 d.f.c
Normal tissue
ABCC3	17q22	20	2.00	0.99	4.00	2.22	1.08	4.55	0.080	0.042
C1orf9	1q24	28	1.82	1.25	2.63	1.19	0.85	1.67	0.012	0.190
C6orf62	6p22.3	28	1.32	0.99	1.75	1.47	1.12	1.92	0.018	0.004
FGFR1OP2	12p11.23	28	1.23	0.95	1.61	1.28	1.01	1.61	0.093	0.036
HEY2	6q21	27	0.91	0.56	1.47	1.61	1.04	2.50	0.058	0.050
NDOR1	9q34.3	28	1.18	0.84	1.67	1.39	1.02	1.92	0.126	0.022
STAT3	17q21.31	28	0.68	0.48	0.94	1.09	0.81	1.47	0.035	0.837
UBA3	3p14.1	28	2.22	1.25	4.00	0.90	0.53	1.54	0.018	0.931
Tumour tissue
ANKRD46	8q22.2	30	1.56	1.11	2.22	0.90	0.65	1.23	0.022	0.702
CD47	11q21	30	0.56	0.38	0.82	0.96	0.68	1.37	0.017	0.479
HNRNPH1	5q35.3	30	1.37	0.97	1.92	1.37	1.01	1.85	0.073	0.036
NDOR1	9q34.3	30	1.45	1.03	2.00	1.37	1.02	1.85	0.045	0.030
SCEL	13q22	30	0.24	0.12	0.50	0.40	0.20	0.78	<0.001	0.006
SERPINA1	14q32.1	30	0.16	0.07	0.37	0.98	0.45	2.13	<0.001	0.686

Fold-change: expressed on a linear scale.

2 d.f. test: models included I-131 dose coded in three separate categories for tumour and normal tissues adjusted for age at surgery, sex, and oblast.

1 d.f. trend test: models included linear terms for dose in tumour and dose in normal tissue adjusted for age at surgery, sex, and oblast; models were fitted to all data combined.

Figure 2

Gene expression is shown relative to the reference dose category (lowest I-131 thyroid dose set to 1, dashed grey line) for selected genes in (A) normal tissue (circles with white fills, first page) and (B) tumour tissue (circles with grey fills, second page). Circles represent mean gene expression values and error bars represent corresponding 95% confidence intervals.

Bioinformatic analysis of whole-genome microarray data

To identify genes that were over- or under-represented among those significantly associated with I-131 dose in microarray data, we conducted PANTHER classification analyses. Genes coding for protein classes such as nucleic acid binding, RNA binding, and ribosomal proteins were significantly over-represented in normal as well as tumour tissue analyses (Table 2). However, genes coding for proteins involved in FGF signalling, p53, or EGF signalling pathways were over-represented in tumour tissue analyses only.

Table 2

PANTHER classification of genes significantly associated with I-131 dose (after Bonferroni correction) in normal and/or tumour tissue

	Normal tissue		Tumour tissue
PANTHER classification	Over/under representation	P-values	Over/under representation	P-values
Cellular component
Cytoplasm	+	1.6E−05	+	2.0E−04
Organelle	+	7.4E−04	+	9.9E−04
Pathway
Ubiquitin proteasome pathway			+	6.5E−05
FGF signalling pathway			+	1.3E−04
p53 pathway feedback loops 2			+	1.3E−04
EGF receptor signalling pathway			+	1.6E−04
p53 pathway			+	2.6E−04
Integrin signalling pathway			+	8.2E−04
DNA replication			+	9.1E−04
Biological process
Primary metabolic process	+	2.0E−16	+	2.1E−21
Metabolic process	+	8.7E−17	+	3.9E−21
Protein metabolic process	+	4.1E−08	+	1.7E−13
Nucleobase, nucleoside, nucleotide, and nucleic acid metabolic process	+	1.5E−13	+	2.7E−13
Unclassified	−	2.6E−10	−	4.0E−09
Establishment or maintenance of chromatin architecture	+	2.3E−07	+	2.5E−05
Organelle organisation	+	2.2E−07	+	3.9E−05
Cell surface receptor linked signal transduction			−	2.9E−05
Cell−cell signalling			−	3.1E−05
Pattern specification process			−	5.9E−05
Cell cycle	+	2.7E−06
Protein transport	+	5.3E−05
Intracellular protein transport	+	5.3E−05
Molecular function
Nucleic acid binding	+	3.2E−16	+	1.4E−21
Binding	+	5.3E−16	+	1.6E−15
Structural constituent of ribosome	+	1.3E−08	+	1.8E−11
Unclassified	−	1.5E−08	−	4.9E−05
DNA binding	+	7.0E−06	+	2.9E−05
Catalytic activity	+	1.4E−05	+	9.0E−07
RNA splicing factor activity, transesterification mechanism			+	7.7E−06
Chromatin binding			+	1.1E−05
G-protein coupled receptor activity	−	8.8E−08
Transferase activity	+	1.9E−05
Receptor activity	−	2.1E−05
Protein class
Nucleic acid binding	+	4.3E−15	+	4.2E−18
RNA binding protein	+	7.0E−10	+	7.9E−17
Ribosomal protein	+	1.1E−07	+	4.0E−11
mRNA processing factor	+	8.9E−05	+	6.2E−08
G-protein coupled receptor	−	5.5E−07
DNA binding protein	+	1.2E−06
Homeobox transcription factor	−	5.8E−06
Helix-turn-helix transcription factor	−	5.8E−06
Histone	+	1.0E−05
Transferase	+	1.4E−05
Signalling molecule			−	4.0E−06

+/− indicates number of observed genes coding for certain biological processes to be over-represented (+) or under-represented (−) relative to the expected number of genes based on the whole-genome microarray annotation used.

To evaluate reproducibility of gene set enrichment analyses and to compare our findings with Dom ), analyses were repeated using DAVID software. In the normal tissue genes coding for proteins involved in the ribosomes, translational elongation, protein modification (phosphorylation or acetylation), and intracellular transport were significantly enriched (P-values between 1 × 10−7 and 5 × 10−35). Genes coding for cell-cycle processes were also significantly enriched (P=0.0003) as well as the genes coding for chronic myeloid leukaemia pathway as defined by KEGG (P=0.04). In tumour tissue, genes involved in those pathways found through PANTHER analyses were enriched, although P-values were slightly higher (data not shown).

Discussion

The relationship between irradiation at a young age and risk of thyroid cancer is strong and strikingly consistent, and thus this tumour provides an excellent model for studying radiation carcinogenesis in humans. We employed measurement-based individual I-131 doses and RNA specimens from fresh frozen thyroid tissue of cases with PTC diagnosed in the UkrAm cohort to evaluate dose−expression relationship in tumour and normal thyroid tissues separately. Using the same case series, we previously identified 11 genes potentially involved in radiation-related thyroid carcinogenesis based on analyses of differential (defined as a difference in gene expression between tumour and normal thyroid tissues of the same individuals) dose-dependent gene expression across the entire genome (Abend ). In the analyses of differential gene expression, we observed evidence of dose-related gene expression not only in tumour but also in corresponding normal thyroid tissue and this motivated us to evaluate dose-dependent gene expression within each tissue type separately. Here, we identified eight genes in normal tissue and six genes in tumour tissue that were significantly associated with I-131 dose, including one gene (NDOR1) associated with I-131 in both tissue types. These findings are important as they suggest that radiation-related changes could occur in histologically normal thyroid tissue and may represent early steps in radiation carcinogenesis.

Using a similar approach, Dom recently compared the gene expression in normal and tumour thyroid tissues of cases exposed and unexposed to I-131. They identified a gene expression signature of 403 genes that discriminated normal tissues of exposed and unexposed individuals and validated seven genes by qRT−PCR. Since we started our study and selected promising targets before the publication of Dom's results, there is no overlap in the set of genes validated by qRT−PCR between the two studies. However, our results of gene set enrichment analyses based on whole-genome microarray data were similar to Dom in that, among genes with expression significantly related to I-131 dose in normal tissue, we found strong over-representation of genes coding for nucleic acid processing, RNA binding, and ribosomal proteins. The finding of I-131-related gene expression in normal thyroid tissue in two independent studies is unlikely to be explained by the presence of microcarcinomas as histological samples in both studies were reviewed by an international pathology panel and, to be detectable in whole-genome microarray, microcarcinomas had to involve a substantial proportion of cells. It is also unlikely that normal tissue findings could be explained by ‘field' effect of thyroid cancer because corresponding normal tissue was obtained from contralateral thyroid lobe and such ‘field' effect had to be related to I-131 dose. Unlike in Dom's study, our cases were ascertained within a well-characterised cohort (Tronko ; Brenner ). To control for potential differences in regional iodine intake and other socioeconomic factors, we adjusted analyses for oblast of residence. One important difference from Dom is that we had individual dose estimates and were able to evaluate dose−response relationship rather than comparing exposed and non-exposed cases. In summary, findings from the two independent studies taken together suggest that I-131-related gene expression in histologically normal thyroid tissue may represent long-lasting consequence of radiation exposure and/or early events in a multistep process of radiation carcinogenesis. However, our data and the data in the Dom study represent single time points in each case, but covering several decades after radiation exposure. It would be more straight forward showing gene expression changes over time on an individual base using several samples per individual. Unfortunately, biological samples such as that were not available for this study, but are currently examined in the context of another study. Nevertheless, the consistency of gene expression changes observed in different individuals with post-exposure times covering two decades and its relation to dose which has been shown in two independent studies comprising two different groups of irradiated individuals supports the interpretation of long-lasting gene expression changes to explain our results.

One potential mechanism through which I-131-related expression in normal tissue could perpetuate is epigenetics (Aypar ; Herceg and Vaissiere, 2011). Epigenetic changes indirectly affect DNA by altering DNA methylation, chromatin remodelling, and microRNA expression rather than DNA structure. The over-representation of genes coding for nucleic acid processing, RNA binding, and ribosomal proteins among those genes significantly associated with I-131 in normal thyroid tissue seems consistent with this idea. In contrast to normal thyroid tissue, genes involved in regulation of FGF signalling, p53, or EGF signalling pathways were over-represented in analyses of tumour tissue. Interestingly, recent evidence suggests that these pathways are deregulated in multiple sporadic cancers (Beroukhim ). Iodine-131-related changes in tumour tissue may result not only from epigenetic changes but also from copy number alterations (Kang ), shown to shape cancerous transcriptome (Ortiz-Estevez ). This idea finds some support in that AUTS2 gene located within 7q11.22–7q11.23 region, previously found to be amplified in post-Chernobyl tumours (13), showed a suggestive dose dependency in the tumour tissue based on 2 d.f. and 1 d.f. trend (P=0.10 and P=0.09, respectively) in our study. Given that dose-dependent expression in genes coding for FGF signalling, p53, or EGF signalling pathways was found only in cancerous thyroid tissue, this may represent later events in a multistep process leading to cancer development.

The 13 genes that we validated in an independent case set by qRT−PCR are all located on different chromosomes and are involved in different biological processes including mRNA processing and environment interaction (HNRNPH1), developmental processes (HEY2 and UBA3), transport of electrolytes, anions and molecules (CD47, ABCC3, and NDOR1), regulation of cell proliferation and cell death (STAT3, ANKRD46, and FGFR1OP2); we also identified less-well-characterised genes (SCEL, C6orf62, C1ORF9, and SERPINA1, Supplementary Table 1). In a recent publication, SERPINE1, which belongs to the same serine proteinase inhibitor (serpin) superfamily as SERPINA1, was among the five genes that discriminated sporadic and radiation-related PTCs (either post-radiotherapy or post-Chernobyl cases, Ory ). Preliminary (unpublished) comparison of our data with gene expression changes examined in the peripheral blood of former Mayak workers with occupational prolonged radiation exposure revealed no similarities which might be caused by exposure differences and other materials (blood, not thyroid tissue) used for gene expression measurements.

Our study has several strengths. First, we used individual I-131 dose estimates based on radioactivity measurements taken shortly after the accident (Likhtarev ; Brenner ). Second, the cancer cases arose within a well-characterised cohort screened for thyroid cancer according to a standardized protocol and irrespective of dose, minimising the impact of unmeasured confounding. The total number of cases with RNA samples (n=63) used for whole-genome microarray analysis (phase I, n=32) and qRT−PCR (phase II, n=28 normal tissue analyses and n=30 tumour tissue analyses) is larger than in previous studies of irradiated populations. Comparison of gene expression measurements performed by whole-genome microarrays (phase I) and qRT−PCR (phase II) together with evaluation of methodological variability of qRT−PCR is all consistent with previous work (Port ; Seidl ) and lessens concerns that our findings are due to methodological artifacts. For instance, it is generally agreed to introduce a two-fold gene expression difference over reference values as an upper and lower limit to address three different sources of variance, namely methodological, intra-individual, and inter-individual variance. From previous work we do know that in 95 from 100 cases these sources of variances can be controlled for by introducing a two-fold difference (Riecke ) and therefore, our findings are likely not to be caused by chance. We also selected microarray gene targets for qRT−PCR validation purposes based on the availability of corresponding inventoried TaqMan chemistry, leaving only coding mRNAs available for gene expression analysis, but missing, for example, long non-coding RNA species which are also covered by the microarrays.

However, there are several limitations to keep in mind when interpreting the results of our study. We did not account for uncertainties in the dose estimates, 95% of which are typically attributable to unknown thyroid gland mass and I-131 content in the thyroid gland in 1986 (Likhtarev ). However, these dose estimates compare favourably to other studies of environmentally exposed populations that exclusively relied on retrospective dose reconstruction and did not have individual measurements of radioactivity. The small sample size limited our ability to accurately quantify the magnitude of dose−response relationship and to characterise its shape.

In summary, our study is among the first to provide direct human data on long-term gene expression in thyroid tissue in relation to individual I-131 doses. Our finding of dose-related gene expression found in normal and tumour thyroid tissues with additional changes occurring in the tumour tissue suggests a multistep process of radiation carcinogenesis which may start in histologically normal tissue.