The demographic decline of a sea lion population followed multi-decadal sea surface warming

The population growth of top predators depends largely on environmental conditions suitable for aggregating sufficient and high-quality prey. We reconstructed numerically the size of a resident population of California sea lions in the Gulf of California during 1978–2019 and its relation with multi-decadal sea surface temperature anomalies. This is the first multi-decadal examination of the sea surface temperature of the Gulf of California and of one of its major predators. A three-decade sustained warming explained the population’s trend accounting for 92% of the variance, including a 65% decline between 1991 and 2019. Long-term warming conditions started in the late 80s, followed by the population’s decline from 43,834 animals (range 34,080–58,274) in 1991 to only 15,291 (range 11,861–20,316) in 2019. The models suggested a century-scale optimum sea surface habitat occurring in mildly temperate waters, from 0.18 to 0.39 °C above the 100-year mean. The mechanistic links of this relation are still untested, but apparent diversification of pelagic fish catches suggests a reduction of high quality prey. We propose this population should be considered vulnerable to any disturbance that could add to the negative effects of the current warm sea surface conditions in the Gulf of California.

The physical structure of oceanic habitats often determines, through a multi-step process, the success of animal populations 1 . Sea surface temperature and its variability is widely known to affect oceanic top predators such as pinnipeds, through bottom-up mechanisms [2][3][4] . Abrupt or sustained changes of sea surface temperature affect the abundance and diversity of plankton communities 5 , pelagic fishes 6 , and ultimately marine mammals 7,8 . The latter typically respond with alterations in foraging habits 9,10 , key physiological processes 11,12 , reproductive success, or survival 4 . Although many marine mammal populations can withstand and recover from short-term sea surface warming conditions (e.g. El Niño) 13 , sustained positive environmental trends at a multi-decadal scale typically cause important shifts in the base of marine ecosystems 14 , and arguably the diversity of potential prey. These shifts could lead to the decline of some marine mammal populations as the conditions move away from the optimum habitat to which the species have adapted 15 .
In the Gulf of California (hereafter "the gulf "), there is a resident population of California sea lions (Zalophus californianus) 16 genetically isolated from the other populations of the Northeast Pacific Ocean 17 . Although a ~ 20% reduction of this population during the 1990s has been proposed based on partial counts spanning 1997-2004 18,19 , its total size and temporal trend, as well as the historical environmental context, remain unknown. A recent review mentioned that the population decreased 44% between 1979 and 2016, based on unpublished data 20 . Negative interaction with fisheries, and a temporal decline of the Pacific sardine (Sardinops sagax) in the early 1990s, as a result of unspecified environmental changes, have been proposed as potential causes of the apparent decline 19 . Nevertheless, interannual events such as El Niño do not have a consistent correlation with the dynamics or health indicators of the gulf 's population 21,22 . Given the uncertainty in both the apparent population decline and its drivers, we hypothesized that a long-term sea surface warming in the Gulf of California could be Scientific RepoRtS | (2020) 10:10499 | https://doi.org/10.1038/s41598-020-67534-0 www.nature.com/scientificreports/ related, since this environmental variable has the potential to trigger progressive changes in the base of marine ecosystems, arguably altering the composition and diversity of prey for top predators. To find support for this hypothesis, we made a numerical reconstruction of the California sea lion total population size and its changes during the last 42 years (1978-2019), using individual trends of the 13 reproductive colonies within the Gulf of California (Fig. 1), inferred from all available animal counts. Then, we explored the multi-decadal variability of the gulf 's sea surface temperature during the last 100 years (1920-2019), from which we were able to predict successfully the population's multi-decadal dynamics. Nevertheless, it is important to point out that the lack of information did not allow us to test for the mechanisms involved in such habitat-species , numbered from north to south. The green asterisks are the locations of the predicted sea surface temperature data used in this study. The map was created with the package "ggplot2" (https ://cran.r-proje ct.org/web/packa ges/ggplo t2/index .html) in R (https ://www.rproje ct.org/).  16,23 . Even though they are often grouped into four ecological regions 23,24 , or three genetic groups 17 , they function as individual units for the most part, because California sea lion individuals commonly exhibit natal philopatry to their colonies 22 , which also reflects distinctive foraging habits 25 . Our results confirmed a dramatic population decline in the last 28 years (1991-2019), after reaching its maximum in the early 1990s. The gulf 's annual sea surface temperature anomalies from the 100-year mean (1920-2019) were able to predict the population size dynamics at a multi-decadal level and suggested an optimum habitat at mildly temperate conditions.

Results
For all the 13 reproductive colonies, a second-order polynomial regression was the best model describing the abundance as a function of the year at a multi-decadal scale, over the simple linear trend. The mean dronebased perception bias correction factor added 41% (range [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48] to the original boat-based counts. The mean percentage of non-pups in the population was 80% (range 77-82), to which a 40.5% (range 4.3-93.2) was added to account for animals likely at sea (Table 1). Almost all reproductive colonies exhibited a decreasing trend during the last three decades (1990s-2010s), except for the southernmost small colony Los Islotes, which showed a sustained population growth in that period, with an apparent stabilization in the last 4 years (2015-2019) (Fig. 2). Some colonies showed decline even since 1978, but most of them, including the largest ones, showed a bell-shaped trend, with an increase during the 1980s followed by a sustained decrease since the early 1990s.
The most populated colonies predicted for 2019 were San Esteban with 4,959 animals (range 3,308-7,633), San Jorge with 2,236 animals (range 1,476-3,476), and Lobos with 1,175 (range 677-2,100). The smallest were El Rasito, Partido, and Consag, with abundances of less than 500 animals. The sum of the colonies' annual predictions was augmented by 6% (range 5.7-6.3), corresponding to the estimated proportion of animals at the 16 non-reproductive colonies during a typical breeding season (Table 1). According to this numerical reconstruction, the population had 35 Table 1. Statistical summary of the posterior distributions of the most relevant parameters estimated by the models. The names in JAGS language are those of the code in "Methods" section. A Gelman-Rubin statistic ( R ) close to 1 indicates good convergence of chains. A high percentage of the effective number of iterations (N eff ) indicates less uncertainty in the parameter's estimation.  (Fig. 4). It was especially steep during the 1990s, and less accelerated during the 2000s and the 2010s, reaching + 1.06 °C in 2019 (range + 0.85 to -1.28). The minimum temperatures occurred during the mid-1930s, reaching anomalies of − 0.57 °C (range − 0.77 to − 0.35). During 1991, when the California sea lion population size reached its maximum, the median SSTa was + 0.11 °C (range − 0.11 to + 0.31) (Fig. 4). According to our model, the population size of California sea lions in the Gulf of California can be predicted successfully from the SSTa at a multi-decadal scale (i.e. 25-year running means). After testing for several degrees of curve complexity, the best model describing this relationship was a mixture between a second-order polynomial and a parabola function. This was accomplished through a mixing parameter, whose median closer to 1 (0.77; range 0.06-0.99) indicated that the parabola function contributed with most of the fit ( Table 1). The proportion of the variance explained by this model (i.e. the Bayesian R-squared) was 0.918 (range 0.880-0.926), suggesting an extremely high predictability. The shape of the resulting curve (Fig. 5) implies that the population size is expected to decrease at both very high and very low SSTa. According to the predictions of this model (median and 95% CI), the range at which the population size would reach its maximum would be during anomalies from + 0.17 to + 0.40 °C (i.e. SST means of between 22.49 and 22.72 °C) (Fig. 5). This could be interpreted as the population's optimum physical habitat conditions. The model predictions also indicated that above this optimum, the population size would have reached half of its maximum at an anomaly of + 0.76 °C from the 100-year mean (i.e. around 23 °C).  (Table 1) confirmed its high predictability. Though with higher uncertainty for the last two decades (i.e. during warmer conditions), the model predicted successfully the median population sizes that resulted from the numerical reconstruction (Fig. 3). Moreover, since this ecological model was based on SST data at a higher temporal resolution than that of the animal counts, it was also capable of showing some details in the expected trend that the numerical reconstruction could not. The most important difference was that the numerical reconstruction tended steadily towards lower values, whereas the abundances predicted from the SSTa suggested some level of stabilization of the trend during the late 2000s and the 2010s, or at least, a less pronounced decrease, compared to that of the 1990s and early 2000s (Fig. 3). Annual medians (black dots) and 95%-credible intervals (CI; black error bars) of the numerical reconstruction are shown. Black inner ticks on the X axis mark the years for which any colony count was available. Medians of decadal changes and the 28-year decline in percentage are within gray and purple boxes, respectively. The population size prediction from the multi-decadal anomalies of sea surface temperature (SST) is shown as a green thick line (median) and green shaded area (95%-CI).

Discussion
The negative relation between the multi-decadal SST warming in the Gulf of California (1990s-2010s) and the California sea lion population trend was evident. This environmental condition predicted accurately the annual population size and trend, accounting for 91.8% of the variability (range 0.882-0.927). Such high explained variance could respond to the simple opposite trends of the SSTa means and of the reconstructed sea lion population size (based on the sum of the colonies' second-order polynomials), which represent the multi-decadal scale, and therefore such relation would be much variable at shorter interannual or multiannual scales. California sea lions reflect long-term changes in the marine environment 4,26 because, as marine mammals and top predators in general, they are conspicuous, with high energetic demands, and have a long lifespan 3 . Therefore, it is possible that other predators within the gulf have also reacted to this scale of environmental variability. Such effects at the multi-decadal scale cannot be interpreted as equivalent to those that have been observed in response to inter-annual warm events such as El Niño, which are well described 21,27,28 , and typically include a strong generalized depletion of prey availability and drastic changes in foraging habits 10,29,30 . Compared to these short-term shifts, the knowledge of the effects of long-term environmental changes is very limited [31][32][33] , yet it is essential to understand and predict population dynamics under climate change scenarios.
While the negative effects of El Niño events tend to last for only 2 or 3 years 21,29 , a multi-decadal warming would suppose a progressive, sustained, habitat change affecting multiple generations. It is challenging to understand the mechanisms involved in the response of a predator population to this scale of sea surface warming, to the point of reducing dramatically its size. One of the possible causes of a long-term decline in a spatially restrained pinniped population could be a diet change from high quality prey to what has been referred to as "junk food" 9,32 . When preferred prey becomes less available, predators must resort to a less nutritious diet. In the California Current Large Marine Ecosystem, California sea lions feed commonly on small pelagic fishes, including several species of sardine (Sardinops sagax), anchovy (Engraulis mordax), and mackerel, which are considered of high quality because of their calorie contents that range from 1.31 to 2.17 cal g -1 , and total fat contents from 0.048 to 0.124 g −1 , in comparison to other prey more common during warming conditions, such as the market squid (Doryteuthis opalescens; 0.92 cal g −1 ; 0.014 total fat g −1 ) or the shortbelly rockfish (Sebastes sp.; 0.94 cal g −1 ; 0.016 total fat g −1 ) 9,34 . In that region, this type of prey replacement has resulted in pups with lower body masses. Nevertheless, in a higher latitude environment off Alaska, a more diverse diet has been related to a less pronounced population decline for Steller sea lions (Eumetopias jubatus) 33 , which suggests that each predator can react differently to prey composition and availability depending on the type of oceanic habitat and the species' particular requirements.
An analogous prey replacement could have happened in the Gulf of California in 1978-2019, although this California sea lion population has different foraging habits than those of the California Current Large Marine Ecosystem, and few studies have described its diet 25,35 . It varies between colonies and consists on a large variety of species, including high-trophic level items such as the midshipman (Porichthys sp.) and squid (Leachia sp.), as well as low-trophic level species like the Pacific jack mackerel (Trachurus symmetricus), the Pacific sardine  Table 1 www.nature.com/scientificreports/ (Sardinops sagax), and the northern anchovy (Engraulis mordax). These prey were found in scats collected in the mid 1990s and early 2000s, so the composition of the diet in the late 1980s, before the population decline, is unknown. However, some trends in the fishery of small pelagic species could support the idea of a shift in the base of the trophic web. Sardine fisheries in the Gulf of California increased slowly until 1989; then, a massive collapse occurred. The catches plummeted from almost 300,000 tons to 7,000 tons in less than 3 years. Since then, several rises and falls have occurred, including the most extreme in 2008, when catches fell from over 500,000 tons to 3,000 tons in 5 years 36 . Nevertheless, as sardine fisheries plummeted and their trend became more unpredictable, other small pelagic fishes such as the thread herring (Opisthonema libertate), the Pacific mackerel (Scomber japonicus), the northern anchovy, and the bigmouth sardine (Cetengraulis mysticetus), gained relative importance in the catches 37 . This could indicate a diversification of the potential prey for California sea lions as the sea surface temperature increased. The effect of ocean surface warming on global fisheries is well documented, particularly in temperate regions, where warmer water species are being captured at higher latitudes 38,39 . Nevertheless, since the coverage of available data of potential prey is not comparable to the multi-decadal scale of this study, nor is spatially representative of the entire gulf, and given the high heterogeneity of fishing effort, it is not possible to validate this hypothetical ecological mechanism. If California sea lions have been foraging on more diverse but lower quality prey in the Gulf of California for several generations, it is expected that the main drivers of the population dynamics, fertility and/or survival rates, would have lowered. Although reductions in both parameters are common consequences of short-term warm anomalies 9,13,31 , their link to multi-decadal environmental shifts is unknown. Nevertheless, the mean proportion of pups during boat-based counts before the decline (1979-1991) was ~ 25%, and ~ 42% for adult females 18 . In 2016, we found a proportion of pups of only 15%, whereas the proportion of adult females was similar (~ 42%) ( Supplementary Information, Table S1), suggesting a reduction in fertility or pup survival that should be further explored. In contrast, in a colony of the species in the California Current Large Marine Ecosystem, with a positive population trend until 2014, the latest proportions of pups and adult females were ~ 41% and ~ 40%, respectively 40 .
The California sea lion population off the U.S. seems to have grown close to its estimated carrying capacity 8 . Although we do not have enough information to estimate directly that parameter for the gulf 's population, it could be close to the maximum population size reached in 1991, according to our model (Fig. 3). This could be addressed more accurately in the future by studying other external sources of pressure regulating fertility and survival besides the environmental conditions suitable for aggregating enough and high-quality prey. The fact that both populations present opposite demographic trends would reflect the differential effects that a multi-decadal ocean warming could pose given particular oceanographic dynamics and the availability of different types of prey. Given the genetic isolation between both populations 17 and the lack of connectivity between them, we do not expect their opposite trends are related. The exception to this could be the gulf 's southernmost small colony Los Islotes, which showed a positive multi-decadal trend (Fig. 2), apparently reaching its carrying capacity 22 , and for which some connectivity with colonies in the southern California Current Large Marine Ecosystem could occur 41,42 .
An increased foraging effort by juvenile and adult animals in warmer conditions should not be discarded as a complementary factor that could explain the decline of California sea lions on land. During warm anomalies, telemetry studies have shown that this species displays offshore foraging trips up to three times longer, presumably as a consequence of reduced prey availability 43,44 , resulting in fewer animals on land when counts are made. Further telemetry data are needed to build a dynamic correction factor of animals likely at sea. Increasing foraging effort can also result in decreased pup survival 13,30 by reducing lactation periods. Nevertheless, a progressive decrease in prey quality rather than quantity would not necessarily imply increase in adult females' foraging duration, but rather a reduction in the quality of their milk throughout several generations.
The 15.5% increase (range 13.7-17.4) of the population size during the 1980s (Fig. 5) was also successfully predicted by the SSTa. Since the shape of the curve of that model was unimodal, the population size would also be lower under extreme cold anomalies at the multi-decadal scale. The optimum sea surface temperatures revealed by this model could hypothetically benefit the occurrence of high quality prey in the Gulf of California, such as sardine, whose fishery also increased during the 1970s and early 1980s 36,45,46 . Though it could be argued that this California sea lion population increased during that period as a result of its recovery from documented legal hunting prior to the 1970s 20,47,48 , there is no evidence that the level of hunting was high enough to produce an appreciable decline of the population from which it would have recovered, at least not enough for a bottleneck effect 49 . Unfortunately, without more information about the hunting numbers, the possibility that the population growth was a co-occurrence of a recovery process and a progressively more suitable environment for the species is untestable at this point.
One particular finding of this study, and the exception to the general pattern, was the abundance trend at the colony Los Islotes, which was the only one that showed a steady increase during the study period, tending to stabilize in the last 5 years. It is the southernmost reproductive colony of the Gulf of California, located in La Paz Bay, a region characterized for its high biological production year round, attributed to a local mesoscale gyre phenomenon occurring during summer 50,51 , which would provide a more seasonally stable prey availability for the colony.
More studies focused on the foraging habits and reproductive health of the population, as well as studies on the abundance and diversity of potential prey, are needed to reach satisfactory conclusions on the mechanisms driving the population's long-term dynamics. It is also probable that other predators in the Gulf of California that depend on similar low-trophic-level prey would show similar population trends, but to our knowledge, this is the first study quantifying such a long-term temporal scale of environmental variability and biological response in the gulf. Based on our results, we propose that the California sea lion population of the Gulf of California should be considered as vulnerable to any disturbance that could add to the apparent negative effects of the current sea Scientific RepoRtS | (2020) 10:10499 | https://doi.org/10.1038/s41598-020-67534-0 www.nature.com/scientificreports/ surface warming conditions. Such conservation status should be maintained until the population size reaches at least half of the maximum observed during the early 1990s, which our models predict would occur at a maximum SST anomaly of + 0.76 °C from the 100-year mean.

Methods
All the analyses described below were based on algebraically explicit Bayesian regression models 52,53 , whose parameters of interest were estimated as samples from their posterior distributions through a Markov Chain Monte Carlo procedure, implemented in the language Just Another Gibbs Sampler (JAGS) 54 in R 55 . This approach avoids the loss of information through the propagation of uncertainty between connected estimations and models, accounting effectively for scarce and/or sparse data, if needed. We ran one million iterations in five independent chains, retaining every 20th value to avoid autocorrelation, and discarding the first 20% as a burn-in phase for each chain (see detailed algebraic explanations and JAGS codes below). The results of the estimations are reported throughout the text as their medians, accompanied by ranges that represent the 95%-credible intervals of the posterior distributions. All the data needed to run the models is provided fully arranged as three R's data list objects in the Supplementary Information File S1.
Population abundance and trend. When estimating abundance of pinnipeds, pup counts are typically used to infer those of other age classes from life charts that take into account life expectancy and fertility rates, among other parameters 8,56,57 . Unfortunately, only part of the available historical counts of California sea lions in the Gulf of California contains information on pups (Supplementary Information Table S1), and there are no life charts specifically estimated for each reproductive colony 20 ; therefore, such traditional estimates were not feasible. Instead, we focused on counts of all age/sex categories (i.e. total counts), as a consistent parameter that can be also used to estimate pinniped abundance 18,19,40 , especially if some perception and availability biases are accounted for 58,59 . For this, we used all counts available in the literature 16,18,19,21,59,60 and new, reported here for the first time.
The main challenge for estimating the population size of California sea lions and its multi-decadal trend in the Gulf of California was that there is only one breeding season, in 2016, for which there are available counts of animals in all of the 13 reproductive colonies, an effort of this study for which we implemented both visual and drone-based counts. The rest of the available counts spanning the 42-year time series (1978-2019), both published and described here for the first time, were made separately for different colonies during different breeding seasons, spanning 1978-2019. Therefore, it was impossible to estimate more than one total abundance based on the simple sum of counts from all colonies. Fortunately, although very sparse in time, these 86 counts followed the same protocols 16 18,19,59 , and therefore were comparable. This method consists in circumnavigating the colonies during morning hours, at 15-45 m from the shoreline, with hand-held binoculars, at a speed of 5-7 km h −1 . The counts included all animals detected on land, as well as those swimming near the surface between the boat and the colony.
To solve perception bias and to correct categorization errors, 16 drone-based counts were made parallel to the boat-based counts during the 2016 breeding season (Supplementary Information Table S1), taking aerial photographs of all the areas with presence of California sea lions (see detailed protocol in Adame et al. 59 ). The categorization of California sea lions followed established guidelines for identifying each individual as adult male, sub-adult male, adult female, juvenile, pup, or undetermined 16,61,62 . The counts reported for the first time in this study were made during 2016-2019 under research permits SGPA/DGVS/00050/16-19 issued by Dirección General de Vida Silvestre-Secretaría de Medio Ambiente y Recursos Naturales, with authorization of Comisión Nacional de Áreas Naturales Protegidas.
The first step of the analysis was fitting individual regressions of animal boat-based counts at each reproductive colony as functions of the year (1978-2019), only during breeding seasons. We tested for a simple linear trend and a second-order polynomial, choosing the best fit according to the lowest value of the Deviance Information Criterion (DIC) 63 . A Poisson likelihood was stated for all counts with a logarithmic link function (see model equations below). We used all published counts, and those reported for the first time in this study ( Supplementary  Information, Table S1). Given that most of the counts were very sparse in time during the 42-year study period (1978-2019), these curves were intended to identify multi-decadal variations only, rather than interannual or more fine-scale dynamics. However, the Bayesian framework of the analyses assured that the uncertainty associated with the estimations captured accurately the number of observations and their sparsity in time, which was the especial case of Consag, Lobos, and Partido.
The perception bias of boat-based counts was estimated for all available surveys as the proportion of animals detected from the boat with respect to those detected from the aerial drone photos at each reproductive colony. All proportions were stated as binomial likelihoods, whose inverse main parameters were used as additive correction factors for the annual abundance predictions. We also estimated an availability bias correction factor to add the proportion of animals likely to be at sea when the surveys were made. For that, we used a mean of those proportions reported in a study that compared on-land and at-sea aerial counts at colonies of the same species in the California Current Large Marine Ecosystem 58 . Unfortunately, there was no available information for our study area for this purpose. Since those proportions were with respect to all categories except pups, we had to estimate first the mean proportion of non-pup animals in the Gulf of California from 76 counts with available information ( Supplementary Information, Table S1). Then, the correction factor was applied only to non-pup animals for each reproductive colony. Although the proportion of animals foraging at sea can vary as a function of prey availability 43,64 , there was no available data or previous information that allowed us to address this dynamically. Therefore, our model assumed this proportion as constant. Since El Partido and El Rasito did not have enough observations during the first 10 years of the time series, and given that both colonies have very Scientific RepoRtS | (2020) 10:10499 | https://doi.org/10.1038/s41598-020-67534-0 www.nature.com/scientificreports/ limited area available for California sea lions, we set upper truncation limits for their predictions at the beginning of the time series to their maximum total counts available.
To obtain the annual posterior distributions of the total population size, we summed the annual predictions of abundance at each reproductive colony from 1978 to 2019, all within the same hierarchical Bayesian structure to propagate the uncertainties of the colonies' estimations 53 . We also added a final availability bias correction factor to account for the animals likely present at the 16 known non-reproductive colonies during a typical breeding season, based on the counts reported by the only study that included such colonies 18 . To better visualize the population size dynamics, we also estimated decadal and maximum percentages of change during the 42-year time series.
100-year anomalies of sea surface temperature. We estimated 25-year running means of sea surface temperature within the Gulf of California and their anomalies from the 100-year mean, from 1920 to 2019 (n = 1,193). This allowed us to filter interannual and decadal signals such as El Niño Southern Oscillation (ENSO) or the Pacific Decadal Oscillation (PDO), respectively, keeping only multi-decadal variability in accordance to the scale of the California sea lion population trend explored in this study. The dataset consisted of predictions from a two-stage reduced-space optimal interpolation procedure, and the superposition of quality-improved gridded sea surface temperature observations onto reconstructions, based on historical ice concentrations at the Earth's poles 65 . The product was developed by the Met Office Hadley Centre (https ://hadob s.metoffi ce.com/) for use in climate monitoring and modeling, and is freely distributed at a monthly-one-degree resolution by the Environmental Research Division's Data Access Program of the National Oceanic and Atmospheric Administration (https ://coast watch .pfeg.noaa.gov/erdda p/gridd ap/erdHa dISST .html), for which 14 SST prediction points corresponded to the Gulf of California (Fig. 1).

Habitat-based population trend.
We fitted a regression model of the 42 annual predictions of California sea lion population size as a function of SSTa at an annual basis, that is, only for the predicted values for July each year, the beginning of the breeding season when the counts were made. For this, we used the means and standard deviations of the posterior distributions of these two parameters, estimated by the two models described above. We tested for incremental polynomial degrees for the regression, but also for pairwise combinations of models through a mixing parameter 66 , which can add flexibility to the curve and increase the explained variability. The best fit was chosen based on the lowest DIC. The results of this function allowed us to calculate geometrically a range of optimum habitat, which was defined as the range of SSTa values at which the population size predictions reached the maximum at the lower and upper 95%-credible interval limits of the curve (Fig. 5). The same approximation was used to calculate the SSTa value at which the population size would reach half of its multidecadal maximum.
Algebraic and JAGS code details. For the numerical reconstruction of the California sea lion population size in the Gulf of California, all colony counts (C) available (i) were assumed to come from a Poisson likelihood with a logarithmic link, whose parameter (λ) was a function of the year (y) at each reproductive colony (j), which had fixed effects on the relation: For each reproductive colony (j), boat-based visual counts (V) in 2016 came from a Binomial likelihood, with a mean proportion of detected animals (P det ) respect to the total (T drone ) (i.e. drone-based). The complement proportion of the former corresponds to un-detected animals (P und ): (1) (2) V j = Binomial(P det, Tdrone j )  The only total count (TC) available encompassing all reproductive and non-reproductive colonies was stated as the second parameter of a binomial likelihood whose number of animals at the 16 non-reproductive colonies (N nrc ) represents a mean proportion P nrc : # Proportion of animals at non-reproductive colonies: n_count_non_rep ~ dbin(prop_other_col, tot_count_all_col) # Bionomial likelihhod prop_other_col~ dunif(0, 1) # Uniformative prior The annual (y) predictions of animal abundance (A) for the 13 reproductive colonies (r) summed were derived from the function in Eq. 1, adding the proportion of undetected animals (P und ) in visual surveys (i.e. drone-based correction) (Eq. 2): (3) Cnp k = Binomial(Pnp, T k ) (4) Pats l = N(µ ats , σ 2 ats ) (5) N nrc = Binomial(P nrc , TC) For each year (y), a number of adult animals (Ad) was estimated at each colony (r), using the mean proportion of non-pup animals (P np ) defined in Eq. 3: # Estimation of adult animals per colony: pred_adults[r,y] <-pred_count_col[r,y] * mu_prop_non_pups The latter estimate was used to include the correction factor for at-sea animals, not available during the counts, and to obtain the completely corrected annual estimates of abundance (CA) for each colony (r): # Add availability bias (non-pup animals at sea) pred_count_col_w_at_sea[r,y] <-pred_count_col[r,y] + (pred_adults[r,y] * mu_at_sea_prop) } The annual population sizes for the 42 years of the series (y) were estimated as the sum of those of the 13 reproductive colonies (i), plus the proportion of animals at non-reproductive colonies (Pnrc): # Predicted annual population size # as sum of the colonies' abundances: pred_sum_gulf[y] <-sum(pred_count_col_w_at_sea[,y]) * (1 + prop_other_col) # Add animals at other colonies } As secondary derived quantities, we estimated the percentage of decadal (n) changes in the population size and the main percentage decrease from the highest to the lowest: Ad y,r = A y,r · Pnp (8) CA y,r = A y,r + Ad y,r · µ ats   (2019) Finally, the ecological model stated that the estimated annual (y) populations sizes (N) followed a Normal likelihood, whose means were a function of the SST anomalies (SSTa) (only for July estimations). The function with the lowest DIC was a mixture between a parabola and a second-order polynomial with θ coefficients, through a mixing parameter (ω). The SSTa estimations were also stated to come from a Normal likelihood with a mean (µ SSTa ) and known standard deviations (σ SSTa ): N y ∼ N(µ N y , σ 2 y ) SSTa y ∼ N(µ SSTa , σ 2 SSTa y ) log(µ N y ) = (1 − ̟ ) · (θ 0 + θ 1 · SSTa y + θ 2 · SSTa 2 y ) + ̟ · θ 0 + θ 2 · SSSSTa 2 y Scientific RepoRtS | (2020) 10 The estimation of the variance explained by this model, the Bayesian R-squared, was based on a vector of errors (i.e. the observed minus the predicted): err_vec_ab[y] <-pred_ab[y] -mu_pred_ab[y] # (for B-R2 estimation) } # Priors: theta0 ~ dunif(min_a0, max_a0) theta1 ~ dunif(min_a1, max_a1) theta2 ~ dunif(min_a2, max_a2) omega ~ dunif(0, 1) mu_sst_anom ~ dunif(min_mu_sst_anom, max_mu_sst_anom) # Bayesian R-squared: r_squ_ecol <-pow(sd(mu_pred_ab), 2) / (pow(sd(mu_pred_ab), 2) + (pow(sd(err_vec_ab), 2))) }", fill=TRUE) sink() # end of JAGS model

Data availability
The authors assure the technical replication of all the analyses of our study by providing the original biological data (i.e. animal counts) in Supplementary Information Table S1, as well as a link within the Methods section granting access to the environmental data (i.e. the SST; https ://coast watch .pfeg.noaa.gov/erdda p/gridd ap/erdHa dISST .html). All the detailed code we built and the equations we used within the analyses, fully replicable in the R language, are given in the Methods section. Also, all the data necessary for running the three models is provided fully arranged as three R data list objects in Supplementary Information File S1. This saves all the programming needed to process and arrange the data.