Skip to main content
Download PDF

Transcriptomic analyses of the radiation response in head and neck squamous cell carcinoma subclones with different radiation sensitivity: time-course gene expression profiles and gene association networks

Radiation Oncology volume 11, Article number: 94 (2016) Cite this article

Abstract

Background

Acquired and inherent radioresistance of tumor cells is related to tumor relapse and poor prognosis – not only in head and neck squamous cell carcinoma (HNSCC). The underlying molecular mechanisms are largely unknown. Therefore, systemic in-depth analyses are needed to identify key regulators of radioresistance. In the present study, subclones of the CAL-33 HNSCC cell line with different radiosensitivity were analyzed to identify signaling pathways related to the different phenotypes.

Methods

Subclones with altered radiosensitivity were generated by fractionated irradiation of the parental CAL-33 cells. Differences in radiosensitivity were confirmed in colony formation assays. Selected subclones were characterized at the genomic and transcriptomic level by SKY, array CGH, and mRNA-microarray analyses. Time-course gene expression analyses upon irradiation using a natural cubic spline regression model identified temporally differentially expressed genes. Moreover, early and late responding genes were identified. Gene association networks were reconstructed using partial correlation. The Reactome pathway database was employed to conduct pathway enrichment analyses.

Results

The characterization of two subclones with enhanced radiation resistance (RP) and enhanced radiosensitivity (SP) revealed distinct genomic and transcriptomic changes compared to the parental cells. Differentially expressed genes after irradiation shared by both subclones pointed to important pathways of the early and late radiation response, including senescence, apoptosis, DNA repair, Wnt, PI3K/AKT, and Rho GTPase signaling. The analysis of the most important nodes of the gene association networks revealed pathways specific to the radiation response in different phenotypes of radiosensitivity. Exemplarily, for the RP subclone the senescence-associated secretory phenotype (SASP) together with GPCR ligand binding were considered as crucial. Also, the expression of endogenous retrovirus ERV3-1in response to irradiation has been observed, and the related gene association networks have been identified.

Conclusions

Our study presents comprehensive gene expression data of CAL-33 subclones with different radiation sensitivity. The resulting networks and pathways associated with the resistant phenotype are of special interest and include the SASP. The radiation-associated expression of ERV3-1 also appears highly attractive for further studies of the molecular mechanisms underlying acquired radioresistance. The identified pathways may represent key players of radioresistance, which could serve as potential targets for molecularly designed, therapeutical intervention.

Background

Head and neck squamous cell carcinoma (HNSCC) develops in approx. 139,000 individuals per year in Europe with a survival rate of approx. 70 % at 1 year and approx. 40 % at 5 years after therapy [1]. More than 90 % of head and neck cancers are classified as HNSCC and originate from the oral cavity, nasopharynx, oropharynx, hypopharynx, or larynx, respectively [2]. The major risk factors for HNSCC are tobacco smoking, alcohol abuse, and poor oral health [3, 4]. For oropharyngeal cancers, infection with high-risk human papilloma viruses is another important risk factor [5]. Thus, HNSCC is a very heterogeneous cancer entity also in terms of therapy response. Surgical resection followed by radio(chemo)therapy is the standard treatment of HNSCC patient [6, 7]. In locally advanced HNSCC, surgery is often limited by the complex anatomy of the affected region and, therefore, definitive radiochemotherapy is an important treatment option. However, acquired and/or inherent radioresistance of tumor cells is a common cause for tumor relapse and poor prognosis. Tumor cells derived from HNSCCs after radiotherapy have been reported to be more radioresistant than cell lines established prior to therapy, thus strengthening the clinical relevance of acquired radioresistance [8]. Along these lines, it was proposed that fractionated irradiation might preferentially eradicate radiosensitive cells, whereas radioresistant cells remain largely untouched. Accordingly, recurrent tumors mostly consist of radioresistant cells [8]. Although different potential mechanisms of radioresistance have been proposed and extensively studied, the underlying molecular details remain largely unknown [9]. Systemic in-depth analyses are needed in order to identify the master regulators of acquired radioresistance, which could serve as potential biomarkers and future therapeutic targets in novel combined modality approaches.

In this study, we characterized two subclones (#303 and #327) derived from the CAL-33 HNSCC cell line, which were generated by fractionated radiation treatment of the parental cells. CAL-33 is an HPV-negative HNSCC cell line that has been established by Gioanni et al. (1988) from a biopsy specimen prior to treatment from a squamous cell carcinoma of the tongue from a male patient [10]. The subclones derived thereof differed in radiosensitivity when compared to the parental CAL-33 cell line. Interestingly, one subclone was significantly more radioresistant, whereas the other one was significantly more radiosensitive. In order to identify potential key regulators of altered radiation sensitivity, the subclones were characterized on the genomic and transcriptomic level. Furthermore, time-course gene expression analyses were performed upon irradiation, and gene association network reconstruction and pathway enrichment analyses were utilized to identify signaling pathways related to the observed radiation phenotypes.

Methods

Cell culture

The human head and neck squamous cell carcinoma (HNSCC) cell line CAL-33 was obtained from the German collection of microorganisms and cell cultures (DSMZ). Cells were maintained in DMEM GlutaMAX I medium supplemented with 10 % FCS and 1 % penicillin/streptomycin and cultured at 37 °C/5 % CO2. Cells were mycoplasma-free as tested by the MycoAlert (Lonza) mycoplasma detection kit.

Generation of clones with altered radiation resistance

In order to generate radioresistant subclones of the parental CAL-33 cell line, exponentially growing CAL-33 cells were exposed to fractionated irradiation (Mueller RT-250 γ-ray, tube Thoraeus Filter, 200 kV, 10 mA) according to a schedule commonly used in radiotherapy: A total dose of 20 Gy was given in daily fractions of 2 Gy 5 times per week. For each week, 2 days of recovery time were included. Afterwards, cells were cloned by limiting dilution procedure and grown for 4 to 12 weeks.

Colony forming assay

Clonogenic survival was determined in colony formation assays as described previously [11]. Briefly, cells were seeded into 6-well plates, allowed to adhere for 4 h, and irradiated at 0, 1, 2, 4, 6 or 8 Gy, respectively (Mueller RT-250 γ-ray, Thoraeus Filter, 200 kV, 10 mA). 14 days after irradiation, colonies were fixed and stained with methylene blue. Only colonies consisting of at least 50 cells were scored. Each assay was carried out in duplicates of three different cell densities per each irradiation dose. Results from three independent experiments were subjected to linear-quadratic regression analyses employing the maximum likelihood approach. Differences between curves were evaluated using F-test [12].

Proliferation and cell cycle analyses

Proliferation rates were determined over a period of three days upon seeding 20,000 cells per well into 24-well plates. Cells were harvested by trypsinization, and total cell numbers were determined by manual counting. Cell numbers were plotted against the growth time, and doubling times were calculated by semi-log regression analyses in the exponential growth phase. Dynamics of cell cycle distribution upon irradiation at 4 Gy was analyzed by flow cytometric phospho-histone H3(S10)/propidium iodide (PI) staining as described in [13]. Briefly, cells were collected by trypsinization and fixed in 70 % ethanol. After extensive washing, cells were stained with anti-phospho-histone-H3-Alexa488 (pH3(S10)) antibody (New England Biolabs, Frankfurt, Germany) and PI/RNase staining solution (BD Biosciences, Heidelberg, Germany). Data of 10,000 cells were recorded on an LSRII flow cytometer (BD Biosciences), and cell cycle analyses were performed by using FlowJo 7.6.5 software (Tree Star Inc., Ashland, OR, USA).

Spectral karyotyping (SKY)

Metaphase chromosome spreads were prepared from untreated CAL-33 parental cell line and generated CAL-33 sublines. Colcemid (Roche) was added at a final concentration of 0.1 μg/ml to the culture medium of exponentially growing cells at a density of 6 × 106 cells per 75 cm2. After 3 h of incubation time, cells were washed with PBS, trypsinized, suspended in fresh culture medium followed by hypotonic KCl treatment (75 mM) at 37 °C for 25 minutes. Following centrifugation, cells were resuspended in 2–3 ml of fixation solution and approximately 40–50 μl of cell suspension was dropped on several microscope slides. After one week of ageing at room temperature, spectral karyotyping was performed as described by Hieber et al. [14]. The karyotype of each cell line was determined based on a minimum of 15 metaphases. Chromosomal aberrations were detectable by color junctions within affected chromosomes. Spectral imaging and image analysis were performed with a SpectraCube system and SkyView imaging software (both from Applied Spectral Imaging).

Genomic copy number typing (array CGH)

In order to characterize copy number changes of parental CAL-33 and generated CAL-33 sublines, array comparative genomic hybridization analysis (array CGH) was performed on high-resolution oligonucleotide-based SurePrint G3 Human 180 k CGH microarrays (AMADID 252206, Agilent Technologies). DNA from non-irradiated samples was isolated using the QIAamp DNA Mini Kit (Qiagen). The DNA concentration was quantified with the NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, Germany). Slight modifications of the original Agilent array CGH protocol were introduced. 250 ng isolated cell line DNA and 250 ng sex-mismatched normal reference DNA (Promega) were labeled with Cy3 and Cy5, respectively, using the CGH labeling kit for oligo arrays (Enzo) following the Enzo’s protocol. Microcon YM-30 columns (Millipore) were used to remove the unincorporated nucleotides. Subsequent labeled DNA hybridization, washing and scanning of the CGH arrays were continued according to the Agilent’s protocol. The fluorescence intensities were extracted as text files with the Feature Extraction software 10.7 (Agilent Technologies). Obtained data were imported into the R statistical platform (version 3.2.2, www.r-project.org) and filtered for quality outliers using the QA measurements generated by the Feature Extraction software. Experimental artifacts were removed from the array CGH data using spatial normalization as suggested and described in MANOR R-package manual [15] and [16]. Array CGH profiles containing a wave bias that appear as waves in plots were removed using ridge regression based algorithms implemented in NoWaves R-package available from http://www.few.vu.nl/~mavdwiel/nowaves.html [17]. Following normalization and/or wave bias removal, data were segmented using circular binary segmentation algorithms as implemented in DNAcopy R-package [18] in order to detect breakpoints and levels in single array CGH profiles [19]. Chromosomal gains and losses were determined using CGHcall algorithm implemented in CGHcall R-package [20]. To reduce data complexity, copy number calls were transformed into regions using the R-package CGHregions [21].

Irradiation and sample preparation

Cells were seeded into 6-well plates and allowed to adhere for 16 h. The numbers of plated cells were adjusted to the incubation times: For samples collected 0.25, 2, 7, 12, and 24 h after irradiation, 3.5 × 105 cells/well were seeded, whereas for those collected after 48, 72, and 96 h, 1.75 × 105 cells/well were used. Cells were irradiated at 0 or 8 Gy of gamma-irradiation (Mueller RT-250, Thoraeus Filter, 200 kV, 10 mA) at a dose rate of 1.3 Gy/min, and samples were collected after 0.25, 2, 7, 12, 24, 48, 72, and 96 h by scraping on ice. The washed, dry cell pellet was snap frozen and stored at −80 °C. Samples from 3 independent experiments were used for subsequent transcriptomic analyses. Total RNA was isolated using the Rneasy Mini Kit (Qiagen) including a DNase digestion step, according to the manufacturer’s protocol. The concentration of RNA was quantified with a Qubit 2.0 Fluorometer (Life Technologies), and RNA integrity was confirmed with a Bioanalyzer 2100 (Agilent Technologies). Samples with RNA integrity number (RIN) >7 were used in subsequent gene expression microarrays analyses.

Global gene expression profiling

Global gene expression profiling of all CAL-33 cell lines was performed on SurePrint G3 Human Gene Expression 8x60k microarrays (Agilent Technologies, AMADID 28004) using 60 ng of total RNA according to the manufacturer’s protocol (one-color Low Input Quick Amp Labeling Kit, Agilent Technologies). Raw gene expression data were extracted as text files with the Feature Extraction software 11.0.1.1 (Agilent Technologies). All data analysis was conducted using the R statistical platform (version 3.2.2, www.r-project.org) [22]. Data quality assessment, filtering, preprocessing, normalization, batch correction based on nucleic acid labeling batches and data analyses were carried out with the Bioconductor R-packages limma, Agi4x44PreProcess and the ComBat function of the sva R-package [23–25]. All quality control, filtering, preprocessing and normalization thresholds were set to the same values as suggested in Agi4x44PreProcess R-package user guide [25]. Only HGNC annotated genes were used in the analysis. For multiple microarray probes representing the same gene the optimal probe was selected according to the Megablast score of probe sequences against the human reference sequence (http://www.ncbi.nlm.nih.gov/refseq/) [26]. If the resulted score was equal for two or more probes, the probe with the lowest differential gene expression FDR value was kept for further analyses since only one expression value per gene was allowed in subsequent gene association network (GAN) reconstruction analysis.

Differential gene expression analysis

The time-course differential gene expression analyses were conducted between irradiated and control cells (sham-irradiated) using a natural cubic spline regression model with three degrees of freedom as described in splineTimeR R-package [27]. Obtained p-values were adjusted by the Benjamini-Hochberg method for false discovery [28]. Genes with an adjusted p-value (FDR, false discovery rate) lower than 0.05 were considered as differentially expressed and associated with radiation response. For the status quo experiment that compares the derived CAL-33 clones with the parental CAL-33 cell line, genes were considered as differentially expressed when a log2 fold-change was higher than 0.5 and a FDR value lower than 0.05.

Identification of early and late responding genes

Temporally differentially expressed genes with fold-change above 2.0 or below 0.5 in any measured time points within the first day after irradiation were considered as early responding genes. Respectively, genes with fold-change above 2.0 or below 0.5 within the second, third or fourth day of irradiation were considered as late responding.

Gene association network reconstruction and identification of important nodes in the reconstructed networks

Temporally differentially expressed genes were subjected to gene association network (GAN) reconstruction using a regularized dynamic partial correlation method [29]. Pairwise relationships between genes over time were inferred based on a dynamic Bayesian network model with shrinkage estimation of covariance matrices as implemented in the GeneNet R-package [30]. Analyses were conducted with a posterior probability of 0.95 for each potential undirected edge. Further, in order to determine the importance of each node in the reconstructed association networks, graph topological analyses based on centrality measures were applied [31]. Three most commonly used centrality measures: degree, shortest path betweenness and closeness describing the importance of gene in a network were combined into one centrality measure [32]. For each gene the three centrality values where ranked and the consensus centrality measure for each node was defined as the mean of the three independent centrality ranks.

Pathway enrichment analysis

The Reactome pathway database was used to conduct the pathway enrichment analysis in order to further investigate the functions of the selected sets of differentially expressed genes [33]. Only pathways containing not more than 600 genes and not less than 20 genes were considered. Thereby, too general and too specific pathways were excluded from the analysis. Statistical significance of enriched pathways was determined by one-sided Fisher’s exact test. The resulting p-values were adjusted for FDR using the Benjamini-Hochberg method.

qRT-PCR technical validation of gene expression data

For technical validation of the gene expression microarray data, RNA samples (500 ng) were reversely transcribed using the QuantiTect Reverse Transcription Kit (Qiagen) according to the manufacturer’s protocol and subjected to qRT-PCR reactions (10 μl) on a ViiA 7 qPCR system (Life Technologies). The following Taqman® Assays were used (Life Technologies): AKT3 (Hs00987350_m1), GADD45A (Hs01077132_m1), MAL (Hs00360838_m1), HOPX (Hs04188695_m1), HYAL3 (Hs00185910_m1), TUBGCP3 (Hs00902139_m1), RGS16 (Hs00892674_m1), TNFAIP3 (Hs00234713_m1). ACTB (Hs01060665_g1) and GAPDH (Hs99999905_m1) served as endogenous reference genes. Relative expression levels were calculated using the ΔΔCt method and Spearman correlation analyses with microarray data were performed. Validation was considered successful for Spearman’s rho > 0.5. Additionally, the fold-change values obtained from microarrays and qRT-PCR were compared.

Results and discussion

Tumor relapse after radiochemotherapy in HNSCC is often linked to intrinsic and/or acquired radioresistance of tumor cells. However, the underlying molecular mechanisms remain largely unknown [9]. To gain knowledge on this fundamental and clinically relevant process, we established an in vitro model of acquired phenotypes of discrepant radiosensitivity in CAL-33 cells. The underlying molecular mechanisms were investigated by static and dynamic global mRNA expression analyses with subsequent network reconstruction and pathway enrichment analyses.

CAL-33 subclones with different phenotypes of radiosensitivity and cytogenetic characteristics

The parental cell line CAL-33 was repeatedly irradiated in order to generate subclones with different phenotypes of radiosensitivity. To analyze acquired alterations in radiosensitivity of the derived CAL-33 subclones, long-term survival upon gamma-irradiation was assessed by colony formation assays (Fig. 1). For further analyses we selected two subclones. Both subclones #303 and #327 showed statistically significant differences (p-values < 0.0001) when compared to the parental CAL-33 cells. Interestingly, subclone #303 showed increased radiosensitivity (sensitive phenotype, SP), whereas subclone #327 was more radioresistant (resistant phenotype, RP) – particularly in the dose range > 4 Gy.

Fig. 1
figure1

Dose-survival curves of parental CAL-33 cell line and derived subclones. The linear-quadratic cell survival curves were fitted to the measured data using maximum-likelihood method. Both subclones (SP and RP) showed statistically significant difference (p-values < 0.0001) in response to ionizing radiation compared to the parental CAL-33 cell line

Full size image

At first glance, the emergence of more radiation sensitive subclones appears unexpected and might be related to clonal evolution that was initiated by the irradiation-induced genomic alterations and that might occur at sublethal doses. This phenomenon has also been observed in previous studies [34, 35]. In comparison to other HNSCC cell lines, CAL-33 is very radioresistant a priori and this might explain why it appears to be very difficult to generate subclones with an even more resistant phenotype.

To analyze structural and numerical chromosomal aberrations in comparison to the parental cell line, the subclones were cytogenetically characterized by SKY analyses. Structural and numerical aberrations identified by SKY in the parental CAL-33 cell line involved chromosomes 3, 7, 8, 9, 16, 18, 20, X, and Y (Fig. 2a). Structural and numerical aberrations of the SP subclone included chromosomes 2, 3, 7, 8, 9, 11, 16, 18, 20, 21, X, and Y (Fig. 2b). Chromosomes 1, 3, 4, 5, 7, 8, 9, 14, 16, 18, 20, X, and Y were affected by aberrations in subclone RP (Fig. 2c).

Fig. 2
figure2

SKY results of the CAL-33 cell lines. The yellow arrows point the common for all CAL-33 cell lines marker chromosomes. The additional chromosomal rearrangements in analyzed clones compared to the parental CAL-33 cells are marked with white arrows. Structural and numerical aberrations in the parental CAL-33 cell line (a) involve chromosomes 3, 7, 8, 9 16, 18, 20, X, and Y. Aberrations of subclone SP (b) include chromosomes 2, 3, 7, 8, 9, 11, 16, 18, 20, 21, X, and Y. Chromosomal aberrations of subclone RP (c) affect chromosomes 1, 3, 4, 5, 7, 8, 9, 14, 16, 18, 20, X, and Y

Full size image

The obtained SKY results were complemented with copy number analysis by array CGH (Fig. 3 and Additional file 1: Table S1). Array CGH analysis identified 173 regions with aberrant copy number status from which 78, 111 and 132 regions were affected by DNA gains or DNA losses in the CAL-33 parental cells, SP, or RP subclones, respectively (Fig. 3 and Additional file 1: Table S1 and Additional file 2: Table S2). 68 copy number alterations (for SP) and 85 copy number alterations (for RP) were different from the parental CAL-33 cell line. Cytogenetic studies of both clones showed distinct genomic changes in comparison to the parental cell line indicating genomic key alterations for irradiation-related phenotypes on chromosomes 1, 2, 3, 4, 5, 8, 11, 14, 16 and 21. In addition, recurrent CAL-33-specific alterations on chromosomes 3, 7, 18, X and Y were observed showing the authenticity of the newly generated cell lines.

Fig. 3
figure3

Array CGH profiles of the CAL-33 cell lines. The array CGH profiles of all of the three cell lines show copy number alterations on several chromosomes: (a) parental CAL-33 cell line, (b) subclone SP, (c) subclone RP. The green bars (starting from the top) represent DNA copy number gains at the corresponding position in the genome, whereas the red bars (starting from the bottom) indicate DNA copy number losses. Bars reaching beyond the middle axis (probability >0.5) were called as gains or losses

Full size image

Analysis of proliferation rates and cell cycle distribution

Proliferation rate and cell cycle distribution are important factors, which can affect radiosensitivity and/or resistance. Accordingly, we performed proliferation and cell cycle analyses. In comparison to CAL-33 parental cells, both subclones displayed prolonged doubling times (29 h for subclone SP, 30 h for subclone RP, and 24 h for the parental cell line, Additional file 3: Figure S1, A). This might be due to the observation that under exponential growth conditions, the percentage of mitotic cells in both subclones was decreased (2.1 % for subclone SP, 1.5 % for subclone RP, and 2.9 % for the parental cells), but fails to explain the differences in radiosensitivity (Additional file 3: Figure S1, B). With regard to irradiation-induced G2-arrest, the sensitive subclone SP revealed virtually identical dynamics as the parental CAL-33 cells, whereas in the resistant subclone RP G2-arrest was initially delayed, but also reached its maximum around 12 h after irradiation, and afterwards appeared to be prolonged (Additional file 3: Figure S1, C). In fact, prolonged cell cycle arrest can contribute to radioresistance as cells have more time to repair irradiation-induced DNA damage. However, extended cell cycle arrest can also be indicative for delayed DNA repair. In order to address the mechanisms underlying radioresistance on a molecular level, next we therefore performed unbiased transcriptome analyses of the CAL-33 subclones.

mRNA gene expression analysis of the CAL-33 subclones

Microarray analyses allowed the identification of differences in basal gene expression levels between the derived subclones and the parental CAL-33 cell line. We identified 523 and 1292 differentially expressed genes (FDR < 0.05) for clones SP and RP, respectively, whereas 361 of the genes were overlapping (Additional file 4: Table S3). It is interesting to note that the RP clone exhibited more pronounced transcriptional differences than the SP clone.

Eight of the differentially expressed genes (SP and/or RP clone versus the parental CAL-33 cell line) were arbitrarily chosen for technical validation of the microarray data. Correlation analysis between qRT-PCR and microarray data showed a strong correlation for seven out of eight validated genes (Additional file 5: Table S4). The microarray and qRT-PCR derived fold-changes were in a good agreement.

Subsequently, the differentially expressed genes were subjected to pathway enrichment analyses in order to determine pathways common or specific to the radioresistant or radiosensitive phenotype. In total, 65 and 455 pathways were significantly enriched (FDR < 0.1) for the SP and RP clone, respectively (Additional file 6: Table S5). The top 100 identified pathways ordered according to the highest matching (percentage of differentially expressed genes of all genes in a pathway) were grouped into major pathways (Table 1).

Table 1 Significantly enriched pathways of genes differentially expressed in subclones SP and RP compared to the parental CAL-33 cells
Full size table

This resulted in a set of commonly deregulated pathways compared to the parental cell line but not specific for a particular phenotype of radiation sensitivity. Moreover, pathways specific to the radiosensitive or radioresistant phenotype were identified. These comprised mainly pathways that are known to be affected by ionizing irradiation in HNSCC [36–38].

Integration of copy number changes with differentially expressed genes

To verify whether the aberrant expression of genes in the SP and RP subclones can be explained by the observed copy number changes, integration of genomic (array CGH data) and transcriptomic data was performed. For subclone SP, we identified 18 genes with DNA gains being up-regulated and 31 genes with DNA losses being down-regulated, whereas for subclone RP, 44 up-regulated genes showed copy number gains and 5 genes with copy number losses were down-regulated (Table 2).

Table 2 Integration of differentially expressed genes with array CGH data
Full size table

This integrative data analysis of DNA copy number and gene expression followed by pathway enrichment analysis allowed us to identify related pathways encompassing signaling by Rho GTPases as one of the deregulated pathways in the SP subclone. At the same time we observed DNA loss and downregulation of the TIAM1 gene that belongs to the Rho GTPases signaling pathway. Yang et al. recently showed that high expression of TIAM1 is associated with poor clinical outcome in patients with HNSCC [39]. Similarly, for the RP subclone a gain and upregulation of RAD1 and RICTOR genes was observed which is in accordance with the deregulation of the homologous recombination and PI3K signaling pathways in HNSCC as previously described in [40]. For both of those clones a gain and upregulation of the TLR4 gene and the deregulation of the toll-like kinases signaling pathway including upregulation of the genes PLCG2, TAB3, RPS6KA2 has been detected. Even though contradictory reports exist, the overexpression of TLR4 and activation of related pathway has been described to promote HNSCC tumor development and to ensure tumor protection from the immune system [41].

Time-dependent gene expression in response to ionizing radiation in CAL-33 clones

It was hypothesized that CAL-33 subclones with different phenotypes of radiosensitivity also show differences in gene expression profiles in response to irradiation. Therefore, we performed differential time-course microarray analyses between irradiated (8 Gy) and sham irradiated control cells. 7299, 6980, and 8111 genes were differentially expressed in the CAL-33 parental, SP, and RP subclones, respectively (Table 3). Although the number of differentially expressed genes after irradiation with 8 Gy was comparable for all cell lines studied, approximately 50 % of the affected genes were not identical. More than 1000 genes were exclusively involved in the radiation response of each of the CAL-33 cell clones (Fig. 4). The entirety of differentially expressed genes in response to ionizing radiation is listed and compared for all CAL-33 cell lines in Additional file 7: Table S6.

Table 3 Comparison of detected and differentially expressed genes after irradiation for analyzed cell sublines
Full size table
Fig. 4
figure4

Venn diagram displaying commonly and exclusively differentially expressed genes of each CAL-33 cell line after irradiation with 8 Gy

Full size image

To identify common and separate pathways of the radiation response, pathway enrichment analyses were performed on early and late responding genes, respectively (Additional file 8: Table S7). The most interesting pathways with regard to their known role in radiation response were grouped in to major pathways (Table 4).

Table 4 Significantly enriched pathways of early and late responding genes after ionizing radiation
Full size table

It is interesting to note that all cell lines share a set of deregulated pathways (Table 4). In addition, for each of the radiosensitivity phenotypes (parental, SP, RP) specific pathways were observed. The early responses in gene expression include interferon signaling and in particular for the resistant phenotype the cellular response to hypoxia which is known for a long time to have a crucial role in radioresistance [42, 43]. Also, some other pathways from an early gene expression response that are involved in DNA repair, cell cycle, cellular response to stress and apoptosis impact on survival after irradiation and therefore contribute to radioresistance [9, 44–46].

Deregulated pathways related to late responding genes include EGFR- and PI3K/AKT signaling that were already linked to poor clinical outcome and therapy response in HNSCC [47–49] in association with ERBB2, ERBB3 and ERBB4 co-expression [50–52]. Interestingly, also involvement of ERBB2 and EERB4 receptor signaling was observed in all cell lines as a late response to ionizing radiation. Further late responding pathways include toll-like receptor cascades, interleukin signaling, NF-kB activation and interferon signaling all of which are frequently detected after treatment with ionizing radiation [53, 54]. However, the role of interleukin signaling and immune response is largely unknown so far. This also applies to cellular senescence and the impact of senescence pathways on radioresistance that were also discovered as late responding pathways in our CAL-33 subclones. Although evidence exists that senescence might be associated with the disruption of the tissue microenvironment leading to the secretion of senescence-associated pro-inflammatory factors and to the development of a pro-oncogenic environment [34] its role in radioresistance in rather unclear so far.

Gene association network reconstruction and network analysis

To go beyond pathway enrichment analyses of differentially expressed genes, gene association networks were reconstructed. Parameters of the obtained networks (provided as igraph R-objects in Additional file 9: File S1) for all three CAL-33 cell lines are presented in Table 3. Subsequently, the combined topological centrality measure was used to characterize the biological importance of genes in the reconstructed association networks. We identified nodes (genes) that are likely to control the network by combining three network centrality measures: degree, closeness and shortest path betweenness [32, 55]. The 5 % of the highest ranked genes (listed in Additional file 10: Table S8) were mapped to the Reactome pathways to further evaluate their biological roles. The top ten pathways according to the FDR values are listed in Table 5. All identified pathways are listed in Additional file 11: Table S9. Two of the detected pathways, signaling by Rho GTPases and signaling by Wnt, were in common for all three subclones exposed to irradiation (Table 5). However, most of the pathways were different for the individual subclones. For the parental cell line, we additionally detected pathways associated with cellular response to stress, signaling by EGFR, and cytochrome P450. The most important pathways in the SP subclone were associated with axon guidance, chromatin organization and semaphorin, whereas for the RP subclone the senescence-associated secretory phenotype (SASP) together with GPCR ligand binding were considered as crucial. The co-occurrence of the SASP and GPCR signaling pathways, that are known to be connected [56], indicate that the identification of those pathways might not be accidental and suggest a high importance for the radiation resistant phenotype. In future experiments more HNSCC cell lines should be analyzed.

Table 5 The top pathways after mapping of 5 % highest ranked genes from the reconstructed gene association networks to the Reactome pathways
Full size table

Expression of endogenous retrovirus and related pathways

The analysis of differentially expressed genes between the RP and SP subclones and the parental cell line revealed the ERVMER34-1 gene as differentially expressed in both clones. Approximately 8 % of the human genome consists of endogenous retroviruses (ERVs) that have been derived from exogenous retroviruses following infection and DNA integration into germ line cells [57, 58]. Although most of the ERV sequences have defective structures, some of ERV genes still have an open reading frame (ORF) and protein expression [59]. ERV genes can promote homologous and non-homologous recombination and therefore, may introduce new mutations [60, 61]. Furthermore, ERVs may lead to genome instability, and contribute to tumor initiation and progression [62]. The expression of ERV genes has been demonstrated in various cancers, including breast, ovarian, prostate and melanoma [63–66]. The differential gene expression analyses of the subclones in comparison to the parental CAL-33 cells revealed ERVMER34-1 gene as differentially expressed in both derived clones (Additional file 4: Table S3). Subsequently, we tested whether any of the retroviral genes were differentially expressed in response to ionizing radiation. Apart from ERVMER34-1, we were able to identify another endogenous retroviral gene, ERV3-1, that was temporally differentially expressed following radiation. The time dependent expression of those genes suggests ERV3-1 expression to be associated with the radiation response, whereas ERVMER34-1 expression exhibits rather random fluctuations (Additional file 12: Figure S2 and Additional file 13: Figure S3). Thus, the ERV3-1 expression clearly implies a possible association with the radiation exposure. The radiation-associated upregulation of ERV3-1 is demonstrated for all subclones starting from the second day of irradiation. To our knowledge, the expression of ERV3-1 following radiation and its influence on radiation resistance has not been addressed in detail so far. A recent study by Lee et al. [67] has demonstrated an increase in the expression of ERV3-1 (HERV-R) env related to a fractionated exposure to γ-radiation in radioresistant A549 lung cancer cells but not in less radioresistant H460 cells. The presented results raise the question whether overexpression of ERV3-1 might be involved in the radiation response of HNSCC cells. To gain knowledge about the potential gene interactions with the ERV3-1 gene we used the gene association networks reconstructed for all CAL-33 subclones and extracted the putative direct or indirect ERV3-1 interaction partners resulting in the first neighborhood genes of ERV3-1 differ between the three analyzed cell lines (Fig. 5). The largest first neighborhood gene association network can be observed for the RP subclone where 29 genes are linked to ERV3-1. For the CAL-33 parental cell line and the SP subclone the first neighborhood gene association network consist only of three (OR2A2, U2AF1L4, C11orf94) and one (FEEH2) potential association partners, respectively. The considerably larger first neighborhood of the ERV3-1 gene for the RP cells suggests a more important role of this gene for acquired radiation resistance. In addition, a Reactome pathway enrichment analysis revealed that the first neighborhood genes of the ERV3-1 gene in RP cells were associated with GPCR signaling (DRD4, OPN1MW, TBXA2R), transmembrane transport of small molecules (ATP1B2, AZGP1, SLC22A17), generic transcription pathway (ZNF419, ZNF550, ZNF782), signaling by Rho GTPases (NCKIPSD), and cell cycle (MAX). However, to our knowledge, the interaction partners of the ERV3-1 gene have not been studied in detail so far, which makes an interpretation difficult and highly speculative at this time. Also further studies have to be performed in order to validate the gene associations with ERV3-1independently.

Fig. 5
figure5

First neighborhood of the ERV3-1 gene extracted from the reconstructed gene association networks

Full size image

Conclusion

In conclusion, the present study presents comprehensive gene expression data of CAL-33 subclones of different radiosensitivity. Based on these data networks have been identified that are linked to the radiation response phenotypes. The pathways associated with the resistant phenotype are of special interest focusing on the senescence-associated secretory phenotype (SASP) together and GPCR ligand binding. Also, the radiation-associated expression of the endogenous retrovirus ERV3-1 appears highly attractive for further studies on the molecular mechanisms of acquired radioresistance.

Abbreviations

CGH, comparative genomic hybridization; CNA, copy number aberration; ERV endogenous retrovirus; FDR, false discovery rate; GAN, gene association network; HNSCC, head and neck squamous cell carcinoma; RP, resistant phenotype; SKY, spectral karyotyping; SP, sensitive phenotype

References

  1. 1.

    Sanderson RJ, Ironside JAD. Squamous cell carcinomas of the head and neck. BMJ. 2002;325(7368):822–7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  2. 2.

    Network NCC. NCCN Clinical Practice Guidelines in Oncology: head and neck cancers. Version 2. 2013.

  3. 3.

    Hashibe M, Brennan P, Chuang SC, Boccia S, Castellsague X, Chen C, Curado MP, Dal Maso L, Daudt AW, Fabianova E, et al. Interaction between tobacco and alcohol use and the risk of head and neck cancer: pooled analysis in the International Head and Neck Cancer Epidemiology Consortium. Cancer Epidemiol Biomarkers Prev. 2009;8(2):541–50.

    Article  Google Scholar 

  4. 4.

    Guha N, Boffetta P, Wünsch Filho V, Eluf Neto J, Shangina O, Zaridze D, Curado MP, Koifman S, Matos E, Menezes A, et al. Oral health and risk of squamous cell carcinoma of the head and neck and esophagus: results of two multicentric case–control studies. Am J Epidemiol. 2007;66(10):1159–73.

    Article  Google Scholar 

  5. 5.

    Pannone G, Santoro A, Papagerakis S, Lo Muzio L, De Rosa G, Bufo P. The role of human papillomavirus in the pathogenesis of head & neck squamous cell carcinoma: an overview. Infect Agents Cancer. 2011;6:4.

    Article  PubMed  PubMed Central  Google Scholar 

  6. 6.

    Seiwert TY, Salama JK, Vokes EE. The chemoradiation paradigm in head and neck cancer. Nat Clin Pract Oncol. 2007;4(3):156–71.

    CAS  Article  PubMed  Google Scholar 

  7. 7.

    Bar-Ad V, Palmer J, Yang H, Cognetti D, Curry J, Luginbuhl A, Tuluc M, Campling B, Axelrod R. Current management of locally advanced head and neck cancer: the combination of chemotherapy with locoregional treatments. Semin Oncol. 2014;41(6):798–806.

    Article  PubMed  Google Scholar 

  8. 8.

    Weichselbaum RR, Beckett MA, Schwartz JL, Dritschilo A. Radioresistant tumor cells are present in head and neck carcinomas that recur after radiotherapy. Int J Radiat Oncol Biol Phys. 1988;15(3):575–9.

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Perri F, Pacelli R, Della Vittoria Scarpati G, Cella L, Giuliano M, Caponigro F, Pepe S: Radioresistance in head and neck squamous cell carcinoma: biological bases and therapeutic implications. Head Neck. 2015;37(5):763–70.

    Article  PubMed  Google Scholar 

  10. 10.

    Gioanni J, Fischel JL, Lambert JC, Demard F, Mazeau C, Zanghellini E, Ettore F, Formento P, Chauvel P, Lalanne CM, et al. Two new human tumor cell lines derived from squamous cell carcinomas of the tongue: Establishment, characterization and response to cytotoxic treatment. Eur J Cancer Clin Oncol. 1988;24(9):1445–55.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Ernst A, Anders H, Kapfhammer H, Orth M, Hennel R, Seidl K, Winssinger N, Belka C, Unkel S, Lauber K. HSP90 inhibition as a means of radiosensitizing resistant, aggressive soft tissue sarcomas. Cancer Lett. 2015;365(2):211–22.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Braselmann H, Michna A, Hess J, Unger K. CFAssay: statistical analysis of the colony formation assay. Radiat Oncol. 2015;10:223.

    Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Kinzel L, Ernst A, Orth M, Albrecht V, Hennel R, Brix N, Frey B, Gaipl US, Zuchtriegel G, Reichel CA, et al. A novel HSP90 inhibitor with reduced hepatotoxicity synergizes with radiotherapy to induce apoptosis, abrogate clonogenic survival, and improve tumor control in models of colorectal cancer. Oncotarget. 2016. doi:10.18632/oncotarget.9774.

    PubMed  Google Scholar 

  14. 14.

    Hieber L, Huber R, Bauer V, Schäffner Q, Braselmann H, Thomas G, Bogdanova T, Zitzelsberger H. Chromosomal rearrangements in post-Chernobyl papillary thyroid carcinomas: evaluation by spectral karyotyping and automated interphase FISH. J Biomed Biotechnol. 2011;2011:693691.

    Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Neuvial P, Hupe P. MANOR: CGH Micro-Array NORmalization. R package version 1.42.0.

  16. 16.

    Neuvial P, Hupe P, Brito I, Liva S, Manie E, Brennetot C, Radvanyi F, Aurias A, Barillot E. Spatial normalization of array-CGH data. BMC Bioinformatics. 2006;7(264):1–20.

    Google Scholar 

  17. 17.

    van de Wiel MA, Brosens R, Eilers PHC, Kumps C, Meijer GA, Menten B, Sistermans E, Speleman F, Timmerman ME, Ylstra B. Smoothing waves in array CGH tumor profiles. Bioinformatics. 2009;25(9):1099–104.

    Article  PubMed  Google Scholar 

  18. 18.

    Seshan VE, Olshen A. DNAcopy: DNA copy number data analysis. R package version 1.44.0.

  19. 19.

    Olshen AB, Venkatraman ES, Lucito R, Wigler M. Circular binary segmentation for the analysis of array-based DNA copy number data. Biostatistics. 2004;5(4):557–72.

    Article  PubMed  Google Scholar 

  20. 20.

    van de Wiel MA, Kim KI, Vosse SJ, van Wieringen WN, Wilting SM, Ylstra B. CGHcall: calling aberrations for array CGH tumor profiles. Bioinformatics. 2007;23(7):892–4.

    Article  PubMed  Google Scholar 

  21. 21.

    van de Wiel MA, van Wieringen WN. CGHregions: dimension reduction for array CGH data with minimal information loss. Cancer Informatics. 2007;3:55–63.

    PubMed  PubMed Central  Google Scholar 

  22. 22.

    R Core Team: R: A language and environment for statistical computing 2013.

  23. 23.

    Smyth GK. Limma: linear models for microarray data. In: Bioinformatics and Computational Biology Solutions Using {R} and Bioconductor. New York: Springer; 2005. p. 397–420.

    Google Scholar 

  24. 24.

    Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, Ellis B, Gautier L, Ge Y, Gentry J, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol. 2004;5(10):R80.

    Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Lopez-Romero P: Agi4x44PreProcess: PreProcessing of Agilent 4x44 array data. R package version 1.16.0.

  26. 26.

    Li Q, Birkbak NJ, Gyorffy B, Szallasi Z, Eklund AC. Jetset: selecting the optimal microarray probe set to represent a gene. BMC Bioinformatics. 2011;12:474.

    Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Michna A. splineTimeR: Time-course differential gene expression data analysis using spline regression models followed by gene association network reconstruction. R package version 0995 2015.

  28. 28.

    Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Royal Stat Soc. 1995;57:289–300.

    Google Scholar 

  29. 29.

    Opgen-Rhein R, Strimmer K. Using regularized dynamic correlation to infer gene dependency networks from time-series microarray data. The 4th International Workshop on Computational Systems Biology, WCSB 2006.

  30. 30.

    Opgen-Rhein R, Strimmer K. From correlation to causation networks: a simple approximate learning algorithm and its application to high-dimensional plant gene expression data. BMC Syst Biol. 2007;1:37.

    Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Koschützki D, Schreiber F. Centrality analysis methods for biological networks and their application to gene regulatory networks. Gene Regulation Syst Biol. 2008;2:193–201.

    Google Scholar 

  32. 32.

    Abbasi A, Hossain L. Hybrid Centrality Measures for Binary and Weighted Networks. In: Menezes R, Evsukoff A, González MC, editors. Complex Networks. Volume 424, edn. Berlin: Springer; 2013. p. 1–7.

    Google Scholar 

  33. 33.

    Matthews L, Gopinath G, Gillespie M, Caudy M, Croft D, de Bono B, Garapati P, Hemish J, Hermjakob H, Jassal B, et al. Reactome knowledgebase of human biological pathways and processes. Nucleic Acids Res. 2009;37(Database issue):D619–22.

    CAS  Article  PubMed  Google Scholar 

  34. 34.

    Sabin RJ, Anderson RM. Cellular Senescence - its role in cancer and the response to ionizing radiation. Genome Integr. 2011;2(1):7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Wang J, Hu L, Gupta N, Shamseldin T, Ozawa T, Klem J, Cardell M, Deen DF. Induction and characterization of human glioma clones with different radiosensitivities. Neoplasia. 1999;1(2):138–44.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Bentzen SM, Atasoy BM, Daley FM, Dische S, Richman PI, Saunders MI, Trott KR, Wilson GD. Epidermal growth factor receptor expression in pretreatment biopsies from head and neck squamous cell carcinoma as a predictive factor for a benefit from accelerated radiation therapy in a randomized controlled trial. J Clin Oncol. 2005;23(24):5560.

    CAS  Article  PubMed  Google Scholar 

  37. 37.

    Kalyankrishna S, Grandis JR. Epidermal growth factor receptor biology in head and neck cancer. J Clin Oncol. 2006;24(17):2666–72.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Stadler ME, Patel MR, Couch ME, Hayes DN. Molecular biology of head and neck cancer: risks and pathways. Hematol Oncol Clin North Am. 2008;22(6):1099–124.

    Article  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Yang H, Cai YC, Cao Y, Song M, An X, Xia Y, Wei J, Jiang WQ, Shi YX. The prognostic value of Tiam1 protein expression in head and neck squamous cell carcinoma: a retrospective study. Chin J Cancer. 2015;34(12):614–21.

    CAS  PubMed  Google Scholar 

  40. 40.

    Lui VW, Hedberg ML, Li H, Vangara BS, Pendleton K, Zeng Y, Lu Y, Zhang Q, Du Y, Gilbert BR, et al. Frequent mutation of the PI3K pathway in head and neck cancer defines predictive biomarkers. Cancer Discov. 2013;3(7):761–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Szczepanski MJ, Czystowska M, Szajnik M, Harasymczuk M, Boyiadzis M, Kruk-Zagajewska A, Szyfter W, Zeromski J, Whiteside TL. Triggering of Toll-like receptor 4 expressed on human head and neck squamous cell carcinoma promotes tumor development and protects the tumor from immune attack. Cancer Res. 2009;69(7):3105–13.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Overgaard J. Hypoxic modification of radiotherapy in squamous cell carcinoma of the head and neck - a systematic review and meta-analysis. Radiother Oncol. 2011;100(1):22–32.

    Article  PubMed  Google Scholar 

  43. 43.

    Li JZ, Gao W, Chan JY, Ho WK, Wong TS. Hypoxia in head and neck squamous cell carcinoma. ISRN Otolaryngol. 2012;2012:708974.

    PubMed  PubMed Central  Google Scholar 

  44. 44.

    Meijer TW, Kaanders JH, Span PN, Bussink J. Targeting hypoxia, HIF-1, and tumor glucose metabolism to improve radiotherapy efficacy. Clin Cancer Res. 2012;18(20):5585–94.

    CAS  Article  PubMed  Google Scholar 

  45. 45.

    Parkin DM, Bray F, Ferlay J, Pisani P. Global cancer statistics, 2002. CA Cancer J Clin. 2005;55(2):74–108.

    Article  PubMed  Google Scholar 

  46. 46.

    Vokes EE, Weichselbaum RR, Lippman SM, Hong WK. Head and neck cancer. N Engl J Med. 1993;328(3):184–94.

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Sharafinski ME, Ferris RL, Ferrone S, Grandis JR. Epidermal growth factor receptor targeted therapy of squamous cell carcinoma of the head and neck. Head Neck. 2010;32(10):1412–21.

    Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Ang KK, Berkey BA, Tu X, Zhang HZ, Katz R, Hammond EH, Fu KK, Milas L. Impact of epidermal growth factor receptor expression on survival and pattern of relapse in patients with advanced head and neck carcinoma. Cancer Res. 2002;62(24):7350–6.

    CAS  PubMed  Google Scholar 

  49. 49.

    Nijkamp MM, Span PN, Bussink J, Kaanders JH. Interaction of EGFR with the tumour microenvironment: implications for radiation treatment. Radiother Oncol. 2013;108(1):17–23.

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Ibrahim SO, Vasstrand EN, Liavaag PG, Johannessen AC, Lillehaug JR. Expression of c-erbB proto-oncogene family members in squamous cell carcinoma of the head and neck. Anticancer Res. 1997;17(6D):4539–46.

    CAS  PubMed  Google Scholar 

  51. 51.

    Xia W, Lau YK, Zhang HZ, Liu AR, Li L, Kiyokawa N, Clayman GL, Katz RL, Hung MC. Strong correlation between c-erbB-2 overexpression and overall survival of patients with oral squamous cell carcinoma. Clin Cancer Res. 1997;3(1):3–9.

    CAS  PubMed  Google Scholar 

  52. 52.

    Shintani S, Funayama T, Yoshihama Y, Alcalde RE, Matsumura T. Prognostic significance of ERBB3 overexpression in oral squamous cell carcinoma. Cancer Lett. 1995;95(1–2):79–83.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Di Maggio FM, Minafra L, Forte GI, Cammarata FP, Lio D, Messa C, Gilardi MC, Bravatà V. Portrait of inflammatory response to ionizing radiation treatment. J Inflamm (Lond). 2015;12:14.

    Article  Google Scholar 

  54. 54.

    Roses RE, Xu M, Koski GK, Czerniecki BJ. Radiation therapy and Toll-like receptor signaling: implications for the treatment of cancer. Oncogene. 2008;27(2):200–7.

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Koschützki D. Network Centralities. In: Junker BH, Schreiber F, editors. Analysis of Biological Networks. Hoboken: Wiley; 2007. p. 65–84.

    Google Scholar 

  56. 56.

    Acosta JC, O'Loghlen A, Banito A, Guijarro MV, Augert A, Raguz S, Fumagalli M, Da Costa M, Brown C, Popov N, et al. Chemokine signaling via the CXCR2 receptor reinforces senescence. Cell. 2008;133(6):1006–18.

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Lander ES, Linton LM, Birren B, Nusbaum C, Zody MC, Baldwin J, Devon K, Dewar K, Doyle M, FitzHugh W, et al. Initial sequencing and analysis of the human genome. Nature. 2001;409(6822):860–921.

    CAS  Article  PubMed  Google Scholar 

  58. 58.

    Löwer R, Löwer J, Kurth R. The viruses in all of us: characteristics and biological significance of human endogenous retrovirus sequences. Proc Natl Acad Sci U S A. 1996;93(11):5177–84.

    Article  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Villesen P, Aagaard L, Wiuf C, Pedersen FS. Identification of endogenous retroviral reading frames in the human genome. Retrovirology. 2004;1:32.

    Article  PubMed  PubMed Central  Google Scholar 

  60. 60.

    Hughes JF, Coffin JM. Evidence for genomic rearrangements mediated by human endogenous retroviruses during primate evolution. Nat Genet. 2001;29(4):487–9.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Whitelaw E, Martin DI. Retrotransposons as epigenetic mediators of phenotypic variation in mammals. Nat Genet. 2001;27(4):361–5.

    CAS  Article  PubMed  Google Scholar 

  62. 62.

    Stoye JP. Studies of endogenous retroviruses reveal a continuing evolutionary saga. Nat Rev Microbiol. 2012;10(6):395–406.

    CAS  PubMed  Google Scholar 

  63. 63.

    Wang-Johanning F, Rycaj K, Plummer JB, Li M, Yin B, Frerich K, Garza JG, Shen J, Lin K, Yan P, et al. Immunotherapeutic potential of anti-human endogenous retrovirus-K envelope protein antibodies in targeting breast tumors. J Natl Cancer Inst. 2012;104(3):189–210.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Wang-Johanning F, Liu J, Rycaj K, Huang M, Tsai K, Rosen DG, Chen DT, Lu DW, Barnhart KF, Johanning GL. Expression of multiple human endogenous retrovirus surface envelope proteins in ovarian cancer. Int J Cancer. 2007;120(1):81–90.

    CAS  Article  PubMed  Google Scholar 

  65. 65.

    Wang-Johanning F, Frost AR, Jian B, Azerou R, Lu DW, Chen DT, Johanning GL. Detecting the expression of human endogenous retrovirus E envelope transcripts in human prostate adenocarcinoma. Cancer. 2003;98(1):187–97.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Li Z, Sheng T, Wan X, Liu T, Wu H, Dong J. Expression of HERV-K correlates with status of MEK-ERK and p16INK4A-CDK4 pathways in melanoma cells. Cancer Invest. 2010;28(10):1031–7.

    CAS  Article  PubMed  Google Scholar 

  67. 67.

    Lee JR, Jung YD, Kim YH, Park SJ, Huh JW, Kim HS. Effects of HERV-R env knockdown in combination with ionizing radiation on apoptosis-related gene expression in A549 lung cancer cells. Biochem Physiol. 2016;5:1.

    Google Scholar 

Download references

Acknowledgement

We thank Aaron Selmaier, Isabella Zagorski and Laura Dajka from the Research Unit Radiation Cytogenetics for their excellent technical support.

Funding

This study was supported by the German Federal Ministry of Education and Research (ZiSS - 02NUK024B and 02NUK024C) and the Clinical Cooperation Group “Personalized Radiotherapy in Head and Neck Cancer”.

Availability of data and material

Gene expression microarray data will be available at ArrayExpress.

Authors’ contributions

AM: gene expression analyses, bioinformatics and biostatistics analysis, manuscript draft; US: generation and characterization of subclones, sample preparation for gene expression analyses; MS: bioinformatics and biostatistics analysis; HZ: spectral karyotyping; support for study design, critical revision of the manuscript; KL: conceived generation of subclones; support for experimental design/concept, critical revision of the manuscript; KU: support for study design, experimental design/concept, bioinformatics and biostatistics analysis; JH: array CGH analyses; support for study design, experimental design/concept; draft and critical revision of manuscript; All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Disclosures

Not applicable.

Author information

Author notes

    Affiliations

    1. Research Unit Radiation Cytogenetics, Helmholtz Zentrum München, German Research Center for Environmental Health GmbH, 85764, Neuherberg, Germany

      Agata Michna, Martin Selmansberger, Horst Zitzelsberger, Kristian Unger & Julia Hess

    2. Department of Radiotherapy and Radiation Oncology, Ludwig-Maximilians-University, 81377, Munich, Germany

      Ulrike Schötz & Kirsten Lauber

    3. Clinical Cooperation Group “Personalized Radiotherapy in Head and Neck Cancer”, Helmholtz Zentrum München, 85764, Neuherberg, Germany

      Horst Zitzelsberger, Kirsten Lauber, Kristian Unger & Julia Hess

    Authors
    1. Agata Michna
      View author publications

      You can also search for this author in PubMed Google Scholar

    2. Ulrike Schötz
      View author publications

      You can also search for this author in PubMed Google Scholar

    3. Martin Selmansberger
      View author publications

      You can also search for this author in PubMed Google Scholar

    4. Horst Zitzelsberger
      View author publications

      You can also search for this author in PubMed Google Scholar

    5. Kirsten Lauber
      View author publications

      You can also search for this author in PubMed Google Scholar

    6. Kristian Unger
      View author publications

      You can also search for this author in PubMed Google Scholar

    7. Julia Hess
      View author publications

      You can also search for this author in PubMed Google Scholar

    Corresponding author

    Correspondence to Julia Hess.

    Additional files

    Additional file 1: Table S1.

    Regions of copy number changes identified by array CGH analysis. Array CGH analysis identified 173 regions with copy number changed. 78, 111 and 132 regions were identified as gains or losses for the parental cell line, the SP and RP subclones, respectively. Regions with normal copy number status are labeled with “0”, DNA-gains are marked with “1” and DNA-losses with “-1”. Regions with different copy number variations between clones SP and RP and the parental cell line are marked as “diff”. (XLSX 56 kb)

    Additional file 2: Table S2.

    Overview on the cytogenetic aberrations in CAL-33 cells. Clonal structural and numerical chromosomal aberrations detected by SKY and copy number changes detected by array CGH are listed. (DOCX 94 kb)

    Additional file 3: Figure S1.

    Analyses of proliferation rates and cell cycle distribution. Proliferation curves of the CAL-33 parental cell line and the subclones SP and RP are shown in panel (A). Doubling times were determined in the exponential growth phase by semi-log regression. Cell cycle analyses were performed under exponential growth conditions (B). The dynamics of G2-arrest upon irradiation at 4 Gy are shown in figure panel (C). (PDF 318 kb)

    Additional file 4: Table S3.

    Differentially expressed genes and the corresponding fold-changes. (XLSX 117 kb)

    Additional file 5: Table S4.

    Technical validation of microarray data. Eight of the differentially expressed genes detected with the microarray technology were arbitrarily chosen for technical validation by qRT-PCR. Correlation analysis results (Spearman’s rho coefficient) and fold-changes (microarray and qRT-PCR) are shown. Asterisks (*) indicate fold-change values for genes detected as differentially expressed by microarray analysis. (DOCX 52 kb)

    Additional file 6: Table S5.

    Pathway enrichment analysis results of differentially expressed genes in clones SP and RP in comparison to parental CAL-33 cells. (XLSX 46 kb)

    Additional file 7: Table S6.

    Differentially expressed genes in response to irradiation with 8 Gy for all CAL-33 cell lines. Differentially expressed genes are listed together with the corresponding fold-changes. Additionally, commonly and exclusively differentially expressed genes have been included. (XLSX 2575 kb)

    Additional file 8: Table S7.

    Significantly enriched pathways of early and late responding genes after irradiation with 8 Gy. Pathways with FDR < 0.1 were considered statistically significant. (XLSX 102 kb)

    Additional file 9: File S1.

    Reconstructed gene association networks. All obtained gene association networks are provided as R-objects of type igraph. (RDATA 32863 kb)

    Additional file 10: Table S8.

    5 % most important genes identified by centrality measures. The 5 % highest ranked genes from the reconstructed gene association networks together with corresponding ranks of centrality measures are shown. (XLSX 88 kb)

    Additional file 11: Table S9.

    Involved pathways after mapping of 5 % highest ranked genes from the gene association networks to the Reactome pathways. The pathways together with the mapped most important (5 % highest ranked) genes are inluded. (XLSX 92 kb)

    Additional file 12: Figure S2.

    Temporal expression and corresponding fold-changes of gene ERVMER34-1 after irradiation with 8 Gy. (PDF 535 kb)

    Additional file 13: Figure S3.

    Temporal expression and corresponding fold-changes of gene ERV3-1 after irradiation with 8 Gy. (PDF 531 kb)

    Rights and permissions

    Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

    Reprints and Permissions

    About this article

    Verify currency and authenticity via CrossMark

    Cite this article

    Michna, A., Schötz, U., Selmansberger, M. et al. Transcriptomic analyses of the radiation response in head and neck squamous cell carcinoma subclones with different radiation sensitivity: time-course gene expression profiles and gene association networks. Radiat Oncol 11, 94 (2016). https://doi.org/10.1186/s13014-016-0672-0

    Download citation

    Keywords

    • Radioresistance
    • HNSCC
    • Head and neck cancer
    • Time-course gene expression
    • Gene association network
    • Signaling pathway
    • Differentially expressed genes
    • Endogenous retrovirus

    Radiation Oncology

    ISSN: 1748-717X

    Contact us