mRNA and small RNA gene expression changes in peripheral blood to detect internal Ra-223 exposure

Abstract

Purpose

Excretion analysis is the established method for detection of incorporated alpha-emitting radionuclides, but it is laborious and time consuming. We sought a simplified method in which changes in gene expression might be measured in human peripheral blood to detect incorporated radionuclides. Such an approach could be used to quickly determine internal exposure in instances of a radiological dispersal device or a radiation accident.

Materials and methods

We evaluated whole blood samples from five patients with castration-resistant prostate cancer and multiple bone metastases (without visceral or nodal involvement), who underwent treatment with the alpha emitting isotope Radium-223 dichloride (Ra-223, Xofigo^®). Patients received about 4 MBq per cycle and, depending on survival and treatment tolerance, were followed for six months. We collected 24 blood samples approximately monthly corresponding to treatment cycle.

Results

Firstly, we conducted whole genome screening of mRNAs (mRNA seq) and small RNAs (small RNA seq) using next generation sequencing in one patient at eight different time points during all six cycles of Ra-223-therapy. We identified 1900 mRNAs and 972 small RNAs (222 miRNAs) that were differentially up- or down-regulated during follow-up after the first treatment with Ra-223. Overall candidate RNA species inclusion criteria were a general (≥|2|-fold) change or with peaking profiles (≥|5|-fold) at specific points in time. Next we chose 72 candidate mRNAs and 101 small RNAs (comprising 29 miRNAs) for methodologic (n = 8 samples, one patient) and independent (n = 16 samples, four patients) validation by qRT-PCR. In total, 15 mRNAs (but no small RNAs) were validated by methodologic and independent testing. However, the deregulation occurred at different time points, showing a large inter-individual variability in response among patients.

Conclusions

This proof of concept provides support for the applicability of gene expression measurements to detect internalized alpha-emitting radionuclides, but further work is needed with a larger sample size. While our approach has merit for internal deposition monitoring, it was complicated by the severe clinical condition of the patients we studied.

Keywords:

• Incorporation
• alpha-emitter
• peripheral blood
• biomarker
• gene expression changes
• mRNA
• small RNA
• Ra-223
• prostate cancer
• Xofigo

Introduction

Radioactive sources are widely used in industry as well as medicine and research (International Atomic Energy Agency 2005; Miller 2012) and obtaining such material appears relatively easy. Finding abandoned radioactive material is concerning and regularly reported (IAEA 2004). In the past, orphan radioactive sources have caused large scale exposure events, such as a medical cesium source in Goiania, Brazil (Eichholz 2005). A more worrisome misuse of radioactive material would be a terrorist attack (Rosoff and Von Winterfeldt 2007; Hall and Giaccia 2012; Confer et al. 2018). A radiologic dispersal device (RDD, also called ‘dirty bomb’) is a weapon that combines a radioactive source with conventional explosives.

In such an RDD scenario, exposure of hundreds or thousands of people may occur. Besides mechanical trauma and direct irradiation, a much larger number of people may be externally contaminated by radioactive material with the concomitant danger of radionuclide incorporation and internal contamination (Rump et al. 2018). Acute radiation syndrome would be unexpected from radionuclide(s) incorporation, except in special cases such as the Litvinenko poisoning (Harrison et al. 2017). However, irradiation accompanied by internal contamination may cause long-term stochastic health effects (e.g. cancer). As an initial therapy, decorporation agents to eliminate internal radionuclides can be administered, and thus absorbed dose can be reduced and the adverse health effects minimized. Therapeutic efficacy decreases if therapy with decorporation agents is initiated too late and in most cases the time window is hours to several days in order to achieve acceptable results. This also depends on the types of radionuclides, the physicochemical properties and the metabolic pathway (Rump et al. 2016, 2017). Considering limited doses of decorporation antidotes in a nuclear mass casualty event, exposed individuals needing urgent treatment must be prioritized (Chaudhry 2008; Rump et al. 2018). Radiation dose from physical dosimetry (like measurements by monitoring portals) are thought to be sufficient for sensitive and high throughput triage in gamma-radionuclide exposure situations, but fall short for inhalation, ingestion or absorption routes.

For detection of incorporated alpha-emitting radionuclides, excretion analysis has been the established method for decades (Karlsruhe 1979; Taylor and Carolina 2000; Jefferson et al. 2009). Because internal dose cannot be measured directly, it has to be modeled and calculated indirectly under several assumptions. Knowledge of the dynamic behavior (biokinetics) of the incorporated material in the human body is essential for radiation dose estimation of individual organs. This is laborious and time consuming. Although gamma emissions comprise only about 1% of the total emitted energy in Ra-223, there are reports that quantitative imaging may be feasible for dosimetry in case of primarily alpha-emitting radionuclides (Hindorf et al. 2012). However, the need for gamma cameras and imaging time renders this approach infeasible for screening large numbers of people.

In a large-scale nuclear exposure scenario, a high-throughput diagnostic system is clearly needed (Kabacik et al. 2011; Boldt et al. 2012; Manning et al. 2013; Matthias Port et al. 2016; Jacobs et al. 2020). One promising approach to consider is differential gene expression. Gene expression analysis has been shown as an alternative tool for high throughput biodosimetry and predictive of clinical outcomes (Meadows et al. 2008; Paul and Amundson 2008; Meadows et al. 2010; Paul et al. 2011; Port et al. 2016, 2018, 2019).

We sought an easier method in which changes in gene expression, as measured in human peripheral blood to detect incorporated radionuclides, could be used as a high-throughput biodosimetric approach. Such an approach could quickly determine internal exposure in a nuclear scenario.

In collaboration with the University Medical Center (Mainz, Germany), we measured gene expression in five patients with castration-resistant prostate cancer and multiple bone metastases (without visceral or nodal involvement). These patients underwent intravenous treatments with the alpha-emitting isotope Radium-223 dichloride (Ra-223, Xofigo^®) for up to six months. Blood samples (n = 24) were collected before and during treatment cycles at approximately monthly intervals. We performed a whole genome screen to evaluate potential candidate genes on a transcriptional (protein coding mRNAs) as well as a post-transcriptional (non-coding small RNAs) level by differential gene expression (DGE) relative to pretreatment. These RNA species were then validated using qRT-PCR methodology.

Materials and methods

Patients, sample collection, Ra-223 (Xofigo^®) treatment, ethical approvals, guidelines

Five advanced-stage castration-resistant prostate cancer patients with multiple bone metastases, but without visceral or nodal involvement, were sequentially enrolled and followed over six months. A two-year time period elapsed for recruiting eligible individuals and approximately monthly subsequent blood collections. Twenty-four peripheral whole blood samples (2.5 mL) were obtained via venipuncture using the PAXgene Blood RNA system (BD Diagnostics, PreAnalytiX GmbH, Hombrechtikon, Switzerland).

All patients underwent treatment with the alpha emitting pharmaceutical Radium-223 dichloride (Ra-223, Xofigo^®). The therapeutic plan was for six intravenous injections of Ra-223 (at a dosage level of 55 kBq per kilogram of body weight) every four weeks (corresponding to one cycle). Since the absorbed dose to the blood and the bone marrow increased for each therapy cycle, the values of the x-axes in Figures 2, 4, 5, 6(A and B) are proportional to the cumulative absorbed dose. Blood samples were taken either before each Ra-223-application or at day 15 between two applications (Table 1). During the screening approach (Phase I), we assessed whole genome next generation sequencing (NGS) screening (mRNA seq and small RNA seq) in one patient with one blood sample collected before the first Ra-223 application (reference) and seven more at specific time points during treatment cycles. To validate the candidate RNAs in Phase II, we used the remaining 16 blood samples (Table 1).

Figure 1. Cumulative absorbed bone marrow dose for Ra-223 administered at the different points in time (therapy cycles, reflecting the days of therapy), representative for patient #1 (screening patient), up to 140 days. Linear regression of dose over time is shown as a dashed line.

Figure 2. Number of potential candidate RNAs (y-axis) showing a differential up- or down-regulation for mRNAs (left part) or small RNAs (right part) over time and treatment cycles (in days, x axis) in patient #1 (screening). Depicted are up- and down-regulated RNAs as well as the sum of both. Corresponding measurements over time are shown with a spline. Gray areas represent the time interval with the highest number of deregulated RNAs.

Figure 3. Venn diagrams show the number of up- (left part) and down-regulated (right part) mRNAs (upper part) and small RNAs (lower part) that were differentially expressed between pretreatment and post-treatment samples in patient #1 (screening) at different time points, e.g. treatment cycles. Overlapping circles represent the number of RNAs which were in common at multiple time points. Numbers outside the overlapping region represent RNAs that were not in common for other time points. Numbers in parenthesis represent the total number of differentially expressed genes per time point. Based on the fold-change and our preference for genes differentially expressed in all time points or peaking at certain time points, 72 candidate mRNAs and 101 small RNAs (29 miRNAs) were selected and forwarded in phase II for validation.

Figure 4. Differential gene expression based on next generation sequencing of representative mRNAs (left part) and small RNAs (right part) deregulated over the whole period of treatment (≥|2|-fold) is shown in separate panels over the time of treatment for patient #1 (screening) and measurements are connected with a straight line. Fold changes are calculated with the unexposed sample before first Ra-223 application used as a reference. Gene IDs of the examples are listed in the right panel of both graphs. Gray areas represent the range of gene expression within the methodologic variance (<|2|-fold). Gene names are listed in all graphs and refer to the last data points (day 140) in descending order, but symbols were omitted for clarity reason.

Table 1. Overview of the included samples, split two-phase study design and baseline characteristics of the patients (age, percentage of skeletal metastases compared to the whole skeleton) are shown.

			Sampling time points
		Skeletal metastases	Cycle I		Cycle II		Cycle III	Cycle IV	Cycle V	Cycle VI	Number of samples
Task/ patient ID	Age (years)	Area of whole skeleton (%)	day 0	day 15	day 28	day 43	day 56	day 84	day 112	day 140	n	Total n
Screening + Methodological Validation
1	73	1	x	x	x	x	x	x	x	x	8	8
Independent Validation
3	71	11	x	x	deceased						2	16
2	69	3	x	x	x		x	deceased			4
4	77	3	x		x		x	x	deceased		4
5	76	6	x		x		x	x	x	x	6

The overview shows the samples for different points in time over the treatment cycles used in patients #1–#5 as well as the number of samples collected from each patient.

Ethical approvals were obtained from the responsible Ethics committee (#837.361.15, Landesärztekammer Rheinland-Pfalz, Mainz, Germany) prior to study initiation. Informed consent was obtained from each participating patient. All methods described in this manuscript (e.g. RNA quantity/quality and TaqMan qRT-PCR) were performed according to the relevant guidelines, regulations and standard operating procedures implemented in our laboratory in 2008 when the Bundeswehr Institute of Radiobiology became DIN-accredited by TÜV Süd München, Germany (DIN EN ISO 9001/2008) and were described in detail elsewhere (Port, Hérodin, Valente, Drouet, Lamkowski, et al. 2017; Port, Hérodin, Valente, Drouet, Ullmann, et al. 2017).

Pharmacokinetics and dosimetry

Ra-223, a bone-seeking alpha-emitter, shows a rapid clearance from the blood compartment, with fecal excretion as the major route of elimination (Carrasquillo et al. 2013). After intravenous injection, Ra-223 is rapidly cleared from the blood by distributing primarily into bone, with preferential uptake on the osteoblastic lesions. Fifteen minutes post-injection, about 20% of the injected radioactivity remains in the blood stream, decreasing to less than 1% at 24 hours after the injection (Information FP 2013). We retrospectively determined the doses absorbed by red bone marrow using two different software tools.

Equivalent doses were calculated according to the ICRP (International Commission on Radiological Protection) model for radium using the DOSAGE software (Lassmann and Nosske 2013). The ICRP biokinetic model assumes that 25% of the administered radium localizes in the bone; daughter products are also taken into account. Additionally, bone marrow doses were calculated using the Integrated Modules for Bioassay Analysis (IMBA) software (Birchall et al. 2007). The calculated committed effective dose (50 years) by internal contamination was calculated and the effective doses for the time points given were estimated based on the fraction of the area under the activity-time curve relative to the total area under the curve. This was performed for bone marrow dose, as well as for the absorbed dose to the blood. In order to estimate the absorbed doses, we used a weighting factor of 0.12 for this tissue and a radiation weighting factor of 20 for alpha radiation as recommended by the ICRP guidelines 103 (2007) for stochastic effects (Harrison et al. 2004).

RNA extraction and quality control

Whole blood samples (2.5 mL) were processed following the PAXgene Blood RNA system. In brief, blood was drawn into a PAXgene Blood RNA tube at the University Medical Center Mainz. The tube was gently inverted (10 times), stored at room temperature for 2–12 h, then at −20 °C. After all samples have been collected, the PAXGene tubes were sent to the Bundeswehr Institute of Radiobiology for further processing according the established workflow (Port et al. 2016). After thawing, washing and centrifugation, cells in the supernatant were lysed (Proteinase K, PAXgene Blood RNA protocol) followed by addition of Lysis/Binding Solution taken from the mirVana kit (Life Technologies, Darmstadt, Germany). With the mirVana kit total RNA, including small RNA species, was isolated by combining a Phenol-Chloroform RNA precipitation with further processing using a silica membrane. After several washing procedures, DNA residuals became digested on the membrane (RNAse free DNAse Set, Qiagen, Hilden, Germany). RNA was eluted in a collection tube and frozen at −20 °C. Quality and quantity of isolated total RNA were measured spectrophotometrically (NanoDrop, PeqLab Biotechnology, Erlangen, Germany). RNA integrity was assessed by the 2100 Agilent Bioanalyzer (Life Science Group, Penzberg, Germany) and DNA contamination was controlled by conventional PCR using an actin primer. We used only RNA specimens with a ratio of A₂₆₀/A₂₈₀ ≥ 2.0 (Nanodrop) and RNA integrity number (RIN) ≥ 7.5 for NGS library preparation (IMGM Laboratories, Martinsried, Germany) or RIN ≥ 7.3 for qRT-PCR analysis.

Stage I screening

Whole genome screening for differentially expressed genes (protein coding mRNAs) and non-coding small RNAs was performed on eight RNA samples (one pretreatment sample and seven samples from treatment cycles, as shown in Table 1). Patient #1 provided a complete set of samples over six months of the therapeutic course.

For mRNAs, library preparation was based on the TruSeq^® Stranded Total RNA Kit (Illumina, San Diego, CA) following the manufacturer’s recommendations (Illumina Inc. TruSeq^® Stranded Total RNA Sample Preparation Guide – Catalog # RS-122-9007DOC – Part # 15031048 Rev. E. 2013). For small RNAs, 300 ng of RNA were used for library preparation, which was based on the NEBNext^® Small RNA Library Prep Set for Illumina (Illumina, San Diego, CA) following the manufacturer’s recommendations (New England Biolabs Inc. NEBNext^® Small RNA Library Prep Set for Illumina (Multiplex Compatible) – Catalog # E7330S/L – Version 4, 9/17).

After library quantification using the Qubit^® dsDNA HS Assay Kit (Thermo Scientific, Waltham, MA) and assessment of length distribution by the Agilent Bioanalyzer system (Agilent 2100 Bioanalyzer, Agilent, Santa Clara, CA) sequencing was performed on the Illumina NextSeq 500 sequencing platform (Illumina, San Diego, CA). For in-depth analysis of differential gene expression and annotation of reads, we used the CLC Genomics Workbench (CLC Bio-Qiagen, Aarhus, Denmark; version 11.0.1). After import of read data into the analysis tool, sequence reads were mapped against the human reference genome (GRCh38.p7, NCBI, and mirRBase Release 21). The tool ‘Differential Expression for RNA-Seq’ was used for identification of DGEs. In this study, only genes were classified as induced in a specific comparison with Log₂ fold change ≥2 among compared samples. Those transcripts were considered to represent a candidate gene for validation in phase II.

Stage II: Validation of stage I candidate genes via qRT-PCR

To validate the candidate genes from stage I (screening), we used the remaining RNA samples on a customized combined low-density array (LDA, high throughput qRT-PCR platform). We quantified 72 candidate mRNAs and 75 candidate small RNAs (other than miRNA) using TaqMan chemistry (order IDs for the examined genes are provided in Supplemental Table 1). A 1 µg RNA aliquot of each RNA sample was reverse transcribed using a two-step PCR protocol (High Capacity Kit). 400 µl cDNA (1 µg/0.5 µg RNA equivalent) was mixed with 400 µl 2 x RT-PCR master mix and pipetted into the 8 fill ports of the LDA. Validation of 29 candidate miRNAs was performed similarly, but cDNA-synthesis was performed using the TaqMan Advanced miRNA cDNA Synthese-Kit and TaqMan Advanced miRNA Assays. Cards were centrifuged twice (1,200 rpm, 1 min, Multifuge3S-R, Heraeus, Germany), sealed, and transferred into the 7900 qRT-PCR instrument. The qRT-PCR was run for two hours in a 384-well LDA format. For the LDA Ct-values (threshold cycles), we used the median expression on each LDA for normalization purposes of mRNAs and small RNAs (except miRNAs) and the internal control spotted on the LDA for miRNAs. All technical procedures were performed in accordance with standard operating procedures implemented in our laboratory since 2008 when the Bundeswehr Institute of Radiobiology was certified according to DIN EN ISO 9001/2008.

Bioinformatics

All genes from the screening approach associated with a ≥2-fold gene expression difference (up or down relative to the reference) underwent gene set enrichment analysis using PANTHER pathway software (http://www.pantherdb.org/). PANTHER groups genes with similar biological function based on their annotation (reference list was the current homo sapiens GO database).

Statistical analysis

During screening for potential candidate RNAs showing DGE over time we, used the unexposed pretreatment sample (reference) and the samples per time point (n = 7) after treatment for seven pairwise comparisons. The filter (criteria) that we used was ≥|2|-fold RNA expression changes over time after first treatment with Ra-223. This applied for those RNA species that were up- or down-regulated in all post-treatment samples. The filter for candidate RNAs at specific points in time after Ra-223 (referred to as peaking profiles) was defined as ≥|5|-fold RNA expression change. Due to missing replicates in the screening phase, we could not correct for multiple comparisons (e.g. FDR) on the screening phase I of the study. We considered multiple comparisons in the bioinformatic approach (FDR) and in the validation part (phase II, Bonferroni) and corrected p-values for multiple comparisons accordingly.

In an effort to address the inter-individual radiation-induced gene expression responses over time, we applied an ordinal scale, hereby shifting the time scale so that the expected gene expression event (e.g. down-regulation based on the screening results) fell into the same dimensionless time point for all patients. Calculating the arithmetic mean per point in time would result in a meaningless interpretation. Afterward, we calculated the geometric mean and the standard deviation for each dimensionless time point (ordinal scale) and performed a curve fit to capture gene expression patterns that were in common for all patients. Excel 2010 (Microsoft) was used for descriptive statistics (n, mean, standard error of mean, standard deviation, min, max) and generating the tables. SigmaPlot (Systat Software, version 14.0/2019, Erkrath, Germany) was used for graphical presentations and curve fits.

Results

Material available for the two phase study design

For this study, 24 blood samples from five patients (mean age 73.2 years) with castration resistant prostate cancer and multiple bone metastases, who underwent treatment with Ra-223, were provided for analysis (Table 1). Three patients died during Ra-223 treatment, one patient suspended Ra-223 therapy.

From 2.5 mL whole blood samples, we isolated an average of 7.1 µg (stdev +/– 3.0 µg) total RNA before the first Ra-223 application and 3.6 µg (stdev +/– 1.6 µg) during the therapy cycles. RNA integrity (RIN) showed a mean value of 8.0 (stdev +/– 0.8) before and 8.4 (stdev +/– 0.8) during therapy, respectively, suggesting high quality RNA sufficient for both methods.

Dose estimation

Depending on the duration of Ra-223 therapy, we derived the activity applied during treatment cycles (dosage level 55 kBq/kg body weight per treatment cycle) and determined the equivalent and effective doses absorbed by the bone marrow as well as the peripheral blood. Calculated doses from the IMBA and DOSAGE software are shown in Table 2. From the effective dose and a tissue weighting factor of 0.12 for bone marrow, the total equivalent bone marrow dose was 33,125 mSv (= 3,975/0.12) and the total absorbed dose was 1,656 mGy (= 3,975/(0.12 × 20)) with a radiation weighting factor of 20 for alpha-radiation. The absorbed doses per administered activity (dose coefficient) for the red bone marrow is 70 mGy/MBq. The committed absorbed dose from treatment to the red bone marrow was 1,656 mGy in patient #1 when cumulating the calculated effective doses from the different cycles. The time course of the cumulated absorbed dose is shown for patient #1 (screening) in Figure 1. Comparing the dose estimates calculated from IMBA with the dose values calculated with the DOSAGE software, the results differ only slightly, but are in the same order of magnitude. The values calculated with DOSAGE software exceed the values determined with IMBA by about 2–38% (Table 2). For example, in patient #1, the absorbed dose was 1,656 mGv compared to 1,668 mGv. The absorbed blood doses according to DOSAGE software were about 25-fold lower compared to the red bone marrow (Table 2).

Figure 5. Differential gene expression based on NGS of exemplary mRNAs (upper part) and small RNAs (lower part) deregulated with peaks (≥ |5|-fold) at specific time points of treatment is shown in separate panels over the time of treatment for patient #1 (screening). Fold changes are calculated with the unexposed sample before first Ra-223 application used as a reference. Gray areas represent the time interval with the peaking RNAs. Gene names are listed in all graphs and refer to the last data points (day 140) in descending order, but symbols were omitted for clarity reason.

Table 2. Using the activity during treatment cycles (dosage level 55 kBq/kg body weight per cycle), the effective doses were calculated using IMBA and DOSAGE software.

Patient ID	Number of therapy cycles	Total Activity/Amount (MBq)	Bone marrow Absorbed dose (mGy, cumulative values)		Blood dose Absorbed dose (mGy, cumulative values)
			according to IMBA	according to DOSAGE	according to DOSAGE
1	6	23	1656	1668	65
2	3	13	919	939	37
3	1	3	223	224	9
4	4	19	1393	1401	55
5	6	19	1329	1357	53

Based on the effective doses, the equivalent and absorbed doses were calculated. A tissue weighting factor of 0.12 for bone marrow and a radiation weighting factor of 20 for alpha-emitters were assumed. The numbers refer to the committed doses by internal contamination for bone marrow in each patient calculated by adding up the calculated single doses for each treatment cycle of Ra-223. Blood doses according DOSAGE software were calculated accordingly. The cumulative values of the activity applied are depicted.

Stage I: screening for candidate RNAs

The quality parameter for a successful sequencing run expressed as the percentage of Q30 bases (translates into a 99.9% accuracy) is supposed to be >80% and was on average 92.3% in RNA Seq and 93.7% in small RNA Seq in our study. The average number of passed filter (PF) reads was 15.8 × 10⁶ (stdev 4.4 × 10⁶) for mRNAs and 23.2 × 10⁶ (stdev 2.8 × 10⁶) for small RNAs, respectively. For mRNAs, 40.3–52.9% of the reads mapped to known gene regions (exons and introns).

From about 20,000 protein coding mRNAs, 1900 genes were differentially up- or down-regulated (≥|2|-fold, examples shown in Figure 4, left part) over the six-month period or had peaking profiles at specific points in time (≥|5|-fold, examples in Figure 5, upper part). With respect to the number of differentially expressed mRNAs, a peak was observed at day 56 after the first Ra-223 injection with 120 up-regulated and 554 down-regulated differentially expressed genes (Figure 2, left part). We saw that down-regulated genes were almost 5-times more common than up-regulated genes. Similarly, we observed only 4 up-regulated mRNAs and 22 down-regulated mRNAs that were evident over all points in time (Figure 3, upper part). Bioinformatic analysis for gene enrichment of related biological processes revealed no significant overlap within the examined time points. This might be caused by the FC of 2 introduced as a very robust selection criterion in the absence of statistical analysis during stage I.

Figure 6. Inter-individual dose/time-to-gene expression patterns are depicted for CXCL5 (upper panel, (A) and RNF11 (lower panel, A) for each patient (#1–5). Corresponding mean values for each time point are provided in (B). No significant association over time was found. The time point, in which expected down-regulation of both genes for each patient was observed, appears above the x-scales. Applying an ordinal scale on the x-axis, hereby shifting individual dose/time-gene-expression pattern so that expected gene expression events align into the same time point (C), resulted in significant exposure-to-gene associations based on the time adjusted mean values for all patients (D). Error bars represent the standard deviation. Fitted functions are depicted in the graphs.

We detected 972 small RNAs (comprising 222 miRNAs) that were differentially up- or down-regulated (≥|2|-fold, examples in Figure 4, right part) over the follow-up time or with peaking profiles at specific time points (≥|5|-fold, examples in Figure 5, lower part). As for the number of differentially expressed small RNAs, we observed a time displacement compared to mRNAs with a peak at day 84 after the first Ra-223 treatment (Figure 2, right part). Here, the deregulation was different to mRNAs described before – we observed an almost 3-fold increase in up- rather than down-regulated genes. In total, 113 up-regulated small RNAs and 34 down-regulated small RNAs were in common for all time points (Figure 3, lower part).

Based on the fold-difference and consistent RNA expression changes over time, we selected 72 candidate mRNAs, 29 miRNAs and 75 small RNAs (other than miRNA) for validation at stage II (Supplemental Table 1). Among these, 16 mRNAs were deregulated over the time period (fold change ≥|2|, 3 up-regulated, 13 down-regulated) as well as 3 miRNAs and 16 small RNAs (other than miRNA). In addition, 56 mRNAs, 26 miRNAs and 59 small RNAs (other than miRNA) appeared deregulated with peaks at various points in time (≥|5|-fold), suggesting further explorative analysis.

Stage II: Validation with qRT-PCR

During stage II validation of the 72 candidate mRNAs from stage I, we found 48 mRNAs that could be validated methodologically when the criterion was a fold change of >|2| for NGS and qRT-PCR measurements in at least one time point after first application of Ra-223 (success rate of about 67%), i.e. fold change of ≥|2|. When considering all time points in the screening patient after first application of Ra-223, 105 from 504 NGS measurements could be methodologically validated via qRT-PCR, resulting in a success rate of about 21%. By means of positive methodological validation as well as independent validation, 15 genes could be validated considering inter-individual points in time. We were unable to validate the remaining mRNAs and all small RNAs for different reasons. In some, there were missing amplification plots in all or amplification plots were observed in only a minority of samples. In addition, we observed inconsistent (increased vs. decreased) as well as inverse DGE in the screening and the validation steps.

For the remaining 15 mRNAs, gene expression values within dose and time showed an exposure-to-gene expression association. The association appeared strongly related to the patient and an up- or down-regulation occurred at different points in time (Supplemental Table 2). For instance, for RNF11, DGE (≥|2|-fold) occurred on day 15 in patient #1, on day 28 in patient #2, on day 56 in patient #4 and on day 112 in patient #5.

Regarding these inter-individual results, a time-shift in response of DGE was observed. These individual response characteristics are shown in Table 3: we focused on the points in time when DGE (≥|2|-fold) was detected for the first time in each patient (#1–#5) and our 15 validated mRNA species. The first time point for DGE of all 15 mRNAs consistently occurred at the same time in each examined patient and this pattern was highly significant (p = .0001). So, DGE of most of the 15 mRNA species was observed as early as 15 days (median) for patients #1 and #3, 28 days (median) in patients #2 and #5 or even further delayed at 56 days of therapy in patient # 4. None of these mRNAs remained significant after Bonferroni correction for multiple comparisons (see also p-values in Table 3).

Table 3. Overview of candidate genes that were validated using qRT-PCR (15 mRNAs) and that appeared predictive for detecting a Ra-223 incorporation based on gene expression differences.

Gene name		First time points (days) when significant gene expression differences were detected in patients
	p-value	#1	#3	#2	#5	#4
CXCL5	.004	15	15	28	–	–
DMTN	.08	15	–	28	28	56
FBXO7	.08	15	–	28	28	56
GNG11	.03	15	–	15	–	56
MAP3K7CL	.01	84	15	15	–	–
SH3BGRL2	.08	15	–	28	28	56
TMEM56-RWDD3	nd	15	–	–	56	–
TUBB2A		15	–	28	28	–
RNF11	.004	15	–	28	112	56
PF4	.015	15	–	28	28	–
BZW1		43	–	28	112	56
PDZK1IP1	nd	15	–	28	28	56
PPBP		–	–	28	–	–
RAB27B		140	–	28	28	56
CXCL8		15	–	15	28	–
median		15	15	28	28	56

Shown are the time points (days) when significant gene expression differences (≥ |2|-fold) were detected for the first time after first Ra-223 injection separately in each patient (#1–#5). Median values for the days are provided below the horizontal double bar. p-values refer to regression models as shown in Supplemental figure 1. None of the p-values survived Bonferroni correction for multiple comparison.

In order to illustrate these inter-individual variability in response and the time-shift in DGE response, we introduced an ordinal time scale accordingly as detailed previously. This is graphically shown for CXCL5 (Figure 6, upper panel) and RNF11 (Figure 6, lower panel), corresponding to individual clinical follow-up (e.g. three patients did not reach therapy cycle 6). We observed several individual radiation-induced gene expression patterns (Figure 6(A)). In each patient, the changes in gene expression occurred at different time points and the averaged gene expression values over dose and time did not show a significant exposure-to-gene expression association that was common among all patients (Figure 6(B)). In an effort to account for inter-individual variability in response, we applied an ordinal scale, shifting individual gene expression pattern on an ordinal time scale and aligned expected gene expression events to the same time point (Figure 6(C)). Subsequently, significant exposure-to-gene expression relationships common to all patients were seen (Figure 6(D)) with significant responses of FC means per adjusted points in time. The corresponding graphic analysis is shown for PF4, PDZK1IP1, TMEM56-RWDD3, SH36GRL2, MAP3K7CL, GNG11, FBXO7 and DMTN in Supplemental Figure 1 and the same features, as described, were demonstrated in 11 of the 15 mRNA species.

Discussion

Measuring gene expression changes in peripheral blood represents a fast, early and high-throughput method in easily accessible peripheral blood samples and has proven to be promising as a biodosimetric tool after external radiation exposure (Chaudhry 2008; Rothkamm et al. 2013; O’Brien et al. 2018; Port et al. 2019). Some cytogenetic techniques have already been used for retrospective individual dose assessment; however, the scant availability of proper studies in this field limits conclusions (Giussani et al. 2020). We wondered whether gene expression changes might also serve as a high-throughput biodosimetric approach after radionuclide incorporation in a radiologic or nuclear scenario such as an accident or terrorist attack. Especially when thinking about incorporated alpha emitting radionuclides, an individual dose assessment is indispensable with regard to decorporation therapy and is superior to general dosimetry.

Getting access to peripheral blood samples from radionuclide treated patients, we examined the whole genome for radiation-related transcriptional (protein coding mRNAs) and post-transcriptional (non-coding small RNAs including miRNAs) gene expression changes in the peripheral blood after controlled incorporation of an alpha-emitting radionuclide over several months and performed independent validation for potential candidate genes of interest. We determined the doses absorbed by red bone marrow using two different software tools (DOSAGE and IMBA), hereby showing, that the results differ slightly, but are in the same order of magnitude with the values calculated with DOSAGE software exceeding the values determined with IMBA by about 21%. This may be due to the calculation methodology: IMBA provides the 50-year committed effective dose, so the values for the different time points of interest had to be derived using the fractions of the area under the activity time curve. A rapid uptake of Ra-223 into the bone combined with a slow release were the necessary assumptions for this calculation. In our calculations, we considered a radiation-weighting factor of 20 for alpha-emitters, in agreement with established ICRP-models for stochastic effect calculations (Harrison et al. 2004). For therapeutic applications and, consequently, deterministic effects, a radiation-weighting factor of 3–5 seems to be more appropriate (Lassmann and Eberlein 2018). As we correlate gene expression with absorbed dose on an ordinal scale, the choice of the radiation-weighting factor is of little significance. The accumulation of Ra-223 in bone remodeling processes together with radiation dose to metastatic bone lesions are the assumed mechanism of the observed beneficial therapeutic index (Blake et al. 2001; Miederer et al. 2015). The collateral radiation induced effects on the surrounding bone marrow are presumed to cause a delayed effect on gene expression that was reflected in the peripheral blood. The models used for dose estimation are not addressing modulations of e.g. individual metabolism, which in the case of cancer will have a major impact on the whole metabolism. Therefore, results just provide a magnitude of absorbed dose and GE measurements were associated with dose on an ordinal scale assuming an increased dose over time.

During the Phase I screening, we identified hundreds of mRNA species that also showed a wave-like increase in number of deregulated mRNAs as was found for small RNAs on the same samples. However, the number of deregulated mRNAs peaked at 56 days after first infusion of Ra-223 and for small RNAs this peak was delayed to 84 days with the number of deregulated genes increased by an order of magnitude (Figure 2). This pattern in time supports a causal mRNA – small RNA association which is well established in the literature (Guo et al. 2010; Andrés-León et al. 2017), whereas the time relation of the peaking species was opposite (Ostheim et al. 2020). However, these gene expression effects are probably overexposed by blood compartment associated biological changes (e.g. cell death) that occurred in the bone marrow as the primary radiation target and was reflected in the peripheral blood compartment with a time delay.

Other than in previous studies of our group, validation of candidate genes from stage I succeeded for 15 mRNA species only. No small RNAs could be validated. Several genes detected during the screening approach were not detected with qRT-PCR in all or several samples so that we excluded them from analysis. However, ‘opposite of expectation’ based on earlier work (screening with microarrays), we have identified many genes that were inversely regulated in the validation samples where we changed from NGS to qRT-PCR methodology. We speculate that this might reflect incompatibilities in NGS methodology, used for screening, and qRT-PCR, used for validation, which have to be considered. For instance, (1) based on recently performed NGS studies of our group, results were only successfully validated with qRT-PCR when identifying the exon-region contributing primarily to the radiation-induced fold differences in gene expression and using a TaqMan assay, which specifically covers this region (described elsewhere ⁴¹). (2) The CLC Genomics Workbench (tool ‘Small RNA Analysis – Annotate and Merge Counts’) was applied. For alignment of a read to the reference sequence, +/– 2 additional or missing nucleotides up- and/or downstream were allowed to be counted. Furthermore, the maximum number during strand-specific alignment of allowed mismatches within a read was set to two. This rule adjusts for methodologic imprecision inherent to NGS, but TaqMan Assays using qRT-PCR recognize only one specific sequence. If missing, even for one nucleotide difference in the nucleotide sequence, no amplification plot will be generated. (3) The reference group contained one sample only and, therefore, any inter-individual variance in gene expression was impossible to assess. This design, while imperfect, was chosen because of the slow pace of sample collection, in that only five patients were enrolled and followed over two years. Only two patients completed six months of therapy. Three patients died during therapy, reflecting the severity and refractory nature of their disease. Further sample collection was curtailed because this therapy has declined in use in favor of a radioligand therapy (De Vincentis et al. 2019). Despite these limitations, we speculated that genes consistently differentially expressed over all time points (eight measurements) in patient # 1 would be sufficient to select meaningful candidate genes. All these limitations may have reduced the number of successfully validated candidate genes. A larger sample size (ideally n = 25 according to power calculations), especially for referent samples, will be key for future gene expression studies. It is quite possible that several more genes in excess of the 15 successfully validated genes exist.

Nevertheless, in the course of this study and considering validation stage, expected DGE of most of the successfully validated 15 mRNA species was observed either at 15 days in two patients, 28 days in another two patients or 56 days of therapy in one patient (Table 3). We did not find that the DGE values were directly proportional to the absorbed alpha radiation dose of Ra-223. In this situation, calculating the mean at the same point in time resulted in meaningless associations (Figure 6(B)). In an effort to interpret inter-individual variability in response in the various follow-up time periods, we shifted individual gene expression pattern to an ordinal scale so that expected gene expression events were aligned accordingly (Figure 6(C)). After this alignment of individual time patterns, significant exposure-to-gene expression association common to all observed patients were shown (Figure 6(D)). These individual patterns of DGE onset were observed for most genes among the patients. We cautiously interpret this finding as support for the application of gene expression measurements to detect internal radionuclide exposure. However, there is considerable inter-patient variability in bone marrow absorbed doses (Chittenden et al. 2015). Being aware that there are genes, which are underlying circadian variations, we checked the gene annotations derived from literature, whether our candidate genes have been previously described to be dependent on the time of the day. Hereby, no association could be found for the 15 candidate mRNAs, although all blood samples were taken in the morning.

To address the issue, whether changes in GE might serve as an easy method to detect incorporated radionuclides in human peripheral blood for improved internalized radionuclide diagnostic with the potential of high throughput diagnosis, the current study is a rare human in-vivo study dealing with this. Previously, internalized radionuclide studies were set up, using mainly murine models for examination of ⁹⁰Sr and ¹³⁷Cs incorporation (Ghandhi et al. 2015, 2020; Shuryak et al. 2020). Dose-rate-independent and dose-rate-dependent GE responses to incorporated radionuclides could be found in these published studies, although the GE patterns appeared very complex in all incorporation studies, most likely resulting from continued irradiation at a decreasing dose rate.

Conclusions drawn from our work are limited given the advanced disease stage of our study group. Each patient entered the study in very different health conditions, with different comorbidities, previous therapies and current concomitant treatment (exclusion criteria: application of another radiopharmaceutical for treatment and other malignancies, both neither previous nor simultaneous). While we chose stringent inclusion criteria to increase similarities in the patients, we admit this was difficult to achieve. Measurements in patients who are heavily diseased may also cause individually delayed effects and may mask the underlying association between incorporated alpha-emitting radionuclides and gene expression changes. It is unclear whether gene expression patterns over time might be more similar in healthy patients, which is a prerequisite for application of this methodology for improved internalized radionuclide diagnostic as a biodosimetry tool.

Due to the small sample size, our study might more rightfully considered as a ‘proof of concept’ effort where we explored the potential of easily identified RNA species in the peripheral blood to determine the incorporation of alpha-emitting radionuclides. Future work will require larger sample sizes and results, especially the 15 promising mRNAs described, should be validated using other in-vivo models and different radionuclides.

In summary, our results support a possible application of gene expression measurements on a transcriptional level (mRNAs) for detection of internalized alpha-emitting radionuclides, but further research is required. These relationships could have been masked by a skewed time-shift caused by inter-individual variability in response, necessitating an adjustment of the analysis strategy. This study highlights the difficulties associated with enrolling patients with severe disease and the quality of inference in such a situation.

Author contributions

Patrick Ostheim and Michael Abend wrote the main manuscript text and prepared the figures, tables and supplement table. Tim Nestler, Matthias Miederer and Manuela A. Hoffmann contributed to the conception of the study. Alexis Rump, Michael Lassmann and Uta Eberlein performed the dose reconstructions. Alexis Rump, Michael Lassmann, Uta Eberlein, Matthias Miederer, Mathias Schreckenberger, Tim Nestler, Manuela A. Hoffmann, Vahe Barsegian, Matthaeus Majewski and Matthias Port revised the work. All authors reviewed the manuscript.

Acknowledgements

We appreciate the sophisticated and carefully performed technical assistance of Sven Doucha-Senf, Thomas Müller and Oliver Wittmann.

Data availability statement

The datasets generated and analyzed in the course of the current study are available from the corresponding author on request.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the German Ministry of Defense.

Notes on contributors

Patrick Ostheim

Patrick Ostheim, MD, is a Post-Doctoral Researcher of Radiobiology and a Resident in Radiology at the Bundeswehr Institute of Radiobiology, Munich, Germany.

Matthias Miederer

Matthias Miederer, MD, is a Senior Physician and Senior Scientist at the Clinic and Polyclinic for Nuclear Medicine, University Medical Center of the Johannes Gutenberg University, Mainz, Germany.

Mathias Schreckenberger

Mathias Schreckenberger, MD, is a Professor of Nuclear Medicine and Head of the Clinic and Polyclinic for Nuclear Medicine, University Medical Center of the Johannes Gutenberg University, Mainz, Germany.

Tim Nestler

Tim Nestler, MD, is an Assistant Professor of Urology and Consultant at the Department of Urology, Federal Armed Services Hospital Koblenz, Koblenz, Germany.

Manuela A. Hoffmann

Manuela A. Hoffmann, MD, is a Senior Scientist of Nuclear Medicine and Occupational physician at the Clinic and Polyclinic for Nuclear Medicine, University Medical Center of the Johannes Gutenberg University, Mainz, Germany and Deputy Branch Chief at the Federal Ministry of Defense, Department of Occupational Health & Safety, Bonn, Germany.

Michael Lassmann

Michael Lassmann, PhD, is a Professor of Medical Physics and senior scientist at the Department of Nuclear Medicine, University of Würzburg, Germany.

Uta Eberlein

Uta Eberlein, PhD, is a Post-Doctoral Researcher and Physicist at the Department of Nuclear Medicine, University of Würzburg, Germany.

Vahe Barsegian

Vahe Barsegian, MD, is a Professor of Nuclear Medicine and Head of the Institute of Nuclear Medicine, Helios Kliniken, Schwerin, Germany.

Alexis Rump

Alexis Rump, MD, is a Senior Scientist of Radiobiology and Anesthesiologist at the Bundeswehr Institute of Radiobiology, Munich, Germany.

Mattháus Majewski

Matthäus Majewski, MD, is a Resident in Urology at the Department of Urology, Armed Services Hospital Ulm, Germany.

Matthias Port

Matthias Port, MD, is a Professor of Radiobiology and Internal Medicine and Head of the Bundeswehr Institute of Radiobiology, Munich, Germany.

Michael Abend

Michael Abend, MD, is a Professor of Radiobiology and Deputy Head of the Bundeswehr Institute of Radiobiology, Munich, Germany.