Teeth, scales, and turtle shell fragments from a single locality (V1075) in Bighorn County, Montana were approved for consumptive sampling from the Oklahoma Museum of Natural History (OMNH). The geographic location and stratigraphic position of the V1075 bonebed is shown in Figure 2. Detailed locality data is available from the OMNH upon request. The V1075 site has yielded an abundant and taxonomically rich vertebrate microfossil assemblage, especially in comparison to other OMNH Cloverly sites. Museum records claim that the assemblage was collected by bulk screenwashing from near the base of Ostrom’s “Unit VII” but the precise stratigraphic position is unknown. Although there are no direct age constraints at the V1075 site, recent radiometric dates from multiple sites across the Bighorn Basin have constrained the depositional age of the Cloverly to the Albian stage (Carrano et al., 2022; D’Emic et al., 2019). Recent work has provided high-precision U/Pb dates from cryptotephra zircons collected at multiple points in a section near Crooked Creek, approximately 50 km south-southeast of the V1075 site (J. Fekete, pers. Comm., December 2024). These Crooked Creek dates have further constrained the Little Sheep-Himes boundary interval to ∼110–109 Ma. The V1075 fauna is typical of the Cloverly, having tooth morphotypes sharing affinities with Deinonychus and Tenontosaurus. The dinosaur fauna also includes sauropods, carcharodontosaurids, troodontids, and nodosaurids. Neosuchian crocodylomorph teeth are especially abundant, including thousands of specimens of tooth morphotypes sharing affinities with goniopholids, atoposaurids, bernissartiids, and pholidosaurs. Family-level designations of crocodylomorph teeth from the Cloverly are tenuous due to the lack of described macrofossils. We sampled all four of the present neosuchian tooth morphotypes. However, the pholidosaur-like teeth were small and had thin enamel, making recovery of viable sample sizes very difficult. As a result, only three of the neosuchian ecomorphotypes were analyzed. We refer to these three ecomorphotypes as “Neosuchian G” (cf. Goniopholididae), “Neosuchian B” (cf. Bernissartiidae), and “Neosuchian A” (cf. Atoposauridae). There are also many mammal teeth in the assemblage, including multituberculates, triconodonts, symmetrodontids, and tribosphenids. We attempted to sample mammal teeth, but we failed to recover viable sample sizes from these very small fossils with such thin enamel. Fish, lizards, and amphibians are also represented. No enamel-bearing fossils of lizards or amphibians were present in the assemblage, and these groups were not sampled. Most of these clades are represented by tens, hundreds, or thousands of specimens in the V1075 assemblage.

Powdered samples were drilled from the specimens using a Brassler NSK UM50TM dental drill with a tungsten-carbide burr. To minimize destruction of fossil resources, only a small number of specimens were selected for duplicate sampling and analysis. In some specimens, we took bulk samples across the long axis of the tooth to homogenize sampling along ontogenetic growth. Others were sampled serially to assess variation along growth axes. The densest material was targeted in each turtle skeletal fragment. For teeth, care was taken to limit sampling to enamel only, with efforts to minimize inclusion of dentine. Like bone, dentine is known to be more susceptible than enamel to diagenetic alteration of isotopic compositions due to the smaller hydroxyapatite crystals and porous texture of bone and dentine. Ganoine scales were crushed whole by mortar and pestle. For conical fish teeth (cf. Lepisosteidae), dentine was milled out from within the tooth, and the remaining enamel was crushed by mortar and pestle. Specimens that were too small to hold by hand while drilling were held in place on glass slides using cyanoacrylate glue. After sampling, the fossils were removed from the slides by submersion in acetone (C₃H₆O) for 24 h. For taxa with enough specimens, at least ten specimens were sampled at least once for each taxonomic group. Drilled hydroxyapatite powders were converted to silver phosphate (Ag₃PO₄) following (Vennemann et al., 2002; Suarez et al., 2012). Details of the silver phosphate method used in this study can be found in Supplementary Material. Briefly, samples were treated with acetone to remove cyanoacrylate and then with bleach (NaOCl) to remove organics. Powders are then dissolved in nitric acid (HNO₃). Phosphate ions are then separated from Ca and OH ions by reacting the sample solution with KOH and KF. A silver amine solution (AgNO₃, NH₄NO₃, and NH₄OH) is added, resulting in the bonding of silver and phosphate ions, which are precipitated as Ag₃PO₄ crystals. The Ag₃PO₄ crystals are rinsed thoroughly with deionized water and dried. Yields ranged from 116 to 793 μg. For each sample, 400–600 µg of crystals were weighed and wrapped in silver capsules. The wrapped samples were combusted at high temperatures to produce CO gas on a Thermo TC/EA connected to a MAT 253 isotope ratio mass spectrometer (IRMS). δ¹⁸Ο values were corrected relative to VSMOW using NIST 80 and 81. Reproducibility was monitored by replicate analyses of Acros Ag₃PO₄ (σ = ∼0.35‰ VSMOW).

δ¹⁸O_phosphate (abbreviated as δ¹⁸O_p) data from multiple taxa were used to estimate the δ¹⁸O of parent waters and the temperature of surface waters. Reconstructing these two parameters from δ¹⁸O_p alone requires analysis of two types of taxa: taxa with constrained body temperatures (T-constrained) and taxa with body temperatures that track environmental temperature (T-unconstrained). Broadly, this multi-taxon approach involves a two-step process that leverages empirically observed relationships between δ¹⁸O_p of bioapatite, δ¹⁸O of parent water, and mineral growth temperature. First, δ¹⁸O_{surface_water} (abbreviated as δ¹⁸O_sw) is estimated from the δ¹⁸O_p of T-constrained taxa. Second, δ¹⁸O_sw is combined with δ¹⁸O_p of T-unconstrained taxa to solve for temperature. Here, we outline our approach to this using δ¹⁸O_p of turtles, crocodylomorphs, and fish following Barrick et al. (1999), Amiot et al. (2007), and Fricke and Wing (2004). R Markdown documentation of the R code used in this analysis is provided in the Supplementary Material. A simple deterministic example of the multi-taxon procedure we followed is outlined below using the equations provided:

1. Applying Equation 1 to the mean δ¹⁸O_p value of V1075 Glyptops fossils (13.8‰ VSMOW) yields a δ¹⁸O_sw value of −8.4‰ (VSMOW).

δ¹⁸O_sw was estimated from the δ¹⁸O_p of two semi-aquatic taxa: the turtle Glyptops and the crocodylomorph morphotype “Neosuchian G.” Previous work has demonstrated linear relationships between δ¹⁸O_p and δ¹⁸O_sw for semi-aquatic turtles and crocodilians in modern settings (Barrick et al., 1999; Amiot et al., 2007). Modern crocodilians and turtles thermoregulate through behavior such as sun basking, maintaining constrained body temperatures when conditions allow. Paired measurements of turtle δ¹⁸O_p and habitat δ¹⁸O_sw suggest that in some turtle taxa, bone growth occurs within a body temperature range of ∼31°C–33°C (Barrick et al., 1999). We reanalyzed Barrick et al. (1999) data and accurately reproduced the linear regression between δ¹⁸O_p and δ¹⁸O_sw (adjusted R² = 0.98). The regression model we reproduced is shown as Equation 1. Amiot et al. (2007) found a similar pattern in modern crocodilians. We reanalyzed the data from Amiot et al. (2007) and accurately reproduced the relationship between δ¹⁸O_p and δ¹⁸O_sw (adjusted R² = 0.84). The regression model we reproduced from data in Amiot et al. (2007) is presented as Equation 2.

This approach to reconstructing δ¹⁸O_sw relies on several assumptions. One assumption is that body water δ¹⁸O of Glyptops and “Neosuchian G” was not significantly influenced by factors such as seasonal variability, evaporation, or heterogeneous water sources, but rather was largely determined by habitat preference and diet (Barrick et al., 1999). Semi-aquatic ectotherms that consume mostly aquatic food sources should have body water δ¹⁸O which correlates with δ¹⁸O_sw. Anatomical evidence strongly suggests a semi-aquatic habit for both Glyptops and “Neosuchian G.” Gracile limbs and jaw elements suggest an aquatic habit and diet for pleurosternid turtles, including Glyptops (Evers, 2023). Cretaceous crocodylomorphs occupied a variety of habitats, including fully aquatic, semi-aquatic, and fully terrestrial (Wilberg et al., 2019). This functional diversity is represented in the V1075 crocodylomorph fauna, as suggested by the disparity in δ¹⁸O_p distributions between morphotypes (referred to hereafter as taxa for simplicity). However, we interpret an aquatic or semi-aquatic habit and diet for “Neosuchian G” based on: a) comparison of δ¹⁸O_p with sympatric taxa, and b) the presence of apicobasal ridges on the teeth. We excluded Naomichelys, the other easily identifiable turtle genus in the V1075 assemblage, from δ¹⁸O_sw reconstructions because its anatomical features indicate a terrestrial lifestyle (Joyce et al., 2015; Joyce, 2017). We excluded “Neosuchian A” and “Neosuchian B” because their reconstructed δ¹⁸O_water values differed significantly from “Neosuchian G” and Glyptops. Also, “Neosuchian A” teeth share affinities with Atoposauridae, a complex group which includes both fully terrestrial and semi-aquatic forms (Tennant et al., 2016). Another key assumption we make is that Glyptops and “Neosuchian G” maintained a similar range of body temperatures to that of the modern taxa that the regressions are based on. Barrick et al. (1999) found that some turtle taxa seem to demonstrate a better correlation between δ¹⁸O_p and δ¹⁸O_sw than other taxa. However, taxa with δ¹⁸O_p that correlated poorly with δ¹⁸O_sw were not sampled from a broad range of latitudes (and therefore temperature regimes), and future expansion on Barrick’s observations may show that the strong correlation between δ¹⁸O_p and δ¹⁸O_sw is true for most or all semi-aquatic turtles. Amiot et al. (2007) found that the regression between δ¹⁸O_p and δ¹⁸O_sw did not vary significantly between crocodilian taxa. The distributions of empirical δ¹⁸O_sw estimates from Glyptops and “Neosuchian G” do not differ significantly from one another, suggesting that our assumptions underlying reconstructions from Glyptops and “Neosuchian G” might be valid.

δ¹⁸O_sw estimated from δ¹⁸O_p of T-constrained taxa (in this case Glyptops and “Neosuchian G”), along with δ¹⁸O_p of T-unconstrained taxa, can be used to estimate temperature. A strong relationship between δ¹⁸O_p-δ¹⁸O_water and temperature has been demonstrated repeatedly in aquatic taxa (Longinelli, 1966; Longinelli and Nutti, 1973a; Longinelli and Nutti, 1973b; Kolodny et al., 1983). Puceat et al. (2010) presented an improved calibration using the measured value of NBS120c to correct for inter-laboratory variability. This calibration by Puceat et al. (2010) is shown as Equation 3. Garfish (Lepisosteidae, referred to here simply as “gar”) are good candidates for a T-unconstrained taxon, as they are ectotherms with limited ability to regulate their body temperature through behavior (Fricke and Wing, 2004). We used δ¹⁸O_p of gar teeth and scales, along with our δ¹⁸O_sw estimate from Glyptops and “Neosuchian G”, to estimate temperature using Equation 3.

Several assumptions underlie the multi-taxon approach to estimating temperature. First, we must assume that the taxa used in the reconstructions all sampled environmental waters with similar δ¹⁸O. If taxa derived their δ¹⁸O values from isotopically distinct water sources, this could lead to inaccurate temperature estimates. It is difficult to determine with high certainty whether fossil animals sampled from isotopically distinct waters, but constraints can be estimated in the case of the V1075 aquatic fauna. Interpreting the depositional environment allows constraints on inhabited water conditions. The V1075 assemblage was collected from gray mudstone facies near the base of Ostrom’s “Unit VII” in the Himes Member. Studies of the Cloverly across the Bighorn Basin suggest that gray mudstones in this interval represent deposition on a poorly drained fluvial floodplain (D’Emic et al., 2019; Orchard, 2024). Vertebrate microfossil bonebeds containing aquatic, semi-aquatic, and fully terrestrial animals within floodplain deposits are typically interpreted as in situ accumulations in floodplain wetlands (Rogers et al., 2007; Rogers and Brady, 2010). The V1075 fossils likely accumulated in such an environment, as there is no evidence in the Cloverly of large, deep lakes, and the V1075 assemblage is hosted in fine-grained deposits inconsistent with channel-lag deposits. Still, we can’t know if these floodplain wetlands were ever connected at the surface with nearby stream channels, or if the V1075 aquatic taxa moved between these wetlands and stream channels. Time spent in a river channel which drained higher elevation uplands could have exposed some of the V1075 aquatic fauna to lower δ¹⁸O than that of the local floodplain wetlands.

To derive an accurate and useful multi-taxon temperature estimate, we must also make assumptions about the seasonality of skeletal growth in the taxa used. The multi-taxon approach relies on the assumption that the two taxa grew their skeletal material during approximately the same seasonal period. This assumption is critical because δ¹⁸O values of precipitation can vary significantly between seasons. If the two taxa incorporated oxygen from meteoric waters into their skeletal phosphates at different times of the year, their body waters would likely reflect distinct δ¹⁸O values during skeletal growth, resulting in inaccurate temperature estimates. Also, for meaningful comparison with other temperature data, it is necessary to estimate which part of the annual temperature cycle the multi-taxon estimate represents.

Our multi-taxon δ¹⁸O-derived temperature estimate likely reflects the mean warm-season temperature of the local floodplain wetlands. However, it is challenging to determine which portion of the annual climate cycle the Cloverly ectothermic taxa recorded during tissue growth. In general, ectothermic metabolic activity peaks when habitat temperatures reach the upper range of a species’ thermal niche (Neubauer and Andersen, 2019). This is demonstrated for turtles and crocodilians by Barrick et al. (1999) and Amiot et al. (2007). Yet, it remains unclear during which months of the year the Cloverly semi-aquatic turtles and crocodylomorphs could sustain such warm, stable body temperatures. This uncertainty also applies to V1075 gar, as the growth regime of modern gar is not well understood. However, general metabolic temperature requirements in ectotherms suggest it likely that Cretaceous gar species also favored warm-season temperatures. While “warm season” cannot be precisely defined in this context, we interpret the multi-taxon temperature estimate as representative of the mean temperature during the warmest half of the year.

To account for all sources of uncertainty, we employed a Monte Carlo approach to estimate the uncertainty surrounding our mean proxy estimates by simulating many potential outcomes based on the variability in our δ¹⁸O_p data and the associated regression model parameters. Specifically, for each proxy measurement, we generated thousands of synthetic datasets by randomly sampling from normal distributions defined by the means and standard deviations of the observed values. Additionally, we incorporated the uncertainty in the regression model coefficients (intercept and slope) by sampling these parameters from their respective distributions, defined by their standard errors. Residual variability, representing unexplained error in the model, was also simulated and added to each estimate. This approach allowed us to propagate all sources of uncertainty (measurement error, regression uncertainty, and residual variability) into the final estimates. By aggregating the results of these simulations, we calculated a mean estimate and 95% confidence intervals, providing a robust quantification of uncertainty around our mean estimates.

We generated a large multi-taxon δ¹⁸O_p dataset and used it to characterize variation across taxa and ecomorphotypes. Single measurements ranged from 9.5‰ to 23.2‰ (VSMOW), with ecomorphotype means spanning 12.0‰–19.1‰. Some taxa showed substantial within-specimen or within-taxon variation. Notably, neosuchian crocodyliforms displayed significant differences among ecomorphotypes, and terrestrial taxa such as sauropods and ornithischians yielded consistently higher δ¹⁸O_p values. These data establish the isotopic framework for interpreting water source and paleoenvironmental differences among taxa.

δ¹⁸Ο_p results are shown in Table 1 and Figure 3A, and are available in Supplementary Material. Single measurements averaged 15.2‰ (VSMOW) and ranged from a minimum of 9.5‰ (villiform fish tooth, cf. Lepisosteidae) to a maximum of 23.2‰ (small theropod tooth, cf. Deinonychus sp.). To minimize destruction of specimens, most specimens were sampled only enough for one analysis (>600 µg). However, we sampled some specimens for duplicate analyses. Variability within individual specimens ranged from σ = <0.1‰ to σ = 1.9‰ (VSMOW). Only two specimens, a gar and an aquatic turtle, showed individual variation greater than σ = 1. Mean δ¹⁸Ο_p values by ecomorphotype range from a minimum of 12.0‰ relative to VSMOW (hybodontiform sharks) to a maximum of 19.1‰ (Sauropoda). Measurements from tooth enamel of the semi-aquatic “Neosuchian G” ecomorphotype (cf. Goniopholididae, σ = 2.9‰ VSMOW) averaged 14.2‰, ranging from a minimum individual measurement of 9.8‰ to a maximum of 17.6‰. In addition to enamel, we also analyzed dentine from two “Neosuchian G” teeth. These resulted in values of 16.8‰ and 15.9‰. On average, we found that molariform crocodilian teeth (cf. Bernissartiidae, “Neosuchian B”) and spatulate crocodilian teeth (cf. Atoposauridae, “Neosuchian A”) yielded significantly more positive δ¹⁸Ο_p values than “Neosuchian G”. “Neosuchian B” δ¹⁸Ο_p values ranged from 13.1‰ to 19.2‰, with an average of 16.6‰ (σ = 1.9‰). “Neosuchian A” averaged 17.0‰ (σ = 2.2‰), ranging from a minimum of 15.1‰ to a maximum of 23.1‰. Ornithischian teeth composition ranged from 13.5‰ to 19.5‰, with an average of 15.5‰ (σ = 2.1‰). Sauropod tooth enamel averaged 17.6‰ (σ = 1.9‰) with individual measurements ranging from 15.1‰ to 22.5‰. Lepisosteid scale (mean = 14.5‰, σ = 1.5‰) and teeth (mean = 12.9, σ = 2.0‰) values range from a minimum individual measurement of 9.5‰ to a maximum of 16.2‰. Carapace fragments of the semi-aquatic turtle Glyptops sp. (mean = 14.1‰, σ = 1.8‰) ranged in δ¹⁸Ο_p from a minimum individual measurement of 11.1‰ to a maximum of 17.1‰. The terrestrial turtle Naomichelys sp. (σ = 0.6‰) showed an average δ¹⁸Ο_p of 14.3‰, with values ranging from 13.3‰ to 15‰. Small theropod specimens (maniraptoran teeth) showed δ¹⁸Ο_p values ranging from a minimum of 12.7‰ to a maximum of 23.2‰, with a mean of 15.8‰. One large theropod (allosauroid) tooth was available for sampling, which yielded a mean δ¹⁸Ο_p of 14.7‰ from four serial samples taken along the growth axis (n = 3, σ = 0.16‰).

Using δ¹⁸O_p values from turtles, crocodylomorphs, gar, and a co-analyzed standard (NBS120c), we conducted Monte Carlo simulations and deterministic calculations to estimate surface water δ¹⁸O (δ¹⁸O_sw) and reconstruct paleotemperature. Simulated δ¹⁸O_sw values ranged from −11.2‰ to −3.6‰ (95% CI), and deterministic estimates were in close agreement. Combining distributions yielded a mean δ¹⁸O_sw of −7.9‰ (95% CI: −10.1, −5.5‰). Applying temperature equations to these values, we estimate a mean surface water temperature of 26°C (95% CI: 8, 43°C). This is the first quantitative paleotemperature reconstruction of the Cloverly environment.

We also simulated probability distributions of coefficients and errors for the Barrick et al. (1999), Amiot et al. (2007), and Puceat et al. (2010) regressions. We simulated δ¹⁸Ο_p distributions for each turtle and crocodylomorph taxon, as well as for gar and NBS120c based on our measured values for each. The mean δ¹⁸O_p (VSMOW) estimates, along with their 95% confidence intervals, were as follows: 13.8‰ (12.7, 14.9‰) for Glyptops, 14.0‰ (12.5, 15.6‰) for “Neosuchian G,” 13.8‰ (13.0, 14.5‰) for gar, and 22.6‰ (22.4, 22.7‰) for NBS120c. We then simulated δ¹⁸Ο_sw reconstructions for each of the turtle and crocodylomorph groups. The estimated mean δ¹⁸Ο_sw values and their 95% confidence intervals are shown in Figure 3B. The mean simulated δ¹⁸Ο_sw (VSMOW) estimate for Glyptops is −8.3‰ (95% CI: −10.5, −6.0‰). For “Neosuchian G,” the mean estimate is −7.5‰ (95% CI: −11.2, −3.6‰). For comparison with our simulated values, we made simple deterministic calculations of δ¹⁸Ο_sw for each empirical δ¹⁸Ο_p value for Glyptops using Equation 1 and for “Neosuchian G” using Equation 2. The mean calculated δ¹⁸Ο_sw value for Glyptops was −8.3‰ (σ = 0.6‰). For “Neosuchian G,” it was −7.6‰ (σ = 0.7‰). The deterministic δ¹⁸Ο_sw distributions for Glyptops and “Neosuchian G” are not significantly different (p = 0.43). We combined the two simulated δ¹⁸Ο_sw distributions and ran a Monte Carlo simulation on them to produce a single probabilistic distribution of δ¹⁸Ο_sw. The mean of simulated δ¹⁸Ο_sw means is −7.9‰ (95% CI: −10.1, −5.5‰). To calculate a mean temperature value and construct a 95% confidence interval, we applied Equation 3 using the simulated distributions of the three variables: δ¹⁸Ο_sw, gar δ¹⁸Ο_p, and NBS120c δ¹⁸Ο_p. The mean temperature estimate is 26°C (95% CI: 8, 43°C).

To contextualize the temperature reconstruction, we converted the estimated mean annual warm season water temperature (26°C) into a mean annual warm season air temperature (24°C; 95% CI: 8, 41°C) using a regression model from modern lake and climate data. While prior work cautions against applying such transforms in warm climates, reanalysis of modern datasets supports a moderate correlation between warm season water temperatures and air temperatures in climates with MAAT >12°C. We also estimated a mean annual water temperature of ∼19°C based on an assumed seasonal range in temperatures.

The conversion of our mean annual warm season water temperature (26°C) to mean annual warm season air temperature (24°C) utilized data from modern lakes. Modern lake and air temperature data from Hren and Sheldon (2012) indicated that while MAWSAT and MAAT are closely related in many settings, there is no significant relationship in climates with mean annual air temperature (MAAT) above 12°C–15°C. It is likely that Albian MAAT in the Cloverly region was greater than 15°C, and so this setting is not appropriate for any of the water-air temperature transform functions defined by Hren and Sheldon (2012). However, we reanalyzed data from Hren and Sheldon (2012) and found that even in climates with MAAT above 12°C, there is a significant relationship (p = 9 × 10⁻¹²), albeit moderately strong (adjusted R² = 0.55), between the mean water temperature and mean air temperature during the warm season (i.e., April-September). Given our MAWSWT estimate of 26°C, and assuming a mean annual range of temperatures (MART) of ∼20°C–25°C based on previous estimates for this region during the Campanian, we estimate a minimum mean annual water temperature (MAWT) of ∼19°C. This simple estimate assumes that MAWT = MAWSWT – (MART/4). Assuming a similar MART to the Campanian of this region is justified because the paleolatitude and global temperature regime in the Albian were comparable to those of the Campanian (Burgener et al., 2019). With an estimated MAWT well above 12°C, we find it justified to apply our MAWSWT-MAWSAT regression model to our V1075 MAWSWT estimate. In doing so, we find that the model predicts a MAWSAT of 24°C (95% CI: 8, 41°C).

Our temperature estimates varied considerably, and our mean estimate is highly uncertain. We performed a series of targeted Monte Carlo simulations to estimate the relative contributions of error from each input and regression coefficient. We conducted a series of Monte Carlo simulations to assess different sources of uncertainty to the total variance in the estimates of δ¹⁸Ο_sw and temperature. Using 1,000 iterations, we varied synthetic δ¹⁸Ο_p, regression slope, intercept, and residual error simultaneously to calculate the total variance of the model outputs. We then performed targeted simulations, varying one input at a time while keeping the others fixed at their mean values, to isolate the variance contributed by each input. The variance contributions were normalized by the total variance to determine their relative importance. The variance contributions to the total uncertainty in δ¹⁸Ο_sw and temperature estimates are distributed across several sources.

For δ¹⁸Ο_sw reconstructions from both Glyptops and “Neosuchian G,” residual error in the underlying regression models from Barrick et al. (1999) and Amiot et al. (2007) was the largest source of uncertainty, contributing approximately 31% and 35% of the total variance, respectively. The regression intercept also contributed significantly, accounting for 26% of the variance for Glyptops and 32% for Neosuchian G, reflecting uncertainty in the baseline offset between δ¹⁸Ο_p and δ¹⁸Ο_sw. Slope error accounted for 17% of the variance for Glyptops and 21% for Neosuchian G. Variability in δ¹⁸Ο_p contributed 25% for Glyptops but only 12% for “Neosuchian G”. While δ¹⁸O_p variability was a notable factor, the dominant sources of uncertainty were the residual error and intercept variability, underscoring the limitations of the regression models in precisely capturing the relationship between δ¹⁸Ο_p and δ¹⁸Ο_sw. These results highlight the importance of addressing residual noise and refining regression parameters to improve the precision of this approach to δ¹⁸Ο_sw reconstructions in the future.

For temperature reconstructions, we were unable to reproduce the pooled regression presented by Puceat et al. (2010). The Puceat et al. (2010) pooled regression combined new data with previously published data from Longinelli and Nuti (1973a) and Kolodny et al. (1983). We compiled the Puceat et al. (2010) data and the Longinelli and Nuti (1973b). We were unable to reproduce the Kolodny et al. (1983) data as presented in Puceat et al. (2010). We performed a regression analysis on the combined Puceat et al. (2010) and Longinelli and Nuti (1973a) data and found that the regression coefficients and errors were very similar to those reported by Puceat et al. (2010) for their pooled model, suggesting that the model we present here is comparable in accuracy and precision. We used Monte Carlo simulations to evaluate the relative error contributions of slope, intercept, residual error, δ¹⁸Ο_sw, and gar δ¹⁸Ο_p. The results indicate that δ¹⁸Ο_sw contributes the most uncertainty (∼40% of the total variance) to the temperature estimates. Intercept (29%) and slope (28%) were also significant contributors of error. Gar δ¹⁸Ο_p (1%) and NBS120c (<1%) contributed comparatively negligible errors. This demonstrates the high sensitivity of the temperature regression to variability of δ¹⁸Ο_sw inputs, as well as the need for additional data to more precisely calibrate the regression model.

Error analysis of the MAWSAT estimate indicates that uncertainty in the MAWSWT estimate is the dominant source of error, significantly outweighing the contributions from the regression coefficients and residual error in the transformation function. While the transform function is only based on a moderately strong relationship (adjusted R² = 0.55), the propagated error from δ¹⁸O_sw through the temperature calculations is the largest contributor to the total uncertainty, far exceeding the error introduced by the regression model itself.

Inter-taxon comparisons of δ¹⁸Ο reveal evidence for resource partitioning among the Cloverly fauna. Significant offsets between reconstructed δ¹⁸Ο_sw values of some crocodylomorph taxa suggest distinct habitat preferences. Although we did not reconstruct δ¹⁸Ο_sw from dinosaur or fish δ¹⁸Ο_p, we can make qualitative inferences based on the offsets between groups. For most V1075 taxa, the observed isotopic niches generally correspond with the paleobiological frameworks inferred previously from other evidence. However, there are exceptions, as discussed below.

Both turtle and crocodylomorph lineages have repeatedly adapted to terrestrial and aquatic habitats throughout their evolutionary histories. Comprehensive anatomical descriptions or a clear phylogenetic framework is required to infer an aquatic or terrestrial habit for a fossil taxon in either of these groups. Many turtle and crocodylomorph fossil occurrences in Cretaceous deposits are based on isolated elements or incomplete skeletons. Thus, contextual information is often missing that would help distinguish aquatic and terrestrial taxa. Atoposaurid-like teeth (“Neosuchian A”) from V1075 also result in δ¹⁸Ο_sw reconstructions consistent with the terrestrial habitat inferred from morphology of body fossils of similar taxa in other Mesozoic assemblages. Atoposaurid-like teeth from mid-Cretaceous deposits have occasionally been referred to the terrestrial genus Theriosuchus sp. (e.g., Winkler et al., 1990; Cifelli et al., 1999). Although no Theriosuchus or other atoposaurid body fossils are known from Aptian-Albian deposits in North America, a small terrestrial paralligatorid (Wanchaampsus kirpachi) has been described from the Aptian of Texas (Adams, 2015). The atoposaurid-like teeth from the V1075 assemblage may have belonged to a similar terrestrial paralligatorid in the Cloverly region during the Albian. Paralligatorids and atoposaurids are thought to have diverged from a common ancestor sometime in the Jurassic, and recent phylogenetic analyses suggest that Atoposaurids were likely extinct by the Early Cretaceous (Tennant et al., 2016). North American Cretaceous teeth like those that we refer to here as the “Neosuchian G” morphotype have traditionally been treated as belonging to Goniopholididae. Although the true phylogeny of the Cloverly goniopholid-like tooth-bearer is unclear, the “Neosuchian G” teeth show marginal carina and apicobasal ridges that suggest aquatic feeding adaptations like that of goniopholids (McCurry et al., 2019). Our data show agreement between functional morphology and isotopic niche in the V1075 “Neosuchian A” and “Neosuchian G” forms. The turtles Naomichelys and Glyptops show very similar δ¹⁸O_p distributions, averaging higher than most aquatic taxa and lower than most terrestrial taxa. This pattern is consistent with a semi-aquatic habitat preference and suggests that their isotopic niche accurately reflects their habitat niche.

Unlike the other V1075 taxa, δ¹⁸O_p of molariform crocodylomorph teeth (cf. Bernissartiidae, referred herein as “Neosuchian B”) are not consistent with morphological and isotopic interpretations of habitat in Bernissartiid-like crocs from other Mesozoic assemblages. Based on a brevirostrine skull shape and durophagous tooth morphology, earliest Cretaceous bernissartids from Europe were assumed to having been semi-aquatic, feeding largely on freshwater molluscs whose fossils are numerous in the Wealden of the Isle of Wight (Buffetaut and Ford, 1979; Sweetman et al., 2014; Martin et al., 2020). This was supported by the interpretation of the sediment facies of the Lower Cretaceous Wealden deposits as representing a seasonal wetland to lacustrine paleoenvironment (Spagna et al., 2012). Likewise, Bernissartid-like teeth from the Cenomanian of Utah have been interpreted as having been from a semi-aquatic animal based on relatively low oxygen isotopes from tooth phosphate (Suarez et al., 2014; Pouech et al., 2014). In contrast, δ¹⁸Ο_p values from V1075 “Neosuchian B” are higher than those of most aquatic and semi-aquatic forms in the assemblage. There are several scenarios that are most likely to have resulted in relative isotopic enrichment in the V1075 “Neosuchian B” teeth. One possibility is that there is more functional and taxonomic diversity represented by isolated mid-Cretaceous molariform crocodylomorph teeth than what has previously been thought. As with the Cloverly “Neosuchian G” teeth, the phylogenetic placement of “Neosuchian B” is somewhat uncertain. Although Aptian-Albian molariform neosuchian teeth in North America have often been attributed to or compared with the Bernissartiidae from the Barremian-Hauterivian of Europe and Africa, there are no known body fossils of bernissartiids from North American deposits (Ullmann et al., 2011; Oreska et al., 2013). It is possible that the V1075 “Neosuchian B” teeth belonged to a paralligatorid crocodylomorph that adapted to terrestrial environments. Although evolutionary habitat shifts between the terrestrial, freshwater, and marine realms are common in the crocodylomorph fossil record, these shifts usually correspond with major phylogenetic divergence rather than at the species or genus level (Wilberg et al., 2019). Bernissartia diverged from Paralligatoridae in the middle Jurassic, and Crocodylia evolved from paralligatorids. Some Crocodilian lineages are thought to have adapted from freshwater to terrestrial environments in the late Cenozoic, so it is plausible that a terrestrial durophagous paralligatorid may have lived in the Cloverly region during the Albian (Montefeltro et al., 2013; Joyce et al., 2015; Joyce, 2017; Wilberg et al., 2019). This scenario cannot be ignored given the coarse diagnostic utility of crocodylomorph teeth and the paucity of bernissartiid body fossils in mid-Cretaceous deposits. In addition, the Bernissartiid-like teeth analyzed by Suarez et al. (2014) are likely millions of years younger than those presented here, and this invokes possible evolutionary differences between these two forms. Another explanation might be that the Cloverly bernissartiid was in fact aquatic like its confamilial relatives, but that its δ¹⁸O_sw was higher due to time spent in high-δ¹⁸O waters. If we assume this is the case, there are two most likely sources of high-δ¹⁸O surface waters that would have been available in the Cloverly environment. Either the water body was evaporatively enriched in ¹⁸O relative to precipitation (lake or pond), or it was influenced by seawater (estuary). With δ¹⁸O alone, it is impossible to discern between the two. If we assume that “Neosuchian B” δ¹⁸Ο_sw was higher due to time spent in an estuarine or marine environment, we must infer that “Neosuchian B” adapted a higher salinity tolerance than other V1075 crocodylomorphs, turtles, and fish sampled in this study. Another possible explanation might be that “Neosuchian B” spent more time sun basking than other amphibious reptile taxa in the V1075 fauna. Any one of the scenarios given above might have resulted in higher body water δ¹⁸Ο in “Neosuchian B” relative to presumed aquatic taxa (“Neosuchian G” and Glyptops). Regardless of what combination of factors resulted in higher δ¹⁸O_sw in “Neosuchian B,” it is evident that this taxon differed from co-occurring aquatic reptiles in habitat or thermoregulatory behavior.

There are no established methods for reconstructing body water δ¹⁸Ο from dinosaurs, although some have attempted to use mammal-derived empirical regressions to do so (Cullen et al., 2020). This is problematic due to assumptions about analogous physiology between mammals and dinosaurs and because these equations require quantitative constraints on atmospheric humidity, which are currently unavailable for the Cloverly (Kohn, 1996). However, qualitative inferences about the δ¹⁸Ο composition of body water in the V1075 dinosaur fauna can be made if we assume that the δ¹⁸Ο_p distributions reflect patterns in body water composition. Sauropods show the highest δ¹⁸Ο_p values of any dinosaur group in the assemblage. This is consistent with δ¹⁸Ο_p datasets from many other sauropod-bearing Cretaceous assemblages and has been attributed to the higher isotopic composition of the water in their leafy diet (Amiot et al., 2011; Suarez et al., 2014; Amiot et al., 2015; Amiot et al., 2021). On average, ornithischians from the V1075 assemblage generally have lower δ¹⁸Ο_p values than the sauropods. This could be explained by food resource partitioning, differences in physiology, differences in habitat preference, or differences in water consumption (ratio of food water and drinking water). Maniraptoran δ¹⁸Ο_p values are more like ornithischians than sauropods, but all the dinosaurian groups show higher average δ¹⁸Ο_p values than the aquatic and semi-aquatic taxa. We analyzed two fish taxa: Lepisosteidae and Hybodontiformes. Both groups show similar δ¹⁸Ο_p distributions with lower means than any other V1075 taxa.

As seen in Figure 4, our reconstructed δ¹⁸O_sw values are generally consistent with proxy-based δ¹⁸O_sw reconstructions from other assemblages in the Western Interior when adjusted for paleolatitude (Suarez et al., 2014; Suarez et al., 2021; Harper et al., 2021). Previous reconstructions of δ¹⁸O_sw in the mid-Cretaceous Western Interior significantly underestimated paleolatitude. Recent tectonic modeling suggests a paleolatitude of ∼52 °N at 110 Ma for the Cloverly, roughly 10 °N higher than previously assumed for the Cloverly (Suarez et al., 2014; van Hinsbergen et al., 2015; Suarez et al., 2021; Vaes et al., 2023). Adjusted for these newer paleolatitudinal reconstructions, δ¹⁸O_sw reconstructions plot near the gradient of modern δ¹⁸O_{precipitation}, and well above the range modeled for the Aptian-Albian by Suarez et al. (2011). However, this modeled δ¹⁸O_{precipitation} range presented originally by Suarez et al. (2011) is almost certainly an underestimate given the following. The Suarez model reconstructed a theoretical δ¹⁸O_{precipitation} gradient by combining a large empirical dataset of δ¹⁸O_carbonate (sideritic and calcitic calcretes) with several latitudinal temperature gradients previously modeled for the mid Cretaceous (Barron, 1983; Ufnar et al., 2002; Zhou et al., 2008). These models were based on underestimated paleolatitude parameters. This could also cause shifts in the latitudinal reconstruction such that what was previously thought to be relatively low δ¹⁸O_{precipitation} values are more in line with expectations at more northerly latitudes. In addition, the Suarez model assumes carbonate δ¹⁸O reflects mean annual surface water values, rather than a seasonally-biased estimate. Evidence increasingly suggests that soil carbonates typically (with some exceptions) form in the warm months and thus record seasonal bias towards summer conditions (Breecker et al., 2009; Peters et al., 2013; Burgener et al., 2016; Huth et al., 2019; Kelson et al., 2020). Seasonal biases in temperature estimates may have led Suarez et al. (2011) and Ufnar et al. (2002) to underestimate temperatures in calculation of δ¹⁸O_{precipitation}. This would have the effect of lowering the δ¹⁸O of water reconstructions. We expect that an updated model will shift the predicted δ¹⁸O_{precipitation} range closer to the modern gradient and corroborate that the Western Interior vertebrate δ¹⁸O_p water isotope reconstructions, including those from V1075, will be consistent with warm-season averages.

Many of the taxonomic groups in our dataset show considerable scatter in δ¹⁸Ο_phosphate values, which we interpret as being due in part to short-term climatic variability. Given the geologically short period of time that our data represents (10⁴–10² years of time averaging), we can rule out the influence of orbitally-forced (10⁵ yrs) climate fluctuations (Rogers et al., 2007; Rogers and Brady, 2010). In highly mobile taxa, some of this scatter may be attributed to movement between isotopically distinct reservoirs. For example, fish may have migrated between floodplain wetlands, tributary streams that were primarily sourced by local precipitation (higher δ¹⁸O), low-order (trunk) rivers with lower δ¹⁸O water influenced by high-elevation catchments in the Sevier highlands, and perhaps even brackish estuaries on the WIS coastal margin. Likewise, it is plausible that some Cloverly dinosaurs may have inhabited broad geographic ranges, or at least occasionally migrated to other parts of the region and ingested waters from a variety of reservoirs. Behavioral influence aside, climatic factors that may have contributed to the scatter would be different between habitat types. Some of the dinosaur taxa represented in our dataset show some degree of positive skew in their δ¹⁸O_p distributions. This positive skew may indicate periods of arid conditions with increased evaporation precipitation ratios and/or a significant decrease in relative humidity, either of which might have increased δ¹⁸O_{body_water} in those taxa. Conversely, there are low outlier δ¹⁸O_p values in the aquatic and semi-aquatic taxa. This may indicate periods of high summer rainout effect or periods of increased runoff from the Sevier Highlands.

The MAWSAT estimate of 24°C (95% CI: 10, 41°C) falls within a range predicted previously by proxy-model data assimilation for this paleolatitude and is on trend with other Aptian-Albian empirical estimates from the WIB (Figure 4; Judd et al., 2024). The Albian saw a transition from Early Cretaceous normal greenhouse conditions to the extreme hothouse conditions of the early Late Cretaceous. Between the beginning and end of the Albian stage, global sea surface temperatures increased by roughly 10°C (O’Brien et al., 2017). Our estimate of 24°C suggests that Albian MAAT might have been ∼10°C higher than today in the Bighorn Basin region. If accurate, this estimate implies that temperatures in the Bighorn Basin region during the mid-late Albian were still ∼10°C cooler during Cloverly deposition than during the following Cretaceous Thermal Maximum. This generally agrees with previous reconstructions of the global temperature trend from Albian to Turonian time (O'Brien et al., 2017; Grossman and Joachimski, 2022; Judd et al., 2024). Our estimated warm-season air temperature of 24°C includes a wide uncertainty range (10°C–41°C) yet falls within the range of values expected during the broader warming trend from the Aptian–Albian Cold Snap to the Cretaceous Thermal Maximum. We interpret our estimate as an annual mean of warm season temperatures, and if this is true, then it follows that our value would be higher than most MAAT estimates. Marine TEX₈₆ and δ¹⁸Ο data suggest that the warming between the Aptian-Albian Cold Snap and the Cretaceous Thermal Maximum occurred in pulses during the early-mid Albian, the late Albian, and then again in the late Cenomanian (O'Brien et al., 2017). Our data suggest that the V1075 fauna lived after the end of the Aptian-Albian Cold Snap but prior to late Albian warming.

Multiple lines of evidence suggest that V1075 enamel preserves primary oxygen isotope compositions. First, we analyzed phosphate-bound oxygen which is more robust than carbonate bound oxygen, unless significant microbial alteration occurs, is the preferred phase to analyzed (Blake et al., 1998; Zazzo et al., 2004). Second, δ¹⁸Ο_p values between the sympatric V1075 taxonomic groups show expected magnitudes of offset (several‰ VSMOW) between aquatic, semi-aquatic, and fully terrestrial taxa. If the phosphate-bound oxygen had exchanged significantly with diagenetic fluids, we would expect homogenized δ¹⁸Ο_p values across the assemblage (Kohn and Cerling, 2002). Although it is possible that diagenesis may have caused a δ¹⁸Ο_p bulk shift without homogenization, our reconstructed δ¹⁸Ο_sw values agree with other latitude-adjusted proxy data (Suarez et al., 2021), suggesting that any diagenetic bias in our data might be small (Suarez et al., 2014).

Third, we evaluated differential diagenetic influence between tissue types (i.e., dentine, bone, and enamel). Bone and dentine are more susceptible to diagenetic overprinting than enamel due to differences in crystallinity and porosity (Kohn and Cerling, 2002). We targeted pure enamel in our tooth samples; however, in teeth with thin enamel, it is likely that some dentine was inadvertently included in the drilled sample. We considered the possibility that variable degrees of dentine/enamel mixing during sampling may have been a significant factor behind the heterogeneous distribution of our δ¹⁸Ο_p data. We analyzed both dentine and enamel from one “Neosuchian G” tooth and found that the mean δ¹⁸Ο_p value of dentine (n = 2) was ∼1‰ lower than that of enamel (n = 2). If δ¹⁸Ο_p variability was strongly influenced by mixing isotopically distinct enamel and dentine during sampling, taxa with thin enamel would likely show greater intra-taxon δ¹⁸O_p variation than those with thick enamel, as samples from thin enamel often include some dentine. However, we observe similar standard errors in δ¹⁸O_p of both thick-enamel (Sauropoda and “Neosuchian B”) and thin-enamel taxa (Ornithischia, Maniraptora, “Neosuchian G,” “Neosuchian A”) in the V1075 assemblage (Table 1). V1075 bone phosphates may be slightly altered like dentine. Turtles, unlike other analyzed V1075 taxa, do not have teeth, and so bone is the only material available. We targeted dense cortical bone in carapace fragments and avoided highly porous cancellous regions. Reconstructed δ¹⁸O_sw values between “Neosuchian G” tooth enamel and Glyptops bone are nearly identical, suggesting that both preserve similar signals.

Alteration of bone and dentine towards deep burial diagenetic values could explain the narrow distribution of turtle values near the overall mean. If dentine was more altered than enamel, dentine-enamel mixing in tooth samples might result in mixed samples clustering near a mean value, while pure enamel samples could deviate further, creating skewed distributions toward the overall mean for each taxon. This process would likely produce an overall normal distribution. However, the data do not fully support this hypothesis: the mean Glyptops value does not align with the overall mean, the distribution includes outliers and is not normal, and only two taxonomic groups exhibit non-normal distributions. This suggests that the data is not homogenized around a constrained δ¹⁸Ο as would be the case if bone oxygen had equilibrated with diagenetic fluids.