# Bedrock radioactivity influences the rate and spectrum of mutation

Nathanaëlle Saclier, Patrick Chardon, Florian Malard, Lara Konecny-Dupré, David Eme, Arnaud Bellec, Vincent Breton, Laurent Duret, Tristan Lefébure, Christophe Jean Douady

<figure><img src="image_1.png" /></figure>

## To cite this version:

Nathanaëlle Saclier, Patrick Chardon, Florian Malard, Lara Konecny-Dupré, David Eme, et al. Bedrock radioactivity influences the rate and spectrum of mutation. eLife, 2020, 9, pp.e56830. 10.7554/eLife.56830. hal-03051267

HAL Id: hal-03051267

https://univ-lyon1.hal.science/hal-03051267

Submitted on 10 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Bedrock radioactivity influences the rate and spectrum of mutation

Nathanaëlle Saclier¹†*, Patrick Chardon², Florian Malard¹, Lara Konecny-Dupré¹, David Eme¹‡, Arnaud Bellec¹,³, Vincent Breton², Laurent Duret⁴, Tristan Lefebure¹*, Christophe J Douady¹,⁵

¹Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5023, ENTPE, Laboratoire d'Ecologie des Hydrosysèmes Naturels et Anthropisés, Villeurbanne, France; ²LPC, Université Clermont Auvergne, CNRS/IN2P3 UMR6533, Clermont-Ferrand, France; ³Univ Lyon, Université Jean Moulin Lyon 3, CNRS UMR 5600 Environnement Ville Société, Lyon, France; ⁴Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR 5558, Laboratoire de Biométrie et Biologie Evolutive, Villeurbanne, France; ⁵Institut Universitaire de France, Paris, France

**Abstract** All organisms on Earth are exposed to low doses of natural radioactivity but some habitats are more radioactive than others. Yet, documenting the influence of natural radioactivity on the evolution of biodiversity is challenging. Here, we addressed whether organisms living in naturally more radioactive habitats accumulate more mutations across generations using 14 species of waterlice living in subterranean habitats with contrasted levels of radioactivity. We found that the mitochondrial and nuclear mutation rates across a waterlouse species' genome increased on average by 60% and 30%, respectively, when radioactivity increased by a factor of three. We also found a positive correlation between the level of radioactivity and the probability of G to T (and complementary C to A) mutations, a hallmark of oxidative stress. We conclude that even low doses of natural bedrock radioactivity influence the mutation rate possibly through the accumulation of oxidative damage, in particular in the mitochondrial genome.

## Introduction

Natural radioactivity is the main natural source of exposure to ionizing radiations on Earth. Natural radioactivity is generated by cosmic radiation or by radionuclides released from the bedrock. While levels of cosmic radiation fluctuates over time due to cosmic events such as supernovae or solar flares, bedrock radioactivity remained mainly stable until 2 billion years ago, when it began to slowly decrease (Karam and Leslie, 2005). Bedrock radioactivity depends on the nature of the rocks which extensively varies spatially (e.g. Ielsch et al., 2017). While few extremely naturally radioactive sites such as the India Kerala and Iranian Ramsar region have been monitored for their impact on the human mutation rate (Forster et al., 2002; Masoomi et al., 2006) or on plant physiology (Saghirzadeh et al., 2008), the influence of regional variation in baseline natural radioactivity on the evolution of biodiversity is still unknown (Møller and Mousseau, 2013).

Radioactivity can impact species' molecular evolution by modulating the rate at which different types of mutation appear and accumulate. Ionizing radiations damage DNA by breaking the DNA sugar-phosphate backbone (Hoeijmakers, 2001; van Gent et al., 2001). Alternatively, ionizing radiations can trigger the formation of reactive oxygen species (ROS) directly inside cells through the radiolysis of water (Wallace, 1998; Ward, 1988). ROS being also mutagenic (Barja, 2002), ionizing radiations are direct and indirect mutagens. The most deleterious and studied mutations generated by radioactivity are DNA double-strand breaks (DSBs). Two repair systems are able to manage double-strand breaks: homologous recombination (HR) and non-homologous end-joining (NHEJ). HR

*For correspondence:
nathanaelle.saclier@gmail.com
(NS);
tristan.lefebure@univ-lyon1.fr (TL)

**Present address:** †ISEM, CNRS, Univ. Montpellier, IRD, EPHE, Montpellier, France; ‡IFREMER, unité Ecologie et Modeles pour l'Halieutique, Nantes, France

**Competing interests:** The authors declare that no competing interests exist.

**Funding:** See page 13

**Received:** 11 March 2020

**Accepted:** 30 November 2020

**Published:** 30 November 2020

**Reviewing editor:** Molly Przeworski, Columbia University, United States

© Copyright Saclier et al. This article is distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use and redistribution provided that the original author and source are credited.
uses an homologous sequence to revert the damage, but is only used at specific cell cycles when an
undamaged sister chromatid is available (Karran, 2000). Otherwise, NHEJ is preferentially used. As
NHEJ does not use a template-strand to reconstruct the missing genetic information, it only restores
the continuity of the DNA molecule to the price of frequent deletions and chromosomal exchanges.
Rothkamm et al., 2001 showed that 50% of the double-strand breaks were misrejoined after a
strong irradiation (80 Gy). Numerous studies demonstrate that an exposure to radiation produces
chromosomal abnormalities (Dikomey et al., 2000; Loucas et al., 2004; Hande et al., 2005), dele-
tions (Jostes et al., 1994; Huo et al., 2001; Adewoye et al., 2015; Allegrucci et al., 2015), and
point mutations (Huo et al., 2001; Winegar et al., 1994; Forster et al., 2002; Barber et al., 2002).
Radiation-induced chromosomal abnormalities and deletions have been thoroughly studied because
of their frequent deleterious impact. The impact of radiation on point mutations has received less
attention.

While the mutational impact of exposure to high doses of radioactivity is well characterized (Dubrova et al., 1996; Ziegler et al., 1993), the exposure to low doses of radioactivity is poorly known. Some authors (Tubiana et al., 2006) propose that DNA repair and apoptosis may completely counteract the effect of ionizing radiations for doses below 0.1 Gy, suggesting that low doses have no biological impact. Indeed, in vitro experiments on mammalian cells show that repair systems are more efficient to repair DSBs at low dose than at higher dose (Boucher et al., 2004, 0.05 vs 3.5 Gy/min). However, exposure to an even lower dose of ionizing radiation (less than 0.1 cGy/min) increases the number of mutants in mammals (Vilenchik and Knudson, 2000; Hooker et al., 2004) suggesting that repair systems are not activated at very low radiation doses. When repair systems are not activated and DNA damage accumulates, cells tend to die by apoptosis, preventing the transmission of radioactivity-induced mutations to the next cell generation (Rothkamm et al., 2001; Collis et al., 2004). Unrepaired DSBs are the main cause of radioactivity-induced cell apoptosis during mitosis (Krueger et al., 2007). As point mutations will not lead to cell death like DSBs, this type of radioactivity-induced mutations could stay unrepaired and be transmitted to the next generation. Exposure to low doses of radioactivity also induces an adaptive response: exposed cells are resistant to a following higher dose of radiation (see Rigaud and Moustacchi, 1996). An early or chronic exposure to low doses of radiation may strengthen the antioxidant defense and reduce the sensitivity to radioactivity-induced ROS. While this effect is well demonstrated for preventing DSBs, its efficiency to reduce point mutations is less certain (Rigaud et al., 1995). Moreover, the adaptive response seems to only protect the nuclear genome while the mitochondrial one may not benefit from it (Jarrett and Boulton, 2005).

The short-term impact of high doses of radiation is indisputable, but the long-term impact of
exposure to low doses is unknown. Some studies found an increase in the number of mutations in
the offspring of exposed people (Dubrova et al., 1996; Forster et al., 2002), while others studies
found the opposite (Satoh et al., 1996; Czeizel et al., 1991), leaving open the question of the trans-
mission of mutations generated by low doses of radioactivity. This lack of knowledge is likely contrib-
uted to by three factors. First, studies are mainly focused on the health effect of low dose of
radiation and not on their long-term mutational impact (e.g. Beir, 2006; Tubiana et al., 2006). Sec-
ond, most of the literature focuses on unrepaired double-strand breaks which are highly deleterious
and are de facto not transmitted to the next generation. Third, studying the long-term mutational
impact of natural radioactivity is challenging because it raises a number of methodological difficul-
ties. On the one hand, experimentally exposing multiple generations of multicellular organisms to
low doses of radiation would require years of experimentation and complex experimental controls.
On the other hand, the main obstacles to in naturae studies are the organisms' mobility, which pre-
vents certainty that a population was exposed to the same natural radioactivity for many genera-
tions, and confounding factors such as ultraviolet radiation from the sun.

Here, we overcome these difficulties by coupling in situ radioactivity characterizations with the
distinctive bio-ecological characteristics of subterranean waterlice within a phylogenetic comparative
framework. Subterranean waterlice are never exposed to UV radiation, live in contrasted bedrock
set-ups and have very limited dispersal capacity (Eme et al., 2018), allowing us to make the assump-
tion that different species have persisted in different but nearly constant radioactive habitats for
numerous generations. To test the impact of radioactivity on the transmission of point mutations, we
estimated the long-term rate of neutral mutations in the nuclear and in the mitochondrial genome
independently.
## Results and discussion

In order to build a robust and powerful comparative design aimed for testing the influence of natural radioactivity on the mutation rate, we first prospected for closely related subterranean species living in contrasted radioactive set-ups. Using the map of bedrock uranium content in France (Ielsch et al., 2017), we prospected areas with low and high radioactivity. From this large survey (58 sites with waterlice), we selected 14 sites with contrasted levels of α radioactivity which were inhabited by closely related groundwater waterlice species. We paid special attention that α radioactivity differed at least by a factor of three between two habitats, each containing a closely related species (Figure 1, Supplementary file 1). On average low level of radioactivity was around 0.357 Bq/g of dry sediment and high level around 1.259 Bq/g of dry sediment. Based on transcriptome sequencing and de novo assembly, we used a phylogenetic approach to estimate nuclear and mitochondrial substitution rates (i.e. the rate of mutations which are fixed). While experimental approaches allow to measure the impact of radioactivity on somatic mutations or on the transmission of mutations across few generations, this phylogenetic approach allows us to measure the impact of natural radiation on the germinal mutation rate over the course of a species history.

Using the 14 selected species and locations, we tested whether there was a significant positive relationship between natural radioactivity and the long-term mutation rate. The latter was estimated using the synonymous substitution rate ($d_S$) calculated on the terminal branches of the phylogenetic tree tracing the history of these 14 species. The $d_S$ is the rate at which silent mutations accumulate in protein coding genes and, when calculated using many different loci and in the absence of a strong synonymous codon usage bias (see methods), is an estimator of the average mutation rate across a

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Species and locations selected to study the impact of bedrock radioactivity on the mutation rate and spectrum. Fourteen species with contrasted bedrock radioactivity exposure were selected (black dot: low exposure, red dot: high exposure). Based on their phylogenetic history (a), we further selected six monophyletic pairs of closely related species to compare their mutational spectrum (Vielvic and Montbar are excluded because of unresolved phylogeny, pairs are indicated using superscript numbers). Received dose of radioactivity was determined from measurements of radioactivity in the sediments of the sampled sites (b). For each site, the areal proportion of low-radioactivity sedimentary rocks and high-radioactivity metamorphic and igneous rocks in a radius of 15 km around the sampling (λ15) site is represented with circles next to the map (c). The online version of this article includes the following figure supplement(s) for figure 1:</figcaption></figure>

**Figure supplement 1.** Phylogeny of the 14 species used to compute synonymous substitution rates.
species' genome (Kimura, 1983). To be comparable across species, these $d_S$ has to be divided by the time over which it is measured. To achieve that, we used the software CoEvol (Lartillot and Poujol, 2011) which models directly a synonymous substitution rate relative to the root age ($d_S$/ra). We computed $d_S$/ra using 769 one-to-one nuclear and 13 mitochondrial orthologous protein-coding genes shared by all 14 species. At each sampling site, we measured the global $\alpha$ radioactivity and the activity of all radio-elements in the sediment. The analysis of the composition in radionuclides at each site reveals that two sites (BRETEMIN and BOREON) show a disruption of the secular equilibrium in the U-238 chain. This suggests that nearby industrial activities (e.g. lead mines) have modified the natural radioactivity at these two sites. As these industrial activities are very recent (since 1950), their impact on the substitution rate, which is measured on a much longer time scale, is unlikely. We therefore did not use these two sites to test the correlation between $d_S$ and any site-specific radioactivity measurement (however, see next paragraph for a regional measurement). The $d_S$/ra is positively correlated with the $\alpha$ radioactivity in the nuclear genome as well as in the mitochondrial genome (Table 1; Figure 2). A linear model predicts a $d_S$/ra increase of 31.8% in the nuclear genome and 56.5% in the mitochondrial genome between species living in low (on average 0.357 Bq/g of dry sediment) and high radioactivity (on average 1.259 Bq/g of dry sediment). We also modeled the biologically effective dose of radioactivity received by each species (Received dose in $\mu$Gy/h, Figure 1). This measure takes into account the transfer coefficient from environment to biota and the radio-toxicity of each radio-element for a crustacean model (ERICA tool V1.2.1 Brown et al., 2016). Again, we found positive correlations between the $d_S$/ra and the received dose of radioactivity (Table 1). These results are robust to the presence of influential cases and to variation in the model of evolution used to perform the test (see Methods, pGLS, p. values < 0.05).

As previously explained, the measured radioactivity at two sites overestimates the radioactivity level to which the organisms have been exposed for many generations because it is influenced by recent human activities. Moreover, while most species collected in highly radioactive habitats were from metamorphic or igneous formations, two species were from sedimentary formations (BOREON and GROTTAZE). Contrary to metamorphic and igneous formations, radioactivity in sedimentary formations is often observed in restricted localities (Ielsch et al., 2017) and can show large variations at the meter scale. A single radiation measurement may not therefore accurately represent the average radiation that a species was exposed to. To account for this variability as well as to include the two human-impacted sites into the regression analysis, we calculated the areal proportion of metamorphic and igneous rock within a 15 km radius around each site (later called $\lambda$15, Figure 1). This proportion was used as a proxy for the long-term regional radioactive exposure because the average linear distribution range of a groundwater crustacean is 30 km (Eme et al., 2018). We found a positive and stronger correlation between the $d_S$/ra and $\lambda$15 in both genomes (Table 1, Figure 2, n = 14 species). The linear model predicts that the nuclear and mitochondrial $d_S$/ra of a species living in a metamorphic formation (>50% of metamorphic and igneous rocks) are on average 34.4% and 61.3% higher, respectively, than those of a species living in a sedimentary formation (<50% of metamorphic or igneous rocks).

<caption>**Table 1.** Phylogenetic generalized least square (pGLS) regressions of the nuclear synonymous substitution rate ($d_S$/ra) and mitochondrial $d_S$/ra against the $\alpha$ radioactivity measured in the sediments ($\alpha$ radio.), the received dose (RD) of radioactivity, and the surface of metamorphic and igneous bedrock within a 15 km radius around the sampling sites ($\lambda$15).</caption>

$\alpha$ radioactivity and RD were log transformed to fit with the linear model assumptions. $R^2$ are Cox-Snell pseudo $R^2$.

<table><thead><tr><th rowspan="2"></th><th colspan="4">Nuclear d<sub>S</sub>/ra</th><th colspan="4">Mitochondrial d<sub>S</sub>/ra</th></tr><tr><th>Slope</th><th>L. ratio</th><th>p. value</th><th>R<sup>2</sup></th><th>Slope</th><th>L. ratio</th><th>p. value</th><th>R<sup>2</sup></th><th>N taxa</th></tr></thead><tbody><tr><td>log(<i>α</i> radio.)</td><td>0.034</td><td>5.995</td><td><strong>0.014</strong></td><td>0.393</td><td>0.506</td><td>7.895</td><td><strong>0.005</strong></td><td>0.482</td><td>12</td></tr><tr><td>log(RD)</td><td>0.038</td><td>6.51</td><td><strong>0.011</strong></td><td>0.419</td><td>0.491</td><td>5.981</td><td><strong>0.015</strong></td><td>0.392</td><td>12</td></tr><tr><td>$\lambda$15</td><td>0.076</td><td>9.039</td><td><strong>0.003</strong></td><td>0.476</td><td>1.097</td><td>11.680</td><td><strong>0.001</strong></td><td>0.566</td><td>14</td></tr></tbody></table>

Each line corresponds to one likelihood ratio test between the models with and without the given explanatory variable. The impact of multiple predictors that are collinear are unreliable in the linear model framework (Quinn and Keough, 2002). As the three independent variables are collinear (R2 > 0.6) a model with a combination of these variable is not shown.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Relationships between synonymous substitution rate relative to the root age ($d_S$/ra) and radioactivity measured either as the sediment $\alpha$ radioactivity (left), the received dose (middle) or the proportion of igneous and metamorphic rock in a radius of 15 km around the sampling sites, $\lambda_{15}$ (right). Each dot represents a species nuclear (top) and mitochondrial (bottom) $d_S$/ra. The fit of the pGLS model is indicated with a solid line and the confidence interval of the correlation is indicated with dashed lines. The two species of a pair are labeled with a number as in Figure 1 (the number seven corresponds to the two species not forming a pair), with in red the species of the pair living in the highest $\alpha$ radioactivity. Species in sedimentary formations ($\lambda_{15} < 50\%$) are depicted with circles and species in igneous/metamorphic formations ($\lambda_{15} > 50\%$) with a diamond.</figcaption></figure>

As bedrock radioactivity is positively correlated with the mutation rate, the underlying question is whether radioactivity also modifies the mutational spectrum, that is, the specific types of mutations that tend to occur. To address this question, we reconstructed the mutational spectrum of 6 independent pairs of species, each composed of two species located in low and high bedrock radiation set-ups, respectively, with a minimum of a 3X increase in the received dose of radioactivity between the two species (Figure 1). Briefly, we first estimated species polymorphism across a set of 2490 one-to-one orthologous genes by sequencing transcriptomes for eight individuals per species. After ancestral sequence reconstruction, we then identified mutations that occurred in each species and computed the relative proportion of each type of mutation (from A to T, A to C, ...), pooling together complementary mutations (e.g. $p(C \to A) + p(G \to T) = p(C:G \to A:T)$). Two types of dependencies are present in testing mutational spectrum variation in a comparative data-set: (i) dependency among mutations – if the proportion of one mutation increases, the proportions of the other mutations will decrease – and (ii) the phylogenetic inertia, two closely related species have more chance to display more similar mutational spectrum. To take into account these two types of dependencies, we first used a forward selection approach as described in Harris and Pritchard, 2017 to pull out mutations at different frequencies across habitats (see Materials and methods). Only the $C:G \to A:T$ mutation was significantly more frequent in radioactive habitats (Table 2, Figure 3). Second, we checked that these results were not induced by the phylogenetic structure of the data-set by using a pGLS regression between the proportion of each type of mutation and different radioactivity proxies. Again, we found a strong positive correlation between the proportion of $C:G \to A:T$
<caption>**Table 2.** Variation of the mutational spectrum as a function of bedrock radioactivity.</caption>

An ordered $\chi^2$ test is first used to test if mutation counts vary between low and high radioactive habitats while accounting for the dependency among mutation. In parallel, the phylogenetic dependency was taken into account using a Phylogenetic Generalized Least Square (pGLS) regression of the proportion of each mutation against the sediment $\alpha$ radioactivity ($\alpha$ radio.), Received Dose (RD), and areal proportion of metamorphic and igneous rock within a 15 km radius ($\lambda$15). $\alpha$ radioactivity and RD are log transformed to fit with the linear model assumptions. $R^2$ are Cox-Snell pseudo $R^2$.

<table><thead><tr><th rowspan="2">Mutation type</th><th>Ordered χ² test</th><th colspan="6">pGLS regression</th></tr><tr><th>Ordered p. value</th><th>Radio.</th><th>Slope</th><th>L.ratio</th><th>P.value</th><th>R²</th><th>N</th></tr></thead><tbody><tr><td rowspan="3">P(C:G→A:T)</td><td rowspan="3">0.000</td><td>log(α radio)</td><td>0.013</td><td>13.010</td><td><strong>0.000</strong></td><td>0.662</td><td>10</td></tr><tr><td>log(RD)</td><td>0.014</td><td>12.079</td><td><strong>0.006</strong></td><td>0.635</td><td>10</td></tr><tr><td>λ15</td><td>0.042</td><td>13.791</td><td><strong>0.000</strong></td><td>0.683</td><td>12</td></tr><tr><td rowspan="3">P(A:T→T:A)</td><td rowspan="3">0.111</td><td>log(α radio)</td><td>-0.009</td><td>6.819</td><td>0.072</td><td>0.433</td><td>10</td></tr><tr><td>log(RD)</td><td>-0.011</td><td>9.348</td><td><strong>0.030</strong></td><td>0.541</td><td>10</td></tr><tr><td>λ15</td><td>-0.025</td><td>6.994</td><td>0.080</td><td>0.442</td><td>12</td></tr><tr><td rowspan="3">P(C:G→T:A)</td><td rowspan="3">0.351</td><td>log(α radio)</td><td>0.012</td><td>3.323</td><td>0.612</td><td>0.242</td><td>10</td></tr><tr><td>log(RD)</td><td>0.010</td><td>1.684</td><td>0.776</td><td>0.131</td><td>10</td></tr><tr><td>λ15</td><td>-0.009</td><td>0.183</td><td>1.000</td><td>0.015</td><td>12</td></tr><tr><td rowspan="3">P(A:T→C:G)</td><td rowspan="3">0.520</td><td>log(α radio)</td><td>0.001</td><td>0.121</td><td>1.000</td><td>0.010</td><td>10</td></tr><tr><td>log(RD)</td><td>0.000</td><td>0.015</td><td>1.000</td><td>0.001</td><td>10</td></tr><tr><td>λ15</td><td>0.001</td><td>0.005</td><td>1.000</td><td>0.000</td><td>12</td></tr><tr><td rowspan="3">P(C:G→G:C)</td><td rowspan="3">0.982</td><td>log(α radio)</td><td>0.005</td><td>2.137</td><td>0.864</td><td>0.163</td><td>10</td></tr><tr><td>log(RD)</td><td>0.006</td><td>2.403</td><td>0.726</td><td>0.181</td><td>10</td></tr><tr><td>λ15</td><td>0.006</td><td>0.360</td><td>1.000</td><td>0.030</td><td>12</td></tr><tr><td rowspan="3">P(A:T→G:C)</td><td rowspan="3">0.982</td><td>log(α radio)</td><td>-0.022</td><td>7.778</td><td>0.075</td><td>0.477</td><td>10</td></tr><tr><td>log(RD)</td><td>-0.019</td><td>4.004</td><td>0.180</td><td>0.284</td><td>10</td></tr><tr><td>λ15</td><td>-0.014</td><td>0.394</td><td>1.000</td><td>0.032</td><td>12</td></tr></tbody></table>

Each line of the pGLS tests corresponds to one likelihood ratio test between the models with and without the given explanatory variable. P values have been corrected following Holm's method (k=18).

mutations and bedrock radioactivity (Table 2, Figure 3). Selection is unlikely to be responsible for the increase of C:G → A:T as this correlation is also observed when the data set is limited to mutations found at the third, usually redundant, codon position (Figure 3—figure supplement 1, Supplementary file 2). A weak negative trend is observed for A:T → T:A mutations; however, this trend disappears when measured on third codon position.

Variations in the frequency of genetic variants among populations can originate from variation in the mutagenic environment or from biases in the fixation of mutations that can be due to demographic factors (Mathieson and Reich, 2017), natural selection (Boussau et al., 2008), or biased gene conversion (Duret and Galtier, 2009). While the impact of fixation biases on sequence evolution have been relatively well described, the long term impact of the mutagenic environment is less well-known. The most extensive demonstration of environmentally induced changes in the mutation spectrum comes from the study of carcinogens (Seo et al., 2000). Many carcinogens induce a specific mutational signature, for instance, UVs increase the CC → TT mutations whereas estrogen treatments increase the frequency of A:T → G:C mutations (Flibotte et al., 2010; Nik-Zainal et al., 2015). Here, we found an increase of the proportion of G → T mutations (and complementary C → A) which is a mutation that increases in a variety of contexts and is not specific to a mutagen in particular. The occurrence of this mutation increases after exposure to carbon black (Jacobsen et al., 2011), to tobacco smoke (Hollstein et al., 1991; Hainaut and Pfeifer, 2001), or to polycyclic aromatic hydrocarbons (Nik-Zainal et al., 2015). Oxidative stress is the likely source of this mutation (Wood et al., 1990; Besaratinia et al., 2004; Jacobsen et al., 2011). Indeed, the characteristic damage linked to oxidative stress is the formation of 8-hydroxyguanine (Shigenaga et al.,
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Contrasts ($\pi$) of the relative proportion of each mutation [p(i:j→k:l)] in each pair of sister species: $\pi_m = \log \frac{p^{(j-k:l)+}}{p^{(j-k:l)-}}$ where + and - refer respectively to the species exposed to the higher and lower level of radioactivity in the pair m. Thus, positive bars represent a higher proportion of the given mutation in the species living in the high radioactivity rock. From left to right, bars are in increasing order of difference ($\Delta$) in $\lambda_{15}$ (the areal proportion of igneous and metamorphic rock in a radius of 15 km around the site) between the two species of each pair. From left to right, mutations are in increasing order of correlation with radioactivity. Numbers below the color scale indicate the species pair number as in **Figure 1**.
The online version of this article includes the following figure supplement(s) for figure 3:</figcaption></figure>

**Figure supplement 1.** Contrasts ($\pi$) of the relative proportion of each mutation [$P(i:j \to k:l)$] computed only on third positions in each pair of sister species: [$\pi_m = \log \frac{P^{(i:j \to k:l)+}}{P^{(i:j \to k:l)-}}$] where + and - refer respectively to the species exposed to the higher and lower level of radioactivity in the pair m.

1989) which leads to $G \to T$ mutations (Shibutani et al., 1991; Cheng et al., 1992). Thus, the mutational spectrum modification observed in radioactive environments is suggestive of an increase of the overall oxidative stress in these environments. Similarly, high artificial doses of ionizing radiation were found to increase oxidative damage (Einor et al., 2016; Haghdoost et al., 2006). As oxidative stress triggers more 8-hydroxyguanine formation in the mitochondrial genome than in the nuclear genome (Richter et al., 1988), this could explain the twice as strong effect observed in this study in the mitochondrial compartment.

Radioactive environments can cause mutations in two intertwined ways. First, the ionisation of molecules in the cells can directly affect the DNA structure by breaking the sugar phosphate backbone or can affect DNA indirectly through the radiolysis of water which decomposes the H₂O molecules and create free radicals (Desouky et al., 2015). These free radicals can damage DNA molecules and create mutations. Second, radioactive decay chains also generate heavy metals (lead, polonium, etc) which are toxic for cells and also cause oxidative stress (Quinlan et al., 1988; Pinto et al., 2003). Due to the physicochemical association between radioactivity and heavy metals, in naturae correlative approaches alone cannot discriminate between the toxicity of heavy metals and the direct toxicity of radioactive rays.

While metamorphic and igneous rocks are more naturally radioactive than sedimentary rocks (Ielsch et al., 2017), they also display characteristic compositional and structural features which may also control species life history traits (Cornu et al., 2013; Eme et al., 2015). Here, we consider whether these characteristic features – specifically their poor calcium content and reduced habitat size and permeability – may confound the effect of radioactivity on the mutation rate. Metamorphic and igneous rocks are calcium poor relative to sedimentary rocks. In isopods, the cuticle is made of
calcium carbonate which is acquired by isopods from their environments (Greenaway, 1985). A lower calcium availability may slow down the growth of isopod species. However, a lower growth rate, if it slows down the generation time, is expected to decrease the mutation rate, while we observed the opposite. Because metamorphic and igneous rocks have smaller size and less permeable habitats, they could support smaller effective isopod population sizes. The $d_S$ is an estimator of the mutation rate that is not directly impacted by the effective population size (Kimura, 1983). However, as effective population size modulates the efficacy of natural selection, Lynch, 2010 proposed that species with small effective population size should evolve higher mutation rate as a result of the accumulation of weakly deleterious mutations in genes involved in repair efficiency and replication fidelity. Efficacy of selection can be estimated by computing a transcriptome-wide ratio $d_N/d_S$. Here, we found no correlation between radioactivity and $d_N/d_S$ (pGLS, p.value = 0.3069) suggesting this effect is not at work in this data-set. Altogether, the influence of confounding factors that would drive the observed mutational variations is unlikely.

If oxidative stress is causing the increase of mutation rate, the radiolysis of water alone hardly explains the much higher impact of radioactivity on the mitochondrial mutation rate compared to the nuclear rate. As there is no reason to argue that radiolysis would not evenly occur within cells, it should impact both genomes similarly. However, differences between the two genomes may explain why the mitochondrial genome is more sensitive to radioactivity. First, the mitochondrial genome lacks some repair systems. For example, in a different yet analogous context, UV damage accumulates in the mitochondrial genome while they are repaired in the nuclear genome (Clayton et al., 1974). Second, the two genomes have very different organizations: while the nuclear genome is compacted into chromatin, the mitochondrial one is organized into nucleoids (Chen and Butow, 2005). While the role of the mitochondrial nucleoids is unclear, the chromatin structure and histone proteins protect the nuclear genome against radiation-induced strand breaks (Ljungman, 1991) and oxidative damage (Ljungman and Hanawalt, 1992). Third, direct radioactivity damage to the mitochondria may increase the activity of the mitochondrial respiratory chain and indirectly increase the production of ROS and DNA damage (Yamamori et al., 2012; Kam and Banati, 2013).

While life history traits are known to be central in controlling the mutation rate in most taxa (Nabholz et al., 2008; Bromham et al., 1996; Nikolaev et al., 2007; Smith and Donoghue, 2008), our results suggest that natural variation of radioactivity can have a comparable effect. Indeed, we found a minimum increase of around 30% percent of the nuclear mutation rate (60% in mitochondria) for species of waterlice living in the more naturally radioactive habitats made of igneous and metamorphic rocks. This increase is of the same magnitude as that observed when waterlice species evolve a 5-fold increase in generation time, a key life history trait controlling mutation rate in waterlice (Saclier et al., 2018) and organisms in general (Thomas et al., 2010; Weller and Wu, 2015; Nabholz et al., 2008). However, the influence of life history traits on the mutation rate varies widely among groups (Allio et al., 2017). Thus, the relative influence of radioactivity compared to life history traits could be different in groups like vertebrates, potentially because of different gametogenesis (Saclier et al., 2018). Moreover, as groundwater waterlice ingest sediments (Francois et al., 2016), they are internally exposed to radioactivity, which may cause more mutations than through external exposure only (Sawada, 2007). The influence of environmental radioactivity on the mutation rate should therefore be explored across a wider range of organisms with contrasted diets and biologies.

Although the literature on the effect of low doses of radiation is far from being conclusive it often suggests a negligible biological impact (Tubiana et al., 2006; Tubiana et al., 2009) and only a handful of isolated studies support an impact of natural radioactivity on the mutation rate. For instance, a higher mutation rate was observed in the human mitochondrial genome in the Kerala region (Forster et al., 2002) and in satellite sequences of crickets inhabiting cave with high radon concentration (Allegrucci et al., 2015). In this study, by combining a large number of genes with the characteristics of the subterranean waterlice, namely the absence of UV confounding effect and limited dispersal, within a statistically powerful comparative framework allowing to work on large time scales and with numerous replicates, we found that a mild variation (≈ 3.5X) in natural bedrock radioactivity substantially alters the mutation rate, in particular the mitochondrial one. One key aspect that remains to be described is the shape of the relationship between the dose of radioactivity and the mutation rate: this study invalidates a model where low doses have no impact but falls short in differentiating between, for example, a linear and a hypersentivity (U-shaped) dependency model.
Altogether, while the universality of this finding warrants corroborative studies in other taxa, it suggests that the influence of natural radioactivity on the evolution of biodiversity has been overlooked.

# Materials and methods

## Sampling

To test the impact of radioactivity on molecular evolution, we focused on subterranean species belonging to the Asellidae family. Subterranean species are never exposed to UV radiation, live in contrasted bedrock set-ups and have very limited dispersal capacity (Eme et al., 2018), allowing us to make the assumption that different species have persisted in different but nearly constant radio-active habitats for numerous generations (but see Statistical analyses paragraph). One of the most interesting feature of subterranean Asellidae is their similarities in terms of morphology, lifestyle and life history traits. This high uniformity is likely due to a low rate of phenotypic evolution and a high level of convergence imposed by the subterranean lifestyle. As a result, distinguishing different species requires a high level of expertise and some species cannot be distinguished without molecular tools (Morvan et al., 2013). The birth of the Asellidae family is estimated at -350 My (Morvan et al., 2013). These characteristics allow us to compare species that are divergent enough to compute an accurate rate of molecular evolution but which conserved with very similar traits.

For 58 sites in France selected on the map of uranium (Ielsch et al., 2017), we collected Asellidae species and sampled about 50 g of sediment to measure global α radioactivity (see the following paragraph). Animals and sediments were collected using the Bou-Rouch pumping methods (Bou and Rouch, 1967). Collected specimens were stored in 96% ethanol at -20°C and were morphologically and molecularly identified. For molecular identification, DNA was extracted using an optimized chloroform DNA extraction protocol for the Aselloidea (Calvignac et al., 2011). We amplified DNA with primers targeting the 16S mitochondrial rDNA gene. PCR reactions were done following Morvan et al., 2013. PCR products were sequenced in both directions using the same primers as for amplification (GATC Biotech, Konstanz; Eurofins MWG Operon, Ebersberg; SeqLab, Göttingen, Germany; BIOFIDAL, Vaulx-en-Velin, France). Chromatograms were visualized and cleaned using Finch v1.5.0 (Geospiza, Seattle, USA). 16S have been deposited on the European Nucleotide Archive and are available under the accession number from LR214526 to LR214880. Using Eme et al., 2018 molecular species delimitation, each sequence has been assigned to a species. Based on this taxonomic assignment and radioactivity measurement, 14 species were retained for further analyses (Supplementary file 1). For these 14 selected species, during a new sampling trip, individuals were flash frozen alive in the field.

## Measures of radioactivity

### α radioactivity

In order to estimate the global radioactivity in sediments, we measured the α radioactivity. An α decay occurs when an atom disintegrates by ejecting an α particle, that is, a particle made of two neutrons and two protons. The α radioactivity should be correlated with the global radioactivity in natural systems. For the 58 prospected sites, three samples of about 50 g of sediment were collected in polyethylene bottles. α radioactivity measurements were made by the LABRADOR service (Institut de Physique Nucléaire de Lyon, France) on proportional counter with the NF ISO 18589-6 standard (Data available on Zenodo, DOI: 10.5281/zenodo.4071754).

### Received dose

In order to estimate the received dose of radiation that is impacting organisms, we collected three samples of 100 g of fine sediments (<100 µm) in each of the selected sites. These sediments were prepared with the NF EN ISO 18589-2 standard and measured by gamma spectrometry in conformity with the NF EN ISO 18589-3 standard using the PRISNA-P analysis platform at the Centre d'Etude Nucléaire de Bordeaux Gradignan (CENBG). This platform is certified by the French Nuclear Security Authority (ASN) for measures of natural radioactivity. Samples were dried in open air, and then dried at 100°C. Matters were packed in a waterproof geometry. Geometries were sealed for one month and then counted for a duration of 86500 s on the same chain of measure. The chain used is an ORTEC chain, presenting an efficiency of about 60% and calibrated in May 2016. This
chain is equipped with a cosmic veto device and located in a half buried laboratory in order to: (i) attenuate the background noise, (ii) improve the detection limits, and (iii) reduce the measure uncertainty. The activity of the main radionuclides were measured in sediment and the activity of the remaining radionuclides was deduced based on the hypothesis of a secular equilibrium of the uranium 238 and thorium 232 chains. As activities of the radionuclides of the uranium 235 decay chain are generally low, only measures higher than the decision threshold (according to the measure variability) were taken into account. When the uranium 235 activity was too low to be measured it has been deduced from the uranium 238 activity, using the natural isotopic ratio of 21.6.

The received dose impacting organisms was estimated using the ERICA tool (V1.2.1, *Brown et al., 2016*) with a 'crustacean' model. We assumed that organisms stay 10% of their time on the surface of sediment and 90% inside sediment. All radionuclides available in the tool were taken into account (i.e, *U*²³⁸, *Th*²³⁴, *U*²³⁴, *Th*²³⁰, *Ra*²²⁶, *Pb*²¹⁰, *Po*²¹⁰, *U*²³⁵, *Th*²³¹, *Pa*²³¹, *Th*²²⁷, *Th*²³², *Ac*²²⁸ and *Th*²²⁸). We used the distribution coefficients proposed by the ERICA tool. Concentration factors proposed by the tool were used when available. If not, we used the concentration factor of the closest biogeochemical element available.

Two sites (BRETEMIN and BOREON) show a disruption of the secular equilibrium in the *U*²³⁸ chain. This suggests that nearby industrial activities (e.g. lead mines) have modified the natural radioactivity of these two sites. As these industrial activities are very recent (since 1950), their impact on the substitution rate which is measured on a much longer time scale is unlikely. These two sites were removed from the correlation between *d*<sub>S</sub>/ra and radioactivity measured with the global α radioactivity or with the received dose.

## Proportion of magmatic and igneous rocks in a 15 km radius

Using the geological map of France (scale : 1/1,000,000, BRGM), the areal proportions of magmatic and igneous rocks in a radius of 15 km around sampling sites were computed (noted λ15), 30 km represents the average distribution range for a subterranean isopod (*Eme et al., 2018*).

## Transcriptome sequencing and assembly

### Sequencing

For each species, we sequenced transcriptomes from eight individuals. For each individual total RNA was isolated using TRI Reagent (Molecular Research Center). Extraction quality was checked on a BioAnalyser RNA chip (Agilent Technologies) and RNA concentrations were estimated using a Qubit fluorometer (Life Technologies). Prior to any additional analysis, species identification was corroborated for each individual by sequencing a fragment of the 16S gene. Illumina libraries were then prepared using the TruSeq RNA Sample Prep Kit v2 (Illumina). For each species one library was paired-end sequenced using 100 cycles, and the seven other libraries were single-end sequenced using 50 cycles on a HiSeq2500 sequencer (Illumina) at the IGBMC GenomEast Platform (Illkirch, France). We obtained around 30 million single-end reads per individual and 118 million paired-end reads per species.

### Assembly

Adapters were clipped from the sequences, low quality read ends were trimmed (phred score <30) and low quality reads were discarded (mean phred score <25 or if remaining length <19 bp) using fastq-mcf of the ea-utils package (*Aronesty, 2013*). Paired-end transcriptomes were de novo assembled using Trinity v2.3.2 (*Grabherr et al., 2011*). Open reading frames (ORFs) were identified with TransDecoder ([http://transdecoder.sourceforge.net](http://transdecoder.sourceforge.net)). For each assembled component, only the most express ORF was retained.

## Families of orthologous genes

Gene families were delimited using an all-against-all BLASTP (*Altschul et al., 1990*) and SiLix (*Miele et al., 2011*) on the ORFs delimited in the previous step. We then kept gene families containing the 14 species, with only one sequence for each species in order to remove paralogs. We obtained 2490 families hereafter considered as one-to-one orthologous genes. These genes were aligned with PRANK (*Löytynoja and Goldman, 2008*) using a codon model and sites ambiguously aligned were removed with Gblocks (*Castresana, 2000*).
## Species tree and gene trees

The 2490 genes were concatenated and a phylogenetic tree (hereafter called the concatenated tree) was built using PhyML v3.0 (Guindon et al., 2010) under a GTR+G+I model with 100 bootstrap replicates and was rooted using the Slavus lineage (Proasellus boui and Proasellus slavus) as an outgroup (Morvan et al., 2013). Most nodes have a bootstrap value of 100% (Figure 1—figure supplement 1). Two nodes have values at 84% and 98% in the clade containing P. nsp VIELVIC; P. nsp HYPOPRAT; P. nsp MONTBAR and P. nsp ROSSFELD. To check the relationship between these four species, we built 2490 individual gene trees with PhyML v3.0 under a GTR+G+I model with 100 bootstrap replicates. Twenty-nine gene trees strongly support (bootstraps >90%) the phylogeny of the concatenated tree for this clade, 208 support other various topologies and the remaining 2253 gene trees do not support any relationship in particular for this clade (bootstraps <90%). Thus, the phylogeny for this clade remained unresolved, possibly as the consequence of a concomitant speciation process of these four species. For approaches with pairs of sister species, as we were unable to resolve the phylogeny for this clade, we selected the species living in the highest level of radioactivity (P. nsp HYPOPRAT) and the species living in the lowest level of radioactivity (P. nsp MONTBAR) among these four species to build a pair, resulting in a total of 6 pairs of sister species (sensu Felsenstein, 2004).

## Mitochondrial genes

Mitochondrial genes were not present amongst the 2490 genes obtained above. Indeed, owing to a different genetic code in invertebrate mitochondria, mitochondrial ORFs were systematically missed by the ORF caller (Transdecoder). We reconstructed mitochondrial genomes using the de novo transcriptome assemblies. Large mitochondrial contigs were built with MITObim (Hahn et al., 2013) by using RNA-seq reads. These contigs were mapped on the assembled mitochondrial genome from the closest possible species (taken from Saclier et al., 2018), allowing us to assemble them. Mitochondrial genomes were annotated using the MITOS web server (Bernt et al., 2013). We recovered the 13 mitochondrial protein-coding genes. Mitochondrial genes were aligned with PRANK (Löytynoja and Goldman, 2008) and sites ambiguously aligned were removed with Gblocks (Castresana, 2000).

## Rate of molecular evolution

We used the synonymous substitution rate ($d_S$) computed on the terminal branches of the tree as a proxy for the long-term species mutation rate (Kimura, 1983). This proxy is valid in absence of selection on codon usage. To check for the absence of biased codon usage, we computed the effective number of codons on the 2490 orthologous genes (ENC, Wright, 1990). This number varies between 20 (only a single codon is used for each amino acid) and 61 (all synonymous codons are used with equal frequency for each amino-acid). ENC ranged between 49.17 and 50.48 (Supplementary file 1), indicating a moderate codon usage bias, more importantly, they do not correlate with alpha radioactivity (pGLS, p.value = 0.6378). Altogether, the $d_S$ estimation does not seem impacted by a strongly biased or variable codon usage.

To compute $d_S$ we first removed some genes showing a conflicting phylogeny. Including genes supporting different phylogenies in a concatenation amounts to constrain a wrong phylogeny for these genes which may biases $d_S$ estimations. Indeed imposing a wrong gene tree will tend to generate convergent mutations in terminal branches of the tree. To avoid such bias in our $d_S$ estimation we used ProfileNJ (Noutahi et al., 2016) with a bootstrap threshold of 90% to compute a cost of reconciliation between the concatenated tree and the gene trees. We kept the gene families with a cost of reconciliation of zero and with sequences long enough for all species (at least a half of the alignment) and removed all other genes, resulting in a set of 769 gene families. $d_S$ were estimated using CoEvol (Lartillot and Poujol, 2011). This software program implements a Muse and Gaut codon model (Muse and Gaut, 1994), with Brownian variation in $d_S$ and $d_N/d_S$ along the tree. Bayesian inference and reconstruction of the history of variation in $d_S$ and $d_N/d_S$ along the tree is conducted by Markov Chain Monte Carlo (MCMC). Two independent chains were run, and were stopped after checking for convergence by eye and with the tracecomp program included in the Coevol package (effective sample size >200 and discrepancy between chains < 0.3). Chains were stopped after 7117 generations (4200 generations excluded as burn-in). The age of the root was
arbitrarily set to 1, resulting in synonymous substitution rate estimates that are relative to the root age (dS/ra) (Supplementary file 1). In order to ensure that assumptions made by CoEvol on the dS evolution along branches don't bias the dS/ra estimation, dS were also computed with CodeML (Yang, 2007) using a free ratio model and with the Bio++ suite (Dutheil and Boussau, 2008). For the last one, a non-homogeneous model (NY98 model) was first applied to the alignment with BppML and then the MapNH program (Version 1.1.1) of the TestNH package (Guéguen and Duret, 2018) was used to reconstruct the ancestral states to estimate the number of synonymous substitutions on each branch. We multiplied the CoEvol dS/ra by the time estimated, to obtain a classical dS. This CoEvol dS was highly correlated with CodeML dS (R² = 0.81) as well as with Mapnh dS (R² = 0.82). Regarding the correlation with radioactivity, by dividing the CodeML dS and the Mapnh dS by the divergence time estimated by CoEvol in order to obtained comparable dS/ra among all species, we obtained similar results whatever the method used to compute the dS/ra (Supplementary file 3).

## Mutational spectrum

To compute the mutational spectrum, we used an approach by pairs of sister species. We determined the polymorphism at the population level for each species by mapping the seven single-end transcriptomes on the assembled paired-end transcriptome with BWA (Li and Durbin, 2009). BAM files were produced with SAMtools (Li et al., 2009), and reads2snps (Gayral et al., 2013) was used to detect polymorphic sites. We then conserved only the 2490 orthologous genes shared by all species to compute the mutational spectrum on the same set of genes.

For the two species of a pair, we reconstructed the ancestral sequence using a parsimonious approach. Namely, for each site in the alignment, if the two species had a single shared allele, this allele was considered as ancestral and the other alleles, if they existed, were considered as derived from the ancestral allele. For each species, we estimated the probability of a mutation in their population, $p(i \to j|f(i)_{anc})$, by counting each type of mutation, either on all positions or on third positions, corrected by the ancestral base frequency:

$$p(i \to j|f(i)_{anc}) = \frac{N(i \to j)}{Ni_{anc}}$$

This probability being dependent on the mutation rate μ, we estimated the mutational spectrum by the proportion, when a mutation occurs, of mutation from the base i to the base j, noted $p(i \to j|\mu, f(i)_{anc})$:

$$p(i \to j|\mu, f(i)_{anc}) = p(i \to j|f(i)_{anc}) * \frac{1}{\sum_{i=\{A,C,G,T\}} \sum_{j=\{A,T,C,G\}} p(i \to j|f(i)_{anc})}$$

This proportion takes into account the mutation rate and is so comparable across species. We pooled complementary mutations (e.g. A to C with T to G) to increase the counts by mutational categories and improve statistical power.

## Statistical analyses

Correlations between dS/ra computed on terminal branch of the tree and the different measures of radioactivity were tested using phylogenetic Generalized Least Squares models (pGLS Martins and Hansen, 1997) with the nlme (Pinheiro et al., 2007) and ape packages (Paradis et al., 2004) in R (R Development Core Team, 2020). As dS/ra is computed on terminal branches of the tree, this test assumes that radioactivity remained stable along the time represented by this branch which is disposable. Natural radioactivity has been slowly decreasing for the last 2 Gy (Karam and Leslie, 2005). However, this decrease is global and we can thus consider that the delta of radioactivity is stable over time. Second, the terminal branches on which the dS/ra is estimated do not exceed 10 My, a time frame in which the level of radioactivity can be considered stable. Thus, the present quantitative estimates should be representative of the radioactivity contrasts that persisted between these habitats in the time-frame of this study. That said, even if subterranean species have low dispersal abilities, there distribution have changed over time, but these movements could blur a signal and are unlikely to generate one.

For the α radioactivity and received dose, the two species showing a disruption in the secular equilibrium were removed. The ultrametric tree built by CoEvol was used to calculate the
phylogenetic variance/covariance matrix under a Brownian motion model to take into account the non-independence among species in the pGLS. Normality of residuals was checked for all models, log transformation was applied when the normality was rejected (Shapiro test).

For mutational spectrum, because the proportions of each mutation are not independent from each other, we proceeded to a selection of variables as proposed by Harris and Pritchard, 2017. This procedure consists of performing iterative chi-square tests. To achieve that, for each mutation we summed the counts of the six species living in 'highly' radioactive environments that we compared to the counts of the six species living in 'weakly' radioactive environments. First, 'normal' chi-square tests are performed for each mutation. Mutations are then ordered following p.values of these tests that we call unordered(p). An ordered p.value is then computed. For the mutation with the lowest unordered(p):m1 the p.value remains unchanged (ordered(p)=unordered(p)). For the following mutation with the second lowest unordered(p): m2, an ordered(p) is computed by doing a chi-square test between the count of this mutation (m2) and the sum of the counts of the other mutations excepted the first one (sum of m3 to m6), and so one for each mutation. For the last mutation (m6), ordered(p) is computed by doing a chi-square test between m5 and m6 counts. To take into account the phylogenetic inertia, we then tested the correlation between the proportion of each mutation and radioactivity (α radioactivity, received dose of radioactivity or proportion of metamorphic or igneous rocks) with a pGLS, this time pruning P. nsp VIELVIC and P. nsp MONTBAR from the chronogram built by CoEvol. P.values from pGLS were corrected with Holm's method to adjust for the multiple tests and to control the false discovery rate.

For each tests, pGLS assumptions, namely normality of residuals, homoscedasticity, absence of influential cases, an evolutionary process (here a brownian motion), have been checked. Normality was checked by plotting residuals of models. Homoscedasticity was checked by plotting the residuals against the fitted values of models. Absence of influential cases was tested by a jackknife approach consisting in removing one by one each species of the data and redoing the test. Robustness to the evolutionary model was tested by performing the pGLS with the four main evolutionary models (Blomberg, Martins, Pagel and Brownian).

## Acknowledgements

This work was supported by the French program STYGOMICS (CNRS Défi Enviromics), the Zone Atelier Territoire Uranifère, the Grottes d'Azé, the Conseil Départemental de Saône-et-Loire, the Association Culturelle du Site d'Azé and the Agence Nationale de la Recherche (ANR-15-CE32-0005 Convergenomix, France, ANR-17-EURE-0018 H2O'Lyon, France). We gratefully acknowledge support from the CNRS/IN2P3 Computing Center (Lyon/Villeurbanne, France) for providing a significant amount of the computing resources needed for this work and the CERESE (Univ. Lyon 1) for the storage of the biological material. We thank the Grottes d'Azé and more specifically Lionel Barriquand, the Mercantour national parc (authorization n°2015-251) and more specifically its scientific manager Marie-France Leccia, the Fédération Rhône-Alpes de Protection de la Nature and the owner of Bout du monde mine for giving us access for sampling. We thank Marcel Meysonnier, Aymeric and Audric Berjoan, Josiane and Bernard Lips, Audrey Brechet, Claude Bou, Benjamin Benti and Léa Dantony for their help in the field. We are grateful to Nicolas Lartillot, Gilles Escarguel, Marie Sémon, Laurent Guéguen, Nicolas Galtier, Benoît Nabholz and Bastien Boussau for helpful discussions. We also thank Laurent Simon and Laura Grice for their suggestions in the latter stages of manuscript preparation.

## Additional information

### Funding

<table><thead><tr><td>Funder</td><td>Grant reference number</td><td>Author</td></tr></thead><tbody><tr><td>Centre National de la Recherche Scientifique</td><td>STYGOMICS - Défi enviromix</td><td>Patrick Chardon<br>Florian Malard<br>Lara Konecny-Dupré<br>Tristan Lefebure<br>Christophe J Douady</td></tr></tbody></table>
<table><tr><td>Agence Nationale de la Recherche</td><td>ANR- 15-CE32-0005<br>Convergenomix</td><td>Lara Konecny-Dupré<br>Laurent Duret<br>Tristan Lefebure<br>Christophe J Douady</td></tr><tr><td>H2O'Lyon, France</td><td>ANR-17-EURE-0018</td><td>Christophe J Douady</td></tr></table>

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

## Author contributions

Nathanaëlle Saclier, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing - original draft, Project administration, Did the field work, defined orthologous gene families, built the phylogeny, computed substitution rate, performed statistical analyses and computed the mutational spectrum; Patrick Chardon, Conceptualization, Data curation, Formal analysis, Funding acquisition, Methodology, Writing - review and editing, Did the field work, performed all radioactivity measurements and the received dose estimation; Florian Malard, Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Project administration, Writing - review and editing, Did the field work, dissected and identified morphologically all sampled species; Lara Konecny-Dupré, Data curation, Formal analysis, Supervision, Writing - review and editing, Extracted DNA and RNA from all samples, made PCR and migrations and prepared library for sequencing; David Eme, Data curation, Validation, Investigation, Writing - review and editing, Did the field work, dissected and identified morphologically all sampled species; Arnaud Bellec, Formal analysis, Writing - review and editing, Extracted the outcrop cover of low-radioactivity sedimentary rocks and high-radioactivity metamorphic and igneous rocks in a radius of 15 km around the sampling, and made the map of uranium for Figure 1; Vincent Breton, Validation, Writing - review and editing; Laurent Duret, Validation, Methodology, Writing - review and editing, Supervised the computation of the mutational spectrum; Tristan Lefebure, Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing, Did the field work, assembled transcriptomes and defined ORFs and orthologous gene families; Christophe J Douady, Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing, Did the field work

## Author ORCIDs

Nathanaëlle Saclier https://orcid.org/0000-0003-1522-9644

David Eme https://orcid.org/0000-0001-8790-0412

Laurent Duret http://orcid.org/0000-0003-2836-3463

Tristan Lefebure https://orcid.org/0000-0003-3923-8166

Christophe J Douady https://orcid.org/0000-0002-4503-8040

## Decision letter and Author response

Decision letter https://doi.org/10.7554/eLife.56830.sa1

Author response https://doi.org/10.7554/eLife.56830.sa2

## Additional files

### Supplementary files

* Supplementary file 1. Fourteen sequenced species with sampling coordinates, synonymous substitution rate relative to the root age for nuclear and mitochondrial genome, α radioactivity measured on each site, effective dose of radioactivity (in µGy/h), and this effective dose corrected for recent human impact (in µGy/h), non-synonymous substitution rate over synonymous substitution rate (dN/dS), the areal proportion of magmatic and metamorphic rocks in a radius of 15 km around the sampling point (λ15), the GC content for all positions or for third positions, the effective number of codon (ENC) and the raw dS computed with the CodeML program and with the MapNH program (not relative to root age).
* Supplementary file 2. Statistical tests for the mutational spectrum analysis. An ordered Chi² test as described in *Harris and Pritchard, 2017* has been performed on mutation counts computed on third positions. A Phylogenetic Generalized Least Square (pGLS) regression of the proportion of each mutations against the α radioactivity measured in sediment (α radio.), the Received Dose (RD) modeled with ERICA tool, and the areal proportion of metamorphic and igneous rock within a 15 km radius (λ15) has also been performed. α radioactivity and RD were log transformed to fit with linear model assumptions. R² are Cox-Snell pseudo R².

* Supplementary file 3. Phylogenetic Least Square (PGLS) regression of *d*<sub>s</sub> computed with CoEvol, CodeMI, or mapNH against radioactivity measured as the α radioactivity measured in sediments, as the effective dose received by organisms or as the areal proportion of metamorphic and magmatic rocks in a radius of 15 km around the sampled point (λ15). Each test corresponds to one likelihood ratio test between the models with and without the given explanatory variable. For α radioactivity and received dose, sampled sites with a break in the secular equilibrium were removed, resulting in tests with only 12 taxa.

* Transparent reporting form

## Data availability

16S sequences have been deposited on the European Nucleotide Archive and are available under the accession numbers from LR214526 to LR214880 (https://www.ebi.ac.uk/ena/data/view/LR214526-LR214880). Alignments and the list of genes used to compute synonymous substitution rate have been deposited on Zenodo (https://zenodo.org/deposit/2563829). Transcriptome reads have been deposited on the European Nucleotide Archive and are available under accession numbers from LR536601 to LR536626 in the study ID PRJEB14193 (https://www.ebi.ac.uk/ena/data/search?query=PRJEB14193). Number of reads and data used for correlations, namely measures of radionuclides and mutations counts have been deposited on Zenodo (https://doi.org/10.5281/zenodo.4071754).

The following datasets were generated:

<table><thead><tr><td>Author(s)</td><td>Year</td><td>Dataset title</td><td>Dataset URL</td><td>Database and Identifier</td></tr></thead><tbody><tr><td>Saclier N, Chardon P, Malard F, Konecny-Dupré L, Emé D, Bellec A, Breton V, Duret L, Lefébure T, Douady CJ</td><td>2020</td><td>Data used in the article</td><td>https://doi.org/10.5281/zenodo.4071754</td><td>Zenodo, 10.5281/zenodo.4071754</td></tr><tr><td>Saclier N, Chardon P, Malard F, Konecny-Dupré L, Emé D, Bellec A, Breton V, Duret L, Lefébure T, Douady CJ</td><td>2020</td><td>Bedrock radioactivity influences the rate and spectrum of mutation - Orthologous genes</td><td>https://doi.org/10.5281/zenodo.2563829</td><td>Zenodo, 10.5281/zenodo.2563829</td></tr><tr><td>Saclier N, Chardon P, Malard F, Konecny-Dupré L, Emé D, Bellec A, Breton V, Duret L, Lefébure T, Douady CJ</td><td>2019</td><td>Aselloidea isopods Sanger sequencing</td><td>https://www.ebi.ac.uk/ena/browser/view/LR214526-LR214880</td><td>ENA, PRJEB30668</td></tr><tr><td>Saclier N, Chardon P, Malard F, Konecny-Dupré L, Emé D, Bellec A, Breton V, Duret L, Lefébure T, Douady CJ</td><td>2019</td><td>Aselloidea isopods transcriptomes sequencing and denovo assembly</td><td>https://www.ebi.ac.uk/ena/data/search?query=PRJEB14193</td><td>ENA, PRJEB14193</td></tr><tr><td>Lefébure T, Saclier N</td><td>2019</td><td>16S sequences</td><td>https://www.ebi.ac.uk/ena/data/view/LR214526-</td><td>EBI European Nucleotide Archive,</td></tr></tbody></table>
<table><tr><td>Lefébure T, Saclier 2019 Transcriptome reads</td><td>LR214880</td><td>LR214526-LR214880</td></tr><tr><td>N</td><td>https://www.ebi.ac.uk/ena/data/view/LR536601-LR536626</td><td>EBI European Nucleotide Archive, LR536601-LR536626</td></tr></table>

## References

* Adewoye AB, Lindsay SJ, Dubrova YE, Hurles ME. 2015. The genome-wide effects of ionizing radiation on mutation induction in the mammalian germline. *Nature Communications* **6**:6684. DOI: https://doi.org/10.1038/ncomms7684, PMID: 25809527

* Allegrucci G, Sbordoni V, Cesaroni D. 2015. Is radon emission in caves causing deletions in satellite DNA sequences of cave-dwelling crickets? *PLOS ONE* **10**:e0122456. DOI: https://doi.org/10.1371/journal.pone.0122456, PMID: 25822625

* Allio R, Donega S, Galtier N, Nabholz B. 2017. Large variation in the ratio of mitochondrial to nuclear mutation rate across animals: implications for genetic diversity and the use of mitochondrial DNA as a molecular marker. *Molecular Biology and Evolution* **34**:2762-2772. DOI: https://doi.org/10.1093/molbev/msx197, PMID: 28981721

* Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. 1990. Basic local alignment search tool. *Journal of Molecular Biology* **215**:403-410. DOI: https://doi.org/10.1016/S0022-2836(05)80360-2, PMID: 2231712

* Aronesty E. 2013. Comparison of sequencing utility programs. *The Open Bioinformatics Journal* **7**:1-8. DOI: https://doi.org/10.2174/1875036201307010001

* Barber R, Plumb MA, Boulton E, Roux I, Dubrova YE. 2002. Elevated mutation rates in the germ line of first- and second-generation offspring of irradiated male mice. *PNAS* **99**:6877-6882. DOI: https://doi.org/10.1073/pnas.102015399, PMID: 11997464

* Barja G. 2002. Endogenous oxidative stress: relationship to aging, longevity and caloric restriction. *Ageing Research Reviews* **1**:397-411. DOI: https://doi.org/10.1016/S1568-1637(02)00008-9, PMID: 12067594

* Beir VII. 2006. Research Council of the National Academies of USA. Health Risks From Exposure to Low Levels of Ionizing Radiation: Beir VII Phase 2. National Academies Press. DOI: https://doi.org/10.17226/11340

* Bernt M, Donath A, Jühling F, Externbrink F, Florentz C, Fritzsch G, Pütz J, Middendorf M, Stadler PF. 2013. MITOS: improved de novo metazoan mitochondrial genome annotation. *Molecular Phylogenetics and Evolution* **69**:313-319. DOI: https://doi.org/10.1016/j.ympev.2012.08.023, PMID: 22982435

* Besaratinia A, Synold TW, Xi B, Pfeifer GP. 2004. G-to-T transversions and small tandem base deletions are the hallmark of mutations induced by ultraviolet a radiation in mammalian cells. *Biochemistry* **43**:8169-8177. DOI: https://doi.org/10.1021/bi049761v, PMID: 15209513

* Bou C, Rouch R. 1967. Un nouveau champ de recherches sur la faune aquatique souterraine. *Comptes Rendus de L'Academie des Sciences* **265**:369-370.

* Boucher D, Hindo J, Averbeck D. 2004. Increased repair of gamma-induced DNA double-strand breaks at lower dose-rate in CHO cells. *Canadian Journal of Physiology and Pharmacology* **82**:125-132. DOI: https://doi.org/10.1139/y04-006, PMID: 15052293

* Boussau B, Blanquart S, Necsulea A, Lartillot N, Gouy M. 2008. Parallel adaptations to high temperatures in the Archaean eon. *Nature* **456**:942-945. DOI: https://doi.org/10.1038/nature07393, PMID: 19037246

* Bromham L, Rambaut A, Harvey PH. 1996. Determinants of rate variation in mammalian DNA sequence evolution. *Journal of Molecular Evolution* **43**:610-621. DOI: https://doi.org/10.1007/BF02202109, PMID: 8995058

* Brown JE, Alfonso B, Avila R, Beresford NA, Copplestone D, Hosseini A. 2016. A new version of the ERICA tool to facilitate impact assessments of radioactivity on wild plants and animals. *Journal of Environmental Radioactivity* **153**:141-148. DOI: https://doi.org/10.1016/j.jenvrad.2015.12.011, PMID: 26773508

* Calvignac S, Konecny L, Malard F, Douady CJ. 2011. Preventing the pollution of mitochondrial datasets with nuclear mitochondrial paralogs (numts). *Mitochondrion* **11**:246-254. DOI: https://doi.org/10.1016/j.mito.2010.10.004, PMID: 21047564

* Castresana J. 2000. Selection of conserved blocks from multiple alignments for their use in phylogenetic analysis. *Molecular Biology and Evolution* **17**:540-552. DOI: https://doi.org/10.1093/oxfordjournals.molbev.a026334, PMID: 10742046

* Chen XJ, Butow RA. 2005. The organization and inheritance of the mitochondrial genome. *Nature Reviews Genetics* **6**:815-825. DOI: https://doi.org/10.1038/nrg1708, PMID: 16304597

* Cheng KC, Cahill DS, Kasai H, Nishimura S, Loeb LA. 1992. 8-Hydroxyguanine, an abundant form of oxidative DNA damage, causes G—T and A—C substitutions. *The Journal of Biological Chemistry* **267**:166-172. DOI: https://doi.org/10.1074/jbc.M409600200, PMID: 1730583

* Clayton DA, Doda JN, Friedberg EC. 1974. The absence of a pyrimidine dimer repair mechanism in mammalian mitochondria. *PNAS* **71**:2777-2781. DOI: https://doi.org/10.1073/pnas.71.7.2777, PMID: 4212385

* Collis SJ, Schwaninger JM, Ntambi AJ, Keller TW, Nelson WG, Dillehay LE, Dewese TL. 2004. Evasion of early cellular response mechanisms following low level radiation-induced DNA damage. *Journal of Biological Chemistry* **279**:49624-49632. DOI: https://doi.org/10.1074/jbc.M409600200, PMID: 15377658

* Cornu J-F, Eme D, Malard F. 2013. The distribution of groundwater habitats in Europe. *Hydrogeology Journal* **21**:949-960. DOI: https://doi.org/10.1007/s10040-013-0984-1
Czeizel AE, Elek C, Susánszky E. 1991. The evaluation of the germinal mutagenic impact of chernobyl radiological contamination in Hungary. *Mutagenesis* **6**:285-288. DOI: https://doi.org/10.1093/mutage/6.4.285, PMID: 1834914

Desouky O, Ding N, Zhou G. 2015. Targeted and non-targeted effects of ionizing radiation. *Journal of Radiation Research and Applied Sciences* **8**:247-254. DOI: https://doi.org/10.1016/j.jrras.2015.03.003

Dikomey E, Brammer I, Johansen J, Bentzen SM, Overgaard J. 2000. Relationship between DNA double-strand breaks, cell killing, and fibrosis studied in confluent skin fibroblasts derived from breast Cancer patients. *International Journal of Radiation Oncology*Biology*Physics **46**:481-490. DOI: https://doi.org/10.1016/S0360-3016(99)00335-1, PMID: 10661357

Dubrova YE, Nesterov VN, Krouchinsky NG, Ostapenko VA, Neumann R, Neil DL, Jeffreys AJ. 1996. Human minisatellite mutation rate after the chernobyl accident. *Nature* **380**:683-686. DOI: https://doi.org/10.1038/380683a0, PMID: 8614461

Duret L, Galtier N. 2009. Biased gene conversion and the evolution of mammalian genomic landscapes. *Annual Review of Genomics and Human Genetics* **10**:285-311. DOI: https://doi.org/10.1146/annurev-genom-082908-150001, PMID: 19630562

Dutheil J, Boussau B. 2008. Non-homogeneous models of sequence evolution in the bio++ suite of libraries and programs. *BMC Evolutionary Biology* **8**:255. DOI: https://doi.org/10.1186/1471-2148-8-255, PMID: 18808672

Einor D, Bonisoli-Alquati A, Costantini D, Mousseau TA, Møller AP. 2016. Ionizing radiation, antioxidant response and oxidative damage: a meta-analysis. *Science of the Total Environment* **548-549**:463-471. DOI: https://doi.org/10.1016/j.scitotenv.2016.01.027

Eme D, Zagmajster M, Fišer C, Galassi D, Marmonier P, Stoch F, Cornu J-F, Oberdorff T, Malard F. 2015. Multi-causality and spatial non-stationarity in the determinants of groundwater crustacean diversity in Europe. *Ecography* **38**:531-540. DOI: https://doi.org/10.1111/ecog.01092

Eme D, Zagmajster M, Delić T, Fišer C, Flot J-F, Konecny-Dupré L, Pálsson S, Stoch F, Zakšek V, Douady CJ, Malard F. 2018. Do cryptic species matter in Macroecology? sequencing european groundwater crustaceans yields smaller ranges but does not challenge biodiversity determinants. *Ecography* **41**:424-436. DOI: https://doi.org/10.1111/ecog.02683

Felsenstein J. 2004. *Inferring Phylogenies*. Sinauer Associates, Inc.

Flibotte S, Edgley ML, Chaudhry I, Taylor J, Neil SE, Rogula A, Zapf R, Hirst M, Butterfield Y, Jones SJ, Marra MA, Barstead RJ, Moerman DG. 2010. Whole-genome profiling of mutagenesis in *Caenorhabditis elegans*. *Genetics* **185**:431-441. DOI: https://doi.org/10.1534/genetics.110.116616, PMID: 20439774

Forster L, Forster P, Lutz-Bonengel S, Willkomm H, Brinkmann B. 2002. Natural radioactivity and human mitochondrial DNA mutations. *PNAS* **99**:13950-13954. DOI: https://doi.org/10.1073/pnas.202400499, PMID: 12370437

Francois CM, Mermillod-Blondin F, Malard F, Fourel F, Lécuyer C, Douady CJ, Simon L. 2016. Trophic ecology of groundwater species reveals specialization in a low-productivity environment. *Functional Ecology* **30**:262-273. DOI: https://doi.org/10.1111/1365-2435.12484

Gayral P, Melo-Ferreira J, Glémin S, Bierre N, Carneiro M, Nabholz B, Lourenco JM, Alves PC, Ballenghien M, Faivre N, Belkhir K, Cahais V, Loire E, Bernard A, Galtier N. 2013. Reference-free population genomics from next-generation transcriptome data and the vertebrate-invertebrate gap. *PLOS Genetics* **9**:e1003457. DOI: https://doi.org/10.1371/journal.pgen.1003457, PMID: 23593039

Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnrirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, et al. 2011. Full-length transcriptome assembly from RNA-Seq data without a reference genome. *Nature Biotechnology* **29**:644-652. DOI: https://doi.org/10.1038/nbt.1883, PMID: 21572440

Greenaway P. 1985. Calcium balance and moulting in the crustacea. *Biological Reviews* **60**:425-454. DOI: https://doi.org/10.1111/j.1469-185X.1985.tb00424.x

Guéguen L, Duret L. 2018. Unbiased estimate of synonymous and nonsynonymous substitution rates with nonstationary base composition. *Molecular Biology and Evolution* **35**:734-742. DOI: https://doi.org/10.1093/molbev/msx308, PMID: 29220511

Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. 2010. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. *Systematic Biology* **59**: 307-321. DOI: https://doi.org/10.1093/sysbio/syq010, PMID: 20525638

Haghdooost S, Sjölander L, Czene S, Harms-Ringdahl M. 2006. The nucleotide pool is a significant target for oxidative stress. *Free Radical Biology and Medicine* **41**:620-626. DOI: https://doi.org/10.1016/j.freeradbiomed.2006.05.003, PMID: 16863995

Hahn C, Bachmann L, Chevreux B. 2013. Reconstructing mitochondrial genomes directly from genomic next-generation sequencing reads-a baiting and iterative mapping approach. *Nucleic Acids Research* **41**:e129. DOI: https://doi.org/10.1093/nar/gkt371, PMID: 23661685

Hainaut P, Pfeifer GP. 2001. Patterns of p53 G→T transversions in lung cancers reflect the primary mutagenic signature of DNA-damage by tobacco smoke. *Carcinogenesis* **22**:367-374. DOI: https://doi.org/10.1093/carcin/22.3.367, PMID: 11238174

Hande MP, Azizova TV, Burak LE, Khokhryakov VF, Geard CR, Brenner DJ. 2005. Complex chromosome aberrations persist in individuals many years after occupational exposure to densely ionizing radiation: an mFISH study. *Genes, Chromosomes and Cancer* **44**:1-9. DOI: https://doi.org/10.1002/gcc.20217

Harris K, Pritchard JK. 2017. Rapid evolution of the human mutation spectrum. *eLife* **6**:e24284. DOI: https://doi.org/10.7554/eLife.24284, PMID: 28440220
Hoeijmakers JH. 2001. Genome maintenance mechanisms for preventing Cancer. *Nature* **411**:366-374.
DOI: https://doi.org/10.1038/35077232, PMID: 11357144

Hollstein M, Sidransky D, Vogelstein B, Harris CC. 1991. p53 mutations in human cancers. *Science* **253**:49-53.
DOI: https://doi.org/10.1126/science.1905840, PMID: 1905840

Hooker AM, Bhat M, Day TK, Lane JM, Swinburne SJ, Morley AA, Sykes PJ. 2004. The linear no-threshold model does not hold for low-dose ionizing radiation. *Radiation Research* **162**:447-452. DOI: https://doi.org/10.1667/RR3228, PMID: 15447037

Huo L, Nagasawa H, Little JB. 2001. *HPRT* mutants induced in bystander cells by very low fluences of alpha particles result primarily from point mutations. *Radiation Research* **156**:521-525. DOI: https://doi.org/10.1667/0033-7587(2001)156[0521:HMIIBC]2.0.CO;2, PMID: 11604065

Ielsch G, Cuney M, Buscail F, Rossi F, Leon A, Cushing ME. 2017. Estimation and mapping of uranium content of geological units in France. *Journal of Environmental Radioactivity* **166**:210-219. DOI: https://doi.org/10.1016/j.jenvrad.2016.05.022, PMID: 27266726

Jacobsen NR, White PA, Gingerich J, Møller P, Saber AT, Douglas GR, Vogel U, Wallin H. 2011. Mutation spectrum in FE1-MUTA(TM) Mouse lung epithelial cells exposed to nanoparticulate carbon black. *Environmental and Molecular Mutagenesis* **52**:331-337. DOI: https://doi.org/10.1002/em.20629, PMID: 20963790

Jarrett SG, Boulton ME. 2005. Antioxidant up-regulation and increased nuclear DNA protection play key roles in adaptation to oxidative stress in epithelial cells. *Free Radical Biology and Medicine* **38**:1382-1391. DOI: https://doi.org/10.1016/j.freeradbiomed.2005.02.003, PMID: 15855056

Jostes RF, Fleck EW, Morgan TL, Stiegler GL, Cross FT. 1994. Southern blot and polymerase chain reaction exon analyses of HPRT- mutations induced by radon and radon progeny. *Radiation Research* **137**:371-379.
DOI: https://doi.org/10.2307/3578712, PMID: 8146281

Kam WW, Banati RB. 2013. Effects of ionizing radiation on mitochondria. *Free Radical Biology & Medicine* **65**: 607-619. DOI: https://doi.org/10.1016/j.freeradbiomed.2013.07.024, PMID: 23892359

Karam PA, Leslie SA. 2005. Changes in terrestrial natural radiation levels over the history of life. In: Karam P. A (Ed). *Radioactivity in the Environment*. Elsevier. p. 107-117. DOI: https://doi.org/10.1016/S1569-4860(04)07011-1

Karran P. 2000. DNA double strand break repair in mammalian cells. *Current Opinion in Genetics & Development* **10**:144-150. DOI: https://doi.org/10.1016/S0959-437X(00)00069-1, PMID: 10753787

Kimura M. 1983. *The Neutral Theory of Molecular Evolution*. Cambridge University Press. DOI: https://doi.org/10.1002/ajpa.1330650413

Krueger SA, Joiner MC, Weinfeld M, Piasentin E, Marples B. 2007. Role of apoptosis in low-dose hyper-radiosensitivity. *Radiation Research* **167**:260-267. DOI: https://doi.org/10.1667/RR0776.1, PMID: 17316076

Lartillot N, Poujol R. 2011. A phylogenetic model for investigating correlated evolution of substitution rates and continuous phenotypic characters. *Molecular Biology and Evolution* **28**:729-744. DOI: https://doi.org/10.1093/molbev/msq244, PMID: 20926596

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, 1000 Genome Project Data Processing Subgroup. 2009. The sequence alignment/Map format and SAMtools. *Bioinformatics* **25**:2078-2079. DOI: https://doi.org/10.1093/bioinformatics/btp352, PMID: 19505943

Li H, Durbin R. 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics* **25**:1754-1760. DOI: https://doi.org/10.1093/bioinformatics/btp324, PMID: 19451168

Ljungman M. 1991. The influence of chromatin structure on the frequency of radiation-induced DNA strand breaks: a study using nuclear and nucleoid monolayers. *Radiation Research* **126**:58-64. DOI: https://doi.org/10.2307/3578171, PMID: 2020739

Ljungman M, Hanawalt PC. 1992. Efficient protection against oxidative DNA damage in chromatin. *Molecular Carcinogenesis* **5**:264-269. DOI: https://doi.org/10.1002/mc.2940050406, PMID: 1323299

Loucas BD, Eberle R, Bailey SM, Cornforth MN. 2004. Influence of dose rate on the induction of simple and complex chromosome exchanges by gamma rays. *Radiation Research* **162**:339-349. DOI: https://doi.org/10.1667/RR3245, PMID: 15447049

Löytynoja A, Goldman N. 2008. Phylogeny-aware gap placement prevents errors in sequence alignment and evolutionary analysis. *Science* **320**:1632-1635. DOI: https://doi.org/10.1126/science.1158395, PMID: 18566285

Lynch M. 2010. Evolution of the mutation rate. *Trends in Genetics* **26**:345-352. DOI: https://doi.org/10.1016/j.tig.2010.05.003, PMID: 20594608

Martins EP, Hansen TF. 1997. Phylogenies and the comparative method: a general approach to incorporating phylogenetic information into the analysis of interspecific data. *The American Naturalist* **149**:646-667.
DOI: https://doi.org/10.1086/286013

Masoomi JR, Mohammadi S, Amini M, Ghiassi-Nejad M. 2006. High background radiation Areas of ramsar in Iran: evaluation of DNA damage by alkaline single cell gel electrophoresis (SCGE). *Journal of Environmental Radioactivity* **86**:176-186. DOI: https://doi.org/10.1016/j.jenvrad.2005.08.005, PMID: 16376699

Mathieson I, Reich D. 2017. Differences in the rare variant spectrum among human populations. *PLOS Genetics* **13**:e1006581. DOI: https://doi.org/10.1371/journal.pgen.1006581, PMID: 28146552

Miele V, Penel S, Duret L. 2011. Ultra-fast sequence clustering from similarity networks with SiLiX. *BMC Bioinformatics* **12**:116. DOI: https://doi.org/10.1186/1471-2105-12-116, PMID: 21513511

Møller AP, Mousseau TA. 2013. The effects of natural variation in background radioactivity on humans, animals and other organisms. *Biological Reviews* **88**:226-254. DOI: https://doi.org/10.1111/j.1469-185X.2012.00249.x, PMID: 23136873
Morvan C, Malard F, Paradis E, Lefébure T, Konecny-Dupré L, Douady CJ. 2013. Timetree of aselloidea reveals species diversification dynamics in groundwater. *Systematic Biology* **62**:512–522. DOI: https://doi.org/10.1093/sysbio/syt015, PMID: 23461879

Muse SV, Gaut BS. 1994. A likelihood approach for comparing synonymous and nonsynonymous nucleotide substitution rates, with application to the chloroplast genome. *Molecular Biology and Evolution* **11**:715–724. DOI: https://doi.org/10.1093/oxfordjournals.molbev.a040152, PMID: 7968485

Nabholz B, Glémin S, Galtier N. 2008. Strong variations of mitochondrial mutation rate across mammals—the longevity hypothesis. *Molecular Biology and Evolution* **25**:120–130. DOI: https://doi.org/10.1093/molbev/msm248, PMID: 17998254

Nik-Zainal S, Kucab JE, Morganella S, Glodzik D, Alexandrov LB, Arlt VM, Weninger A, Hollstein M, Stratton MR, Phillips DH. 2015. The genome as a record of environmental exposure. *Mutagenesis* **30**:763–770. DOI: https://doi.org/10.1093/mutage/gev073, PMID: 26443852

Nikolaev SI, Montoya-Burgos JI, Popadin K, Parand L, Margulies EH, Antonarakis SE, National Institutes of Health Intramural Sequencing Center Comparative Sequencing Program. 2007. Life-history traits drive the evolutionary rates of mammalian coding and noncoding genomic elements. *PNAS* **104**:20443–20448. DOI: https://doi.org/10.1073/pnas.0705658104, PMID: 18077382

Noutahi E, Semeria M, Lafond M, Seguin J, Boussau B, Guéguen L, El-Mabrouk N, Tannier E. 2016. Efficient gene tree correction guided by genome evolution. *PLOS ONE* **11**:e0159559. DOI: https://doi.org/10.1371/journal.pone.0159559, PMID: 27513924

Paradis E, Claude J, Strimmer K. 2004. APE: analyses of phylogenetics and evolution in R language. *Bioinformatics* **20**:289–290. DOI: https://doi.org/10.1093/bioinformatics/btg412, PMID: 14734327

Pinheiro J, Bates D, DebRoy S, Sarkar D, R Core Team. 2007. Linear and nonlinear mixed effects models. *R Package Version* **3**:1–89. DOI: https://doi.org/10.4148/2475-7772.1273

Pinto E, Sigaud-kutner TCS, Leitao MAS, Okamoto OK, Morse D, Colepicolo P. 2003. Heavy metal-induced oxidative stress in algae1. *Journal of Phycology* **39**:1008–1018. DOI: https://doi.org/10.1111/j.0022-3646.2003.02-193.x

Quinlan GJ, Halliwell B, Moorhouse CP, Gutteridge JM. 1988. Action of lead(II) and aluminium (III) ions on iron-stimulated lipid peroxidation in liposomes, erythrocytes and rat liver microsomal fractions. *Biochimica Et Biophysica Acta (BBA) - Lipids and Lipid Metabolism* **962**:196–200. DOI: https://doi.org/10.1016/0005-2760(88)90159-2, PMID: 3167077

Quinn G, Keough M. 2002. Multiple testing. *Experimental Design and Data Analysis for Biologists* **1**:48–50. DOI: https://doi.org/10.1017/cbo9780511806384

R Development Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria, R Foundation for Statistical Computing. https://www.R-project.org/

Richter C, Park JW, Ames BN. 1988. Normal oxidative damage to mitochondrial and nuclear DNA is extensive. *PNAS* **85**:6465–6467. DOI: https://doi.org/10.1073/pnas.85.17.6465, PMID: 3413108

Rigaud O, Laquerbe A, Moustacchi E. 1995. DNA sequence analysis of HPRT- mutants induced in human lymphoblastoid cells adapted to ionizing radiation. *Radiation Research* **144**:181–189. DOI: https://doi.org/10.2307/3579257, PMID: 7480644

Rigaud O, Moustacchi E. 1996. Radioadaptation for gene mutation and the possible molecular mechanisms of the adaptive response. *Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis* **358**:127–134. DOI: https://doi.org/10.1016/S0027-5107(96)00113-3, PMID: 8946017

Rothkamm K, Kühne M, Jeggo PA, Löbrich M. 2001. Radiation-Induced genomic rearrangements formed by nonhomologous End-Joining of DNA Double-Strand breaks. *Cancer Research* **61**:3886–3893. PMID: 11358801

Saclier N, François CM, Konecny-Dupré L, Lartillot N, Guéguen L, Duret L, Malard F, Douady CJ, Lefébure T. 2018. Life history traits impact the nuclear rate of substitution but not the mitochondrial rate in isopods. *Molecular Biology and Evolution* **35**:2900–2912. DOI: https://doi.org/10.1093/molbev/msy184, PMID: 30247705

Saghirzadeh M, Gharaati MR, Mohammadi S, Ghiassi-Nejad M. 2008. Evaluation of DNA damage in the root cells of Allium cepa seeds growing in soil of high background radiation Areas of Ramsar-Iran. *Journal of Environmental Radioactivity* **99**:1698–1702. DOI: https://doi.org/10.1016/j.jenvrad.2008.03.013, PMID: 18524436

Satoh C, Takahashi N, Asakawa J, Kodaira M, Kuick R, Hanash SM, Neel JV. 1996. Genetic analysis of children of atomic bomb survivors. *Environmental Health Perspectives* **104 Suppl 3**:511–519. DOI: https://doi.org/10.1289/ehp.96104s3511, PMID: 8781374

Sawada S. 2007. Cover-up of the effects of internal exposure by residual radiation from the atomic bombing of hiroshima and Nagasaki. *Medicine, Conflict and Survival* **23**:58–74. DOI: https://doi.org/10.1080/13623690601084617, PMID: 17370859

Seo KY, Jelinsky SA, Loechler EL. 2000. Factors that influence the mutagenic patterns of DNA adducts from chemical carcinogens. *Mutation Research/Reviews in Mutation Research* **463**:215–246. DOI: https://doi.org/10.1016/S1383-5742(00)00047-8, PMID: 11018743

Shibutani S, Takeshita M, Grollman AP. 1991. Insertion of specific bases during DNA synthesis past the oxidation-damaged base 8-oxodG. *Nature* **349**:431–434. DOI: https://doi.org/10.1038/349431a0, PMID: 1992344

Shigenaga MK, Gimeno CJ, Ames BN. 1989. Urinary 8-hydroxy-2'-deoxyguanosine as a biological marker of in vivo oxidative DNA damage. *PNAS* **86**:9697–9701. DOI: https://doi.org/10.1073/pnas.86.24.9697, PMID: 2602371
Smith SA, Donoghue MJ. 2008. Rates of molecular evolution are linked to life history in flowering plants. *Science* **322**:86-89. DOI: https://doi.org/10.1126/science.1163197, PMID: 18832643

Thomas JA, Welch JJ, Lanfear R, Bromham L. 2010. A generation time effect on the rate of molecular evolution in invertebrates. *Molecular Biology and Evolution* **27**:1173–1180. DOI: https://doi.org/10.1093/molbev/msq009, PMID: 20083649

Tubiana M, Aurengo A, Averbeck D, Masse R. 2006. The debate on the use of linear no threshold for assessing the effects of low doses. *Journal of Radiological Protection* **26**:317–324. DOI: https://doi.org/10.1088/0952-4746/26/3/N01, PMID: 16926474

Tubiana M, Feinendegen LE, Yang C, Kaminski JM. 2009. The linear no-threshold relationship is inconsistent with radiation biologic and experimental data. *Radiology* **251**:13–22. DOI: https://doi.org/10.1148/radiol.2511080671, PMID: 19332842

van Gent DC, Hoeijmakers JH, Kanaar R. 2001. Chromosomal stability and the DNA double-stranded break connection. *Nature Reviews. Genetics* **2**:196–206. DOI: https://doi.org/10.1038/35056049, PMID: 11256071

Vilenchik MM, Knudson AG. 2000. Inverse radiation dose-rate effects on somatic and germ-line mutations and DNA damage rates. *PNAS* **97**:5381–5386. DOI: https://doi.org/10.1073/pnas.090099497, PMID: 10792040

Wallace SS. 1998. Enzymatic processing of radiation-induced free radical damage in DNA. *Radiation Research* **150**:S60. DOI: https://doi.org/10.2307/3579809, PMID: 9806610

Ward JF. 1988. DNA Damage Produced by Ionizing Radiation in Mammalian Cells: Identities, Mechanisms of Formation, and Reparability. In: Cohn W. E, Moldave K (Eds). *Progress in Nucleic Acid Research and Molecular Biology*. Academic Press. p. 95–125. DOI: https://doi.org/10.1016/s0079-6603(08)60611-x

Weller C, Wu M. 2015. A generation-time effect on the rate of molecular evolution in Bacteria. *Evolution* **69**:643–652. DOI: https://doi.org/10.1111/evo.12597

Winegar RA, Lutze LH, Hamer JD, O'Loughlin KG, Mirsalis JC. 1994. Radiation-induced point mutations, deletions and micronuclei in lacI transgenic mice. *Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis* **307**:479–487. DOI: https://doi.org/10.1016/0027-5107(94)90258-5, PMID: 7514721

Wood ML, Dizdaroglu M, Gajewski E, Essigmann JM. 1990. Mechanistic studies of ionizing radiation and oxidative mutagenesis: genetic effects of a single 8-hydroxyguanine (7-hydro-8-oxoguanine) residue inserted at a unique site in a viral genome. *Biochemistry* **29**:7024–7032. DOI: https://doi.org/10.1021/bi00482a011, PMID: 2223758

Wright F. 1990. The 'effective number of codons' used in a gene. *Gene* **87**:23–29. DOI: https://doi.org/10.1016/0378-1119(90)90491-9, PMID: 2110097

Yamamori T, Yasui H, Yamazumi M, Wada Y, Nakamura Y, Nakamura H, Inanami O. 2012. Ionizing radiation induces mitochondrial reactive oxygen species production accompanied by upregulation of mitochondrial electron transport chain function and mitochondrial content under control of the cell cycle checkpoint. *Free Radical Biology and Medicine* **53**:260–270. DOI: https://doi.org/10.1016/j.freeradbiomed.2012.04.033, PMID: 22580337

Yang Z. 2007. PAML 4: phylogenetic analysis by maximum likelihood. *Molecular Biology and Evolution* **24**:1586–1591. DOI: https://doi.org/10.1093/molbev/msm088, PMID: 17483113

Ziegler A, Leffell DJ, Kunala S, Sharma HW, Gailani M, Simon JA, Halperin AJ, Baden HP, Shapiro PE, Bale AE. 1993. Mutation hotspots due to sunlight in the p53 gene of nonmelanoma skin cancers. *PNAS* **90**:4216–4220. DOI: https://doi.org/10.1073/pnas.90.9.4216, PMID: 8483937