Procedures Common for All Assays

Blood sampling in heparinized vials (2–3 ml whole blood was provided from one healthy 30-year-old female donor for calibration samples and one healthy 32-year-old male donor for blinded, coded samples), radiation exposure (X-ray source, Yxlon, Hamburg, Germany), distribution of calibration samples (to all teams except the team performing DNA microarrays) and blinded, coded samples to participating laboratories followed procedures as described in detail in the inter-assay comparison manuscript, of our coordinated article published in this special issue. Comparison between teams was done descriptively using three levels of granularity: (1) precise dose estimates, (2) dose estimates falling within the ±0.5 Gy uncertainty interval for reference doses <2.5 Gy and ±1.0 Gy for reference doses >3 Gy, as introduced for triage biodosimetry (25), and (3) allocation of reported dose estimates to reference doses representative for clinically relevant groups.

Sample Preparation and Gene Expression Assays Using qRT-PCR or Microarrays

The Bundeswehr Institute of Radiobiology (BIR) team incubated the irradiated blood samples for 12 h at 37°C, 5% CO₂ (Heraeus BBD 6220) using whole blood, and an equal volume of antibiotic free Roswell Park Memorial Institute (RPMI) 1640 medium containing 10% heat-inactivated fetal calf serum. After washing the blood mixture with EL buffer, cells were lysed in β-mercaptoethanol-RLT buffer (QIAamp RNA Blood Mini Kit, Qiagen) according to the manufacturer's instructions and stored at –20°C overnight. Frozen samples were shipped by overnight courier service on wet ice under defined conditions according to United Nations Regulation 650.

The teams ran assays according to locally established protocols. Quantitative RT-PCR was run by seven teams using different (TaqMan, SYBR Green) chemistries and different sets of genes (Table 1). One team used a DNA microarray platform and either a combination of four genes or only one gene (FDXR) for dose estimation. Table 2 provides details for each team regarding RNA isolation, cDNA synthesis and PCR parameters, and details for performing microarrays following the MIQE criteria for publishing gene expression studies (26, 27).

TABLE 1

Overview of Participating Teams, Utilized Platforms, Number and Names of Genes or Gene Combinations Used, the Origin of Calibration Samples, and Further Details

TABLE 2

Overview of Methodological Details of Either qRT-PCR (Quantitative Reverse Transcription Polymerase Chain Reaction) or Microarrays Used by the Contributing Teams

TABLE 2

Extended.

After sample receipt, each team isolated total RNA and converted RNA to cDNA according to their protocols as shown in Table 2. Quantitative RT-PCR teams employed different kits for cDNA synthesis, different genes for normalization purposes (e.g., HPRT or 18S rRNA), different radiation-responsive genes either alone or combined (e.g., FDXR, DDB2, BAX, CDKN1A, GADD45A) and different amplification protocols and instruments (Table 2). For analyses, either normalized threshold cycles (Ct-values), relative threshold cycles (Cp-values) or gene expression fold-changes relative to the unexposed control of the calibration curve were employed. Ct-values and Cp-values represent different algorithms to determine a quantitative threshold cycle (Cq-value) along the exponential part of the amplification plot. For convenience, we use only the term Cq-value in this manuscript. Cq-values represent a surrogate for gene expression changes. They reflect inverse log₂-transformed RNA copy number changes so that, e.g., increased Cq-values (e.g., increase of Cq-values from 16 to 17) represent decreased RNA copy numbers and each additional Cq-value refers to a halving of RNA copy numbers. Seven teams used dose estimates based on 1–4 radiation-responsive genes. Five teams used a combination of 2-8 genes for dose estimation.

Some teams generated calibration samples in addition to the calibration samples sent by BIR. For instance, Team 4 (Table 3) obtained additional blood samples at two time points from a healthy 61-year-old male donor and generated two additional calibration sample sets, resulting in three reported dose estimates based on different calibration curves. These calibration samples were irradiated using a ¹³⁷Cs γ-radiation source with a dose rate of 0.356 Gy/min. Team 2 combined their own (X-ray irradiation, dose rate: 1.14 Gy/min) and BIR calibration samples into one calibration curve and, thus, one dose estimate was reported per blinded, coded sample (Table 3). Finally, different software tools such as Statistica 9.0 software (StatSoft), GraphPad Prism, version 6.04 (GraphPad Software, La Jolla, CA, USA), R 3.6.2, or Sigma Plot 14.4 (Jandel Scientific, Erkrath, Germany) were employed to fit the calibration curves and to perform dose estimation.

TABLE 3

The Table Depicts Team Contributions (from Left to Right) Regarding Employed Genes, Reported Dose Estimates per Reference Sample 1–3, Differences among Reported and Reference Dose-Values as well as the Summed Absolute Difference over all Reference Samples (SAD), a Correct (Yes) or Incorrect (No) Order of Dose Estimates (from Lowest to Highest) Corresponding to Three Dose Categories [Unexposed, Low (1.2 Gy) and Highly Exposed (3.5 Gy)], the Use of FDXR Gene Expression Changes for dose estimation, as well as the Report Time

TABLE 3

Extended.

Gene expression analysis using DNA microarrays was performed on the Agilent platform as described previously (22). 200 ng total RNA was transcribed into cDNA with an oligo-dT primer, followed by transcription to cRNA labeled with cyanine 3-CTP (Quick-Amp Labeling Kit, One-color, Agilent). cRNA purification was performed with the RNeasy Mini Kit (Qiagen) and dye incorporation and cRNA yields were measured with the NanoDrop-1000 spectrophotometer (Thermofisher). Labeled cRNA samples were applied on the DNA microarray slides (44k whole human genome, G4112F, Agilent). For hybridization, DNA microarrays were placed into a hybridization oven (Agilent) at 65°C for 17 h. After hybridization, DNA microarrays were washed, and slides were immediately scanned with the Microarray Scanner (G2505 B, Agilent) as recommended by Agilent. The preprocessing procedure and subsequent statistical analysis were applied separately using Agilent Feature Extraction Software Version 9.5.1. The processed signals (Agilent gProcessedSignals including background subtraction) were subsequently log₂-transformed and normalized via GAPDH. Dose estimation was performed using a 4-gene signature as well as one gene (FDXR) employing internal in-house calibration curves (11, 22, 23).

Further details on each method are provided in the  Supplementary Materials¹ (598-615_rare-199-06-07_s01.docx) ( https://doi.org/10.1667/RADE-22-00206.1.S1).

Statistical Methods

Differences (for examining over- or underestimation) as well as absolute differences (AD, allowing for quantitative comparisons per reference sample) of estimated doses to their corresponding reference doses were calculated considering the uncertainty intervals introduced for triage dosimetry as outlined above. Reported dose estimates lying within corresponding reference dose intervals are named “correctly” reported dose estimates within this manuscript. For a more aggregated comparison of the teams’ contributions, the differences between estimated and reference doses over the three reference exposures were summed and are called summed absolute difference (SAD). In addition, allocation of reported dose estimates to reference dose groups corresponding to clinically relevant groups was examined. These groups comprised unexposed (no clinical resources required), low dose exposed (1.2 Gy group, no hospitalization required) and high dose exposed individuals (3.5 Gy group, immediate and intensive care required). Group comparisons were performed using either the t-test or Mann-Whitney Rank sum tests, where applicable. All calculations were performed using either SAS (release 9.4, Cary NC) or Excel (Microsoft Corporation). Graphs were created using Sigma Plot 14.5 (Jandel Scientific, Erkrath, Germany). Linear or linear-quadratic curve fits were performed where applicable (Sigma Plot 14.5).  Supplemental File 1 (598-615_rare-199-06-07_s01.docx) includes more statistical details for each team where required.

Participating Institutions, Contributions and Reported Dose Estimates

Table 1 provides an overview on the platforms, genes and calibration samples used by the eight teams. Further methodological details are summarized in Table 2. Team numbers in Table 1 are the same throughout all figures and tables to better follow results and workflow characteristics. Eight teams (Table 1) provided a total of 105 dose estimates (Table 3). For instance, Team 4 provided dose estimates based on four genes (GADD45A, FDXR, CDKN1A, MDM2) and their geometric mean, using three different calibration curves. This generates five dose estimates per blinded, coded sample and calibration curve and a total of 45 (5x3x3) reported dose estimates (Table 3).

Team 8 using a DNA microarray platform estimated the dose based on a previously developed algorithm using a gene set consisting of 4 genes (TNFSF4, FDXR, PHLDA3, THC2705989) (expressed sequence tag with a yet unknown function) and FDXR only (22, 23). No further calibration samples were required. Instead, internal calibration curves of the four signature genes were used. The derived expression values of the four signature genes in sample 1 (0 Gy) of the blinded, coded RNA samples were similar to or slightly below the values of the non-irradiated controls of the internal in-house calibration curves. The two contributions (4-gene signature and FDXR alone) resulted in two dose estimates per blinded, coded sample and a total of 6 dose estimates (Table 3).

The earliest report times of categorical and precise dose estimates were around 5 h and 9 h, respectively (Table 3). Further dose estimates were reported within 3 days and five teams provided dose estimates later. The median report time was about 10 days (244 h) and included groups that did not process these samples with priority.

Methodological Accuracy of qRT-PCR

The accuracy of qRT-PCR gene expression measurements was examined using the technical replicates run in duplicate/triplicates and/or performed by two experts. In the case of microarrays, no replicate measurements of calibration samples were performed. The coefficient of variation (CV, standard deviation relative to the mean RNA copy number) for 84.6% and 58.9% of all 436 measurements did not exceed 10% and 5%, respectively (see horizontal dotted lines in Fig. 1). Lower CVs (<5%) were observed for most measurements produced by two of the teams (left side, Fig. 1). Three other teams (middle, Fig. 1) produced CVs for most genes ranging between 5–10%, while CVs exceeding 10% were calculated for one team (right side, Fig. 1). FDXR was the only gene used either alone or in combination by all teams (Fig. 1, Table 3).

FIG. 1

The coefficient of variation (CV) was calculated based on all available technical replicates and genes. The distribution of CVs per gene is reflected by a box plot showing the 10th (lower whisker), 25th (lower end of the box), 50th (median, straight line), 75th (upper end of the box) and 90th (upper whisker) percentiles. Genes are ordered with increasing raw Cq-values. Teams are numbered as outlined in Table 1.

Calibration Curves and Reported vs. Reference Dose Estimates

Teams plotted gene expression values versus dose in different ways. For instance, two teams plotted normalized Cq-values of a combination of genes and aggregated values either by applying a geometrical mean or building a sum (Fig. 2A and B). In most cases, teams plotted normalized Cq-values of single genes (e.g., FDXR) vs. dose either using linear scales (Fig. 2C and D) or log-scales (Fig. 2E). Plotting Cq-values and dose on log-scales facilitated a linear regression model covering most of the variance (r² = 0.99) as opposed to other approaches where a linear regression model diverged from the calibration values (e.g., r² = 0.59, Fig. 2C). For normalized Cq-values increased copy numbers of radiation-induced genes resulted in a negative slope of the calibration curves. When using unexposed calibration samples as the reference (value = 1), differential gene expression increased with increasing radiation dose, resulting in a positive slope of the calibration curve (Fig. 2F and G). Unexposed dose estimates (referred to as “1” in Fig. 2) in some plots appear above the calibration curve (e.g., Fig. 2A and E). They are referring to negative dose estimates. To avoid extension of the x-axis into the negative dose range, these samples are plotted above the calibration curve and are referring to the 0 Gy (linear scale) or 0–0.01 Gy dose range (log scale). A saturation in gene expression values was observed between 3–4 Gy for most calibration curves or even at 2 Gy (Fig. 2A).

FIG. 2

Calibration curves were generated before the exercise by seven teams and those shown herein represent typical examples. Either multiple (panels A and B) or single gene (panels C–G) expression changes are plotted versus the known radiation doses. Gene expression values are often given as raw or normalized Cq-values or as a fold change (panels F and G) using unexposed samples as the reference. The presented formula in panel 1 (Team 1) refers to the gene expression response after 0–2 Gy and excludes higher doses where the response saturates. Further details on fitting of data are presented in  Supplementary Materials S1 (598-615_rare-199-06-07_s01.docx) ( https://doi.org/10.1667/RADE-22-00206.1.S1). Circles represent mean gene expression values from technical replicates. Error bars represent standard deviation of duplicate measurements (and standard error of mean in panel F) and are visible when larger than the symbols. Reported dose estimates are plotted as squares and numbers 1, 2 and 3 refer to unexposed, 1.2 Gy and 3.5 Gy irradiated blinded, coded samples, respectively. Unexposed (no. 1) dose estimates are sometimes not plotted on but above the calibration curve, referring to negative dose estimates. To avoid extension of the x-scale into the negative dose range, these samples are plotted above the calibration curve and are referring to the 0 Gy (linear scale) or 0–0.01 Gy dose range (log scale). Teams are numbered as outlined in Table 1.

The unexposed control sample (sample 1) was correctly estimated to be in the range between 0 and 0.5 Gy by all teams except Team 4, where 2 out of 12 dose estimates were not in this range (Fig. 3, Table 3). Differences to the reference dose were zero in all except two dose estimates which belonged to Team 4 (Fig. 3A and B). Reported dose estimates of 1.2 Gy (Fig. 3C and D) and 3.5 Gy (Fig. 3E and F) irradiated samples were systematically overestimated by Team 4, while other teams provided dose estimates close to the reference dose as well as over- and underestimations with a trend to an underestimation observed at 3.5 Gy exposures in three teams (Fig. 3 F).

FIG. 3

Reported dose estimates are provided per team (teams are numbered as outlined in Table 1) on the left panels and corresponding differences of reported versus reference doses are shown in the right panels. The panels A, C, and E reflect reported dose estimates corresponding to the unexposed, 1.2 Gy and 3.5 Gy irradiated samples, respectively. Panels B, D, and F reflect corresponding differences of reported versus reference doses. Dotted lines refer to the ±0.5 Gy (reference doses below 2.5 Gy) and ±1 Gy (reference doses >3 Gy) uncertainty interval for dose estimation of triage dosimetry. Solid lines reflect values of the reference dose. Circles represent reported dose estimates and squares reflect mean values with error bars showing the standard deviation. Results are provided in descending order of mean reported dose estimates. Stars correspond to dose estimates based on gene expression changes using FDXR.

Altogether, 50% of the reported dose estimates were within the accepted uncertainty interval for triage dosimetry regarding the 1.2 Gy and 3.5 Gy irradiated blinded, coded samples (see bold dose estimates and calculations in lower part of Table 3). The dose estimates from all teams correctly identified the ordering with regard to the dose effect relationship of the blinded, coded samples (lowest, middle, highest exposures), using all reported and semiquantitative dose estimates (Table 3).

Reported dose estimates using contributions from all teams revealed overlapping dose estimates among unexposed, low, and highly exposed groups, which hindered their discrimination (Fig. 4 A). Overlapping reported dose estimates between exposure groups were reduced after exclusion of Team 4's contributions (Fig. 4 B) and using reevaluated dose estimates generated from newly synthesized cDNA of BIR calibration samples (Fig. 4 C). This resulted in reported dose estimates that identified unexposed and highly exposed samples with 100% certainty, meaning that if the dose estimate predicted 0 Gy or 3.5 Gy, in all cases the corresponding samples belonged to the unexposed or highly exposed clinical group, respectively (Fig. 4 B and C, arrows). Reported dose estimates below 3.5 Gy could not discriminate between the low (1.2 Gy) and the highly (3.5 Gy) exposed group in most cases (Fig. 4 B and C).

FIG. 4

Reported doses estimates of all teams and corresponding dose categories (x-scale) are shown in panel A. Panel B reveals the same data, excluding dose estimates from Team 4. Panel C comprises reported dose estimates of all teams as well as reevaluated dose estimates from Team 4 generated on normalized Cq-values of newly synthesized cDNA from BIR calibration samples as provided in Fig. 5B. Circles represent reported dose estimates and squares reflect mean values with error bars showing the standard deviation. Arrows in panels B and C are related to unexposed and highly exposed groups. They reflect reported dose estimates (0 Gy and ≥3.5 Gy, represented by horizontal dashed lines), where dose groups are identified correctly and without misclassification.

Team 4 Contributions

Team 4 experienced a systematic shift of dose estimates towards higher values (Table 3, Fig. 5A, B), so these results were examined more closely. Dose estimates were derived using either BIR calibration samples or two additional sets of calibration samples generated by Team 4, which differed from the BIR calibration samples regarding the radiation exposure, as well as donors age and sex (Fig. 5C). Reevaluation of FDXR-based dose estimates performed after the exercise by the same team were 0 Gy, 2.2 Gy and 2.8 Gy (using the raw Cq-for the reported BIR calibration samples), 0 Gy, 2.0 Gy and 3.0 Gy (using new cDNA of BIR calibration samples), and 0 Gy, 2.3 Gy and 3.6 Gy (using normalized Cq-values from re-run reported BIR calibration samples) corresponding to the 0, 1.2 and 3.5 Gy irradiated blinded, coded samples (Fig. 5).

FIG. 5

Here we show Team 4 data. Panel A reflects FDXR gene expression changes plotting raw (left graphs) and normalized Cq-values (middle graphs) of reported dose estimates using either BIR (Bundeswehr Institute of Radiobiology) calibration samples (top panels), or two sets of calibration samples generated by Team 4, termed Team A and Team B (bottom panels), respectively. Panels to the right represent reevaluated data generated on newly synthesized cDNA of RNA extracted from calibration samples as indicated. Dashed horizontal lines refer to Cq-values corresponding to blinded, coded samples no. 1 (0 Gy), no. 2 (1.2 Gy), and no. 3 (3.5 Gy) as indicated in the panels. Panel B shows dose estimates generated both during (left) and after the exercise (right) based on the different calibration curves used by Team 4 for FDXR, GADD45A, CDKN1A and MDM2 genes. “Reported dose estimates” corresponds to dose estimates determined during the exercise, and dose estimates based on raw and normalized Cq-values are shown. “Reevaluated dose estimates” include dose estimates generated after the exercise. It comprises raw and normalized Cq-values generated from remaining cDNA as used for the reported dose estimates as well as data generated from newly synthesized cDNA. The table in panel C reflects exposure, age and sex differences between the different sets of calibration samples employed as well as a calculation of fold-differences among the four genes examined, which were generated at the dose which created the maximal fold change from 0–4 Gy irradiated calibration samples employing the 0 Gy Cq-values as the reference.

Reasons for this discrepancy were most likely technical issues, either pipetting error, primer dilution and/or use of defective mastermix for certain samples, which impaired the accuracy of the 18S rRNA Cq-values obtained during the reporting of dose estimates, and, in turn, had a negative impact on reported normalized Cq-values. This was concluded since dose estimates using raw Cq-values from reported dose estimates provided better dose estimates than after normalization (Fig. 5A, B), and that 18S rRNA Cq-values obtained within the exercise for blind samples no. 2 and no. 3 were higher compared to additional runs performed after the exercise. This in turn resulted in lower reported normalized Cq-values, and consequently higher dose estimates as compared to those observed using newly synthesized cDNA, where both blind coded samples no. 2 and no. 3 did intersect the BIR calibration curve (Fig. 5A). Additionally, re-run of qPCR using the same cDNA used for the reported FDXR-based dose estimates using BIR calibration samples and Team 4 calibration A samples (insufficient material to re-run calibration B samples), resulted in similar dose estimates to those derived from the newly synthesized cDNA after the exercise for all calibration samples, and also from raw Cq-values from reported data in the case of BIR samples (Fig. 5A). Moreover, 1:2 dilution of 18S rRNA leading to raw Cq-values of 9 seem to represent Cq-values within the linear dynamic range of the method and do not explain the observed discrepancies ( Supplementary Materials (598-615_rare-199-06-07_s01.docx);  https://doi.org/10.1667/RADE-22-00206.1.S1). Also, dose estimates based on the GADD45A, CDKN1A and MDM2 genes improved after the exercise using either raw Cq-values for reported or new cDNA (Fig. 5B). A lower magnitude of gene expression, i.e., higher Cq-values, with early saturation at approximately 1–2 Gy was observed for the additionally generated two calibration sample sets of Team 4 resulting in higher dose estimates as compared to those derived from BIR calibration samples, which presented a later saturation at approximately 4 Gy (Fig. 5C). Additionally, fitting the curvi-linear dose-to-Cq-value-association with a linear regression model introduced deviations from the reference doses (Fig. 2C), yet the use of a linear quadratic model to fit the reported normalized Cq-values (Fig. 5A) would have made the dose estimation of blind sample no. 3 impossible. Finally, the Cq-value difference calculated within the 0–4 Gy dose band of irradiated calibration samples converted into fold changes in gene expression of 3–4 fold for GADD45A and MDM2, 3–14 fold for CDKN1A and 13–23 fold for FDXR, making FDXR the most promising candidate gene, based on increased robustness in the detection of high fold changes (Fig. 5C). Also, BIR calibration samples provided significantly higher fold changes for all genes compared with the two calibration sets generated by Team 4 (P = 0.01).

Using Single Genes, FDXR, and Combined Gene Signatures for Dose Estimation

FDXR was employed by all teams and 17 out of 35 dose assessment approaches employed FDXR either alone or combined with other genes (Table 3). Altogether 78.1% of correctly reported dose estimates made use of FDXR gene expression measurements (see calculations in the lower part of Table 3). Also, reported dose estimates (excluding Team 4) based on a gene combination provided more precise dose estimates (median SAD of 0.8 Gy, ±0.8 Gy) than those employing a single gene, e.g., CDKN1A or DDB2 (median SAD of 1.7 Gy, ±0.8 Gy) or FDXR only (median SAD of 1.8 Gy, ±0.8 Gy, Fig. 6). Interestingly, when reported separately, difference of FDXR based reported versus reference doses (marked as stars in Fig. 3F) was consistently above 0 for the 3.5 Gy reference sample. The lowest FDXR based SAD of 0.5 Gy was accomplished by Team 7 with an accuracy comparable to teams combining several genes for dose estimation (Table 3, Fig. 6).

FIG. 6

The summed absolute difference (Gy) was plotted depending on the application of single genes for dose estimation (FDXR only, other single genes, left and middle part of the graph) or the use of a combined set of genes (right side). Circles represent reported dose estimates and squares are reflecting median values with error bars showing the standard deviation.

Testing for Improvements in Dose Estimations Using Remaining RNA and cDNA Samples when Employing the Same Workflow (After-Exercise Experiments)

Six teams entrusted remaining isolated RNA (n = 5) or synthesized cDNA (n = 3) of calibration curves and blinded, coded samples for further study following the workflow of Team 7 (Fig. 7) to elucidate the impact of a similar workflow on the accuracy of dose estimates. Employing the same workflow (from Team 7) to the 1.2 Gy and 3.5 Gy reference samples decreased the absolute differences (ADs) relative to the employment of different workflows in three of the entrusted RNA/cDNA samples, while increasing the ADs in three other RNA/cDNA samples (Fig. 8, underlying details are shown in  Supplementary Fig. S1 (598-615_rare-199-06-07_s01.pptx);  https://doi.org/10.1667/RADE-22-00206.1.S2). Improved ADs were more pronounced at 3.5 Gy. Here, ADs employing the same workflow were roughly halved over RNA/cDNA samples processed using different workflows. In contrast, the 1.2 Gy and 3.5 Gy blinded, coded RNA samples originating from Teams 3 and 6 did not gain, i.e., SADs became larger, from the employment of the same workflow (Fig. 8).

FIG. 7

Overview of samples prepared at BIR (Bundeswehr Institute of Radiobiology), sent to the participants, and processed by the Teams no. 1–7 starting with isolation of RNA after the arrival of lysed peripheral blood aliquots (upper row). The reported dose estimates were generated based on these different workflows of the teams. To further elucidate the origin of differences in reported dose estimates among teams, either cDNA (second row) or isolated RNA (third row) was entrusted by six teams to the team showing reported dose estimates closest to the reference doses. The same workflow and its impact on dose estimates was applied to entrusted cDNA (middle row) and RNA (lower row).

FIG. 8

The absolute difference (Gy) was compared per team employing either the different workflows established at Team 1–7 laboratories (white bars) or applying the same workflow for entrusted RNA or cDNA samples by one team (gray filled bars). FDXR gene expression changes after a reference dose of 1.2 Gy (upper panel) and 3.5 Gy (lower panel) were examined. Dose estimates that improved after applying the same workflow are shown on the left (labeled by original team number as outlined in Table 1), and those where dose estimates worsened are shown on the right side of the graph.