Peripheral blood samples were collected from healthy donors in 9 ml EDTA vacutainer tubes. Blood collection was approved by the local SCK CEN Ethics Committee and was carried out in accordance with the ethical standards of the Helsinki Declaration of 1975, as revised in 2000. Prior to blood donation all the donors involved in the present study signed an informed consent form. Within 30-60 min of blood drawing, PBMCs were isolated by centrifugation on Histopaque-1077 density gradient (Sigma-Aldrich, Bornem, Belgium) according to the manufacturer’s instructions. Isolated cells were suspended at a density of 10⁶ cells/ml in LGM-3 culture medium (Lonza, Walkersville, MD, USA) and were allowed to equilibrate to culture conditions at 37°C in a humidified 5% CO₂ atmosphere.

Three human lymphoblastoid cell lines TK6, WTK1 and NH32 were used in this study. TK6 cell line was purchased from DSMZ, Leibniz-Institute, Braunschweig, Germany, while WTK1 and NH32 cells were a generous gift from the laboratory of Prof Schwartz, University of Washington. All three cell lines have the same progenitor, WIL2 cell line (29). The TK6 cells express wild-type p53, while WTK1 cells overexpress a mutant form of p53 due to methionine to isoleucine substitution at codon 237 and no wild-type p53 protein (30). The NH32 cell line is genetically homologous to TK6, but its TP53 gene has been genetically inactivated by a homozygous knockout (31). Cell lines were maintained as exponentially growing stationary cultures in RPMI-1640 medium supplemented with 10% heat-inactivated foetal bovine serum at 37°C in 5% CO₂ at densities of 4–10 × 10⁵ cells/ml. Cell line irradiations were performed in three replicates (1 × 10⁶ cells per sample) as described below. Following irradiations cells were incubated at standard conditions for 24 h prior to RNA extraction.

X-ray, carbon and iron ion irradiations were performed independently, on different days. X-ray irradiation experiments were performed at the irradiation facility at SCK CEN, Belgian Nuclear Research Centre, Mol, Belgium (for microarrays) and at GSI, Darmstadt, Germany (for qRT-PCR validation and γH2AX staining). At SCK CEN, PBMCs were exposed to 1.00 Gy of X-rays, using a Pantak HF420RX machine (250 kV, 15 mA, dose rate of 0.26 Gy/min). Cells were irradiated “free-in-air” at 21°C in a horizontal position with single doses of 0.1 and 1.0 Gy of X-rays from a Pantak HF420 RX generator at an air kerma (K_air) rate of 0.26 Gy/min or were sham-irradiated. The beam quality can be approximated to H-250 (ISO4037): 250 kV, 15 mA, 1.2 mm Al equivalent inherent filtration and 1 mm Cu additional filtration. First HVL was 2.43 mm Cu and the second HVL was 3.52 mm Cu. The K_air at the reference position was measured using a NE2571 ionization chamber (SN309) connected to a Farmer 2500 electrometer. The chamber, together with the electrometer, was calibrated in terms of K_air and the traceability to the international standards was assured. The reference point of the ionization chamber was placed at the same distance with the reference position of the samples. The ionization chamber was always placed in the beam, next to the samples, for a precise measurement of the time integrated K_air. The stability of the X-ray generator during the irradiation was verified in this way using a monitor chamber. Cell line samples were exposed to 1.00 and 3.00 Gy at SCK CEN following the same procedures. For samples that were irradiated at GSI, freshly isolated PBMCs were transported for 4 h by car to GSI using a transportable incubator. X-ray exposures at GSI were performed using an IV320-13 X-ray tube (250 keV, 16 mA, dose rate of 0.5 Gy/min; Seifert, Germany) at 0.25 and 1.00 Gy. Irradiation with heavy ions was performed at GSI on the heavy-ion synchrotron SIS. Carbon ion exposure (0.25 and 1.00 Gy) was performed in the middle of a 25 mm spread-out Bragg peak (center depth 42.5 mm, realized with a PMMA bolus), obtained by active energy variation of the beam in the range of 114.6 – 158.4 MeV/nucleon. Accordingly, the dose averaged LET at the proximal and distal part of the samples (5-ml plastic tube, inside diameter 10 mm) was 60-80 keV/µm. Irradiation with iron ions (0.25 and 1.00 Gy) was performed with a monoenergetic beam (1 GeV/nucleon; LET 155 keV/µm). The beam monitor calibration was performed according to the procedure described by Luoni et al. (32). For verification of the applied dose, additional absolute dosimetry was performed by measuring the absorbed dose to water using a PTW TM30013 Farmer ionization chamber positioned at the sample depth. The dose-averaged LET was calculated using the TRiP98 treatment planning system (33). Sham-irradiated samples were always subjected to the same procedures as the irradiated ones, except for the radiation exposure itself. After in vitro irradiation, cells were incubated at 37°C in a humidified 5% CO₂ atmosphere for the indicated time until further processing.

For RNA isolation, a combination of the TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA) extraction method and the clean-up on Qiagen RNeasy columns (Qiagen, Venlo, The Netherlands) was used. Briefly, 5 x 10⁶ cells were lysed in 1 ml of TRIzol^® reagent and further processed following the manufacturer’s recommendations. Following the RNA precipitation with isopropanol, the obtained pellet was re-suspended in 1 ml of ethanol and transferred to the RNeasy column. Further purification was done according to the manufacturer’s instructions. RNA concentration was measured on a NanoDrop-2000 spectrophotometer (Thermo Scientific, Erembodegem, Belgium) and the quality of total RNA samples was assessed using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Only samples with an RNA integrity number >7 were considered as suitable for further microarray hybridization. RNA extraction from cell line samples was performed following the same procedures.

Gene expression profiling was performed using the GeneChip^® Human Gene 1.0 ST Array (Affymetrix, Santa Clara, CA, USA) which interrogates 28,536 well-annotated genes with 253,002 distinct probe sets, allowing expression analysis at both gene and exon level. Ten µg of cRNA, synthesized and purified from 0.25 µg of total RNA using the Ambion^® WT Expression kit (Ambion, USA) was used for cDNA synthesis, followed by cDNA fragmentation and labeling with the GeneChip^® Terminal Labeling kit (Affymetrix). Fragmented and labeled cDNA was hybridized to Human Gene 1.0 ST arrays (Affymetrix) using the GeneChip^® Hybridization, Wash and Stain kit (Affymetrix) (hybridization module) and hybridization controls (Affymetrix) with rotation at 45°C for 16 hours. After hybridization, arrays were washed and stained using the GeneChip^® Hybridization, Wash and Stain kit (stain module) after which the arrays were immediately scanned using an Affymetrix GeneChip^® Scanner. Raw data of X-ray and heavy ion experiments have been submitted to the ArrayExpress database under accession numbers E-MTAB-3463 and E-MTAB-5761, respectively.

The obtained microarray data were imported into Partek Genomics Suite, version 6.6 (Partek Inc., St Louis, MO, USA) as.CEL-files. The probe summarization and probe set normalization were done using the Robust Multichip Analysis (RMA) algorithm (34) which includes background correction, quantile normalization and log₂ transformation. Microarray data were analyzed using ANOVA with dose, donor and time point (whenever applicable) as factors. To correct for multiple testing, we used the false discovery rate (FDR) as described by Benjamini and Hochberg (35) to adjust p-values (FDR < 0.05). Genes were considered significantly differentially expressed between the two groups if adjusted p-values were < 0.05. In some cases, a more stringent additional cut-off of fold-change ≥|2| was used, as explained in the text.

We also performed Alternative Splicing ANOVA in Partek to detect genes which were alternatively spliced in response to different radiation types. An FDR-corrected p-value of < 0.05 was considered significant for alternative splicing events. To further reduce the number of false positives, the probe sets with log₂ values below the noise level in all samples were excluded from analysis, except for the cases where there was a significant difference in expression of a single exon between the groups (p < 0.05).

The Venny on-line tool (36) was used to compare gene lists and create Venn diagrams: http://bioinfogp.cnb.csic.es/tools/venny/index.html.

The following genes were selected for qRT-PCR validation: PCNA, GADD45A, RPS27L, ASTN2, NDUFAF6, FDXR, MAMDC4. The same RNA samples as those used for microarray hybridization plus two additional samples irradiated on the same day but not selected for microarray hybridization (n=6), were used for cDNA synthesis with the GoScript™ Reverse Transcription System (Promega, Leiden, The Netherlands) with random hexamer primers. For each gene, qRT-PCR reactions were run in duplicate using the MESA GREEN^® qRT-PCR kit (Eurogentec, Seraing, Belgium) on an Applied Biosystems^® 7500 Real-Time PCR instrument following the manufacturer’s instructions. To determine the efficiency and specificity of the designed primers, a standard curve experiment with melt curve was run for every primer pair. qRT-PCR data were analyzed by 7500 Software v2.0.6 and Microsoft Excel using the Pfaffl method (37). The relative amount of transcript of the selected genes was normalized to PGK1 and HPRT1 using the geometric mean of these reference genes (38). cDNA synthesis and qRT-PCR of the cell line samples was performed following identical procedures for the following genes: EDA2R, NDUFAF6, PTPN14 and VWCE. All primer sequences can be found in Supplementary Table 1 .

The RRHO algorithm allows for the comparison of two microarray datasets. Each dataset is processed as a ranked list based on expression differences between two classes of samples (0 Gy and 1 Gy, in our case). RRHO analysis (39) was performed using the online tool (http://systems.crump.ucla.edu/rankrank/index.php). As this algorithm only allows the comparison of two gene lists at a time, the following comparisons were performed: X-rays vs carbon ions, X-rays vs iron ions and carbon ions vs iron ions using a step size of 100.

Transcription factor and Gene Ontology terms enrichment analysis was performed using the Enrichr online tool (http://amp.pharm.mssm.edu/Enrichr/) (40, 41) which uses input gene lists to calculate enrichment of genes based on different databases of chromatin immunoprecipitation experiments and Ontologies. We used the “ENCODE and ChEA Consensus TFs from ChIP-X” and “GO Biological Process 2015” databases to calculate enrichment of transcription factor binding and biological processes, respectively, in significantly differentially expressed genes. The same tool was used to calculate enrichment of GO biological processes terms in significantly alternatively spliced genes.

Gene set enrichment analysis (42) was performed using default settings: the significance of the normalized enrichment score for each gene set was assessed through 1000 gene set permutations. Gene sets with an FDR q-value < 0.25 were considered significant, as suggested by the GSEA tutorial. For each radiation type, 1-Gy and sham-irradiated samples analyzed at 8 h after exposure were used for comparison. To have a general view of response to each radiation type, Hallmark Gene Sets collection of the Molecular Signatures Database was used. This collection consists of 50 gene sets representing specific well-defined biological states and processes, which helps to reduce noise and redundancy in different available databases and provides a better delineated biological space for GSEA.

PBMCs from four donors were fixed in 4% paraformaldehyde (Merck KGaA, Darmstadt, Germany) at several time points (see Figure 1 ). Following the fixation step, cells were cytospun on glass slides using Shandon™ EZ Double Cytofunnels™ and permeabilized with 0.25% Triton X 100 (Sigma Aldrich, Belgium) for 5 min, blocked with 3% bovine serum albumin (Sigma Aldrich, Belgium) for 30 min and incubated overnight at room temperature with monoclonal mouse anti-γH2AX (phospho S139) antibody [3F2] (ab22551, Abcam, Cambridge, MA, USA) at 4°C. Cells were then incubated for 1 h with polyclonal goat anti-mouse secondary antibody coupled to FITC (F2012, Sigma Aldrich, Belgium) at 37°C and then mounted in Vectashield mounting medium containing DAPI (Vector Laboratories, Burlingame, CA, USA). Between each of the previous steps, the slides were washed with phosphate-buffered saline.

An automated inverted fluorescence microscope (Eclipse Ti, Nikon, Tokyo, Japan), equipped with a motorized XYZ stage was used for image acquisition of immunostained slides. Images were acquired with a 40X Plan Fluor oil objective (Numerical aperture 1.3) and an Andor iXon3 camera (Andor Technology, South Windsor, CT, USA), providing images with a lateral resolution of 0.2 µm/pixel. For each sample, 25 fields were acquired on 7 Z-planes (separated by 1 μm). The obtained images were analyzed with the CellBlocks.ijm script (43), written for FIJI image analysis freeware (44), essentially as described before (45). In brief, the image analysis workflow starts by segmenting each nucleus in the DAPI channel, using an automatic thresholding algorithm, after noise reduction and flat field correction. Subsequently, γH2AX foci signals are selectively enhanced by means of a multiscale Laplacian and segmented by means of automatic thresholding. Within each nucleus, the intensity of the γH2AX channel is measured along with the number of γH2AX foci and the foci occupancy, i.e., the total projected area of the nucleus that is occupied by spots (total spot area divided by the nucleus area). On average, 500 nuclei were analyzed per sample.

To compare the effects of high- and low-LET radiation exposure on gene expression in human PBMCs, microarray analysis was performed at 8 h after exposure to 1 Gy of X-rays, carbon ions or iron ions. This time point was chosen because we observed a prominent gene expression response after 8 h in our previous studies (46, 47), moreover we found that prolonged culturing times trigger activation of apoptosis-related genes also in control samples, which might complicate data interpretation (48). We also observed dose-dependent gene expression up-regulation for doses ranging from 0.025 to 2.00 Gy, however, following 2.00 Gy exposure a clear saturation of the effect was observed (48), therefore in this study we opted for using 1 Gy for all PBMCs exposures. Following the exposure to X-rays, carbon ions or iron ions, 69, 95 and 78 genes, respectively, were detected as differentially expressed (FDR-corrected p < 0.05) compared to control samples ( Figures 2A – E and Supplementary Tables 2 – 4 ). The majority of these genes were induced after irradiation ( Figures 2A – E ), including 30 genes that were differentially expressed in response to all radiation types. Of these, 14 genes were up-regulated more than 2-fold ( Figure 2E ). The lists of genes differentially expressed exclusively following the exposure to X-rays, carbon or iron ions can be found in Supplementary Tables 2 – 4 , respectively.

When comparing two independent high-throughput gene expression experiments with different sample numbers, threshold-free methods outperform threshold-based ones in providing reliable results (39). Thus, to obtain a better impression of the similarity in gene expression after exposure to different radiation types, the Rank-rank Hypergeometric Overlap algorithm was used. This revealed a very significant degree of overlap among the genes up-regulated in both conditions (shown in the top right corner of the heatmap) for the comparisons between X-rays and carbon ions (47 overlapping genes, lowest p≈10^-56) ( Figure 2F ) as well as X-rays and iron ions (59 overlapping genes, lowest p≈10^-65) ( Figure 2G ). However, for the comparison between the two high-LET ions ( Figure 2H ) the degree of overlap was much more significant (1715 overlapping genes, lowest p≈10^-146). Together, our data show that irrespective of the radiation type, the majority of the affected genes are up-regulated after exposure, and that the identity of these genes is highly similar, although some radiation type-specific genes do seem to exist.

Transcription factor enrichment analysis suggested that, for all radiation types, the affected genes were transcriptionally regulated by p53 ( Figure 3 , left panel), and they were enriched in functions related to canonical p53-dependent pathways such as response to (ultra violet) radiation, negative regulation of the cell cycle, DNA repair and apoptosis ( Figure 3 , right panel).

The list of p53 targets is constantly growing and it remains probably the most studied transcription factor. There are at least 350 confirmed p53 targets and over 3500 potential targets (49). To confirm the p53-dependent activation of some of the significantly up-regulated genes, three lymphoblastoid cell lines with different p53 status were irradiated with different doses of X-rays (0 Gy, 1 Gy and 3 Gy). qRT-PCR was performed to measure gene expression of selected hits EDA2R, NDUFAF6, PTPN14 and VWCE, all of which were induced after X-ray, carbon and iron irradiation. EDA2R is a well-known bona fide p53 target (50), while PTPN14 has been recently identified as such (51). VWCE is often found in genome-wide data sets of activated p53 targets (49). NDUFAF6 has not yet been validated as a direct p53 target gene, but was identified to be radiation-responsive in our previous study (47). In TK6 cells, which have wild-type p53, we observed a significant, dose-dependent induction of all four genes. In contrast, none of these genes were induced after irradiation of WTK1 (p53 mutated) or NH32 (p53 null) cells ( Figure S1 ).

Next, to detect modest but coordinated changes, we performed GSEA (52). A classical DNA damage response, with p53-pathway, apoptosis and DNA damage repair-related gene sets being very significantly enriched in up-regulated genes was observed after exposure to 1 Gy of all radiation types ( Figures 4 and 5A ). Interestingly, especially after exposure to heavy ions, also several immune response and inflammation-related gene sets were identified as significantly enriched in irradiated samples ( Figures 4 and 5B, C ). For instance, genes related to the inflammatory response showed no preferential enrichment in either sham- or X-irradiated PBMCs. In contrast, exposure to heavy ions, especially iron ions, resulted in a significant up-regulation of these genes ( Figure 5B ). Similarly, the radiation effect on genes involved in TNFα signaling via NF-κB was more pronounced after heavy ion irradiation as compared to X-irradiation ( Figure 5C ). Together, these results corroborate the observation that exposure of PBMCs to heavy ion irradiation induces similar pathways but with more pronounced changes in gene expression when compared to X-rays, while certain pathways, especially those related to inflammation are particularly triggered by heavy ions.

Seven genes (PCNA, GADD45A, RPS27L, ASTN2, NDUFAF6, FDXR, MAMDC4) were selected for qRT-PCR validation of microarray data and temporal follow-up. The rationale of gene selection for qRT-PCR validation is explained below. PCNA, FDXR, GADD45A and RPS27L were significantly up-regulated 8 h after exposure to all radiation types, while ASTN2 showed significant up-regulation after exposure to X-rays and iron ions, but not to carbon ions according to the microarray data. NDUFAF6 and MAMDC4 were alternatively spliced (see next section for detailed description) in response to exposure to all radiation types. To obtain better insight in the dose- and time-dependence of expression of these genes, a lower dose (0.25 Gy) as well as an additional time point (24 h) were included ( Figure 6 ). All selected genes showed significant up-regulation at 8 h after exposure to 1 Gy of all radiation types. In many cases their induction was significant even at a lower dose of 0.25 Gy. FDXR showed the highest degree of induction following the exposure to all radiation types. The expression of ASTN2, in contrast to microarray results, was also significantly induced by carbon ions. Importantly, while all genes, except MAMDC4, reduced in expression with time in X-irradiated cells, their up-regulation was in general retained, or often even further induced in cells exposed to heavy ions, especially in the case of carbon ions.

Exposure to low-LET radiation not only changes gene expression but also triggers the production of alternative transcripts (due to alternative splicing or transcription) (47, 53–55). A core signature of genes that become alternatively spliced in response to all radiation types was identified ( Figure 7A ). The majority of these genes were also differentially expressed at the gene level (36 out of 46), which aligned well with our previous results (47). More overlap was observed between the iron- and carbon-ion irradiated groups – 33% of the genes were in common, while between the X-ray irradiated cells and heavy-ion irradiated cells the overlap was only about 15%. We also compared the number of differentially expressed exons between different radiation types to assess the levels of induction of transcript variants. Exposure to 1 Gy of X-rays resulted in significant (FDR < 0.05) up-regulation of 724 exons, to carbon ions – of 511 exons and to iron ions – of 708 exons ( Supplementary Table 6 ). In this case, more overlap was observed between iron ions and X-rays – 32.8% of exons were in common ( Figure 7B ). When comparing the fold-changes in expression of the overlapping 246 exons ( Supplementary Table 6 ), the highest induction levels were shown by carbon ions ( Figure 7C ). In addition, changes in expression of the 20-exon signature identified earlier (47) were compared between different radiation types ( Supplementary Table 6 ). These 20 exons are particularly responsive to X-rays and important for the sample classification according to the exposure dose. The comparison revealed that most of the above-mentioned 20 exons are in general less responsive to heavy ions compared to X-irradiation ( Figure 7D ). The detailed results for four genes overlapping for all radiation types (PCNA, VWCE, FDXR and MAMDC4) are shown in Figures 7E – H . In this case, the most pronounced alternative splicing response was observed after carbon ions exposure. This was especially the case for MAMDC4 and VWCE. The Gene Ontology Biological Processes terms enriched in alternatively spliced genes common for all radiation types were predominantly related to apoptosis and DNA damage repair ( Supplementary Table 5 ). Carbon ion exposure resulted in alternative splicing of several genes coding for classical HLA class I molecules (HLA-A, HLA-B, and HLA-H) and class II molecules (HLA-DMB) and histone-coding genes (HIST2H3A, HIST2H3PS2, HIST2H3C, and HIST2H3D), which was not observed following the exposure to iron ions or X-rays.

Next to a transcriptional profiling, we evaluated the genotoxic impact of radiation exposure in the PBMCs over time for all radiation types. While quantification of DNA damage is commonly performed by counting the number of nuclear γH2AX foci (56, 57), high-LET radiation induces strongly clustered breaks along the track of the beam that may result in few microscopic foci, but with large relative size when a cell is visualized perpendicular to the orientation of the beam track ( Figure 8A ). γH2AX foci nuclear occupancy may therefore be a better measure for damage severity (45). Indeed, when calculated as the number of foci per nucleus, the absolute number of unrepaired breaks after 24 h was similar for all radiation types ( Figure 8B ). However, when considering foci occupancy, the amount of unrepaired DNA double-strand breaks (DSBs) 24 h after irradiation was 23% for X-rays, 42% for carbon ions and 31% for iron ions. When considering the foci occupancy per nucleus, the amount of damage still present 24 h after exposure to iron ions was comparable to the amount of damage observed in X-irradiated cells at 0.5 h ( Figure 8C ). For X-rays and iron ions the maximal foci occupancy was detected at 0.5 h post-irradiation, while for carbon ions this peak was observed at 2 h post-irradiation. The severity of DNA damage and kinetics of the repair were therefore clearly LET-dependent.

To compare the changes in gene expression with DNA repair kinetics at the level of individual donors, samples of four individuals that were irradiated with iron ions were used. All four individuals showed a clear reduction in the number of DSBs with time ( Figure 9A ). However, for Donor 1 the percentage of unrepaired DNA DSBs after 24 h was 43.6% and the difference between the damage detected at 0.5 h and 24 h time point was not statistically significant. In contrast, for the other donors this difference was significant, and the percentages of unrepaired damage were lower, ranging between 25.2-28.7%.

Hierarchical clustering of the gene expression profiles of DNA repair-related genes showed time- and subject-dependent expression. This resulted in two major clusters of samples depending on the time point, with 24-h samples segregating from 8-h and 12-h samples ( Figure 9B ). Within each time cluster, expression profiles of Donor 1 clustered separately from those of the other three subjects ( Figure 9B ). Several radiation-induced DNA repair genes (e.g., PCNA, DDB2, RBM14) showed an enhanced radiation response in Donor 1 compared to other donors, especially after a high dose ( Figures 9C – E ). This donor also showed elevated levels of expression of several DNA damage response-related genes, (e.g., ATM, ATR, RAD51D, MRE11A) independent of the irradiation dose and time point ( Figures 9F – I ). This indicates that individual differences in the overall and radiation-induced expression levels of DNA repair genes exist, which may possibly explain individual differences in DNA repair kinetics.