# Mitigation of the diagenesis risk in biological apatite
## δ18O interpretation

Christophe Lécuyer, Jean-Pierre Flandrois

► To cite this version:

Christophe Lécuyer, Jean-Pierre Flandrois. Mitigation of the diagenesis risk in biological apatite δ18O interpretation. Palaeogeography, Palaeoclimatology, Palaeoecology, 2023, 630, pp.111812. 10.1016/j.palaeo.2023.111812. hal-04616065

HAL Id: hal-04616065

https://hal.science/hal-04616065v1

Submitted on 24 Jul 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Mitigation of the diagenesis risk in biological apatite δ¹⁸O interpretation

2

3

4

Christophe Lécuyer¹ and Jean-Pierre Flandrois²

5

6

7 ¹LGL-TPE, CNRS UMR5276, ENSL, Univ Lyon, Univ Lyon 1, Institut Universitaire de

8 France, 43 bd du 11 Novembre 1918, 69622 Villeurbanne, France.

9

10 ²LBBE, CNRS UMR5558, Univ Lyon, Univ Lyon 1, 43 bd du 11 Novembre 1918, 69622

11 Villeurbanne, France. E-mail : jean-pierre.flandrois@univ-lyon1.fr.

12

13 Corresponding author :

14 Christophe Lécuyer. E-mail : christophe.lecuyer@univ-lyon1.fr. ORCID : 0000-0001-6462-

15 8727

16

17
18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

**Abstract** – The assessment of the potential blurring effect of diagenetic processes represents the greatest challenge for geochemists whose purpose is to reconstruct environmental parameters based on isotopic proxies. Numerous studies performed on phosphatic remains from terrestrial and marine sedimentary deposits have shown post-depositional isotopic modifications due to exchange reactions with surrounding waters. If diagenesis may be diagnosed through several lines of evidence such as mineralogical or chemical modifications of the pristine *in vivo* apatite after burial in soils and sediments, diagenesis does not prove that bioapatite phosphate exchanged oxygen isotopes with aqueous fluids. Therefore, we propose a new method of detection of diagenetic alteration of bioapatite phosphate based on a statistical approach designed to analyze sample collections. The approach requires that phosphatic remains belonged to animals living during a short time interval and within an environment characterized by steady air or water temperature and the δ¹⁸O of ambient water.

We show that the distribution of bioapatite phosphate δ¹⁸O follows a normal distribution resulting from the combination of independent factors such as those related to the genetics, environment and analytical measurements. In the case of oxygen isotope exchange between the bioapatite phosphate and various amounts of circulating aqueous fluids, our numerical simulations reveal that the distributions of bioapatite phosphate δ¹⁸O are asymmetric with shapes depending on the temperature of reaction and the composition of the reacting aqueous fluid. This new statistical tool based on mass balance equations and Monte Carlo simulations is considered valuable to mitigate the diagenesis risk in the interpretation of bioapatite phosphate δ¹⁸O as paleoecological and paleoclimatic proxies.

**Keywords:** apatite; vertebrate; diagenesis; oxygen isotope; Monte Carlo simulation; paleoclimate
44

45

# 1. Introduction

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

Since the pioneering studies published devoted to biogenic apatites (Longinelli and Nuti, 1973a,b; Kolodny et al., 1983; Longinelli, 1984), the oxygen isotope compositions of the phosphate group in biogenic apatites have gradually become essential for reconstructing paleoenvironments and paleoclimates in both terrestrial and marine environments (Kolodny and Raab, 1988; Fricke and O'Neil, 1996; Lécuyer et al., 1996; Amiot et al., 2004; Tütken et al., 2007; Joachimski et al., 2012; Goedert et al., 2017; Rey et al., 2018). The thermo-dependence of the isotopic fractionation between phosphate and water (Kolodny et al., 1983; Blake et al., 1997; Lécuyer et al., 1999; 2013) has also allowed the thermophysiology and ecology of large groups of extinct vertebrates such as dinosaurs or giant marine reptiles of the Phanerozoic (Barrick and Showers, 1994; Fricke and Rogers, 2000; Amiot et al., 2006; Tütken and Vennemann, 2009; Bernard et al., 2010; Faure-Brac et al., 2022). In the same way, Rey et al. (2017) also identified Permian mammalian reptiles as probably the first terrestrial vertebrates able to regulate their body temperature. Phosphate oxygen isotope ratios have also helped elucidate the ecology of certain dinosaurs, such as spinosaurs, which had to spend much of their time partially submerged in freshwater, lowering their body temperature compared to strictly terrestrial theropods (Amiot et al., 2010). This natural isotopic tracer is also used in archaeology, often in conjunction with Sr radiogenic isotopes, to identify the origins of human populations and their migratory routes (Fricke et al., 1995; White et al., 1998; Evans et al., 2006; Touzeau et al., 2013; Lightfoot and O'Connell, 2016; Pellegrini et al., 2016; Clauzel et al., 2023).

67

68

However, as with carbonates of biogenic origin, bioapatites may be subject to physical (ultrastructure, texture, decomposition of organic tissues) and chemical (mineralogy,
geochemistry) modification processes that occur post-mortem when skeletal elements are buried in soils or sediments. This process is known as early diagenesis, and can be accompanied by changes in the original isotopic composition of the vertebrate bone or tooth in the presence of an aqueous solution and possibly microorganisms involved in dissolving the mineral phase. These processes of alteration or diagenesis of bioapatites have been the subject of numerous studies (Reynard et al., 1999; Trueman and Tuross, 2002; Zazzo et al., 2004a,b; Pfretzschner, 2006; Tütken et al., 2008; White and Booth, 2014; Keenan, 2016; Turner-Walker et al., 2023). They continue to mobilize the scientific community involved, as understanding the mechanisms of diagenesis is essential for the correct interpretation of the fossil record of Phanerozoic terrestrial and marine vertebrates.

However, assessing the potential blurring effect of diagenetic processes on the pristine oxygen isotope compositions of apatite phosphate represents the greatest challenge for geochemists to reconstruct environmental parameters. Indeed, numerous studies performed on phosphatic remains from terrestrial and marine sedimentary deposits have shown post-depositional isotopic modifications due to exchange reactions with surrounding waters (e.g. Longinelli, 1966; McArthur et al., 1980; McArthur and Herczeg, 1990; Shemesh, 1990; Sánchez Chillón et al., 1994; Kolodny et al., 1996; Sharp et al., 2000). The loss of original oxygen isotope compositions of phosphates does not appear to be really age-dependent but more related to the post-depositional environment and taphonomy (Shemesh, 1990; Kolodny and Luz, 1992).

Various experimental approaches have been performed to simulate the alteration processes in order to predict the direction and magnitude of the oxygen isotope exchange that may take place between the biomineral and the water circulating after burial in soils or sediments (Blake et al., 1997; 1998; 2005; Paytan et al., 2002; O'Neil et al., 2003; Geisler et al., 2009; Jaisi et al., 2011; Sperber et al., 2014; Chang and Blake, 2015; Aufort et al., 2019).
On the one hand, under inorganic conditions and at ambient temperature, it has been shown that the oxygen isotope exchange between dissolved phosphate and water is very slow (the rate constant k = 10⁻¹³ to 10⁻¹¹ s⁻¹), which precludes any significant change in the isotopic composition of a Ca-phosphate mineral at the 10³–10⁴ year- scale (Lécuyer et al., 1999). On the other hand, in the course of biologically-mediated reactions, the oxygen isotope compositions of apatite can be shifted by several δ units as it was observed by using modern bone and enamel powders, which were reacted at ambient temperatures with ¹³C- and ¹⁸O-labelled waters (Zazzo et al., 2004a). These experiments revealed that any interpretation of apatite δ¹⁸O values must remain careful.

Empirical approaches are also commonly performed by researchers to assess the preservation state of the initial isotopic signatures of bioapatite phosphate as listed below:

- Trace element ratios including REE that may reveal dissolution and recrystallization of the apatite crystals (Shemesh, 1990; Reynard et al., 1999; Trueman and Tuross, 2002). However, changes in REE contents do not necessarily change the oxygen isotope composition of apatite phosphate as they either adsorb on crystal surfaces or substitute to calcium in the apatite crystal lattice (Trueman et al., 2011),

- Raman or FTIR spectroscopy to calculate the size of apatite crystals. Indeed, the nanocrystals of apatite are thermodynamically unstable and are prone to recrystallization after animal's death. Bioapatite crystals tend to recrystallize by Ostwald ripening into a phosphate phase characterized by a more ordered crystal structure with larger crystal dimensions (Shemesh, 1990; Pucéat et al., 2004; Hedges, 2002; Berna et al., 2004; Dal Sasso et al., 2018),

- The amount of carbonate in substitution to phosphate in the crystal lattice of the apatite (Lécuyer et al., 1998; Keenan and Engel, 2017; Végh et al., 2022). However, Vennemann et
118

al. (2001) have shown that bioapatite phosphate may keep its original oxygen isotope composition for a $CO_3^{2-}$ content ranging from 1.8 wt% to 13.4 wt%.

119

120

- The $\delta^{18}O$ of the carbonate compound compared to that of phosphate in the case of terrestrial mammals (Iacumin et al., 1996; Zazzo et al., 2004b). However, Vennemann et al. (2001) have shown that the difference between the $\delta^{18}O$ values of the carbonate and phosphate compounds range from 6.5‰ to 14.7‰ in modern shark teeth.

121

122

123

124

However, it must be kept in mind that all the mineralogical and geochemical tests mentioned above do not provide definitive evidence of preservation of the *in vivo* oxygen isotope composition of apatite phosphate noted here $\delta^{18}O(A)$. Therefore, we propose a new method of detection of diagenetic alteration of bioapatites based on a statistical approach designed to analyze phosphatic remains that belonged to animals living during a short time interval and within an environment characterized by similar parameters such as the temperature (T) and the oxygen isotope composition of ambient water noted $\delta^{18}O(W)$. This new tool, developed in this study, is useful to detect scrambled or partly modified post-depositional oxygen isotope compositions of bioapatite (skeletal remains) in association with other taphonomical, mineralogical and geochemical criteria.

125

126

127

128

129

130

131

132

133

134

135

## 2. Methods

136

### 2.1. Principles of oxygen isotope exchange between water and apatite phosphate based on mass balance considerations

137

138

139

Of specific interest, both direction and magnitude of changes in the $\delta^{18}O(A)$ values of bioapatite can also be explored by using simple mass balance calculations that involve a temperature-dependent isotopic fractionation. Firstly, we consider a law of mass conservation between the initial and final stages of the water-apatite isotopic exchange reaction. Secondly,

140

141

142
we also assume that the hydration of the mineral during its interaction with an external water source is negligible compared to the total amount of circulating water. Knowing that apatite tissues contain about 10-20 vol% of water according to Welborn (2021) and Surowiec et al. (2022), our model only considers the isotope exchange between the phosphate molecule and an external water source. Thus, the following relationship (eq. 1) is established as follows:

$$W.[O]_w \cdot \delta^{18}O(W)_i + A.[O]_A \cdot \delta^{18}O(A)_i = W.[O]_w \cdot \delta^{18}O(W)_f + A.[O]_A \cdot \delta^{18}O(A)_f$$

(1)

with W = water mass and [O]$_w$ the oxygen molar ratio; A = apatite mass and [O]$_A$ the oxygen molar ratio; $\delta^{18}O(W)_i$ is the oxygen isotope composition of water before reaction and $\delta^{18}O(W)_f$ is the oxygen isotope composition of water after reaction; $\delta^{18}O(A)_i$ is the oxygen isotope composition of apatite phosphate before reaction and $\delta^{18}O(A)_f$ is the oxygen isotope composition of apatite phosphate after reaction. This formulation describes the oxygen isotope exchange that takes place between apatite phosphate and water in closed system, which means that the apatite reacts continuously with a given amount of water until the isotopic equilibrium is reached.

Re-arranging the terms in eq. 1 allows the water-apatite ratio (W/A) to be expressed (eq. 2) as follows:

$$\frac{W}{A(\text{closed})} = \frac{[O]_A}{[O]_W} \times \left[ \frac{\delta^{18}O(A)_f - \delta^{18}O(A)_i}{\delta^{18}O(W)_i - \delta^{18}O(W)_f} \right] \quad (2)$$
In the fossil record, generally, $\delta^{18}O(W)_f$ is not known but it can be replaced by an expression that involves the temperature-dependent oxygen isotope fractionation between apatite and water. Indeed, $1000\ln\alpha_{\text{apatite-water}}$ can be classically approximated by the difference between $\delta^{18}O(A)_f$ and $\delta^{18}O(W)_f$, also noted $\Delta$, for which temperature dependency is expressed (eq. 3) according to the equation determined by Lécuyer et al. (2013):

$$ \Delta = \frac{117.4-T}{4.50} \qquad (3) $$

Consequently, eq. 2 is re-arranged to the following form (eq. 4):

$$ \frac{W}{A(\text{closed})} = \frac{[O]_A}{[O]_W} \times \left[ \frac{\delta^{18}O(A)_f - \delta^{18}O(A)_i}{\delta^{18}O(W)_i - \delta^{18}O(A)_f + \left\{ \frac{117.4-T}{4.50} \right\}} \right] \qquad (4) $$

Here, we consider the hypothetical case where a fish lived in a water mass characterized by a homogenous temperature of 18°C and a $\delta^{18}O$ value of 0‰ (SMOW), its teeth should have a $\delta^{18}O(A)_i$ value close to 22‰ according to equation (3). After the burial of teeth in sediments, let us consider that the tooth apatite starts reacting with meteoric waters having a $\delta^{18}O(W)_i$ value of -8‰ (SMOW), the mass balance equation developed above predicts that the oxygen isotope compositions of the teeth ($\delta^{18}O(A)_f$) decrease by several $\delta$ units depending on the amount and temperature of the reacting water (Figure 1).

Simon et al. (2017) proposed a general dimensionless mass balance equation for single-pass open systems that describes the equilibrium elemental distribution and the isotopic
composition of reacting rocks as a function of the amount of circulating water. The isotopic ratio k₂/k₁ of the altered rocks is calculated considering a dynamic single-pass open system where infinitesimal amounts of water dW, of fixed elemental concentration C<sub>w</sub><sup>i</sup> and isotopic ratio R<sub>w</sub><sup>i</sup>, react completely with the rock, then leaves the system. The mass fraction of water (f<sub>w</sub>) and rock (f<sub>r</sub>) in the bulk system are assumed to be constant. Unlike Taylor (1977; 1978) and Criss and Taylor (1986), Simon et al. (2017) considered that the bulk concentration of the element k may vary during the process of water-rock interaction. Simon et al. (2017) also assumed that the atomic weight of the element k is identical in both rock and water. Except for unrealistic differences in isotopic ratios, this approximation is generally valid. The evolution of the ratio of k₂ over k₁ is thus described (eq. 5) by the following first-order differential equation:

$$d[M_T(f_w \cdot C_w^{k_1} \cdot R_w + f_r \cdot C_r^{k_1} \cdot R_r)] = (C_w^{k_1})^i \cdot R_w^i \cdot dW - C_w^{k_1} \cdot R_w \cdot dW_i \quad (5)$$

where C<sub>w</sub><sup>k1</sup> and C<sub>r</sub><sup>k1</sup> are the concentrations of isotope k₁, k₂, R<sub>w</sub> and R<sub>r</sub> the isotopic ratio k₂/k₁, in water and rock, respectively. The i exponent notes initial values, and M<sub>T</sub> is the total mass of the system with f<sub>r</sub> + f<sub>w</sub> = 1.

In the case of exchanging <sup>18</sup>O and <sup>16</sup>O between apatite and water in an open system, the general mass balance equation (eq. 5) is simplified and has the following analytical solution (eq. 6), which has been already reported, for example, in Criss and Taylor (1986):

$$\frac{W}{A(\text{open})} = \frac{[O]_A}{[O]_W} \ln \left( \frac{[O]_W}{[O]_A} \times \frac{W}{A(\text{closed})} + 1 \right) \quad (6)$$
Equation (6) was applied to the hypothetical fish case mentioned above and the model output is shown in Figure 2. In open water/apatite systems, less water is required to achieve the same shift in oxygen isotope composition from the pristine to the altered sample. For high W/A ratios (> 1) both expressions are reduced to equation 7 as follows:

$$1000\ln\alpha(PO_4-H_2O) \approx \delta^{18}O(A)f - \delta^{18}O(W)f = f(T) \quad (7)$$

With $\alpha(PO_4-H_2O)$ representing the oxygen isotope fractionation factor between apatite and water. Note that $\delta^{18}O(W)_f \approx \delta^{18}O(W)_i$ when W/A >> 1 as the apatite oxygen reservoirs becomes negligible compared to that of water (Criss and Taylor, 1986). We mention here that the international reference for stable oxygen isotopes in apatite phosphate is VSMOW; the ‰ unit being implicit when a numerical value is assigned to the $\delta$ notation.

## 2.2. Why normal distributions are expected for samples collected from a sedimentary layer deposited over a short time interval?

A normal distribution, also called Gaussian distribution (eq. 8), is a probability density function for a random variable $x$ expressed as follows:

$$f(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2} \quad (8)$$

The parameter $\mu$ is the mean of the distribution and the parameter $\sigma$ is the standard deviation that quantifies the amplitude of the spread of the distribution. One of the remarkable
233

features of a normal distribution is the equality between the mean and the median, reflecting
the symmetric shape of the distribution around the mean value.

234

235

As a canonical example to illustrate a normal distribution for modern terrestrial vertebrates, consider the case of the height of males or females from a given geographical area. As shown by A'Hearn et al. (2009), the presence of adolescents in the data base of height distributions for Italian provinces, drawn from conscription records, has been shown to systematically distort the distribution of heights whatever the geographic province considered. If height only resulted from a genetic background, only two possibilities would be observed; short and tall people. In fact, there are many independent genetic, cultural and environmental factors that influence the height of an individual, including uncertainties associated with the measurement of this variable. The combination of several independent factors explains why the data distribution is normal, in agreement with the Central Limit Theorem. Indeed, this theorem states that the average of many observations of a random variable is itself a random variable, of which its distribution tends to a normal distribution with the increasing number of samples.

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

Two additional examples presented below are of interest for our study as they concern
the body temperature and the diet of human populations; two variables that contribute to the
oxygen isotope composition of vertebrate skeletal remains. Firstly, Waalen and Buxbaum
(2011) have shown that men and women body temperatures follow distinct but normal
distributions for individuals whose age is comprised between 20 years and 59 years. In fact,
the study was applied to a little more than 18,000 individuals whose age was comprised
between 20 years and 98 years. Above, 59 years, ageing is accompanied by decreasing body
temperatures for both women and men categories. Secondly, McCarron et al. (2013)
investigated a 24-hour urinary sodium excretion in a human population consisting of more
than 50,000 individuals over at least 5 years and across 45 countries. Based on their survey,
258

they concluded that human sodium intake is characterized by a narrow range of data
following a normal distribution.

259

260

Consequently, we expect that the δ¹⁸O(A)ₜ distribution of a population of vertebrates
that lived in a given area under similar climatic conditions will follow a normal distribution.
Indeed, the main controlling variables, such as air or water temperature, body temperature,
and δ¹⁸O(W)ᵢ of the diet, must follow a normal distribution considered apart from each other.

261

262

263

264

265

## 2.3. How to perform simulations of diagenetic apatite δ¹⁸O(A)ₜ?

266

267

268

269

270

271

272

273

274

275

276

277

278

279

The simulation tool ApaOxIMOD (https://github.com/jpflandrs/ApaOxIRA) is
basically a Monte Carlo process simulation, based on the creation of a huge array of A/W
values taken from a uniform distribution between given threshold values. A moving A/W
window may be set within the threshold values to study the influence of modifications of the
A/W range. This represents the diversity of the deposit situation and physical state of the
biological apatite. Typically, 2x10⁶ situations are simulated. For each situation, equation (4) is
used to compute the δ¹⁸O in biological apatite δ¹⁸O(A)ₜ at the equilibrium on the whole array
given the temperature T, initial δ¹⁸O of biological apatite (δ¹⁸O(A)ᵢ) and the δ¹⁸O of reacting
water δ¹⁸O(W)ᵢ. The three parameters (T, δ¹⁸O(W)ᵢ, δ¹⁸O(A)ᵢ) may be set and their variations
may be explored by steps for each simulated A/W range to analyze the relative influence of
each parameter. The output is a family of histograms summarizing the possible fate of
δ¹⁸O(A)ᵢ in the population after the diagenetic episode, either related to early or late
diagenesis as the model is dimensionless.
The risk-analysis tool ApaOxIRA uses a similar approach except that the initial array of W/A situations is composed of subarrays of the size n representing the desirable sampling strategy. The three parameters (T, δ¹⁸O(W)ᵢ, δ¹⁸O(A)ᵢ) are fixed but n may vary on a defined step basis. The analytical uncertainty may also be considered. Finally, each subarray representing a sample set of length n of δ¹⁸O(A)ᵢf values is submitted to the normality tests. Shapiro-Wilk and Anderson-Darling are combined. Both have demonstrated a high potential in normality detection (Razali and Wah, 2011; Yap and Sim, 2011) but they are using different approaches and have different sensitivity to skew and tails. Normality is rejected at the given risk α when one or both tests does reject the normality. The number of sets in the array where normality cannot be rejected is computed. The ratio of Normal to total simulated samples gives the "Normal distribution by chance" ratio. The standard deviation of this ratio is estimated by a bootstrap process. The output is a table of the "Normal distribution by chance" ratio for each value of n. When real data are submitted, aside from basic statistics and skew and kurtosis tests, both Shapiro-Wilk and Anderson-Darling tests are computed, their results are combined as previously to decide on the normality of the submitted set. A set of graphs helps the analysis, including the fitting to a normal law and the qq-plot. ApaOxIRA is available at GitHub (https://github.com/jpflandrs/ApaOxIRA).

297

## 2.4. *Structure, functions and properties of the Python™ program*

299

Two Python3 programs ApaOxIMOD and ApaOxIRA have been designed to perform all the tasks. Both are entirely relying on the array structure of NumPy to take advantage of the speed and efficiency of NumPy in manipulating large arrays. The statistics are also computed with NumPy function and python3 is minimally used to connect Numpy functions and for reporting. The 15 possible parameters and options are taken from a pilot file using the
YAML format that requires only a plain-text editor to be quickly changed. For ApaOxIRA, the computation is speeded-up by multithreading. and the speed depends on the parameters, especially the number of sets and the increasing complexity of the MC simulation. Typically, using an iMacPro (2017) with 3 GHz Intel Xeon W by using 18 threads and less than 256 Mo of RAM, when asking for a MC with 250,000 sets, the mean individual time to result is 2.3 seconds when testing repetitively sets of length 33-75. In the range of array size 100,000-1,000,000, completion occurs in near linear time.

2.5. Simulations of diagenetic apatite δ¹⁸O(A)ₜ values obtained from the Python program

In order to clarify expectations, the effects of both variations in T and δ¹⁸O(W)ᵢ on the resulting diagenetic bioapatites were investigated in distinct W/A ranges comprised between 0.05 and 0.95 with the choice of range widths and their boundary values being chosen arbitrary in this case (Figures 3 and 4). At any T, δ¹⁸O(A)ₜ decreases with decreasing δ¹⁸O(W)ᵢ. We can also observe that with increasing temperature of reaction in the range 10°C to 50°C, the range of δ¹⁸O(A)ₜ tends to increase. In addition, for a given temperature T, the normality departure of the distribution increases with the decreasing δ¹⁸O(W)ᵢ (Figures 3 and 4). This normality departure can be quantified by using the skewness parameter (eq. 9) that is calculated by using the Fisher-Pearson standardized moment coefficient (Ziegel, 2001):

$skewness = \frac{\sum_{i}^{n}(X_i - \mu)^3}{(n-1)\sigma^3} \qquad (9)$
where n is the number of variables of the distribution and Xᵢ equals the iᵗʰ random variable.

Moreover, the shape of either tail of distribution may be described by the kurtosis parameter (eq. 10) that is expressed as follows (Ziegel, 2001):

$$kurtosis = n \left( \frac{\sum_{i}^{n} (X_i - \mu)^4}{\sum_{i}^{n} (X_i - \mu^2)^2} \right) \quad (10)$$

It is worthy to note that δ¹⁸O(A)ₜ values are characterized by increasing ranges when

W/A ranges approach 1. However, those δ¹⁸O(A)ₜ range widths become narrower with increasing T because the isotopic fractionation between the mineral and water phases decreases with increasing T as predicted by thermodynamics (Figure 3).

## 2.6. From normal distributions to non-gaussian distributions of vertebrate δ¹⁸O(A)ₜ values

According to specific oxygen isotope fractionation equations that were determined for fish (Kolodny et al., 1983; Lécuyer et al., 2013), reptiles (Barrick et al., 1999; Amiot et al., 2007), marine mammals (Séon et al., 2022), herbivorous (Iacumin and Longinelli, 2002; Amiot et al., 2004), omnivorous (Longinelli, 1984) and carnivorous mammals (Iacumin and Longinelli, 2002), only animal δ¹⁸O(A) values that may be treated by isotopic fractionation equations of similar properties can be pooled together to be statistically tested for normal distributions. In the case of a sample collection that suffered diagenetic alteration, the oxygen isotope compositions are expected to follow non-gaussian distributions as a result of temperature-dependent equation (3) described above.
350 3. Results: model outputs

351 3.1. *Examples of normal distributions in the modern and fossil isotopic records*

352

353 Please note that Q-Q plots and frequency histograms are provided for each case study
354 in the Supplementary material (Figures 1 to 4). Here we present four examples of normal
355 distributions observed for modern and fossil fauna. Firstly, bioapatites from Neogene sharks
356 and rays from the Brazilian "Blue Amazon" (Aguilera et al., 2017) display a normal
357 distribution with a *p*-value of 0.68 > α = 0.05 (H₀ = normal distribution is assumed; n = 62)
358 and a Shapiro-Wilk test yielding W = 0.9834 comprised in the 95% region of acceptance
359 [0.9555 – 1]. Secondly, deeper in time, a sample collection of Jurassic sharks retrieved from
360 Swiss marine sediments (Leuzinger et al., 2015) also display a normal distribution as attested
361 by a *p*-value of 0.43 > α = 0.05 (H₀ = normal distribution is assumed; n = 50) and a Shapiro-
362 Wilk test yielding W = 0.9768 comprised in the 95% region of acceptance [0.954 – 1].
363 Thirdly, normal distributions also apply to terrestrial vertebrates as shown by the oxygen
364 isotope study of bioapatite remains from a faunal assemblage dated to the Last Glacial
365 Maximum (≈ 20,000 years old) of southwestern France and northern Spain (Lécuyer et al.,
366 2021). The normal distribution observed for Cervidae, Caprinae and Equidae is described by a
367 *p*-value of 0.43 > α = 0.05 (H₀ = normal distribution is assumed; n = 63) and a Shapiro-Wilk
368 test yielding W = 0.9809 comprised in the 95% region of acceptance [0.9621 – 1].

369

370 3.2. *Example of non-gaussian distribution in the fossil isotopic record: witness of post-
depositional alteration processes*

372

373 Both Q-Q plot and frequency histogram are provided for the case study presented below
374 (see the Supplementary material, Figure 4). Late Cretaceous (≈ 70 Ma) fish and reptiles from
Laño, Spain (Lécuyer et al., 2003) are characterized by a low *p*-value of 0.036 (n = 27) rejecting the hypothesis H₀ that observations follow a normal distribution. The density distribution of δ¹⁸O(A)ₜ values is highly skewed (-0.758) and not drawn from a normal population as shown by the Shapiro-Wilk test (W = 0.9182 below the 95% region of acceptance [0.9242 – 1]).

380

## 4. Discussion

### 4.1. Fossilization of vertebrate skeletons and post-burial diagenesis

383

Fossilization is an exceptional process that requires specific conditions such as a rapid burial of the dead organisms in sediments that isolates the skeleton from the external oxidizing fluids. Biogenic skeletal remains may have been altered or transformed in response to the circulation of aqueous fluids of meteoric origin, having temperature, pH and Eh resulting from water-rock interactions, which make them able to dissolve the fossil biominerals (Nielsen-Marsh and Hedges, 2000; Berna et al., 2004; Pfretzschner, 2004; Trueman et al., 2004; Kohn, 2008). This process of diagenesis is usually considered to occur in two main stages, termed ‘early diagenesis’ and ‘late diagenesis’ (Pfretzschner, 2004; see also Keenan, 2016 for a review of apatite diagenesis). Early diagenesis is characterized by inorganic and/or microbially-mediated degradation of soft tissues, partial to complete loss of collagen, and changes in the chemical and structural composition of bioapatite (Trueman and Martill, 2002; Trueman and Tuross, 2002; Pfretzschner, 2006; Jans, 2008; White and Booth, 2014; Turner-Walker, 2019; Procopio et al., 2021). Late diagenesis is characterised by the filling of tooth or bone porosity by authigenic minerals (Pfretzschner, 2001a,b; Pfretzschner and Tütken, 2011).

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398
The recrystallized mineral may acquire a new isotopic composition if the interacting aqueous fluids have a temperature and a δ¹⁸O value distinct from those that originally formed the biominerals during animal's lives (Allan and Matthews, 1990). Minerals are more or less sensitive to dissolution processes depending on their solubility products. For example, at 25°C, the constant of solubility of calcite equals 3.4x10⁻⁹ while that of apatite is only 2.10x10⁻³³ (Chase, 1998) according to JANAF thermodynamic tables, making the latter more resistant to dissolution than calcite. Vertebrates are characterized by hard tissues of various robustness. Tooth enamel is more resistant to dissolution processes and may be preserved for tens and even hundreds of million years (Kolodny and Raab, 1988; Kohn et al., 1999; Passey and Cerling, 2002; Zazzo et al., 2004a). Consequently, tooth enamel has been privileged to measure δ¹⁸O phosphate values in order to retrieve parameters of the living environment of vertebrates. In the case of small mammals such as rodents, it is difficult to mechanically separate tooth enamel from dentin (Royer et al., 2013); the latter is characterized by greater porosity than enamel. When bones are weak and often deeply altered in their host sediments, many vertebrates provide only more or less mineralised bone remains, such as birds and turtles. However, they have been the subject of stable isotope studies that have provided invaluable information on the palaeoenvironments (Barrick, 1998; Matson and Fox, 2008; Angst et al., 2013; 2014; Suarez et al., 2021).

417

418

## 4.2. *The Laño case study: estimation of T and δ¹⁸O(W) during diagenesis*

419

420

421

422

423

The examination of supplementary Figure 4 suggests that pristine δ¹⁸O(A)ᵢ values could have been comprised between 20‰ and 21‰ while the lowest value of 16.5‰ should reflect the largest amplitude of alteration (W/A) suffered by this sample. The surficial conditions in which those vertebrate fossils remained buried suggest T values that did not
exceed 20°C and were not below 10°C regarding the latitude of the sedimentary deposit (Lécuyer et al., 2003). The shape of the density distribution and the lowest δ¹⁸O(A)ₜ values suggest small W/A and negative δ¹⁸O(W)ᵢ values that are compatible with meteoric waters as the main agent (possibly along with microorganisms of the soil) of alteration of bioapatites. If meteoric waters (MW) are considered to be the reacting aqueous fluid, i.e. δ¹⁸O(MW)ᵢ, in deposits located under mid- to high-latitudes, T and δ¹⁸O(MW)ᵢ are dependent variables according to equation (11) expressed as follows:

$$MAAT (°C) = 1.40(±0.08)δ^{18}O(MW)_{i} + 21.86(±0.64) \quad (11)$$

where MAAT represents the mean annual air temperature and δ¹⁸O(MW)ᵢ, the weighted mean isotopic composition of annual precipitation (Lécuyer et al., 2021).

By combining equations (7) and (11), both variables can be solved. Such a technique has already been proposed by Prud'homme et al. (2016) for the reconstruction of climatic conditions during the Last Glacial Maximum in Western Europe. In the present analysis of the Laño isotope data, a δ¹⁸O(MW)ᵢ range from -8.5 to -1.3 corresponds to a T range comprised between 10°C and 20°C according to equation (10). Consequently, for a δ¹⁸O(A)ᵢ of 20.5‰, a δ¹⁸O(A)ₜ value of 16.5‰ may have resulted from W/A values in the range 0.3 to 9.3. However, the shape of the observed distribution is not compatible with W/A ≥ 1, thus the T range is restricted to 10°C–16°C for a δ¹⁸O(W)ᵢ range of -8.5 to -4.2 and W/A from 0.3 to 0.8.
446

4.3. How can the Monte Carlo approach be used to interpret the δ¹⁸O(A)ₜ of fossil biogenic apatite phosphate?

447

448

449

We are aware that our model presents several limits to detect the alteration of the pristine oxygen isotope composition of bioapatite. However, following some rules presented below should allow researchers to detect whether or not their bioapatite samples suffered post-depositional oxygen isotope exchange with aqueous solutions:

450

451

452

453

- The risk of incorrectly assuming a normal distribution is a function of the number of samples of skeletal phosphate remains analysed. For example, if sixty bioapatite samples are analysed, a normality test can be performed with a 5% risk (Figure 5).

454

455

456

- The normality test also requires that the phosphate δ¹⁸O(A) of the sample collection can be analysed by similar isotopic fractionation equations that can be pooled, such as different species of herbivorous mammals. Other equations must be used when applied to phosphatic skeletal remains of carnivorous mammals, reptiles, fish or birds.

457

458

459

460

- The fossil vertebrate population must have lived in an area where environmental conditions (temperature, oxygen isotope composition of drinking water and food water) remained fairly stable during the time span corresponding to the sedimentary deposit. The outcrop scale is recommended for paleoenvironments such as those represented by Lagerstätten sedimentary deposits. See for example, the Late Jurassic deposits of Western Europe (Solnhofen, Crayssac, Canjuers, Cerin, Chassiron), which contain numerous marine and terrestrial vertebrates (Billon-Bruyat et al., 2005).

461

462

463

464

465

466

467

- In the case of a sample collection that has undergone diagenetic alteration, the oxygen isotope compositions of phosphatic skeletal remains can be interpreted using a unique oxygen isotope phosphate-water temperature-dependent equation (Kolodny et al., 1983; Lécuyer et al., 2013)

468

469

470
471 - Note that a non-normal distribution of δ¹⁸O(A)ₜ values is only expected to be observed when W/A ratios remained variable and lower than 1. Above the threshold value of 1, the δ¹⁸O(A)ₜ reflects the δ¹⁸O(W)ᵢ at a given T that fixes the magnitude of the isotope fractionation. So, the distribution of data will result in a cluster that may resemble a normal distribution because of the stochastic nature of analytical errors generated by the Isotope Ratio Mass Spectrometers (IRMS).

477

478 **4. Conclusions**

479 On the basis of some theoretical aspects as well as modern and past examples of 480 faunal associations, we show that the distribution of δ¹⁸O(A) within a sample collection 481 follows a normal distribution whose density distribution results from the combination of 482 independent factors such as those related to the genetics, environment and analytical 483 measurements. In the case of geochemical isotopic interactions between various amounts of 484 circulating aqueous fluid and the bioapatite (W/A molar ratio for oxygen) in the Earth's 485 surface layers (soils, sediments), our numerical simulations show that δ¹⁸O(A) distributions 486 become strongly asymmetric, with shapes that depend on the reaction temperature and the 487 δ¹⁸O(W) of the reacting aqueous fluid.

488 For water-apatite mass ratios comprised between 0.05 and 0.95, the bioapatite δ¹⁸O 489 values decrease with decreasing δ¹⁸O values of reacting water at any temperature, while with 490 increasing temperature of reaction in the range 10°C to 35°C, the range of bioapatite δ¹⁸O 491 tends to slightly increase. For a given temperature T, the normality departure of the 492 distribution increases with the decreasing δ¹⁸O value of reacting water; this normality
493 departure and the shape of either tail of distribution being estimated by the skewness and
494 kurtosis parameters, respectively.

495 The statistical approach developed in this study requires some conditions to be applied
496 in the strategy for sampling fossils in the field. Those conditions imply, firstly, that ideally at
497 least sixty bioapatite samples should be collected in order to perform a reliable normality test.
498 Secondly, the skeletal remains must belong to animals that lived within a short time interval
499 and in an environment characterised by similar air or water temperatures and a steady oxygen
500 isotope composition of the surrounding water, whether marine or fresh. Our advice is
501 therefore to work at the outcrop scale and ideally in Lagerstätten sedimentary deposits. Third
502 and finally, the normality test also requires that the phosphate δ¹⁸O(A) of the sample
503 collection can be treated by similar isotopic fractionation equations that can be pooled, such
504 as various species of herbivorous mammals. In the case of a sample assemblage that has
505 undergone diagenetic alteration, the oxygen isotope compositions of phosphate skeletal
506 remains can be interpreted using a unique oxygen isotope phosphate-water temperature-
507 dependent equation.

508 Interestingly, the asymmetric distribution of δ¹⁸O(A) offers the possibility of
509 estimating some variables of the depositional environment of skeletal remains as exemplified
510 by the Laño case study. Indeed, knowing the ambient temperature as well as the oxygen
511 isotope composition of the diagenetic aqueous phase (δ¹⁸O(W)) open a “frosted glass
512 window” on the climatic conditions that prevailed during the diagenetic alteration of the
513 bioapatites.

514

515 Authors contributions

516 C.L. and J.-P.F. designed the study and wrote the paper; both authors reviewed the
517 manuscript.
518

519 **Declaration of Competing Interest**

520 The authors declare that they have no conflict interest.

521

522 **Data availability**

523 Data used to test normality distributions and Python™ code are available as supplementary
524 files.

525

526 **Acknowledgements**

527 The authors are grateful to R. Amiot and Th. Clauzel for their careful reading of a draft of this
528 manuscript that contributed to improve the scientific content of this study. This scientific
529 content of this study has also been largely improved by the comments and suggestions
530 provided by two anonymous reviewers.

531
533 **References**

534 A'Hearn, B., Peracchi, F., Vecchi, G., 2009. Height and the normal distribution: Evidence from Italian military data. Demography 46, 1-25. https://doi.org/10.1353/dem.0.0049

535 Aguilera, O., Luz, Z., Carrillo-Briceño, J.D., Kocsis, L., Vennemann, T.W., de Toledo, P.M., 536 Nogueira, A., Amorim, K.B., Moraes-Santos, H., Polck, M.R., de Lourdes Ruivo, M., 537 Linhares, A.P., Monteiro-Neto, C., 2017. Neogene sharks and rays from the Brazilian 538 'Blue Amazon'. PLoS One 12, e0182740. https://doi.org/10.1371/journal.pone.0182740

540 Allan, J.R., Matthews, R.K., 1990. Isotope signatures associated with early meteoric diagenesis, in: Tucker, M.E., Barthurst, R.G.C. (Eds.), Carbonate Diagenesis. John Wiley 541 & Sons, pp. 197-217. https://doi.org/10.1002/9781444304510.ch16

543 Amiot, R., Buffetaut, E., Lécuyer, C., Wang, X., Boudad, L., Ding, Z., Fourel, F., Hutt, S., 544 Martineau, F., Medeiros, M.A., Mo, J., Simon, L., Sutheethorn, V., Sweetman, S., Tong, 545 H., Zhang, F., Zhou, Z., 2010. Oxygen isotope evidence for semi-aquatic habits among 546 spinosaurid theropods. Geology 38, 139-142.

547 Amiot, R., Lécuyer, C., Buffetaut, E., Escarguel, G., Fluteau, F., Martineau, F., 2006. Oxygen isotopes from biogenic apatites suggest widespread endothermy in Cretaceous dinosaurs. Earth and Planetary Science Letters 246, 41-54.

550 Amiot, R., Lécuyer, C., Buffetaut, E., Fluteau, F., Legendre, S., Martineau, F., 2004. 551 Latitudinal temperature gradient during the Cretaceous Upper Campanian-Middle 552 Maastrichtian: δ¹⁸O record of continental vertebrates. Earth Planet. Sci. Lett. 226, 255- 553 272. https://doi.org/10.1016/j.epsl.2004.07.015

554 Amiot, R., Lécuyer, C., Escarguel, G., Billon-Bruyat, J.-P., Buffetaut, E., Langlois, C., 555 Martin, S., Martineau, F., Mazin, J.-M., 2007. Oxygen isotope fractionation between 556 crocodilian phosphate and water. Palaeogeogr. Palaeoclimatol. Palaeoecol. 243, 412-420. 557 https://doi.org/10.1016/j.palaeo.2006.08.013.

558 Angst, D., Buffetaut, E., Lécuyer, C., Amiot, R., 2013. "Terror Birds" (Phorusrhacidae) from 559 the Eocene of Europe imply trans-Tethys dispersal. PLoS One, 8, e80357. 560 https://doi.org/10.1371/journal.pone.0080357

561 Angst, D., Lécuyer, C., Amiot, R., Buffetaut, E., Fourel, F., Martineau, F., Legendre, S., 562 Abourachid, A., Herrel, A., 2014. Isotopic and anatomical evidence of an herbivorous diet 563 in the Early Tertiary giant bird *Gastornis*. Implications for the structure of Paleocene 564 terrestrial ecosystems. Naturwissenschaften 101, 313-322.
565 Aufort, J., Gervais, C., Segalen, L., Labourdette, N., Coelho-Diogo, C., Baptiste, B., Beyssac, 566 O., Amiot, R., Lécuyer, C., Balan, E. (2019). Atomic scale transformation of bone in 567 controlled aqueous alteration experiments. Palaeogeography, Palaeoclimatology, 568 Palaeoecology 526, 80–95.

569 Barrick, R.E., 1998. Isotope paleobiology of the vertebrates: ecology, physiology, and 570 diagenesis. The Paleontological Society Papers 4, 101–137.

571 Barrick, R.E., Fischer, A.G., Showers, W.J., 1999. Oxygen isotopes from turtle bone: 572 applications for terrestrial paleoclimates? Palaios 14, 186–191.

573 Barrick, R.E., Showers, W.J., 1994. Thermophysiology of *Tyrannosaurus rex*: evidence from 575 oxygen isotopes. Science 265, 222–224. DOI: 10.1126/science.265.5169.222

576 Berna, F., Matthews, A., Weiner, S., 2004. Solubilities of bone mineral from archaeological 577 sites: the recrystallization window. Journal of Archaeological Science 31, 867–882.

578 Bernard, A., Lécuyer, C., Vincent, P., Amiot, R., Bardet, N., Buffetaut, E., Cuny, G., Fourel, 579 F., Martineau, F., Mazin, J.-M., Prieur, A., 2010. Regulation of body temperature by some 580 Mesozoic marine reptiles. Science 328, 1379–1382. DOI: 10.1126/science.1187443

581 Billon-Bruyat, J.P., Lécuyer, C., Martineau, F., & Mazin, J.M., 2005. Oxygen isotope 582 compositions of Late Jurassic vertebrate remains from lithographic limestones of western 583 Europe: implications for the ecology of fish, turtles, and crocodilians. Palaeogeography 584 Palaeoclimatology Palaeoecology 216, 359–375.

585 Blake, R.E., O'Neil, J.R., Garcia, G.A., 1997. Oxygen isotope systematics of biologically 586 mediated reactions of phosphate: I. Microbial degradation of organophosphorus 587 compounds. Geochim. Cosmochim. Acta 6, 4411–4422. https://doi.org/10.1016/S0016- 588 7037(97)00272-X

589 Blake, R.E., O'Neil, J.R., Garcia, G.A., 1998. Effects of microbial activity on the δ¹⁸O of 590 dissolved inorganic phosphate and textural features of synthetic apatites. American 591 Mineralogist 83, 1516–1531.

592 Blake, R.E., O'Neil, J.R., Surkov, A.V., 2005. Biogeochemical cycling of phosphorus: 593 insights from oxygen isotope effects of phosphoenzymes. American Journal of 594 Science 305, 596–620.

595 Chang, S.J., Blake, R.E., 2015. Precise calibration of equilibrium oxygen isotope 596 fractionations between dissolved phosphate and water from 3 to 37°C. Geochimica et 597 Cosmochimica Acta 150, 314–329.
598 Clauzel, T., Richardin, P., Ricard, J., Prilaux, G., Martin, J. E., Amiot, R., Fourel, F., Nguyen
599 Anh, V., Arnaud-Godet, F., Vinçon-Laugier, A., Flandrois J.-P., Lécuyer, C., 2023.
600 Geographic origin and social status of the Gallic warriors from Ribemont-sur-Ancre
601 (France) studied through isotope systematics of bone remains. International Journal of
602 Osteoarchaeology 33, 39–50.

603 Criss, R.E., Taylor, H.P., 1986. Meteoric-hydrothermal systems, in: Valley, J.W., Taylor, H.P., Jr., O'Neil, J.R. (Eds.), Stable isotopes in high temperature geological processes, vol. 16. Min. Soc. Amer., Washington, DC, pp 373–424.

606 Dal Sasso, G., Asscher, Y., Angelini, I., Nodari, L., Artioli, G., 2018. A universal curve of
607 apatite crystallinity for the assessment of bone integrity and preservation. Scientific
608 Reports 8, 1–13.

609 Evans, J., Stoodley, N., Chenery, C., 2006. A strontium and oxygen isotope assessment of a
610 possible fourth century immigrant population in a Hampshire cemetery, southern England.
611 J. Archaeol. Sci. 33, 265–272. https://doi.org/10.1016/j.jas.2005.07.011

612 Faure-Brac, M.G., Amiot, R., de Muizon, C., Cubo, J., Lécuyer, C., 2022. Combined
613 paleohistological and isotopic inferences of thermometabolism in extinct Neosuchia, using
614 *Goniopholis* and *Dyrosaurus* (Pseudosuchia: Crocodylomorpha) as case
615 studies. Paleobiology 48, 302–323. https://doi.org/10.1017/pab.2021.34.

616 Fricke, H.C., O'Neil, J.R., 1996. Inter-and intra-tooth variation in the oxygen isotope
617 composition of mammalian tooth enamel phosphate: implications for palaeoclimatological
618 and palaeobiological research. Palaeogeogr. Palaeoclimatol. Palaeoecol. 126, 91–99.
619 https://doi.org/10.1016/S0031-0182(96)00072-7.

620 Fricke, H.C., O'Neil, J.R., Lynnerup, N., 1995. SPECIAL REPORT: Oxygen isotope
621 composition of human tooth enamel from medieval Greenland: linking climate and
622 society. Geology 23, 869–872.

623 Fricke, H.C., Rogers, R.R., 2000. Multiple taxon–multiple locality approach to providing
624 oxygen isotope evidence for warm-blooded theropod dinosaurs. Geology 28, 799–802.
625 https://doi.org/10.1130/0091-7613(2000)28<799:MTLATP>2.0.CO;2

626 Geisler, T., Kasioptas, A., Menneken, M., Perdikouri, C., Putnis, A., 2009. A preliminary in
627 situ Raman spectroscopic study of the oxygen isotope exchange kinetics between H₂O and
628 (PO₄)aq. Journal of Geochemical Exploration 1, 37–37.

629 Goedert, J., Amiot, R., Arnaud-Godet, F., Cuny, G., Fourel, F., Hernandez, J. A., Pedreira-
630 Segade, U., Lécuyer, C., 2017. Miocene (Burdigalian) seawater and air temperatures
631 estimated from the geochemistry of fossil remains from the Aquitaine Basin,
632

France. Palaeogeogr.

Palaeoclimatol.

Palaeoecol. 481,

14–28.

633

https://doi.org/10.1016/j.palaeo.2017.04.024

634

Hedges, R.E., 2002. Bone diagenesis: an overview of processes. Archaeometry 44, 319–328.

635

Iacumin, P., Bocherens, H., Mariotti, A., Longinelli, A., 1996. Oxygen isotope analyses of co-existing carbonate and phosphate in biogenic apatite: a way to monitor diagenetic alteration of bone phosphate? Earth and Planetary Science Letters 142, 1–6.

636

637

638

Iacumin, P., Longinelli, A., 2002. Relationship between δ¹⁸O values for skeletal apatite from reindeer and foxes and yearly mean δ¹⁸O values of environmental water. Earth Planet. Sci. Lett. 201, 213–219. https://doi.org/10.1016/S0012-821X(02)00635-0

639

640

641

Jaisi, D.P., Kukkadapu, R.K., Stout, L.M., Varga, T., Blake, R.E., 2011. Biotic and abiotic pathways of phosphorus cycling in minerals and sediments: insights from oxygen isotope ratios in phosphate. Environmental Science & Technology 45, 6254–6261.

642

643

644

Jans, M.M. (2008). Microbial bioerosion of bone – a review, in: Wisshak, M., Tapanila, L. (Eds.) Current Developments in Bioerosion. Erlangen Earth Conference Series. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-77598-0_20

645

646

647

Joachimski, M.M., Lai, X., Shen, S., Jiang, H., Luo, G., Chen, B., Chen, J., Sun, Y., 2012. Climate warming in the latest Permian and the Permian–Triassic mass extinction. Geology 40, 195–198. https://doi.org/10.1130/G32707.1

647

648

649

650

Keenan, S.W., 2016. From bone to fossil: A review of the diagenesis of bioapatite. Am. Mineral. 101, 1943–1951. https://doi.org/10.2138/am-2016-5737

650

651

652

Keenan, S.W., Engel, A.S., 2017. Early diagenesis and recrystallization of bone. Geochimica et Cosmochimica Acta 196, 209–223.

653

654

Kohn, M.J., 2008. Models of diffusion-limited uptake of trace elements in fossils and rates of fossilization. Geochimica et Cosmochimica Acta 72, 3758–3770.

655

656

Kohn, M.J., Schoeninger, M.J., Barker, W.W., 1999. Altered states: effects of diagenesis on fossil tooth chemistry. Geochimica et Cosmochimica Acta 63, 2737–2747.

657

658

Kolodny, Y., Luz, B., 1992. Isotope signatures in phosphate deposits: Formation and diagenetic history, in: Clauer, N., Chaudhuri, S. (Eds.), Isotopic Signatures and Sedimentary Records. Lecture Notes in Earth Sciences, vol 43. Springer, Berlin, Heidelberg. https://doi.org/10.1007/BFb0009862

659

660

661

662

Kolodny, Y., Luz, B., Navon, O., 1983. Oxygen isotope variations in phosphate of biogenic apatites, I. Fish bone apatite—rechecking the rules of the game. Earth Planet. Sci. Lett. 64, 398–404. https://doi.org/10.1016/0012-821X(83)90100-0

663

664
665 Kolodny, Y., Luz, B., Sander, M., Clemens, W.A., 1996. Dinosaur bones: fossils or pseudomorphs? The pitfalls of physiology reconstruction from apatitic 666 fossils. Palaeogeogr. Palaeoclimatol. Palaeoecol. 126, 161–171. 667 668 https://doi.org/10.1016/S0031-0182(96)00112-5

669 Kolodny, Y., Raab, M., 1988. Oxygen isotopes in phosphatic fish remains from Israel: 670 paleothermometry of tropical Cretaceous and Tertiary shelf waters. Palaeogeogr. 671 Palaeoclimatol. Palaeoecol. 64, 59–67. https://doi.org/10.1016/0031-0182(88)90142-3

672 Lécuyer C., Grandjean P., Sheppard, S.M.F., 1999. Oxygen isotope exchange between 673 dissolved phosphate and water at temperatures ≤135°C: Inorganic versus biological 674 fractionations. Geochim. Cosmochim. Acta 63, 855–862. https://doi.org/10.1016/S0016- 675 7037(99)00096-4

676 Lécuyer, C., Amiot, R., Trotter, J., Touzeau, A., 2013. Calibration of the phosphate δ¹⁸O 677 paleothermometer with the calcium carbonate-water oxygen isotope fractionation 678 equations. Chem. Geol. 347, 217–226. https://doi.org/10.1016/j.chemgeo.2013.03.008

679 Lécuyer, C., Bogey, C., Garcia, J.-P., Grandjean, P., Barrat, J.-A., Floquet, M., Bardet, N., 680 Pereda-Superbiola, X., 2003. Stable isotope compositions and rare earth element content of 681 vertebrate remains from the late Cretaceous of northern Spain (Laño): Did the 682 environmental record survive? Palaeogeogr. Palaeoclimatol. Palaeoecol. 193, 457–471. 683 https://doi.org/10.1016/S0031-0182(03)00261-X

684 Lécuyer, C., Grandjean, P., Barrat, J. A., Nolvak, J., Emig, C., Paris, F., Robardet, M., 1998. 685 δ¹⁸O and REE contents of phosphatic brachiopods: a comparison between modern and 686 lower Paleozoic populations. Geochimica et Cosmochimica Acta 62, 2429–2436.

687 Lécuyer, C., Grandjean, P., Emig, C.C., 1996. Determination of oxygen isotope fractionation 688 between water and phosphate from living lingulids: potential application to 689 palaeoenvironmental studies. Palaeogeogr. Palaeoclimatol. Palaeoecol. 126, 101–108. 690 https://doi.org/10.1016/S0031-0182(96)00073-9

691 Lécuyer, C., Hillaire-Marcel, C., Burke, A., Julien, M. A., Hélie, J.F., 2021. Temperature and 692 precipitation regime in LGM human refugia of southwestern Europe inferred from δ¹³C 693 and δ¹⁸O of large mammal remains. Quat. Sci. Rev. 255, 106796. 694 https://doi.org/10.1016/j.quascirev.2021.106796

695 Leuzinger, L., Kocsis, L., Billon-Bruyat, J. P., Spezzaferri, S., Vennemann, T., 2015. Stable 696 isotope study of a new chondrichthyan fauna (Kimmeridgian, Porrentruy, Swiss Jura): an 697 unusual freshwater-influenced isotopic composition for the hybodont shark 698 *Asteracanthus*. Biogeosciences 12, 6945–6954. https://doi.org/10.5194/bg-12-6945-2015
699 Lightfoot, E., O'Connell, T.C., 2016. On the use of biomineral oxygen isotope data to identify human migrants in the archaeological record: intra-sample variation, statistical methods and geographical considerations. *PloS one*, 11(4), e0153850.
700
701
702 https://doi.org/10.1371/journal.pone.0153850.

703 Longinelli, A., 1966. Ratios of oxygen-18: oxygen-16 in phosphate and carbonate from living and fossil marine organisms. *Nature* 211, 923–927. https://doi.org/10.1038/211923a0

705 Longinelli, A., 1984. Oxygen isotopes in mammal bone phosphate: a new tool for paleohydrological and paleoclimatological research? *Geochim. Cosmochim. Acta* 48, 385–390. https://doi.org/10.1016/0016-7037(84)90259-X

708 Longinelli, A., Nuti, S., 1973a. Oxygen isotope measurements of phosphate from fish teeth and bones. *Earth and Planetary Science Letters* 20, 337–340.

710 Longinelli, A., Nuti, S., 1973b. Revised phosphate-water isotopic temperature scale. *Earth Planet. Sci. Lett.* 19, 373–376. https://doi.org/10.1016/0012-821X(73)90088-5

712 Matson, S.D., Fox, D.L., 2008. Can oxygen isotopes from turtle bone be used to reconstruct paleoclimates. *Palaios* 23, 24–34.

714 McArthur, J.M., Coleman, M.L., Bremner, J.M., 1980. Carbon and oxygen isotopic composition of structural carbonate in sedimentary francolite. *J. Geol. Soc. 137*, 669–673.
715
716 https://doi.org/10.1144/gsjgs.137.6.0669

717 McArthur, J.M., Herczeg, A., 1990. Diagenetic stability of the isotopic composition of phosphate-oxygen: palaeoenvironmental implications. Geological Society, London, Special Publications 52, 119–124.

720 McCarron, D.A., Kazaks, A.G., Geerling, J.C., Stern, J.S., Graudal, N.A., 2013. Normal range of human dietary sodium intake: a perspective based on 24-hour urinary sodium excretion worldwide. *Am. J. Hypertens.* 26, 1218–1223. https://doi.org/10.1093/ajh/hpt139

723 Nielsen-Marsh, C.M., Hedges, R.E., 2000. Patterns of diagenesis in bone I: the effects of site environments. *Journal of Archaeological Science* 27, 1139–1150.

725 O'Neil, J.R., Vennemann, T.W., McKenzie, W.F., 2003. Effects of speciation on equilibrium fractionations and rates of oxygen isotope exchange between (PO₄)aq and H₂O. *Geochimica et Cosmochimica Acta* 67, 3135–3144.

728 Passey, B.H., Cerling, T.E., 2002. Tooth enamel mineralization in ungulates: implications for recovering a primary isotopic time-series. *Geochimica et Cosmochimica Acta* 66, 3225–3234.

731 Paytan, A., Kolodny, Y., Neori, A., Luz, B., 2002. Rapid biologically mediated oxygen isotope exchange between water and phosphate. *Global Biogeochemical Cycles* 16, 13–1.
733 Pellegrini, M., Pouncett, J., Jay, M., Pearson, M.P., Richards, M.P., 2016. Tooth enamel oxygen "isoscapes" show a high degree of human mobility in prehistoric Britain. Sci.
734
735 Rep. 6, 1-9. DOI: 10.1038/srep34986

736 Pfretzschner, H.U., 2001a. Iron oxides in fossil bone. Neues Jahrbuch für Geologie und Paläontologie-Abhandlungen, 417–429. DOI: 10.1127/njgpa/220/2001/417

738 Pfretzschner, H.U., 2001b. Pyrite in fossil bone. Neues Jahrbuch für Geologie und Paläontologie-Abhandlungen, 1–23. DOI: 10.1127/njgpa/220/2001/1

740 Pfretzschner, H.U., 2004. Fossilization of Haversian bone in aquatic environments. Comptes Rendus Palevol. 3, 605–616.

742 Pfretzschner, H.U., 2006. Collagen gelatinization: the key to understand early bone-diagenesis. Palaeontographica Abteilung A, 135-148. DOI: 10.1127/pala/278/2006/135.

744 Pfretzschner, H.U., Tütken, T., 2011. Rolling bones–Taphonomy of Jurassic dinosaur bones inferred from diagenetic microcracks and mineral infillings. Palaeogeography
745
746 Palaeoclimatology Palaeoecology 310, 117–123.

747 Procopio, N., Mein, C. A., Starace, S., Bonicelli, A., Williams, A., 2021. Bone diagenesis in short timescales: insights from an exploratory proteomic analysis. Biology 10, 460. https://doi.org/10.3390/biology10060460.

750 Prud'Homme, C., Lécuyer, C., Antoine, P., Moine, O., Hatté, C., Fourel, F., Martineau, F., Rousseau, D.D., 2016. Palaeotemperature reconstruction during the Last Glacial from δ¹⁸O of earthworm calcite granules from Nussloch loess sequence, Germany. Earth and Planetary Science Letters 442, 13–20.

754 Pucéat, E., Reynard, B., Lécuyer, C., 2004. Can crystallinity be used to determine the degree of chemical alteration of biogenic apatites? Chemical Geology 205, 83–97.

756 Razali, N.M., Wah, Y.B., 2011. Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests. J. Stat. Model. Anal. 2, 21–33.

758 Rey, K., Amiot, R., Fourel, F., Abdala, F., Fluteau, F., Jalil, N. E., Liu, J., Rubidge, B.S., Smith, R.M.H., Steyer, J.S., Viglietti, P.A., Wang, X., Lecuyer, C., 2017. Oxygen isotopes suggest elevated thermometabolism within multiple Permo-Triassic therapsid clades. *Elife*, 6, e28589.

762 Rey, K., Day, M. O., Amiot, R., Goedert, J., Lécuyer, C., Sealy, J., Rubidge, B.S., 2018. Stable isotope record implicates aridification without warming during the late Capitanian mass extinction. Gondwana Res. 59, 1–8. https://doi.org/10.1016/j.gr.2018.02.017.
765 Reynard, B., Lécuyer, C., Grandjean, P., 1999. Crystal-chemical controls on rare-earth element concentrations in fossil biogenic apatites and implications for paleoenvironmental reconstructions. *Chemical Geology* 155, 233–241.

768 Royer, A., Lécuyer, C., Montuire, S., Amiot, R., Legendre, S., Cuenca-Bescós, G., Jeannet, M., Martineau, F., 2013. What does the oxygen isotope composition of rodent teeth record? *Earth Planet. Sci. Lett.* 361, 258–271. https://doi.org/10.1016/j.epsl.2012.09.058.

771 Sànchez Chillón, B., Alberdi, M.T., Leone, G., Bonadonna, F.P., Stenni, B., Longinelli, A., 1994. Oxygen isotopic composition of fossil equid tooth and bone phosphate: an archive of difficult interpretation. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* 107, 317–328. https://doi.org/10.1016/0031-0182(94)90103-1

775 Séon, N., Amiot, R., Suan, G., Lécuyer, C., Fourel, F., Demaret, F., Vinçon-Laugier, A., Charbonnier, S., Vincent, P. (2022). Intra-skeletal variability in phosphate oxygen isotope composition reveals regional heterothermies in marine vertebrates. *Biogeosciences* 19, 2671–2681.

779 Sharp, Z.D., Atudorei, V., Furrer, H., 2000. The effect of diagenesis on oxygen isotope ratios of biogenic phosphates. *Am. J. Sci.* 300, 222–237. https://doi.org/10.2475/ajs.300.3.222

781 Shemesh, A., 1990. Crystallinity and diagenesis of sedimentary apatites. *Geochimica et Cosmochimica Acta* 54, 2433–2438.

783 Shemesh, A., 1990. Crystallinity and diagenesis of sedimentary apatites. *Geochim. Cosmochim. Acta* 54, 2433–2438. https://doi.org/10.1016/0016-7037(90)90230-I

785 Simon, L., Lécuyer, C., Putelat, T., 2017. The calculation of water-rock ratios using trace element (Li, B) stable isotopes. *Ann. Univ. Mariae Curie-Skłodowska, sectio AAA – Physica* 71, 79. http://dx.doi.org/10.17951/aaa.2016.71.79

788 Sperber (von), C., Kries, H., Tamburini, F., Bernasconi, S.M., Frossard, E., 2014. The effect of phosphomonoesterases on the oxygen isotope composition of phosphate. *Geochimica et Cosmochimica Acta* 125, 519–527.

791 Suarez, C.A., Frucci, M.N., Tompkins, T.B., Suarez, M.B., 2021. Quantification of a North American greenhouse hydrological cycle: using oxygen isotopic composition of phosphate from Early Cretaceous (Aptian–Albian) turtles. *Geological Society, London, Special Publications* 507, 335–345.

795 Surowiec, R.K., Allen, M.R., Wallace, J.M., 2022. Bone hydration: How we can evaluate it, what can it tell us, and is it an effective therapeutic target? *Bone Reports* 16, 101161. https://doi.org/10.1016/j.bonr.2021.101161
798 Taylor, H.P., 1977. Water/rock interactions and the origin of H₂O in granitic batholiths:

799 Thirtieth William Smith lecture. J. Geol. Soc. 133, 509–558.

800 https://doi.org/10.1144/gsjgs.133.6.0509

801 Taylor, Jr, H.P., 1978. Oxygen and hydrogen isotope studies of plutonic granitic rocks. Earth

802 Planet. Sci. Lett. 38, 177–210. https://doi.org/10.1016/0012-821X(78)90131-0

803 Touzeau, A., Blichert-Toft, J., Amiot, R., Fourel, F., Martineau, F., Cockitt, J., Hall, K.,

804 Flandrois, J.-P., Lécuyer, C., 2013. Egyptian mummies record increasing aridity in the Nile

805 valley from 5500 to 1500 yr before present. Earth and Planetary Science Letters 375, 92–

806 100.

807 Trueman, C.N., Behrensmeyer, A.K., Tuross, N., Weiner, S., 2004. Mineralogical and

808 compositional changes in bones exposed on soil surfaces in Amboseli National Park,

809 Kenya: diagenetic mechanisms and the role of sediment pore fluids. Journal of

810 Archaeological Science 31, 721–739.

811 Trueman, C.N., Kocsis, L., Palmer, M. R., Dewdney, C., 2011. Fractionation of rare earth

812 elements within bone mineral: a natural cation exchange system. Palaeogeography

813 Palaeoclimatology Palaeoecology 310, 124–132.

814 Trueman, C.N., Martill, D.M., 2002. The long-term survival of bone: the role of

815 bioerosion. Archaeometry 44, 371–382.

816 Trueman, C.N., Tuross, N., 2002. Trace elements in recent and fossil bone apatite. Reviews in

817 Mineralogy and Geochemistry 48, 489–521.

818 Turner-Walker, G., 2019. Light at the end of the tunnels? The origins of microbial bioerosion

819 in mineralised collagen. Palaeogeography Palaeoclimatology Palaeoecology 529, 24–38.

820 Turner-Walker, G., Galiacho, A.G., Armentano, N., Hsu, C.Q., 2023. Bacterial bioerosion of

821 bones is a post-skeletonisation phenomenon and appears contingent on soil

822 burial. Quaternary International 660, 75–83.

823 Tütken, T., Furrer, H., Vennemann, T.W., 2007. Stable isotope compositions of mammoth

824 teeth from Niederweningen, Switzerland: implications for the Late Pleistocene climate,

825 environment, and diet. Quaternary International 164, 139–150.

826 Tütken, T., Vennemann, T.W., 2009. Stable isotope ecology of Miocene large mammals from

827 Sandelzhausen, southern Germany. Paläontol. Z. 83, 207–226.

828 https://doi.org/10.1007/s12542-009-0011-y
829

Tütken, T., Vennemann, T.W., Pfretzschner, H.U., 2008. Early diagenesis of bone and tooth apatite in fluvial and marine settings: constraints from combined oxygen isotope, nitrogen and REE analysis. Palaeogeography Palaeoclimatology Palaeoecology 266, 254–268.

830

831

832

Végh, E.I., Márquez-Grant, N., Schulting, R.J., 2022. Physicochemical changes in bone bioapatite during the late postmortem interval pre-and post-burning. Applied Spectroscopy 76, 1080–1099.

833

834

835

Vennemann, T.W., Hegner, E., Cliff, G., Benz, G.W., 2001. Isotopic composition of recent shark teeth as a proxy for environmental conditions. Geochimica et Cosmochimica Acta 65, 1583–1599.

836

837

838

Waalen, J., Buxbaum, J.N., 2011. Is older colder or colder older? The association of age with body temperature in 18,630 individuals. J. Gerontol. – Biomed. Sci. Med. Sci. 66, 487–492. https://doi.org/10.1093/gerona/glr001

839

840

841

Welborn, V.V., 2021. Environment-controlled water adsorption at hydroxyapatite/collagen interfaces. Physical Chemistry Chemical Physics 23, 13789–13796.

842

843

White, C.D., Spence, M.W., Stuart-Williams, HL. Q., Schwarcz, H.P., 1998). Oxygen isotopes and the identification of geographical origins: The Valley of Oaxaca versus the Valley of Mexico. J. Archaeol. Sci. 25, 643–655. https://doi.org/10.1006/jasc.1997.0259.

844

845

846

White, L., Booth, T.J., 2014. The origin of bacteria responsible for bioerosion to the internal bone microstructure: results from experimentally-deposited pig carcasses. Forensic Science International 239, 92–102.

847

848

849

Yap, B.W., Sim, C.H., 2011. Comparisons of various types of normality tests. J. Stat. Comput. Simul. 81, 2141–2155. https://doi.org/10.1080/00949655.2010.520163

850

851

Zazzo, A., Lécuyer, C., Mariotti, A., 2004a. Experimentally-controlled carbon and oxygen isotope exchange between bioapatites and water under inorganic and microbially-mediated conditions. Geochim. Cosmochim. Acta 68, 1–12. https://doi.org/10.1016/S0016-7037(03)00278-3

852

853

854

855

Zazzo, A., Lécuyer, C., Sheppard, S.M., Grandjean, P., Mariotti, A., 2004b. Diagenesis and the reconstruction of paleoenvironments: a method to restore original δ¹⁸O values of carbonate and phosphate from fossil tooth enamel. Geochimica et Cosmochimica Acta 68, 2245–2258.

856

857

858

859

Ziegel, E.R., 2001. Standard probability and statistics tables and formulae. Technometrics 43, 249.

860

861

862
863

864 **Figure captions:**

865

866 Figure 1: Predicted changes in the δ<sup>18</sup>O(A)ₜ values of a fish tooth that reacts in closed system with a post-depositional aqueous solution of meteoric origin as a function of its relative amount (W/A) and temperature (20°C and 50°C chosen as arbitrary values).

869

870 Figure 2: Predicted changes in the δ<sup>18</sup>O(A)ₜ values of a fish tooth that reacts in open system with a post-depositional aqueous solution of meteoric origin as a function of its relative amount (W/A) and temperature (20°C and 50°C chosen as arbitrary values).

873

874 Figure 3: Simulations of δ<sup>18</sup>O(A)ₜ values were computed by using the ApaOxIMOD program. The Monte-Carlo simulations were performed using an array of 2,000,000 events. The following parameters were used: temperature T varies from 10°C to 35 °C; δ<sup>18</sup>O(W)ᵢ varies from -8.0 to -2.0, and δ<sup>18</sup>O(A)ᵢ was always fixed to 20.0. A/W was considered in the 0.05-0.15 range. An average analytical uncertainty was simulated (1σ = 0.1). The histograms were built with a bin number of 500.

880

881 Figure 4: Simulations of δ<sup>18</sup>O(A)ₜ values were computed by using the ApaOxIMOD program. The Monte-Carlo simulations were performed using an array of 2,000,000 events. The following parameters were used: temperature T varies from 10°C to 60°C; δ<sup>18</sup>O(W)ᵢ and δ<sup>18</sup>O(A)ᵢ were always fixed to -8.0 and 20.0, respectively. A/W was considered in the 0.05-0.95 range (left column) with five subranges highlighted in the right column. An average
886 analytical uncertainty was simulated (1σ = 0.1). The histograms were built with a bin number
887 of 500.

888
889 Figure 5: Risk (%) of incorrectly assuming a normal distribution as a function of the number
890 of samples of skeletal phosphate remains analysed. For example, if sixty bioapatite samples
891 are analysed, a normality test can be performed with a 5% risk, and fifty samples with a 10%
892 risk.

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

$\delta^{18}\text{O(A)}_i = 22\text{‰} \text{ (VSMOW)}$

$\delta^{18}\text{O(H}_2\text{O)}_i = -8\text{‰} \text{ (VSMOW)}$
<figure><img src="image_3.png" /></figure>

$\delta^{18}\text{O(A)}_i = 22\text{‰} \text{ (VSMOW)}$

$\delta^{18}\text{O(H}_2\text{O)}_i = -8\text{‰} \text{ (VSMOW)}$
Effect of Variation of T and δ18OWi
within a range of W/A from 0.05 to 0.15
C. Lecuyer and JP. Flandrois 2023

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

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

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

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

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

<figure><img src="image_9.png" /></figure>
Effect of Variation of T in multiple W/A ratio ranges and for multiple Temperature values

C. Lecuyer and JP. Flandrois 2023

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

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>