δ ¹³C and δ ¹⁸O values were obtained from vaquita bone samples from the inner dorsal section of the skull of 13 specimens from the Vertebrate Collection at the Centro de Investigación en Alimentación y Desarrollo (CIAD-Guaymas, Sonora) and the Mammal Collection at the Instituto de Biología-Universidad Nacional Autónoma de México (IB-UNAM). These specimens were collected from the Gulf of Santa Clara, Son., to Puerto Peñasco, Son., from 1983 to 1993 ( Figure 1B ; Table 1 ).

The bone fragments were ground, and no pretreatment was applied because many studies show that pretreatment can affect δ ¹⁸O values (Garvie-Lok et al., 2004; Pellegrini and Snoeck, 2016; Snoeck and Pellegrini, 2015). A minimum weight of 1 mg of each bone sample was collected using a high-precision drill Micro-Mill System at the Marine Mammals Laboratory, IPN-CICIMAR, La Paz, Baja California Sur, México.

The stable isotope analysis in bone was conducted at the University of California Davis Stable Isotope Laboratory (https://stableisotopefacility.ucdavis.edu/) using a Thermo Scientific GasBench II coupled to a Thermo Finnigan Delta Plus XL isotope-ratio mass spectrometer. These were equilibrated for 24 h at 30°C. The resulting data were normalized using commercially available laboratory reference materials (USGS-44, NBS-18, NBS-19, and LSVEC) calibrated to international standard calcium carbonate of Vienna Pee Dee Belemnite (V-PDB), with a standard deviation of ± 0.10‰ for δ ¹³C and ± 0.15‰ for δ ¹⁸O. Results are expressed as δ-values using the equation δ = [(¹⁸O/¹⁶O or ¹³C/¹²C) _sample/(¹⁸O/¹⁶O or ¹³C/¹²C)_standard) - 1] × 1000. Delta (δ) units are reported as parts per thousand (‰).

In marine mammal studies, the δ ¹⁸O values are more commonly presented relative to the Vienna Standard Mean Oceanic Water (V-SMOW) index. In the present study, for comparative purposes, the δ ¹⁸O values in vaquitas were converted from PDB to SMOW according to the following equation (Koch et al., 1997):

δ ¹⁸O (SMOW) = [δ ¹⁸O (PDB) × 1.03086] + 30.86

In order to eliminate the seasonal signal in the variables, we estimated the monthly anomalies of the time-series data for sea surface temperature (SST) and sea surface salinity (SSS) utilizing daily composites from 1983 to 1993. We obtained SST and SSS measurements for the Upper Gulf of California (UGC), delineated to the south by San Felipe (Baja California) and Puerto Peñasco (Sonora) ( Figure 1B ). An anomaly was the deviation between the recorded environmental variable at a specific period and the long-term average of the same environmental variable for the same period and region. A high-resolution SST analysis was developed by NOAA using an optimal interpolation (OISST), with a spatial grid resolution of 0.25 degrees and a temporal resolution of 1 day (Huang et al., 2020). SST values were obtained from the Met Office Marine Data Bank (MDB), while SSS values were extracted from ESA MIRAS SMOS Level-2 swath Network Common Data Format (NetCDF) files (currently V6.62). These environmental data were obtained via ERDDAP (Environmental Research Division’s Data Access Program) (Simons, 2019) using the “rerddap” v1.0.4 and “rerddapXtracto” v1.1.7 packages in R (Chamberlain and Boettiger, 2017; R Core Team, 2023). Additionally, we obtained a record of the average monthly inflow of the Colorado River (FLX) from the International Boundary and Water Commission (http://www.cila.gob.mx/).

To evaluate possible changes in the carbon source of primary producers and the vaquita’s diet changes in relation to the variability in the flow of the Colorado River ( Figure 2C ), two periods were previously defined: (i) from 1982 to 1988, when freshwater arrived in the UGC and seven vaquita bone samples were collected, and (ii) from 1989 to 1993, when the river water stopped entering and six vaquita bone samples were collected.

In order to determine the possible correlation between the environmental variables, the normality of the data was assessed with a Shapiro test based on annual average values ​​derived from monthly composites spanning from 1980 to 2010. Then we conducted a Spearman correlation test, with p-values ​​below the significance threshold of 0.05, suggesting a non-normal distribution.

To test the existence of significant differences between the two periods and age categories, the datasets of δ¹³C and δ¹⁸O bone values were tested separately for normality using Kolmogorov–Smirnov (K–S) tests and assessed for homogeneity of variances with a Levene’s test. Because the datasets were not normally distributed, a Mann–Whitney U test was used. Statistical analyses were conducted using IBM SPSS Statistics v22.

To explore the pattern of δ¹³C and δ¹⁸O data and its relationship with the groups defined based on the flow of the Colorado River and the changes in salinity and temperature that they present, a principal component analysis (PCA) was carried out. Overall, to statistically confirm the link between the dependent variable (stable isotope) and the independent variables (SST, SST, FLX), a generalized linear model (GLM) was used. This model was the most robust for small sample sizes and non-normally distributed variables. To enhance the stability of the model results, we removed the random effect and employed a fixed-effects model. The FLX variable was standardized because it had a much higher significance value than the other variables. Also, we employed Akaike information criterion (AIC) weights to select the optimal model based on reduced AIC scores (Akaike, 1998; Burnham and Anderson, 2002). The delta-AIC (ΔAIC) and AIC weight (AICw) were used as evaluation metrics to rank the models. The AIC quantifies the difference in AIC scores between the best model and the model under assessment. On the other hand, the AIC weight, also known as the AICc weight, represents the proportion of predictive power that the evaluated model contributes to the complete collection of models (Akaike 1998). We generated the models using the glm function from the “Ime4” version 1.1-35 package (Bates et al., 2015) in the R programming language and evaluated the models using the “AICcmodavg” version 2.3-2 package in R (Mazerolle, 2023; R Core Team, 2023).

From satellite data, an SST/SSS diagram was built for each dataset (from 1982 to 1988 and from 1989 to 1993) to determine signs of change in the UGC’s water mass.

Finally, boxplots of the original values of each stable isotope (δ ¹³C and δ ¹⁸O) were built to compare the average values and ranges of bone (Drago et al., 2020) and tooth (Clementz and Koch, 2001) values of marine mammals from different ecosystems, compiled from the literature with respect to the vaquita as a hypersaline habitat species.

The Spearman correlation showed a non-normal distribution (p < 0.05). The p-values between SST and SSS were not significant ( Appendix 1 ). The association between Colorado River flow and SSS was significant, although the correlation coefficient was quite low.

To determine whether the El Niño event influenced the pulses of water released into the UGC and, consequently, the environmental variability, a red polygon was incorporated into Figure 2 to denote the months with documented warm periods, with a threshold of ± 0.5°C for the Oceanic Nino Index (ONI; https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). The ONI did not directly influence the environmental variables. The SST had negative anomalies during warm events, indicating a temperature decrease relative to previous years.

Temperature anomalies exhibited interannual variations during the study period (from 1983 to 1993). However, a range of negative anomalies (from -0.5 to -1.5°C) was observed between 1983 and 1991 ( Figure 2A ). During the later years of the time series (from 1992 to 1993), the highest positive anomalies (from 0.5 to 1.5°C) were observed.

Salinity anomalies showed a clear regime shift from negative to positive anomalies ( Figure 2B ). Negative salinity anomalies (from -0.1 to -1 psu) dominated from 1982 to 1988, with the highest negative values (-1 psu) in 1983. After 1988, a predominance of positive anomalies (from 0.1 to 0.5 psu) was observed until 1993, with a maximum value of +0.5 psu in 1989.

The Colorado River’s flux displayed the inverse of the salinity pattern ( Figure 2C ), with high positive anomalies (from 130 to 730 m³ s^-1) between 1982 and 1988. The highest value (730 m³ s^-1) was observed in 1983. From 1989, negative anomalies (< -70 m³ s^-1) were observed until 1993. A positive anomaly peak (430 m³ s^-1) was present in 1993.

The low salinity values (~< 34.5 psu) in the period from the end of 1983 to 1987 were aligned with high flow (~> 250 m³ s^-1) from the Colorado River into the UGC. In contrast, the increase in salinity in the 1990s (~> 34.5 psu) corresponded with a decrease in the flow of the Colorado River (~< 250 m³ s^-1), except in 1993, when an increment of water from the river was present ( Figure 3 ).

Based on the two datasets (from 1982 to 1988 and from 1989 to 1993), it was tested whether there were significant differences in each variable described above ( Table 2 ). There were no significant differences between the two datasets for SST (t = -1.947, gl = 70, p = 0.056), but there were significant differences for SSS (t = -6.940, gl = 70, p < 0.05) and Colorado River inflow (U = 174.000, Z = -4.741, p < 0.05).

Each dataset’s SST/SSS diagrams showed different water masses, based mainly on the salinity values ( Figure 4 ). From 1982 to 1988, the water mass (blue area) displayed temperature values from ~16 to 24°C, with an average value of 19.4°C, and salinity values from 33.8 to 35 psu, with an average value of 34.4 psu. From 1989 to 1993, the water mass (yellow area) had shorter intervals than in the previous period. It showed temperature values between ~15 and 24°C, with an average of 19.7°C, and salinity values between 34.8 and 35.4 psu, with an average of 34.8 psu.

The hypothesis test showed significant differences for δ ¹³C (U = 144.000, Z = -5.233, p < 0.05) but not for δ ¹⁸O (U = 432.000, Z = -1.744, p = 0.081) between the two datasets (1983–1988 vs. 1989–1993). Richer δ ¹³C values were found for the earlier period (δ ¹³C = -9.1 ± 1.3‰) compared with the later period (δ ¹³C = -10.8 ± 2.1‰). In the case of oxygen, the values of the two groups were very similar, from 1983 to 1988 (δ ¹⁸O = 30.5 ± 0.4‰) and from 1989 to 1993 (δ ¹⁸O = 30.4 ± 0.6‰) ( Table 2 ).

Also, there were significant differences for δ¹³C (U = 1.000, Z = -2.205, p < 0.032) but not for δ¹⁸O (U = 8.000, Z = -0.490, p = 0.730) between the two age categories (mature and immature). Richer δ ¹³C values were found for mature animals (δ ¹³C = -8.7 ± 0.8‰) and lower values for immature animals (δ ¹³C = -11.7 ± 2.0‰). In the case of oxygen, the values of the two groups were also similar, but the richest values were found in immature animals (δ ¹⁸O = 30.6 ± 0.7‰) compared with mature animals (δ ¹⁸O = 30.3 ± 0.2‰).

The PCA based on the average values of each variable per dataset showed that the determinate values (0.047) had no collinearity among the variables, reinforcing the accuracy of our analysis.

According to the PCA, the first two components explained 81% of the variance ( Table 3 ). The two datasets are superimposed over the PCA biplot ( Figure 5 ). The dataset from 1982–1988 (blue circles and ellipse) is located over the right quadrants (superior and inferior), but most of the data are close to the X-axis, as the centromere indicates (the largest circle of the group). This dataset is associated with the Colorado River inflow and δ ¹³C vectors and has most values in the right quadrants.

The dataset from 1989–1993 (yellow ellipse and circles) is located in the left quadrants, and most of the data are close to the X-axis, as indicated by the centromere (the largest circle of the group). This dataset is related to the SST and SSS vectors and has most values in the left quadrants.

One data point from each set is associated with the δ ¹⁸O vector. These two data points are located outside the ellipses in the upper part of the quadrants.

The GLM results broadly support the PCA predictions. The best model produced by the GLM to explain the behavior of δ ¹³C was the Colorado River flow, with an intercept value of 1.17e^-14 and a p-value of 0.0209. The best model to explain the behavior of δ ¹⁸O was SSS, with an intercept value of 2.0e^-16 and a p-value of 0.593 ( Table 4 ).

The bone and tooth δ ¹³C and δ ¹⁸O values of marine mammals in different habitats, as obtained from the literature, were compared with those obtained from the vaquita in this study ( Table 5 ). The δ ¹³C values in coastal mammal species ranged from ~-9 to -11‰. The vaquita presented a mean value of δ ¹³C = -9.1 ± 1.3‰ at the end of the 1980s and a mean value of -10.8 ± 2.1‰ at the beginning of the 1990s. Both values are inside the range of the other coastal mammals. Extreme values were observed in otter species: Enhydra lutris, an otter species that inhabits kelp fields, showed a clear enrichment (-6.1 ± 0.9‰), and Lontra canadensis, an otter species of rivers, displayed an impoverishment (–17.3 ± 4.3‰) (Clementz and Koch, 2001; Drago et al., 2020).

In the case of δ ¹⁸O, the most enriched values were found in estuarine species (~27–28‰), followed by coastal species (~26–28‰), and the most impoverished values were detected in otter species that inhabit rivers (~23–26‰). The average values of δ ¹⁸O in vaquitas were the most enriched, with mean values of 30.5 ± 0.4‰ and 30.4 ± 0.6‰, respectively, for each period defined here To complement the analysis above, the vaquita isotopic values and the values of the coastal and estuarine species considered in Table 4 were compared in a plot ( Figure 6 ). It shows the distance between the vaquita δ ¹³C values and those of the other species in both periods (1.7‰), which was higher than that between species such as Phocoena phocoena (white square), Tursiops truncatus (black circle), and Globicephala macrorhynchus (black circle) (< 0.03‰). On the other hand, even though there was no difference in vaquita δ ¹⁸O values between the two periods, the distance from the values of the rest of the plotted species was > 2‰.

Considering that the vaquita, P. sinus, an endemic species of the UGC, is classified as critically endangered by the Mexican government and is on the IUCN Red List (NOM-059-SEMARNAT-2010, Diario Oficial de la Federación (2010), https://iucn-csg.org/wp-content/uploads/2023/06/Vaquita-Survey-2023-Main-Report.pdf), this study analyzed signals of the damming of the Colorado River across U.S. territory and the Morelos Dam in Mexico in the vaquita’s habitat, using satellite data and stable isotopes of carbon (δ ¹³C) and oxygen (δ ¹⁸O) from vaquita bones.

In the study period (from 1983 to 1993), in the context of the evolution of the UGC from an estuarine system to an inverse estuary (Álvarez-Borrego et al., 1975; Lavín et al., 1998; Lavín and Sánchez, 1999), the UGC was in the latter condition. However, a big influx of freshwater arrived in the UGC between ~1980 and 1988 because of abnormal snow melts in the Upper Colorado Basin, causing the dams to reach their maximum capacity (Lavín et al., 1998) and allowing us to compare two different conditions.

Based on the drastic change in the Colorado River’s flow toward the UGC, both periods were statistically identified using SSS and SST variation. These results were supported by SST/SSS diagrams showing two water masses inside a salinity gradient from 33.8 to 35.2 psu ( Figure 4 ). The water mass from the end of the 1980s had an average salinity of 34.4 psu, while the water mass detected at the beginning of the 1990s had an average salinity of 34.8 psu. Although the diagrams were built with satellite data, and without in situ data, a temporal change was observed in the water mass, demonstrating that variations in the Colorado River’s flow could have modified the vaquita’s habitat. According to Wright and Colling (1995), an average difference of 0.4 psu in salinity is a significative change in water mass properties.

Sometimes, marked hypersaline conditions develop in estuaries, exposing the primary producers to high osmotic stress (Kinne 1964), and the change in productivity can reduce an ecosystem’s diversity and alter trophic levels (Reist et al., 2006; Sora et al., 2024). Thus, when hypersalinity increases, the diversity decreases, sometimes resulting in mass mortalities and even the local extinction of flora and fauna (Molony and Parry, 2006; Kim et al., 2013; Tweedley et al., 2019). In the case of the UGC, it has been proven in various faunal groups such as fish, shrimp, and clams that the flow of the Colorado River modulates both their breeding areas and their life cycles, negatively impacting the species as salinity and temperature increase (Kowalewski et al., 2000; Rodriguez et al., 2001; Rowell et al., 2005; Rowell et al., 2008; Lozano-Montes et al., 2008; Galindo-Bect et al., 2021). Even a decrease in biomass and diversity of fish larvae has been observed in the most saline larval habitats within the UGC (Sánchez-Velasco et al., 2012). In this context, we could assume that the vaquita has been affected by changes in environmental conditions caused by fluctuations in river flow, at least indirectly through a decrease in prey availability.

Carbon isotope values (δ ¹³C) significantly differed between periods, with a mean value of -9.1‰ at the end of the 1980s and a mean value of -10.8‰ at the beginning of the 1990s ( Table 2 ). These were interpreted as indicators of variable carbon sources between both periods. When the river’s flow is blocked, the nutrient source also decreases, giving rise to primary producers with depleted ¹³C values. This is likely because the low growth rate of phytoplankton allows greater discrimination of the heavy isotope as they do not reproduce fast enough to exhaust the ¹²C in the medium. In contrast, when the flow toward the delta region is active, organic matter would be expected to enter from the marshes, leading to higher concentrations of nutrients and giving rise to rapid phytoplanktonic reproduction in a short period of time (phytoplanktonic bloom). Due to its high reproductive rate and different carbon sources, the phytoplankton would be enriched in ¹³C (Fry and Sherr, 1989), leading to bioaccumulation along the trophic web (Miller et al., 2008). This agrees with the results obtained by the PCA ( Figure 5 ), in which a positive association was observed between the δ ¹³C isotopes and the flow of the Colorado River during the first period analyzed.

On the other hand, it has been observed that isotopic values ​​in primary producers in estuarine environments are influenced by saline concentration, with higher values ​​in areas of lower salinity (Fogel et al., 1992; Brutemark et al., 2009), as reported in this study ( Table 2 ). In addition to salinity differences, this is probably related to the entry of freshwater that carries organic matter enriched in ¹³C from primary producers with higher isotopic values, ​​such as macrophytes in the marsh area. This phenomenon has been observed in sedimentary organic carbon in nearshore environments with inputs of organic matter from macrophytes to the system and values ​​between -9.8‰ and -14.0‰ (Fry et al., 1977; France, 1995). Another explanation for the differences in the δ ¹³C in both periods may be a change in the available prey. Suppose the δ ¹³C isotopic values ​​show only a small difference between diet and consumer due to isotopic fractionation (~1‰ enrichment with each trophic step) (Clementz and Koch, 2001). In that case, it is reasonable to assume that during the period of greater input from the Colorado River, the vaquita had access to prey at a higher trophic level than during the period of no input from the Colorado River.

It is essential to mention that differences between age categories may affect δ¹³C isotopes since a difference of 3.0‰ was found between both categories. Previous studies have reported differences, although not significant, in the prey captured by adult animals (predominance of demersal and benthic organisms) and immature animals, which tend to feed on smaller pelagic prey (Findley et al., 1994). In addition to a greater niche breadth, a greater variability assessed by δ¹⁵N is observed in immature animals than in mature animals, suggesting that mature animals had a more selective feeding spectrum with probably higher trophic levels (Rodríguez-Pérez et al., 2021). Given that most mature specimens belonged to the 1980s, it is also essential to consider that the differences in isotopic signals between periods were possibly due to this factor.

The oxygen (δ ¹⁸O) values ​​were not significantly different between periods or age categories. This similarity may be due to the turnover rate of bone tissue, which stores information from long periods (Hedges et al., 2004), possibly making it impossible to separate the periods in the case of δ ¹⁸O.

The comparisons of marine mammals were made in different tissues. Those taken from Clementz and Koch (2001) came from dental enamel, which includes the isotopic composition of the mother’s diet and body water, plus any fractionations associated with diffusion across the placenta, nursing, and the rapid growth of the young animal, and bone, which also offers a long-term record of the life history of an animal and includes the last years or probably months in immature specimens. Despite the above, since both tissues come from apatite, the isotope fractionations for water to carbonate are roughly equal for the different materials (Clementz and Koch, 2001), so we consider the presented comparison viable.

Conversely, the δ ¹⁸O values ​​in vaquitas differed from those of other aquatic mammals ( Figure 6 ) in previous studies (Clementz and Koch, 2001; Drago et al., 2020). It is important to note that our study focused on coastal/estuarine or freshwater specimens ( Table 4 ). These high values are in line with Drago et al. (2020), who found that isotopic values increase in environments where aquatic mammals live in more saline conditions. As a result, freshwater or estuarine organisms record lower values ​​( Table 4 ) than oceanic organisms (~> 30.0‰; Drago et al., 2020). The δ ¹⁸O isotope was not associated with any of the environmental variables analyzed in this study, as was observed in the PCA.

Although it was not possible to demonstrate changes in the vaquita’s habitat use by δ ¹⁸O values, changes in the environment and in the trophic source of the primary producers were evident. This is significant because our hypothesis is that the vaquita’s enrichment, which is similar to that of other groups in the same habitat referred to as the VPDB, may result from the high evaporation rate and salinity in the UGC. This enrichment in clams and fish has been related to a population decline in species such as Mulinia coloradoensis, Cynoscion othonopterus, and Totoaba macdonaldi and has been attributed to high temperatures due to environmental change caused by the blocking of the flow of the Colorado River toward the UGC (Rodriguez et al., 2001; Dettman et al., 2004; Rowell et al., 2005, 2008). However, in the vaquita, an endothermic organism, temperature is not a primary determining factor (Clementz and Koch, 2001), as the impact on its diet may be due to an increase in salinity and temperature. Studies focusing on vaquita structures with longer temporal records, such as teeth and bullae, will help answer whether the damming of the Colorado River affects the vaquita’s habitat use.