Transcriptional and Secretomic Profiling of Epidermal Cells Exposed to Alpha Particle Radiation

Cell Exposure and Harvesting

Primary human derived keratinocytes (HEKn) obtained from the Invitrogen (ATCC, Manassas, VA, USA), were maintained in a humidified incubator at 37°C, 5% CO₂, 95% air in 75 cm² flasks (T-75). The cells were grown in EpiLife media (Invitrogen) containing Human Epidermal Growth Supplements (Invitrogen) as per manufacturer’s instructions. Cells were seeded onto 35 mm culture dishes containing 2 layers of 2.5 μm-thick mylar films at a density of 2.5 x 10⁵ cells with 5 ml of culture media and allowed to grow to 90% confluency over 3 days. Cells were then washed using Phosphate Buffered Saline (PBS) and an additional 5 ml of fresh media was added prior to α-particle radiation exposure. These exposures were performed at doses of 0 (control), 0.5, 1.0 or 1.5 Gy using americium (²⁴¹Am) electroplated discs with an activity level of 66.0 kBq ± 3% (dose rate of 0.98 ± 0.01 Gy/h, LET of 127.4±0.4 keV/µm). The absorbed dose of α-particle radiation to which the cells were exposed was calculated using the GEANT4 v.9.1 Monte Carlo toolkit [19]. For the α-particle exposures, cells were cultured in thin mylar based plastic dishes (MD) (Chemplex Industries, Palm City, FL, USA), which allowed for penetration of α-particles. Cells destined for X-radiation at doses of 0 Gy to 1.5 Gy were exposed using the X-RAD 320 X-ray irradiation system at a low dose rate of 0.98 ± 0.05 Gy/h, 120 keV (Precision X-ray, Inc., North Branford, CT, USA ). Following irradiation, cells were either washed and harvested 24 h after exposure in 350 μl RLT buffer (Qiagen Inc, Mississauga, ON), for transcriptional profiling or 200 μl of media was repeatedly sampled at 24, 48, 72 and 96 h for secretomic analysis. All samples were flash frozen and stored at -80°C until ready for use. For each exposure group and end point, a total of 5 independent experiments were conducted. Cell viability was assessed 96 h post-exposure via Trypan Blue (Bio-Rad) viability assay. After 96 h both the control and 1.5 Gy α- and X-ray exposed cells remained ~90% viable (data not shown).

Microarrays

Frozen cell lysates were pipetted onto a QIAshredder spin column, and total RNA was extracted using the RNeasy Mini kit according to the manufacturer’s instructions (Qiagen Inc.). Additionally, Qiagen’s On-Column RNase-free DNase was used to eliminate possible DNA contamination. The concentration and quality of the RNA sample isolation was determined through spectrophotometric means (optical density (OD) ratio of A260:A280), and using an Agilent 2100 Bioanalyzer, following the manufacturer’s instructions. All RNA Integrity Numbers (RINs) were above 9.0. An input of 200 ng of total RNA was used for whole genome analysis following the Illumina(r) Whole Genome Expression Profiling Assay Guide (11317302 Rev. A). Samples were hybridized on Illumina human-12 v4 RNA BeadChips. BeadChips were imaged and quantified with the Illumina iScan scanner and data were processed with Illumina GenomeStudio v2010.2 software.

The preprocessing includes three steps: background correction (performed in GenomeStudio software: www.illumina. com), quantile normalization [20] and log2 transformation of normalized data. The last two steps were performed in lumi R package [21]. LIMMA (linear models for microarray data) [22, 23] was used to identify differentially expressed gene signatures under different conditions. Briefly, it starts by fitting a linear model for each gene in the data, and then an empirical Bayes (EB) method is used to moderate the standard errors for estimating the moderated t-statistics/F-statistics for each gene, which shrinks the standard errors towards a common value. This test is similar to an ANOVA method for each gene except that the residual standard deviations are moderated across genes to ensure more stable inference for each gene. The moderated standard deviations are a compromise between the individual genewise standard deviations and an overall pooled standard deviation. The false discovery rate (FDR) was evaluated using the Benjamini and Hochberg multiple testing procedure [24].

Quantitative Real-time Polymerase Chain Reaction (qRT-PCR)

Selected genes identified by microarray analysis as displaying statistical significance, with fold changes of 1.5 or higher and for which primers were validated were further assessed by qRT-PCR. Total RNA (100 ng) isolated from cells were reverse transcribed into complementary DNA using the RT2 First Strand Kit (SABiosciences Corp., Frederick, Maryland, USA). Gene profiling was performed according to the manufacturer instructions using custom RT2-profiler PCR arrays (SABioSciences Corp.). Reactions were prepared in 96-well plates and performed on a spectrofluorometric thermal cycler (Biorad iCycler; Hercules, CA). The relative expression of each gene was determined by using the comparative threshold (Ct) method. Analysis of qRT-PCR expression profiles and statistical analysis of data were performed using the super array biosciences web portal for data analysis of their products.

(SABiosciences http://www.sabiosciences.com/pcr/arrayanaly-sis.php).

Bio-Plex Assay

Supernatants from exposed (0.5, 1.0 and 1.5 Gy α-particle and X-rays) and control samples were analyzed for secretion levels of 27 cytokines starting at 24 h post exposure until 96 hr. Cytokines examined included interleukin (IL)-1b, IL-1ra, IL-2, IL-4, IL-5, IL-6, IL-7, IL-8, IL-9, IL-10, IL-12, IL-13, IL-15, IL-17, eotaxin, fibroblast growth factor (FGF), Granulocyte-Colony Stimulating Factor (G-CSF), granulocyte/ macrophage colony-stimulating factor (GM-CSF), interferon-γ (IFN-γ), interferon gamma-induced protein 10 kDa (IP-10), monocyte chemotactic protein-1 (MCP-1), macrophage inflammatory protein (MIP)-1a, MIP-1b, platelet derived growth factor (PDGF)-bb, regulated upon activation, normal T-cell expressed, and secreted (RANTES), tumor necrosis factor-α (TNF-α) and vascular endothelial growth factor (VEGF). Analysis was performed using a multiplex suspension array assay according to the manufacturer’s instructions (Bio-Rad). Briefly, conjugated beads were allowed to react with a sample containing a known (standard) or unknown amount of cytokines for 30 min. Conjugated beads with bound target were then washed and incubated with biotinylated detection antibodies that were directed against specific cytokine epitopes. The resulting complexes were then incubated for a further 10 min with streptavidin-phycoerythrin and excess reagent was washed off and assessed for bound cytokine using a microtiter plate reader (Bio-Rad). The concentration of cytokines in supernatants was then assessed from the generated standard curves for each individual cytokine using Bio-Plex software (Bio-Rad). To mine for differentially expressed cytokines the approach of Kleynhans et al. [25] was employed. Briefly, for each radiation type and cytokine, a two-way ANOVA (TIME *DOSE) linear mixed effects model was fit to take the repeated measures into account. Compound symmetry correlation structure was assumed, meaning for all samples the slopes of the lines regressing the dependent variable (intensity of cytokine) were the same over time. DOSE * TIME were treated as a factorial design and used to generate the main effects (DOSE * TIME) and their interaction (DOSE * TIME). Repeated measures at different time-points were treated as random factors.

Pathway Analysis

For the two gene lists (67-gene list and 5-gene list), the significance of predefined pathways, functions and network using Ingenuity Pathway Analysis software (IPA, version 11904312, http://www.ingenuity.com Core Analysis) was conducted. Functional analysis identified the biological functions and/or diseases that were significant to the given gene list. Molecules from the list that were associated with biological functions and/or diseases in the Ingenuity Knowledge Base were considered for the analysis. Right‐tailed Fisher’s exact test was used to calculate a p-value determining the probability that each biological function and/or disease assigned to that list was due to chance alone. Canonical pathway analysis identified the pathways from the IPA library that were significant to the gene list. Molecules from the gene list that were associated with a canonical pathway in the Ingenuity Knowledge Base were considered for the analysis. The significance of the association between the gene list and the canonical pathway was measured in two ways: 1) A ratio of the number of molecules from the gene list that map to the pathway divided by the total number of molecules that map to the canonical pathway is displayed. 2) Fisher’s exact test was used to calculate a p-value determining the probability that the association between the genes in the list and the canonical pathway is explained by chance alone. The functional analysis of a network identified the biological functions and/or diseases that were significant to the molecules in the network. A network is a graphical representation of the molecular relationships between molecules. Molecules are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). All edges are supported by at least one reference from the literature, from a textbook, or from canonical information stored in the Ingenuity KnowledgeBase. The intensity of the node color indicates the degree of up (red) or down (green) regulation. Nodes are displayed using various shapes that represent the functional class of the gene product. Edges are displayed with various labels that describe the nature of the relationship between the nodes (e.g., P for phosphorylation, T for transcription).The network molecules associated with biological functions and/or diseases in the Ingenuity Knowledge Base were considered for the analysis. Similar to the functionality analysis a right-tailed Fisher’s exact test was used to calculate a p-value determining the probability that each biological function and/or disease assigned to that network was due to chance alone.

Differential Gene Expression Profiles and Pathway Analysis

To mine for differentially expressed genes, LIMMA (22, 23) was employed. In order to reduce the size of the gene lists and to identify more reliable radiation specific genes, significant genes identified by LIMMA (FDR≤0.05) were further filtered on a |1.5| fold change (FC) cut-off. Cells exposed to low dose-rate X-rays at all three doses tested (0.5, 1.0, 1.5 Gy) did not express differentially modulated transcripts that were shown to be statistically significant using Benjamini-Hochberg (BH) corrected p-values relative to the control treatment groups (data not shown). Using less stringent statistical methods, the X-ray treatment did result in changes in gene expression, however these transcripts warrant further validations using alternate methods in order to ensure the responses are robust.

On the contrary, keratinocytes exposed to α-particles showed modulated expression in a number of genes (Table 1). At the lowest dose of α-particle radiation (0.5 Gy) no statistically significant genes were shown to be expressed however the medium dose (1.0 Gy) of radiation caused changes in expression of a total of 5 genes (CDC20, CCNB2, AURKA, BUB1 and NEK2) and all were shown to be upregulated relative to the control treatment group. These 5 genes were also observed to be expressed at the high dose of α-particle radiation (Fig. 1). All 5 of these genes were shown to be upregulated with fold changes approximately in the range of 1.5 fold of the unirradiated control treatment group (Table 2). Pathway analysis of these 5 differentially expressed genes showed that the top high-level functional categories were in relation to DNA replication, recombination and repair, and tissue morphology (FDR≤0.05). Analysis of the low-level functions showed involvement in cell cycle regulation and mitotic repair (Table 3). The top scoring network associated with these genes was centered around TP-53 signaling (Fig. 2).

A total of 68 genes were identified that were differentially expressed at the high dose (1.5 Gy) (Table 4a,b) of α-particle radiation when compared to control. Of these genes, 75% were upregulated and 25% were downregulated. Ninety three percent (93%) of the genes were exclusive to this dose of radiation and 7% of the total responding genes were also observed to be expressed at the medium (1.0 Gy) dose of radiation. The top expressing upregulated genes with fold changes greater than 1.8 fold were KIF20A, NEFM, CSN3, C11orf75, SLC16A9. Strongly repressed transcripts with FC≤-2 included C7orf10, HIST1H2BD, BMP6 and HIST1H2AC. A detailed functional analysis of the high expressing genes using a pathway analysis tool indicated that the top canonical pathways associated with this gene set were mitotic roles of polo-like kinase, cyclins and cell cycle regulation, chromosomal replication, ATM signaling and p53 signaling (Table 5). Nine genes were associated with these canonical pathways including KIF23, CDC20, CCNB2, CCNB1, CCNA2, MCM6, MCM2, TP53INP1 and BIRC5. Two top scoring pathways were obtained for this gene set. The first network was centered on TP-53 (Fig. 3A) signaling and the second top scoring network was clustered around ERK1 signaling (Fig. 3B).

Validation of Altered Expression of Radiation-Inducible Genes

For a selected few genes, for which validated primers were available, the observed microarrays results were verified using qRT-PCR. As shown in Table 6, all genes that exhibited a significant response at the high dose (1.5 Gy) of α-particle radiation 24 h post-exposure using microarray analysis were also observed to exhibit a similar trend using qRT-PCR. A comparable pattern of expression was also observed at the low dose (0.5 Gy) using the two methodologies. At the medium dose, although the fold changes were similar, the significance level did not always correspond. The five genes CDC20, CCNB2, AURKA, BUB1 and NEK2 that were shown to be significantly modulated at both the medium and high dose of radiation using microarray technology did express similar fold changes using qRT-PCR. However, these genes were not shown to be statistically significant with the exception of CDC20. Overall, these results suggest a comparable response between the two technologies, particularly at the high dose level.

Bio-Plex Analysis

Supernatants from cells exposed to α-particle radiation and X-rays were assessed for the relative expression of 27 secreted cytokines over a time-course of 96 h. Analysis using two-way ANOVA (TIME*DOSE) revealed dose main effects (p<0.05) at each of the doses. Cytokines that were shown to be statistically significant at all three doses of α-particle radiation were plotted (PDGF-bb, IL-13 and RANTES) over the time-course of the experiment (Fig. 4). A similar analysis conducted for X-irradiated cells also revealed dose main effects (p<0.05). The expression of 6 cytokines (IP-10, IL-12, IL-10, RANTES, IL-5 and MIP-1a) was observed to be modulated for all three doses (Fig. 5). For both radiation types, the DOSE *TIME interaction did not reveal statistical significance at all interactions, however a clear biological trend was observed across the doses and time-points (Table 7).