Adineta vaga under fire: simulating the impact of radiation

Abstract

Previous studies have demonstrated the remarkable resistance of bdelloid rotifers to ionizing radiation, making them an interesting model system for studying radiation effects on living organisms. In this study, we use simulations, instead of direct experimental exposures, to examine whether all bdelloids are affected equally by radiation exposure and to explore the relationship between biological data and energy deposition patterns induced by low and high linear energy transfer (LET) radiation. To this end, tool for particle simulation (TOPAS) a simulation tool, widely used in the field of medical physics and radiation therapy, was utilized. Using simulations for proton, iron ions, and X-ray exposure, our findings showed that all individuals, cells, and nuclei were effectively hit by the administered doses of 4 MeV protons, 0.5 GeV/n ⁵⁶Fe, and X-ray radiation. The results support that the impact on survival and fertility rate measured in Adineta vaga is caused by radiation-induced damage rather than the absence of hits in certain individuals or germinal cells. Notably, simulations revealed significant differences between low- and high-LET radiation concerning irradiated individuals' nuclei. Specifically, for an equivalent dose, high-LET radiation requires fewer incident particles compared to low-LET radiation, resulting in a sparser distribution of radiation hits on the nucleus surface. In conclusion, the study supports the idea that reduced fertility described in high-LET exposed samples is associated with complex DNA damage caused by the condensed energy deposition pattern of high-LET radiation compared to low-LET.

Introduction

With a worldwide distribution, bdelloids are among the tiniest metazoans on the planet. One remarkable feature of bdelloid rotifers, being even more unusual among animals, is the absence of males, vestigial male structures, or hermaphrodites in any of the populations studied within the 460 described morphospecies. Bdelloid rotifers are females apparently cloning themselves since millions of years (Maynard Smith, 1986; Judson & Normark, 1996; Segers, 2007). Approximately 90% of bdelloid species have been documented to inhabit semi-terrestrial habitats such as mosses and lichens (Ricci, 1987, 2017; Melone & Fontaneto, 2005). These habitats are known for sudden and sometimes rapid changes in physical and chemical conditions like temperature, food, chemical composition, and water availability. Adapted to temporary environments, most, but not all, bdelloids are able to enter and survive desiccation at any stage of their life cycle and resume reproduction after rehydration without negative effects. This desiccation resistance is shared with other small animals associated with limno-terrestrial environments like some nematodes and several tardigrade species (Ricci & Pagani, 1997; Rebecchi et al., 2006). It remains however unclear how desiccation resistant organisms deal with this unpredicted decrease in water content, a process that induces cellular damages (e.g., protein aggregation and denaturation), DNA damages, leakage of cellular membranes and generation of reactive oxygen species (ROS) (Potts, 1994; Mattimore & Battista, 1996; Shirkey, 2003; Alpert, 2005; Gusev et al., 2010; Crowe, 2014; Hespeels et al., 2014). Moreover, in order to remain viable, dried organisms must acquire mechanisms that maintain the function of their macromolecules (e.g., proteins, DNA) and membranes, and the capacity to restart their metabolism when water becomes available (Potts et al., 2005; Leprince & Buitink, 2015).

The notorious desiccation tolerance of bdelloid rotifers seems to provide them with a high tolerance to a variety of stresses including high pressure, vacuum, UV or freezing (Ricci et al., 2005; Fischer et al., 2013). Gladyshev and Meselson (2008) were the first to document the extraordinarily resistance of two bdelloid rotifer species to gamma rays. This finding was later confirmed by Hespeels et al. (2014), who demonstrated their high resistance to prolonged desiccation and high doses of different types of ionizing radiation (IR).

Exposure of cells to ionizing radiation triggers a complex chain of physical, chemical, and biological processes, varying in terms of timing, spatial distribution, and energy scale, making the cellular response challenging to predict (Joiner & van der Kogel, 2018). The ionization density in particle tracks is typically characterized by Linear Energy Transfer (LET), which refers to the average energy (in keV) deposited by a charged particle per unit distance traveled (1 μm). High-Z charged particles, such as Fe ions, are considered high-LET due to their high ionization density along their path, while photons, such as X-rays, produce sparse ionizations along their path, and are therefore considered low-LET ionizing radiation (Joiner & van der Kogel, 2018; Hagiwara et al., 2019). While the same radiation dose produces a similar quantity of ionizations, the differences in their spatial distribution result in diverse biological damages. It has been shown that low-LET radiation, with their sparse ionizations, primarily cause DNA base oxidation and single-stranded breaks (Lehnert, 2007). In contrast, exposure to high-LET radiation leads to the formation of complex DNA damages, including clustered double-stranded breaks (DSBs), which results in a reduced survival of the irradiated organisms (Semenenko et al., 2004). To document the energy deposited by ionizing radiation in the matter or cell, the absorbed dose is expressed in Gray (Gy). Mathematically, it was defined as the ratio of the energy deposited by the radiation to the mass of the matter traversed.

Desiccated Adineta vaga (Davis, 1873) individuals were able to survive doses higher than 5,000 Gy of X-ray and proton radiation, which was striking compared to the Lethal Dose 50 (i.e., the dose required to kill 50% of the irradiated population) of mammalian cells, which ranges from 2 to 6 Gy after X-ray exposure (Hespeels et al., 2020). Other studies have demonstrated the resistance of different bdelloid species to various types of radiation. Mniobia russeola (Zelinka, 1891) has been found to withstand heavy ions such as iron (Fe) and helium (He) (Jönsson & Wojcik, 2017), while Philodina roseola Ehrenberg, 1832 exhibits high resistance against gamma radiation (Gladyshev & Meselson, 2008). This radioresistance may appear obscure since none of the exposed radioresistant species were ever confronted with such levels of radiation in their natural environment. By analogy with the extreme desiccation and radiation resistant bacterium Deinococcus radiodurans (Anderson, 1956) (Mattimore & Battista, 1996), it has been hypothesized that the extraordinary radiation resistance of bdelloid rotifers reflects an adaptation to their desiccation resistance (Gladyshev & Meselson, 2008; Gladyshev & Arkhipova, 2010). Indeed, the extreme radiation resistance of D. radiodurans, surviving gamma radiation doses higher than 20,000 Gy, was attributed to their ability to deal with massive amount of DNA DSB accumulated during desiccation. Indeed, a similar impact on genomic integrity was recorded in samples of D. radiodurans desiccated for 6 weeks or irradiated with 5,200 Gy of gamma radiation. Moreover, D. radiodurans was able to restore genomic integrity post rehydration or irradiation (Mattimore & Battista, 1996). Similarly, it has been demonstrated that prolonged desiccation or exposure to high dose of IR generates hundreds of DNA DSB in nuclei of bdelloids that are fully or partially repaired post rehydration or irradiation (Hespeels et al., 2014, 2020). Indeed, A. vaga individuals were able to restore the genomic integrity of their germinal cells after high doses to IR (> 500 Gy) (Terwagne et al., 2022). On the other hand, an incomplete DNA DSB repair was observed in their somatic cells suggesting that genetic alterations at the sequence level may be tolerated in somatic cells due to their non-cycling, postmitotic status or the absence of signaling for unrepaired DNA DSBs (Hespeels et al., 2020; Terwagne et al., 2022).

The radio-resistance of A. vaga individuals was differentially interpreted, in term of magnitude, when focusing on the survival rate or on their capacity to produce viable offspring, and whether it was low- or high-LET radiation exposure. In a previous study (Hespeels et al., 2020) we reported a factor of approximately 3 between doses required to sterilize a population or to kill all exposed animals. For example, dose required to sterilize 100% A. vaga population with X-ray is between 1,500 and 2,000 Gy while it is between 5,000 and 7,500 Gy to kill the entire population (Hespeels et al., 2020). The median sterilizing doses (SD50) for A. vaga individuals were 1,035 ± 20 Gy, 453 ± 23 Gy, and 461 ± 1 Gy for X-ray, protons, and iron radiation exposure, respectively. A 2.3 factor was reported when comparing SD50 of X-ray versus protons and Fe particles. Indeed, it would have taken twice as much dose deposited in the case of X-rays to sterilize individuals in comparison with protons and iron particles (Hespeels et al., 2020). These observations, also reported for tardigrades exposed to IR, were hypothesized to be attributed to the nature of the cells evaluated through survival or fecundity rate assays (Pagani et al., 1993; Beltran-Pardo et al., 2015). The survival rate in those metazoans is indeed related to non-dividing somatic cells, while fertility data is linked to germinal cells that will ensure cellular multiplication during embryonic development and as consequence are more sensitive to damages induced by radiation. The precise causes of these differences remain to be investigated but may be informative to the documentation of damage type induced by low- and high-LET radiation on a radioresistant metazoan.

Until recently, DNA damage quantification of irradiated bdelloid individuals relied mainly on the semiquantitative approach PFGE (Pulsed Field Gel Electrophoresis) done on thousands of individuals. However, this technique was not able to provide information related to each type of cell and was biased toward damages affecting somatic nuclei of bdelloids (Gladyshev & Meselson, 2008; Hespeels et al., 2014, 2020). In a new study, Terwagne et al., reported, using the TUNEL assay approach, the presence of DNA DSBs in primary oocytes of A. vaga individuals exposed to the sterilization dose of 1,280 Gy of proton radiation (Terwagne et al., 2022). It remains plausible that the survival or fertility curve reported previously may be affected by the presence of untargeted animals or germinal cells during irradiation. A simulation approach would be best suited to investigate this hypothesis.

TOPAS (TOol for PArticle Simulation) is a powerful and versatile simulation tool widely used in the field of medical physics and radiation therapy research. It features a user-friendly interface and robust capabilities, enabling users to model and simulate complex scenarios, accurately predicting the behavior of particles as they interact with various materials and biological tissues. TOPAS leverages the Geant4 Monte Carlo toolkit, a widely adopted simulation package in high-energy physics, to perform detailed and accurate particle transport simulations, including electromagnetic, hadronic, and nuclear processes, with high fidelity (Agostinelli et al., 2003; Allison et al., 2006, 2016; Perl et al., 2012; Faddegon et al., 2020). The versatility of TOPAS extends beyond radiotherapy applications, as it can also be utilized in other fields, such as radiation biology, detector development, and space radiation research. It provides a flexible framework for customization, allowing users to incorporate their own models, materials, and scoring algorithms, thus adapting the simulation environment to their specific research needs. For instance, Geant4 has been recently utilized to simulate proton imaging of microscopic samples (Michelet et al., 2019). A numerical phantom of Caenorhabditis elegans (Maupas, 1900) was derived from experimental data and used to reconstruct scanning transmission ion images to determine the minimum number of protons to use per shot. In another study, Geant4 enabled the accurate dosimetry of charged particle microbeams for the irradiation of 2 cell-stage embryos of C. elegans (Torfeh et al., 2019). The energy deposit in various chromatin condensation states was investigated during the cell division phases, allowing subsequent irradiation experiments to be reproducible in term of dose. Although a few simulations studies have been performed on small invertebrates and crustaceans (EL Bakkali et al., 2022), this study is the first to use TOPAS to determine the irradiation pattern from low- and high-LET particles on bdelloid rotifers.

Here, we used TOPAS, to reevaluate how radiation sources target and hit A. vaga individuals at population level, but also at nucleus level. This reevaluation is based on data we previously obtained (Hespeels et al., 2020). In other terms, we sought to determine which fraction of individuals and cells are effectively targeted by radiation. These data are critical to ensure a good comprehension of bdelloid radiation data and exclude the possibility that the high radiation resistance reported for bdelloid rotifers is linked with untargeted individuals/cell types. Additionally, using simulations, we investigated whether biological data could be associated with specific energy deposition patterns at DNA level induced by low- and high-LET radiation. To evaluate the validity of these hypotheses, we modeled all three radiation types (i.e., proton, ⁵⁶Fe ions, and X-rays) previously applied on desiccated A.vaga in Hespeels et al., 2020.

Materials and methods

TOPAS and modeling of bdelloid samples

The software used is TOPAS MC Inc, version 3.7. TOPAS is an extension of Geant4, a Monte Carlo simulation software, designed for the transport and interactions of charged and uncharged particles in matter (Agostinelli et al., 2003; Allison et al., 2006, 2016; Perl et al., 2012; Faddegon et al., 2020). Whole sets of electromagnetic, hadronic, and optic processes called physics lists can be simulated over a wide range of particles and energies. The Livermore list is chosen as it is specialized in weak energy interactions, from 100 GeV down to 250 eV.

The next step involves defining the geometry, materials, and radiation sources. Simulation data were provided in Supplementary Material 1. Three levels of granularity were utilized to replicate experimental setups from previous research (Hespeels et al., 2020). At the first level, a layer of desiccated bdelloids is modeled and subjected to irradiation to determine the dose distribution across individual bdelloids. Figure 1a compares a SEM image of the layer and the modeled geometry. At the second level on Fig. 1b, a single desiccated bdelloid is modeled and exposed to radiation to determine the dose distribution across all its nuclei. To provide a detailed view of a bdelloid, Fig. 1b displays a close-up SEM image of the bdelloid layer, where the outline of a single bdelloid is highlighted in red. Figure 1c corresponds to the third level, the nucleus. a confocal imagery shows somatic and germinal cells of a hydrated A. vaga individual.

Fig. 1

View of the three modeled levels of granularity simulated in TOPAS: monolayer of desiccated bdelloid, bdelloid individual and nucleus level. The modeling is based on SEM reference images. The layer level is a 14 × 14 mm² large and 30 µm thick cube with a 14 mm diameter and 30 µm thickness on average. An overview of a layer of desiccated bdelloid was illustrated with a SEM image. It is divided into 140 × 140 bdelloid voxels that correspond to the desiccated individuals. The bdelloid is 100 × 100 µm² large and 30 µm thick, from the dimensions of a bdelloid individual highlighted in red on the SEM image. It is subdivided into 50 × 50 nuclei voxels. The nucleus level is 2 × 2 µm² large, 30 µm deep and divided into 50 × 50 voxels to provide a DNA level resolution of damage location. The third level confocal imagery shows somatic and germinal cells of hydrated A. vaga individuals. Nuclei were counterstained with DAPI (blue). Arrows indicate the pools of primary oocytes. Scale bar: 10 µm. Picture from Terwagne et al., (2022)

Bdelloid layers are aggregates of desiccated bdelloids as used in Hespeels et al. 2020. Experimentally, the desiccation process takes place on an agarose support of approximately 100 µm thickness. A bdelloid layer has on average a 14 mm diameter and 30 µm thickness (based on the SEM images as in Fig. 1a). The thickness is subject to biological variability such as impurities or overlap of bdelloid individuals. Bdelloid layer is nevertheless modeled as a cube of 14 × 14 mm² with 30 µm thickness. Indeed, a cube can be segmented into cubic volumes. Thus, the layer can be segmented in 140 × 140 voxels of 100 × 100 µm². This corresponds roughly to the dimensions of a desiccated bdelloid (Terwagne et al., 2022). Similarly, for the second level of simulation, a single bdelloid of 100 × 100 µm² surface is segmented in 50 × 50 voxels of 2 × 2 µm², which corresponds to the order of magnitude of a nucleus surface (Terwagne et al., 2022), thus enabling nuclei level resolution. The third level models a single nucleus with 2 × 2 µm² surface. It is segmented in 50 × 50 voxels of 40 × 40 nm², which is convenient to investigate damage repartition at the nucleus level. TOPAS uses user-defined materials. The exact composition of bdelloids is unknown, but the chemical composition of human skin (H₁₀C₂₀N₄O₆₅S_0.2Cl_0.3Na_0.2K_0.1) was selected as giving the best approximations of the experimental LET.

Experimental setups and modeled geometries

The geometries modeled in TOPAS accurately replicate the experimental setups used in reference (Hespeels et al., 2020). All geometry elements are contained in a 1 × 1 × 1 m³ virtual cube of vacuum called world. The world cube is a computational boundary beyond which physical processes are not simulated. Three experimental setups are modeled: 4 MeV protons, 0.5 GeV/n ⁵⁶Fe ions, and X-rays irradiation.

4 MeV proton

A layer of bdelloids is irradiated by a 1 cm² beam of 4 MeV protons, placed at 1 m ahead of a Si3N4 window. The beam is produced in vacuum and extracted in air through the 1 µm Si3N4 window. It then travels through 3 mm of air to reach the bdelloid layer on its agarose support. The doses received by the layer ranged experimentally from 100 to 10,000 Gy (Hespeels et al., 2020). The surface of the beam is larger than the sample to provide uniform irradiation. Preliminary simulation has shown that the presence or absence of the agarose support (located after the bdelloid layer) has no impact on simulation outputs. No backscattered protons have been observed, and the contribution of backward emitted secondary electrons created within the agarose layer is negligible (data not shown). It has therefore been removed to speed up the calculations. The simulated geometry is represented on Fig. 2a.

Fig. 2

Proton (a) and ⁵⁶Fe (b) irradiation setups and PhaseSpace scorer (c). Figure 2a depicts the proton irradiation setup. A 14 × 14 mm² proton beam is placed in vacuum, 1 m away from the setup. Protons are extracted from the vacuum through a 1 µm thick window of Si₃N₄. After traversing 3 mm of air, the protons reach the bdelloid layer. This setup is identical for X-ray, with the only difference being the type of particle used. In b, a 14 × 14 mm^{2 56}Fe beam is placed in vacuum, 1 m ahead of the setup. The beam passes through a 0.25 mm thick layer of polyethylene (PCR tube) before crossing 1 mm of air and reaching the bdelloid layer. Figure 2c illustrates the use of the PhaseSpace scorer. During the layer-level simulation (left), a virtual surface with no interactive properties is created upstream of the selected bdelloid voxel. The scorer is attached to this surface and registers all particles that pass through it. The scorer can be used as source in subsequent simulations to irradiate a unique bdelloid (right). The figures are not drawn to scale for ease of reading

X-rays

In Hespeels et al., 2020, bdelloids were irradiated with a X-Rad 225-XL (PXi, USA), with spectrum energy ranging from 10 to 225 keV. The desiccated samples received doses ranging from 100 to 7,500 Gy. The modeled geometry used for the X-ray is identical to the one used for proton exposure. The only difference is the source particles, changed from protons to high energy photon particle.

0.5 GeV/n ⁵⁶Fe ions

In Hespeels et al., 2020, a bdelloid layer was placed in a 0.25 mm thick polyethylene plastic PCR testing tube. The iron beam was wider than the sample to irradiate its surface uniformly. The system is modeled using squared polyethylene PCR thickness and bdelloid layer for the sake of even surface subdivision. Heavy ions simulation can be computationally costly. Therefore, the following approximations are made. The beam is set just as wide as the system, to avoid simulating useless ions. The simulated geometry is represented in Fig. 2b.

Sources

Particles sources are implemented in TOPAS as beams, with user-defined geometries. All simulations use a 14 × 14 mm² wide squared beam to cover the bdelloid layers. The sources are placed upstream, at the limit of the world cube. Particles are homogeneously distributed within the beam. The energy of proton and ⁵⁶Fe particle sources is unique: 4 MeV and 0.5 GeV/n, respectively, as used in Hespeels et al., 2020. For X rays, the energy spectrum of a 225 kVp Comet MXR 161 XR tube, as used in Hespeels et al., 2020 is implemented using the SpekCalc software (Poludniowski, 2007; Poludniowski & Evans, 2007; Poludniowski et al., 2009). The output of SpekCalc is implemented in TOPAS using a custom python script. Particles sources of continuous energy are implemented in TOPAS as two vectors: one containing the energy values, and one containing their probability of occurrence (as normalized relative intensity). TOPAS being a MonteCarlo simulation software, it simulates each particle individually. The user sets a pre-defined number of particles to simulate.

Scoring: Dose, and other physical measurements, are retrieved with scorer elements. Three scorers are used: DoseToMedium, SurfaceTrackCount, and PhaseSpace. The DoseToMedium scorer is assigned to a geometry element and computes the deposited dose in the volume. It can segment the geometry element to retrieve the dose in every voxel. This scorer can be used in all simulation levels to find the dose distribution over the bdelloid layer, the dose distribution over the nuclei of a single bdelloid or the dose repartition inside a nucleus. The SurfaceTrackCount scorer is assigned to a geometry element and delivers the number of particles crossing its surface. This scorer can also segment the geometry element in voxels and is used to retrieve the distribution of particles over the voxels. Finally, The PhaseSpace scorer is assigned to a surface. It stores information about the particles crossing its surface in an output file: X, Y, and Z positions, momentum direction, type of particle, charge, energy in MeV, weight, Particle ID, a Boolean flag if the direction of cosine is negative and finally if the scored particle is a primary or secondary particle. A PhaseSpace can be used as source in another simulation. Here, the PhaseSpace is recorded at the first level of granularity (bdelloid layer) to save particles crossing one voxel, i.e., a single bdelloid. This PhaseSpace can then be used as a source in the single bdelloid level (second level of granularity) to study the dose repartition over its nuclei. Figure 2c illustrates this principle, which is less resource consuming than scoring both levels of granularity in one simulation. Simulations on a full layer and on a single bdelloid level were conducted for the three types of particle sources described above.

Poisson distribution as a tool to describe the spatial repartition of particles

Modeling bdelloid layers as cubes is convenient, as cubes are easily subdivided into voxels of identical surfaces. As described before, a 14 × 14 mm² bdelloid layer surface can be subdivided in voxels of 0.1 × 0.1 mm², i.e., the dimensions of a single bdelloid. By associating a DoseToMedium scorer to each voxel, the dose distribution over the bdelloids can be retrieved. For high doses (i.e., high number of particles) and for a uniform beam, the distribution follows a normal law, but for lower doses (less particles), the distribution is more adequately described by the Poisson distribution. The Poisson distribution is a discrete probability distribution, defined by a single parameter µ $\in {\mathbb{R}}$. Its probability mass function formula is the following:

$$\mathrm{P}\left(X=k\right)=\frac{{\mu }^{k}{{e}^{-}}^{\mu } }{k !}$$

With P(X = k) the probability to measure X of valor k and µ the mean value of the distribution. Because X is discrete, it cannot directly represent the dose, which is continuous. The latter, fortunately, is proportional to the total number of particles. Thus, scoring with SurfaceTrackCount on every bdelloid provides the spatial Poisson distribution of the particles. In this case, P(X = k) is the probability for a bdelloid to be traversed by k particles. µ represents the mean number of particles passing through a bdelloid, which can be related to the mean dose over the bdelloid layer.

The Poisson description was used on the second level of granularity as well. By subdividing a 100 × 100 µm² bdelloid surface in 50 × 50 voxels corresponding to 2 × 2 µm² nuclei, the Poisson distribution can describe the spatial repartition of the particles over the nuclei. P(X = k) represents here the probability for a nucleus to be traversed by k particles. µ represents the mean number of particles passing through a nucleus, related to the mean dose over the single bdelloid. The exact proportionality relation is, for particles, the broad beam formula:

$$D=\frac{1.6\times {10}^{-9}LET\phi }{\rho }$$

With D the dose (in Gy), LET the linear energy transfer of the particle, $\phi$ the fluence (in particles/cm²) and $\rho$ the density of the target environment (in g/cm³). As $\mu =\frac{{N}_{\mathrm{tot}}}{{N}_{\mathrm{voxels}}}$, the ratio between the total number of particles and the number of voxels, for a given total number of particles, it is possible to relate µ and the expected dose output through the beam fluence. The key advantage of the Poisson description is that it can be used on every level of simulation: either to describe the repartition of particles over the nuclei of a bdelloid or the repartition of particles over the individuals composing a layer. The only difference will be in the values of N_voxels = 19,600 for bdelloid voxels and 2,500 for nuclei voxels.

Results

Simulation of desiccated bdelloids exposed to proton beam

The first aim was to determine the minimal dose required on a layer of dried bdelloid rotifers to damage all nuclei of all individuals. Since the dose on a layer is the average of the dose over its constituent bdelloids, there will be a bdelloid dose distribution centered on the layer dose. The least damaged bdelloids are expected to be located at the lower end of the distribution and are more likely to remain unaffected at the individual or germinal cells level. Therefore, the initial step in this investigation was to establish the minimal dose required to damage all nuclei of the least damaged bdelloids in the layer.

First, a single bdelloid was modeled and exposed to an increasing number of incident 4 MeV proton particles. Figure 3a illustrates the dose distribution across the nuclei in the bdelloid model. When the number of particles is below 23 × 10³, some nuclei receive a dose of 0 Gy, indicating that a nucleus can remain undamaged. For example, at 2 × 10³ particles (the blue curve in Fig. 3a), 1,125 nuclei (45% of the total) score a null dose. At 13 × 10³ particles (the orange curve in Fig. 3a), 16 nuclei (0.64%) score a null dose. From 23 × 10³ particles onwards (green curve in Fig. 3a), every single nucleus of the bdelloid is damaged and the corresponding dose distribution on the nuclei ranges from 0.30 Gy to 7.14 Gy. On average, the dose delivered to a bdelloid is 2.99 Gy. Therefore, based on our modeling results, a minimal dose of 2.99 Gy is required to damage every nucleus in desiccated bdelloid individuals using 4 MeV protons.

Fig. 3

Analysis of Dose Distribution, Particle Count, and Irradiation Variability in bdelloid model exposed to Proton 4 MeV. a Presents the dose distribution for three different total particle counts (the dose for each voxel is obtained by DoseToMedium scorer). The histograms, represented in blue, orange, and green, count the number of nuclei that have received a given dose. b Presents the particles distribution for the three same total particle count (using SurfaceTrackCount scorer to count particles crossing each nucleus). The histograms, represented in blue, orange, and green, represent the probability for a nucleus to be crossed by a given number of particles. For each number of total particles, the corresponding Poisson distribution is drawn in red (dashed lines). µ is the mean number of particles crossing a voxel. Histogram c is obtained by simulating a layer of bdelloids and attaching a DoseToMedium scorer to every bdelloid voxel. PhaseSpace scorers are attached to three selected bdelloids to investigate the dose disparities at the nuclear level, leading to the histograms in d: a less irradiated bdelloid (in blue at 3.22 Gy), an average irradiated bdelloid (in green at 3.27 Gy) and a highly irradiated bdelloid (in red at 3.32 Gy). e Uses a colormap to represent the dose distribution on an average irradiated bdelloid. Each color dot corresponds to a nucleus. The dose ranges from 0.40 to 7.48 Gy, with no nucleus receiving 0 Gy

The second simulation output is obtained by maintaining the same geometry and employing the SurfaceTrackCount scorer to count particles crossing each nucleus's surface. Figure 3b presents the particle distribution for an increasing number of incident particles. Regarding the dose distribution, at 2 × 10³ particles (blue curve in Fig. 3b), 1,125 voxels (45%) register 0 particles, with a mean value per voxel of 0.80 particles. At 13 × 10³ particles (orange curve in Fig. 3b), 16 voxels (0.64%) register 0 particles with a mean value per voxel of 5.2 particles, and at 23 × 10³ particles (green curve in Fig. 3b), all voxels register between 1 and 22 particles, with a mean value µ = 9.2 particles. The 45%, 0.64%, and 0% proportions of untouched voxels remain consistent for both dose and particle distributions. Since protons continuously interact with matter, a particle passing through a voxel is guaranteed to deposit energy. Overall, this simulation level indicates that 23 × 10³ incoming protons will deposit around 2.99 Gy in a single rotifer, and damage all its nuclei. The Poisson description adequately describes the particle count distribution and indicates that all the nuclei are hit by 9 particles on average at this threshold dose (see Fig. 3b). To simulate at the rotifer layer level, a minimum particle count per rotifer of 23 × 10³ will ensure that no rotifer or nucleus is spared by radiation. To account for the width of the Poisson distribution and push threshold-irradiated individuals to the lower end of the dose distribution, the total dose delivered to the full layer should be superior to 2.99 Gy. For instance, using µ = 23.7 × 10³ (as mean particle count per voxel, corresponding to 3 Gy) at the layer level satisfies to the previous requirement.

In the third part of the simulation, a layer of desiccated bdelloids measuring 14 × 14 mm² is modeled. This simulation aims to demonstrate the variability of exposure across the population of bdelloids in a layer. A mean dose of 3.27 Gy was delivered, slightly exceeding the 3 Gy threshold determined above. The results are presented in Fig. 3c and 3d. Histogram c in Fig. 3 illustrates the dose distribution across the bdelloids in the layer, ranging from 3.17 to 3.77 Gy. As expected, no bdelloid scores 0 Gy at the individual or nucleus level. To further highlight the dose disparities, we assign three PhaseSpace scorers to three different bdelloids, which are highlighted on histogram c: a bdelloid that receives a lower dose (irradiated at 3.22 Gy, shown in blue in histograms c and d), an average irradiated bdelloid (irradiated at 3.27 Gy, shown in green in histograms c and d), and a highly irradiated bdelloid (irradiated at 3.32 Gy, shown in red in histograms c and d). Figure 3d highlights the dose distribution among the nuclei of the three bdelloids. The dose received by the nuclei of the lesser, average, and highly irradiated rotifers ranges 0.40–7.97 Gy, 0.35–6.98 Gy, and 0.41–7.48 Gy, respectively, with mean doses of 3.22, 3.27, and 3.32 Gy. In all cases, no nucleus scored 0 Gy. Figure 3e illustrates in a different way the dose repartition over the average irradiated bdelloid using a color map. The doses range from 0.41 to 7.48 Gy, but no nucleus scored 0 Gy.

Simulation of desiccated bdelloids exposed to ⁵⁶Fe radiation

A similar approach was applied to simulate the exposure of desiccated bdelloid to 0.5 GeV/n ⁵⁶Fe radiation. For the first part of the simulation, a single modeled bdelloid was exposed to an increasing number of ⁵⁶Fe particles. The DoseToMedium scorers attached to the nuclei retrieved the dose distribution over the bdelloid, as presented in Fig. 4a. At 2 × 10³ particles (the blue curve on Fig. 4a), 1,125 nuclei (45% of the total) scored 0 Gy. At 13 × 10³ (the orange curve on Fig. 4a), 16 nuclei (0.64%) scored 0 Gy, and at 23 × 10³ particles (the green curve on Fig. 4a) the dose ranges from 5.90 to 120.65 Gy, with a mean dose of 51.69 Gy. From these modeling results it was concluded that a 51.69 Gy dose is required to damage every nucleus in bdelloid individuals using 0.5 GeV/n ⁵⁶Fe particles.

Fig. 4

Evaluating dose and particle distributions in modeled bdelloids at different irradiation levels of ⁵⁶Fe particles. a Presents the dose distribution for three different total particle counts. The histograms, represented in blue, orange, and green, report the number of nuclei that have received a given dose (the dose for each voxel is obtained by DoseToMedium scorer). b Presents the particles distribution for the three same total particle counts. The histograms, represented in blue, orange, and green, represent the probability for a nucleus to be crossed by a given number of particles (using SurfaceTrackCount scorer to count particles crossing each nucleus). For each number of total particles, the corresponding Poisson distribution is drawn in red (dashed lines). µ is the mean number of particles crossing a voxel. Histogram c is obtained by simulating a layer of bdelloids and attaching a DoseToMedium scorer to every bdelloid voxel. PhaseSpace scorers are attached to three selected bdelloids to investigate the dose disparities at the nuclear level, leading to the histograms in d: a less irradiated bdelloid (in blue at 55.32 Gy), an average irradiated bdelloid (in green at 56.12 Gy), and a highly irradiated bdelloid (in red at 56.98 Gy). e Uses a colormap to represent the dose distribution on an average irradiated bdelloid. Each color dot corresponds to a nucleus. The doses range from 2.50 to 129.88 Gy, with no nucleus receiving 0 Gy

Proton and ⁵⁶Fe particles both interact continuously with matter. Due to this property, the second readout of this simulation level yields identical results to the proton simulation. Particle distributions in this case are purely ballistic and do not depend on the nature of the ion. Similarly to the proton case, it is found that at 23 × 10³ particles all voxels register a mean value µ = 9.2. At the rotifer layer level, 23 × 10³ Fe particles ensures an average dose of 56.12 Gy and push 51.31 Gy threshold-irradiated individuals at the lower end of the dose distribution.

In the third part of the simulation, a layer of desiccated bdelloids was modeled. As for the proton simulation, the results on Fig. 4c and 4d illustrate the diversity of irradiation among the layer of desiccated individuals. The simulation reached the minimum 51.31 Gy threshold for all voxels. Figure 4c represents the dose distribution over the bdelloids of a layer and ranges from 54.78 to 57.60 Gy. No bdelloid scores 0 Gy. Histograms on Fig. 4c are obtained by assigning three PhaseSpace scorers to three different bdelloids: a less irradiated bdelloid (irradiated at 55.32 Gy, in blue), an average irradiated bdelloid (irradiated at 56.12 Gy, in green), and finally a highly irradiated bdelloid (irradiated at 56.98 Gy, in red). The PhaseSpace scorers serve as sources in bdelloid-level simulations to retrieve the dose distribution over the nuclei. The nuclei dose distributions are represented in blue (for less irradiated bdelloid in Fig. 4d), green (for average irradiated bdelloid in Fig. 4d), and red (highly irradiated bdelloid in Fig. 4d). No nucleus scored 0 Gy. The dose distribution over the nuclei is wider than the dose distribution over the bdelloids: 7.91 to 113.02 Gy for the nuclei of the less irradiated bdelloid, 2.50 Gy to 129.88 Gy for the average irradiated bdelloid and 9.35 Gy to 129.19 Gy for the highly irradiated bdelloid. Figure 4e illustrates in a different way the dose repartition on the average irradiated bdelloid using a color map.

Simulation of desiccated bdelloids exposed to 225 kVp X-rays

In contrast to proton and iron ions, X-ray particles do not interact continuously with matter. As a result, the presence or absence of damage in a nucleus is not solely determined by the crossing of a particle but also whether an interaction has occurred. Consequently, it is expected that more than 23 × 10³ particles will be needed to damage all nuclei.

Figure 5a illustrates the dose distribution across the nuclei at the rotifer level. At 0.01 Gy (corresponding to 1 × 10⁸ particles), 8 nuclei (0.36% of the total) still score 0 Gy. 1.55 × 10⁸ particles is the threshold above which all nuclei are damaged. The corresponding dose distribution ranges from 4.68 × 10^–4 to 8.67 × 10^–2 Gy, with an average dose over the bdelloid of 0.02 Gy. This dose of 0.02 Gy is required to damage every nucleus in desiccated bdelloid individuals using 225 kVp X-rays. Similarly, at 0.03 Gy (corresponding to 2 × 10⁸ particles), no nucleus scores 0 Gy, with a dose ranging from 1.62 × 10^–3 to 1.74 × 10^–1 Gy. Figure 5b shows the particle distributions per voxel for an increasing number of incident particles. The three histograms correspond to 1.00 × 10⁸ (in blue), 1.55 × 10⁸ (in orange), and 2.00 × 10⁸ incident particles (in green). The corresponding Poisson distributions are drawn in red, with µ values of 40,000.27, 62,003.57, and 80,004.62, respectively.

Fig. 5

Evaluation of dose and particles distribution of X-ray irradiated bdelloid level for three different total particle counts. The histograms, represented in blue, orange, and green, count the number of nuclei that have received a given dose. b Presents the particles distribution for the same three total particle counts. The histograms, represented in blue, orange, and green, represent the probability for a nucleus to be crossed by a given number of particles. For each number of total particles, the corresponding Poisson distribution is drawn in red (dashed lines). µ is the mean number of particles crossing a voxel. c Uses a colormap to represent the dose distribution on an average irradiated bdelloid. Each color dot corresponds to a nucleus. The dose ranges from 24.79 × 10^–5 to 11.53 × 10^–2 Gy, with no nucleus receiving 0 Gy

To ensure a minimum dose of 0.02 Gy at the lower end of the layer dose distribution, all rotifers should be hit by a minimum of 6.2 × 10⁴ X-rays. Although conducting a full layer irradiation simulation would require 3 to 4 months of computation time, we can confidently estimate the dose received by the bdelloid layer based on the reliable results obtained from simulations with proton and ⁵⁶Fe particles. For instance, with an average layer dose of 2.67 × 10^–2 Gy it is expected that all bdelloids are exposed to a non-zero dose of radiation and that 0.02 Gy threshold-irradiated individuals are pushed to the lower end of the dose distribution. Figure 5c illustrates the dose repartition on the average irradiated (2.67 × 10^–2 Gy) bdelloid using a color map. The dose ranges from 24.79 × 10^–5 to 11.53 × 10^–2 Gy, but no nucleus scored 0 Gy.

Dose distribution at the nucleus level according to particle type

In addition to the above determined threshold doses ensuring irradiation of all nuclei of all rotifers in a layer, each one associated to a given particle type, the distribution patterns of energy deposition at the nucleus level are also valuable, as it can result in distinct biological responses. Figure 6 presents dose deposition map for protons, ⁵⁶Fe and X-ray particles, at comparable doses of 55.10, 55.44, and 54.32 Gy, respectively. The dose distributions exhibit striking differences: good homogeneity for X-rays, sparsity for protons with unaffected DNA regions, and highly localized energy deposits for ⁵⁶Fe ions. In the case of charged particles, the primary particles are responsible for most of the damage, while secondary particles, such as electrons, contribute negligibly to the dose accumulation (as visualized by the light-yellow spots corresponding to secondary electrons for ⁵⁶Fe irradiation in Fig. 6).

Fig. 6

Views of the dose distribution for proton, ⁵⁶Fe and X-ray radiation at nucleus level. Comparable doses of 55.10, 55.44, and 54.32 Gy, respectively, are attained. The squares represent bdelloid nuclei of 2 × 2 µm², segmented in 40 × 40 nm² voxels. Voxels share a common color scale, representing the amount of dose deposited within it, ranging from 0 Gy (white), to 1.45 × 10⁴ Gy (dark brown)

Discussion

In 2020, Hespeels et al., characterized the biological response of desiccated A. vaga individuals to high doses of low- and high-LET radiation including X-rays, 4 MeV protons, and 0.5 GeV/n ⁵⁶Fe particles. Experiments showed that desiccated A. vaga individuals have a lethal dose range of 5,000 to 7,500 Gy when exposed to low- and high-LET radiation. Specifically, the survival rate at a dose of 5,000 Gy was significantly higher for animals exposed to X-ray radiation (82.7% SD ± 11.4) than to proton radiation (19.1% SD ± 34.1) (Hespeels et al., 2020). This finding was consistent with the concept of relative biological effectiveness, which suggests that high-LET radiation is more biologically damaging than low-LET radiation (Goodhead, 1999). Second, results showed a clear difference between low- and high-LET radiation at the reproductive level. A 2.3 factor was reported when comparing SD50 (sterilization dose 50, i.e., the dose of radiation required to sterilize 50% of the population) of X-rays versus protons and ⁵⁶Fe particles.

The progressive decline in fertility of irradiated bdelloids has been previously attributed to the accumulation of cellular damage, particularly DNA damage. When exposed to radiation beyond a certain threshold, the ability of bdelloids to efficiently repair their damaged cells is compromised, leading to sterility. Alternatively, two hypotheses have been proposed to explain the decline with increasing dose observed in previous studies. The first hypothesis suggests that some bdelloids could be partially or completely spared by radiation due to statistical disparities across the bdelloid layer. Secondly, germinal cells of certain bdelloids may remain unaffected by the radiation, preventing sterilization. However, as the radiation dose increases, it may be hypothesized that every nucleus of every bdelloid individual would gradually be impacted, ultimately leading to total sterilization. Simulations in TOPAS are powerful tools to study this phenomenon and, more generally, to provide insights on which dose is needed for each particle source to target all nuclei of all bdelloid individuals, sterilizing the population. To evaluate this, the exposure of desiccated bdelloids to all three radiation sources (proton, ⁵⁶Fe ions, and X-rays) was modeled. The DoseToMedium scorer was used on the entire layer of bdelloids to calculate the global dose, and voxel subdivision was used to determine the dose distribution over individual bdelloids. The PhaseSpace scorer was then used on a bdelloid from the layer as a source to simulate the individual level, which enabled subdivision of the bdelloid into nuclei voxels to determine the dose distribution over all the nuclei. This method enables the verification of whether germline nuclei are spared or not. By iteratively adjusting the number of particles and evaluating the nuclear damage (reported as energy deposition events, leading ultimately to damage on biological molecules, such as DNA) on selected bdelloid individuals, the threshold dose above which all nuclei (both somatic and germline) are exposed can be determined. If these dose thresholds are significantly lower than the SD50 and the doses used in experiments, it would rule out both hypotheses. Thus, simulations were used to determine the minimum macroscopic dose (threshold dose) required on a bdelloid layer to ensure that every nucleus of all bdelloid individuals receives a minimum dose greater than 0 Gy. The same method was used for every particle. Data were summarized in Table 1.

Table 1 Summary of TOPAS simulation data

Comparison of particles and dose thresholds associated with Proton, Fe and X-rays radiation.

Here, TOPAS simulations highlight that exposure to protons 4 MeV, 0.5 GeV/n ⁵⁶Fe, and X-ray radiation at doses of 3.00 Gy, 53.00 Gy, and 2.67 × 10^–2 Gy, respectively, effectively hit all individuals, cells, and nuclei within the samples. All these values were far under the level of dose required to observe any biological response in our previous results supporting the idea that the impact of the radiation on the survival or fertility rate was indeed associated with damages induced by the radiation, and not by the absence of hits in some individuals or germline nuclei. For example, no impact on the reproductive rate was reported at doses below 250 Gy for proton and ⁵⁶Fe exposure and under 750 Gy for X-ray (Hespeels et al., 2020). Doses of at least 1 mGy have been previously shown to induce DNA double-strand breaks (DSBs). Furthermore, a linear relationship between DSB induction and dose has been reported in the past (Rothkamm & Löbrich, 2003). Therefore, based on our simulations, we can extrapolate that it's unlikely that some bdelloids were not hit or some germline nuclei were not affected by radiation in previously published experiments.

Particle radiation generates a unique form of DNA damage through the creation of numerous ionizations localized along the path of the particle through the cell. These damages are not evenly distributed and vary greatly depending on the linear energy transfer (LET) of the radiation. Low-LET radiation results in diffuse damage, while high-LET radiation generates clusters of different types of damage, such as base lesions, abasic sites, single-strand breaks (SSBs), and double-strand breaks (DSBs), which are collectively referred to as complex DNA damage (Ward, 1994; Lomax et al., 2013; Mavragani et al., 2019). Understanding the nature of radiation is essential in predicting the type and extent of DNA damage. Here TOPAS software was used to visualize the distribution of dose deposition within a single bdelloid nucleus. Our results confirm that high-LET particles, such as ⁵⁶Fe and protons, generate more localized dose deposition along the particle track than low-LET X-ray radiation (Fig. 6). Indeed, the pattern generated by TOPAS simulations is consistent with the biological response observed in desiccated bdelloids exposed to low- and high-LET radiation. Previous attempts to experimentally observe differences between DNA damages induced by proton and X-ray radiation on desiccated bdelloids show an increased amount of short-sized DNA fragments after proton exposure, compared to X-ray irradiated samples (Hespeels et al., 2020). This difference was previously attributed to the clustered damages generated by proton exposure, resulting in a higher relative biological effectiveness of proton radiation versus X-rays. However, quantifying short DNA fragments induced by radiation remains a challenging task (Höglund et al., 2000; Alloni et al., 2013; Barbieri et al., 2019). The complexity of damage caused by high-LET radiation is a crucial factor for understanding the effects of radiation exposure, as it can reduce survival rates and the capacity to produce viable offspring in desiccated bdelloids.

In the absence of techniques that enable the examination of specific damages induced by each type of radiation, employing simulation approaches offers valuable understanding of the physical processes at play in the methods applied to desiccated bdelloid specimens. This approach can also be expanded to include other metazoans like tardigrades, which have been previously utilized as a reference for radiotolerance. Incorporating a model into simulations inherently creates some divergence between our experimental design, as simulated by TOPAS, and the actual samples. However, it is impossible to replicate every real-world constraint. For example, our model constrains the nucleus size to a 2 µm square, which may not precisely represent real-world conditions, given that genomic material compaction degree and volume extent can differ between cell types and developmental stages. Estimation based on A. vaga microscopy captures revealed that germline nuclei were approx. 2 µm diameter, as included in our model. However, most of the somatic nuclei have an approx. 3 µm size and germovitellaria nuclei were close to 4 µm (Terwagne et al., 2022). Therefore, despite these limitations, we have confidence in the model's ability to closely approximate the average structural features previously documented in desiccated bdelloids. Finally, it is essential to emphasize that our TOPAS model focused on the description of energy deposition in simulated rotifers, without providing specific data on the resulting damage from dose deposition. Ionizing radiation damage to biological samples can occur through two distinct pathways. Based on "direct damage," we inferred a connection between dose deposition and DNA damage, but additional factors might significantly impact the biological response. The second pathway, called "indirect action," involves generating effectors such as reactive oxygen species (ROS), which are often produced through water radiolysis due to incident radiation, subsequently damaging cellular components (Dartnell, 2011). The proportion of water molecules was found to be 94% in hydrated A. vaga individuals and approximately 6.5% in desiccated samples, indicating significant targets for interaction with ionizing radiation (Hespeels et al., 2014). While it was not feasible to simulate and quantify here the contribution of these elements to the biological response of irradiated bdelloids, the localization of indirect damage is anticipated to follow the pattern of Low-LET, sparse localization” vs. “High-LET, clustered damage." In conclusion, investigating the complexity of damage induced at the cellular level, encompassing proteins and membranes, remains a critical aspect of future research, alongside DNA damage characterization.

By showcasing the differential impact of low- and high-LET radiation on desiccated bdelloid rotifer samples, the current study opens up new avenues for exploration in subsequent research. Key questions to explore include the specific repair mechanisms triggered by different types of radiation-induced damage, as well as how these mechanisms are regulated. Furthermore, it remains unknown how radiation resistance pathways and specific molecular actors themself are impacted by increasing doses of radiation and may be ultimately inactivated. Additionally, it is important to discriminate between DNA repair mechanisms occurring in somatic cells versus germ cells, as this distinction may have significant implications for the understanding of survival and reproduction data. How are the various types of damage repaired? Can we observe any differences in the genomes of the offspring of irradiated organisms? Several ongoing studies employing transcriptomic and comparative genomic approaches should provide more insights into these matters. Among ongoing experiments, these questions are also being tackled in the unique environment of space, on board of the International Space Station (ISS) and will also evaluate the impact of microgravity on the DNA repair system of A. vaga species. Furthermore, examining the repair of such damage across different bdelloid species, both resistant and sensitive to desiccation and irradiation, may help refine our understanding of the mechanisms that enable these animals to withstand massive doses of ionizing radiation. This unique ability sets bdelloid rotifers apart in the animal kingdom.

Conclusion

In this study, the impact of radiation exposure on bdelloid rotifers was simulated, examining whether every nucleus is equally affected or if some rotifers are partially or entirely spared based on parameters used in previously published studies (Hespeels 2014, Hespeels 2020). Additionally, the research aimed to correlate biological data with energy deposition patterns induced by low and high linear energy transfer (LET) radiation, using simulations of protons, iron ions, and X-rays. Our simulations indicated that exposure to 4 MeV protons, 0.5 GeV/n ⁵⁶Fe, and X-ray radiation from doses of 3.00 Gy, 53.00 Gy, and 2.67 × 10^–2 Gy, respectively, would likely uniformly affect all individuals, cells, and nuclei within the samples. Simulated data suggest that the impact of radiation on survival and fertility rates in bdelloid rotifers stem from radiation-induced damage, ruling out the possibility of unaffected rotifers or germinal cells in previous studies, as all administered doses were above the biological response threshold. Furthermore, our simulations unveiled substantial differences between low- and high-LET radiation in terms of irradiated individuals' nuclei. Simulations verified the sparse versus diffuse distribution of radiation hits among cell nuclei when exposed to high- or low-LET radiation, respectively. In summary, our simulations support the notion that the diminished fertility observed in high-LET exposed samples is associated with complex DNA damage due to the concentrated energy deposition pattern of high-LET radiation compared to low-LET.