Reconstructing population dynamics of a threatened marine mammal using multiple data sets

Models of marine mammal population dynamics have been used extensively to predict abundance. A less common application of these models is to reconstruct historical population dynamics, filling in gaps in observation data by integrating information from multiple sources. We developed an integrated population model for the Florida manatee (Trichechus manatus latirostris) to reconstruct its population dynamics in the southwest region of the state over the past 20 years. Our model improved precision of key parameter estimates and permitted inference on poorly known parameters. Population growth was slow (averaging 1.02; 95% credible interval 1.01–1.03) but not steady, and an unusual mortality event in 2013 led to an estimated net loss of 332 (217–466) manatees. Our analyses showed that precise estimates of abundance could be derived from estimates of vital rates and a few input estimates of abundance, which may mean costly surveys to estimate abundance don’t need to be conducted as frequently. Our study also shows that retrospective analyses can be useful to: (1) model the transient dynamics of age distribution; (2) assess and communicate the conservation status of wild populations; and (3) improve our understanding of environmental effects on population dynamics and thus enhance our ability to forecast.

). Point estimates of realized population growth rates were below 1 in only 3 years (2003, 2005, and 2013); in more than half the years the lower credible limit was also above 1. We estimated that the population declined in 2013 by 331 animals (217-459) (Fig. 5). By contrast, in an average year the population increased by an estimated 50 manatees . The simulation-based hindcast models also provided estimates of abundance, but with less precision than the IPM (Supplementary Note and Supplementary Figs. S1 and S2, online). The estimates of the number of dead manatees by coarse stage and year from the hindcast models were often inconsistent with the carcass recovery data ( Supplementary Fig. S3, online). Estimated age class structure varied over time, but the proportion of adults in the population was always more than 65% (Fig. 6). Cohorts can be seen progressing through the other stages over time; for example, a peak in relative abundance of first-year calves in 2007 translates into a peak in secondyear calves in 2008, third-year subadults in 2009, and fourth-year subadults in 2010.
Posterior estimates of adult survival probability were similar in most years to prior estimates ( Fig. 7) but more precise. The last 2 years of the analysis (2014 and 2015) stand out as having low estimates of survival in the prior but close to average estimates of survival in the posterior. Posterior estimates of average calf and subadult survival probability were less than prior estimates and showed a pattern of increase with age (Fig. 8). Except for first-year calves, however, the posterior estimates were not noticeably more precise than prior estimates. Posterior estimates of adult female reproductive probability were mostly very similar to the prior estimates, with little gain in precision ( Supplementary Fig. S4, online). The low posterior estimate of reproductive probability for fourthyear subadults (0.00083; 0.000-0.108) was similar to the prior estimate (0.00100; 0.000-0.285).
Estimated average carcass recovery rates were high for subadults (0.948; 0.870-0.982) and adults (0.971; 0.899-0.995) but lower for calves (0.670; 0.477-0.844; Supplementary Fig. S5, online). Estimated recovery  11,23 Stage-specific survival (s) and reproductive (γ) probabilities govern the transitions between stages. Females (♀) and males (♂) are both modeled. Calves enter the population model at age 1.5; until that age they are tracked with their mother (0.5-year-old calves are shown in the diagram for completeness). Color indicates the primary sources of information for parameters within the integrated population model: the parameters in green (s 1 -s 4 ) are informed primarily by carcass data and surrogate mortality ratio estimates; the parameters in orange (s p and s a ) are informed primarily by carcass data and mark-recapture survival estimates; and the parameters in purple (γ 4 , γ p , and γ b ) are informed primarily by mark-recapture reproductive estimates. Other sources of information, including abundance estimates, informed all parameters.  Fig. S7, online) and in age class structure for 1997-2001 but little variation in total abundance for any year or in age class structure for 2002-2016 (Supplementary Figs. S8 and S9, online). Sensitivity analysis for different starting abundances showed little variation in total abundance for any year and no variation in stage structure or age class structure (Supplementary Figs. S10-S12, online).

Discussion
In this study we were able to obtain estimates of parameters that had been missing for the southwest subpopulation, including survival probabilities of younger stages of manatees, recovery rates of manatee carcasses, and abundance in years before and between abundance surveys.
Survival probabilities of younger animals are key parameters in population viability analyses of Florida manatees 11,23,25 . But these probabilities have long been extrapolated from one study of manatees in a small management unit on Florida's east coast 26 . The average probabilities of juvenile survival estimated here are lower than those obtained from that extrapolation (Fig. 8). Independent estimates of the younger manatee survival probabilities for the southwest management unit will soon be available from genetic mark-recapture-recovery modeling, but such data are not forthcoming for the other three Florida manatee management units, making the approach used here for estimating these probabilities more readily applicable.
In addition, our model provided estimates of the effects of red tide and cold events on the population. The red tide event of 2013, during which 353 carcasses were recovered in the southwest (of which at least 268 were killed by red tide), contributed to an estimated net drop in the population of 331 (217-459) manatees (Fig. 5) for an annual population growth rate of 0.89 (0.85-0.93; Fig. 4). Our results support the finding that such red tide events (classified as intense) affect calves particularly ( Supplementary Fig. S13, online) 11 . In contrast, the cold event of 2010, which led to 247 recovered carcasses in the southwest region, did not appear to lead to a net drop in population, according to our model. This may be in part because our prior estimate of adult survival that year was relatively high (Fig. 7), and the model assumes (and estimates) a fixed ratio between age-class survival rates across years (Supplementary Table S1, online). These new estimates can be helpful in communicating the impact these disturbance events had on the population. Unusual mortality events that lead to high carcass counts often attract a lot of attention from the press and the public. The IPM provides a way to put such mortality events in perspective and to answer questions such as "What was the impact of a particular mortality event on the population?" In addition, the average population growth rate (1.02, 1.01-1.03) estimated from our data supports the hypothesis that the manatee population was increasing from 1997 to 2016 (Figs. 3 and 4). This is the first rigorous estimate of historical (realized) population growth rate for this population. This information is complementary to and consistent with the projected population growth rate obtained from the CBM projections 11 . www.nature.com/scientificreports/ Our model also provided more precise estimates of many parameters estimated earlier, such as adult survival and abundance for years in which abundance surveys were carried out (Figs. 3 and 7). In some cases, our approach may reduce bias, although it is also possible for IPMs to introduce or increase bias 27 . Possible biases in some input estimates to our model, such as abundance 28,29 and end-of-time-series survival probabilities 30 , have been noted [28][29][30] . In some cases, the median estimates obtained from the IPM were substantially different from the original estimates (compare prior abundance survey and posterior estimates in Figs. 3 and 7). The IPM might correct for biases in abundance and end-of-time-series survival estimates, although this idea needs to be further www.nature.com/scientificreports/ evaluated. Because it includes a recovery model for carcass data, the IPM does not hindcast impossible numbers of deaths, unlike the simulation-based hindcast model ( Supplementary Fig. S3, online). The IPM results suggest that these results from the simulation-based hindcast model were off both because the 2011 abundance estimate input was too low and because the survival estimate inputs for juveniles (s 1 -s 4 ) were too high. By integrating multiple sources of information, we are synthesizing the best available information but also hedging our bets by not relying on just one source of data in estimating critical demographic parameters. Many of our posterior estimates are consistent with other published results for Florida manatees. Our estimates of realized population growth rates (Fig. 4) are similar to the projected population growth estimates from the CBM and consistent with general trends of growth in synoptic and carcass counts. Our estimates of age structure (Fig. 6), although variable over time, are consistent with the asymptotic stable age structure that projecting from a simple matrix model would provide. Our estimates of the mortality effects of the 2013 red tide ( Supplementary Fig. S13, online) are similar to those from the CBM. The pattern of our estimated recovery probabilities by coarse stage (Supplementary Fig. S5, online) is consistent with an earlier estimate of age-specific recovery rates relative to (unknown) adult recovery probability 31 , although our estimates of subadult and adult recovery probabilities are closer to 1 than we expected. The high estimates of recovery probability may be due to the IPM attempting to harmonize partially incompatible model components ( Supplementary Fig. S3, online). When model components generate incompatible results, either due to model misspecification or not referencing exactly the same populations, an IPM must reconcile those results. This reconciliation can generate bias in some estimates, although the generally higher precision of IPM estimates may still mean higher accuracy. Ground truthing or other research may be needed to determine whether FWC is actually recovering such a high proportion of manatee carcasses.
The results of this study are relevant to the management of Florida manatee populations. The manatee recovery plan used by the USFWS under the Endangered Species Act relies on several metrics that can be obtained from the IPM, such as realized population growth rates and population size. The IPM provides one of the most rigorous assessments to date for these quantities and may be used by natural resource managers in assessing the status of the manatee population. It can also be used to update key model parameters of the CBM, which at present is the primary population assessment tool for managers.
Another important regulatory framework relevant to marine mammal conservation in the United States is the Marine Mammal Protection Act. Here again, an IPM can help in addressing some of the act's requirements. Indeed, the act specifies a formula for computing potential biological removal (PBR; the maximum number of animals that can be removed from a stock while allowing it to reach or remain at its optimum sustainable population) [32][33][34][35] where N min is the minimum population abundance estimate (20th percentile of abundance estimate distribution), R max is the theoretical maximum rate of increase for the stock, F r is a recovery factor (generally 0.5 for threatened species, but see Moore et al. 35 ), and N is the point estimate of population abundance. Based on our estimate from the last year of the analysis (2016), N min for the southwest population of Florida manatees is about 2780. This is lower than N min would be based on the abundance survey (prior) estimate from the same year (about 3140); CV N from the IPM posterior was lower than from the prior ( Supplementary Fig. S2, online) but N was as well (Fig. 3). Estimation of R max requires extrapolating growth rates to conditions of low population density and absence of anthropogenic mortality; our IPM is not designed for that purpose, but future extensions could be developed to address this need. A merging of our IPM, or other matrix model approach, with an allometric approach to estimating R max would allow a more accurate estimate of this parameter 36 . Both matrix model (individual population) and allometric (cross population) approaches to estimating R max are strongly affected by biases caused by using empirical estimates of adult survival instead of what adult survival would be under ideal conditions; however, these biases run in opposite directions, so an integration of these approaches greatly reduces any bias in R max 36 . Another benefit of the IPM is its usefulness for planning monitoring activities, including how to allocate resources to various aspects of the monitoring program, such as aerial surveys, photo-identification, genetic sampling, and carcass recovery. Various sampling scenarios (e.g., 40% of carcasses recovered; 200 genetic samples per year; one aerial survey every 5 years) can be combined with simulated data generated under those scenarios to see how the accuracy of model parameter estimates differs among scenarios. Trade-offs between parameter accuracy/precision and budget allocation can then be examined to improve monitoring efficiency. Optimizing the sampling with an IPM also makes sense in the context of targeted monitoring for adaptive management 37 . In such applications, the IPM can be used to estimate state variables (e.g., abundance) that keep track of system changes, allow managers to implement state-dependent decisions, and update beliefs about which model is the best approximation of reality (through Bayes theorem) 37,38 . A now classic example of an implementation of this adaptive management process is for the sustainable harvesting of waterfowl in North America 37 , where the optimal state-dependent harvest policies are driven, at least partially, by waterfowl abundance. IPMs are now being used to increase precision of abundance and other state variables in adaptive management of waterfowl 39,40 .
(1) www.nature.com/scientificreports/ A monitoring component that could be streamlined is the carcass-recovery and necropsy program. The present protocol is that almost all carcasses reported must be recovered and necropsied, which, along with the growth in the manatee population, is making this program increasingly labor-intensive and expensive. The IPM gives us the first true estimates of carcass recovery probabilities for Florida manatees. These estimates are now being used by FWC in evaluating and improving the efficiency of these programs.
Monitoring populations of marine mammals involves special challenges, such as the difficulty, cost, and risk to researchers involved in counting the population, often through aerial surveys. Several other studies that involved the development of IPM for marine mammals 16,[41][42][43] had at least one thing in common with ours: population surveys were not conducted every year, which differs from most IPMs used for terrestrial birds and mammals. Our approach, like those applied to other marine mammals, could be valuable for filling in abundance estimates for other sirenians and small cetaceans, where estimating survival and reproductive probabilities from mark-recapture data is often easier than obtaining abundance estimates. As explained earlier, the IPM can then be used to determine the optimal frequency of surveys and optimal spatial sampling effort (e.g., how much area to survey and how many survey visits at each location to estimate detection) 28 .
Studies of other marine mammals 16,41-43 collected explicit data on age or stage structure, while for manatees, reliable data were not available for these parameters. We were able to estimate age class structure for the years 2002-2016 using neither stage structure data nor particularly informed priors (Fig. 6). This is likely because of the weak ergodic theorem of demography, which shows that the initial stage structure becomes less relevant with more years of known (or, in our case, estimated) survival and reproductive probabilities 3,44 . Our approach may be useful for other marine species without reliable stage structure information. Modeling stage structure and transient dynamics can be important to improving understanding of the dynamics of wild populations and can have important management implications. For instance, Johnson et al. 45 found that the initial stage structure could have substantial policy consequences for the management of an invasive species.
Our IPM and the associated input models are based on a series of assumptions (Supplementary Table S1, online). One of the assumptions of the IPM is the independence of the data sources for the input analyses. This assumption is violated in our case; the adult survival analysis shares carcass data with the recovery analysis and mark-recapture data with the reproductive analysis. Two simulation studies 17,46 found that violating this assumption had little effect, but as their analyses were not identical to ours, this assumption violation still might diminish the accuracy of our estimates. Simulations by Rieke et al. 47 show that assumption violations in one of the model components can dramatically reduce the accuracy of estimates of latent parameters. Therefore, in our case, the estimates of juvenile survival, recovery probabilities, and abundance in years without abundance surveys should all be interpreted cautiously.
There are several possible extensions of this model, for example for use in the other three Florida manatee management units (Fig. 1). Because we are uncertain about winter within-coast manatee distribution 29 , two coast-wide IPMs that each jointly model the two management units on that coast might be most appropriate. With an initial abundance distribution and yearly vital rate estimates for each management unit (possibly including movement rates between regions, if they become available), subsequent coast-wide abundance estimates could be shared between them. This would allow relaxation of the assumption that the proportion of the winter population in each of the two management units remains fixed over time.
Possible extensions could demonstrate whether and to what extent the IPM decreases bias in input estimates, through simulating estimates with known biases and carcass data, running the IPM with the simulated data, and repeating this process many times. One could similarly test the model's robustness to different assumption violations.
Preliminary analyses suggest that our use of earlier analyses as priors in the integrated model does not bias results but that it might reduce precision. Therefore, it may be useful to estimate more parameters from data directly within a future version of this IPM. In addition, incorporating additional data sources (such as genetic mark-recapture and age estimates using tympanoperiotic ear bones) could improve parameter estimation. Since each parameter can have only one prior, this too requires performing more of the data analysis within the IPM.
Despite these limitations, we believe that this manatee IPM is the most rigorous means of retrospective assessment of the population dynamics of the Florida manatee. Because the model is modular (e.g., abundance module, survival module), as each module is improved, the model as a whole is improved. This offers a compelling framework within which to synthesize and update information about population dynamics. We have shown here that an IPM can be used: (1) to infer historical trends in abundance, improving our understanding of population dynamics and therefore our ability to forecast; (2) to model the transient dynamics of stage distribution, which can be important to some populations; (3) to assess the conservation status of wild populations and to communicate that information to stakeholders (e.g., we can now quantify the impact of the 2013 red tide event on the manatee population); and (4) to improve allocation of effort in complex monitoring programs.
Our modeling frameworks are relevant to population status assessment protocols for management and conservation, such as recovery plans under the Endangered Species Act and potential biological removal under the Marine Mammal Protection Act. Other marine mammal conservation programs, such as that of the Hawaiian monk seal, also have complex monitoring components 48 . We hope that our ideas can inform other programs that focus on the conservation of marine mammals.

Methods
Study population. The Florida manatee is a highly mobile, long-lived marine mammal. Its population is divided into four management units or regions 49 , based on where individuals aggregate in the winter: northwest, southwest, Atlantic coast, and Upper St. Johns River (Fig. 1). We focus here on the southwest management unit.  11,23 . We considered a population with 10 stages, assuming a December census averaging about half a year after birth, labeled for the subsequent year. Like Runge et al. we used a two-sex model (Fig. 2). This population model can be expressed as a series of stochastic binomial draws and deterministic difference equations: where N corresponds to the number of manatees in a stage; s indicates survival probability; subscripts f and m indicate sex; and subscripts 1, 2, 3, 4, p, c, b, and a denote first, second, third, fourth, prebreeder, mothers-withcalf, breeder, and adult stage classes, respectively. S and G correspond to results of binomial draws, numbers of manatees in a stage that survived or bred, respectively, in a time step, tracked to model the demographic stochasticity processes of survival and reproduction: where γ represents reproductive probability. Because first-year calves are not tracked separately in the model but are included in the mothers-with-calf stage, the total abundance in year t can be calculated as: We use the term stage to refer to one of the 10 stages in Eq. (4). We use the term coarse stage to refer to calves (first-and second-year manatees), subadults (third-and fourth-year manatees), or adults (all other stages), without regard to sex. We use the term age class to divide manatees into first-year, second-year, third-year, fourth-year, and adult stages, without regard to sex.
We developed several simple simulation-based hindcast models as our initial approach to applying this population model retrospectively (Supplementary Methods online).
Input data and estimates. Before 2011, the only statewide manatee counts were done using aerial synoptic surveys, carried out most years during the coldest part of winter when manatees tend to aggregate at warmwater sites 28,50,51 , starting in 1991. The synoptic surveys are each flown over a short time (one to a few days), generally with a single observer per airplane. Although these surveys are intended to be comprehensive, they do not account for the number of manatees that are missed because they are absent from surveyed sites, present but underwater or otherwise not available to be detected by the observer, or available yet not detected by the observer. In the retrospective analyses, we used synoptic survey results from southwest Florida as counts that represent the lower bounds for abundance.
In 2011, a new aerial survey method was implemented that accounts for manatee presence at survey locations, availability, and detection, using a stratified random plot design, independent estimates of manatee availability, and a double observer protocol, respectively 28 . The population abundance estimate for southwest Florida in 2011 was obtained from Hostetler et al. 29 , which was an update of the estimate made by Martin et al. 28 We approximated this estimate and its uncertainty using a lognormal distribution, with mean 7.72 and SD 0.166 on the log scale. In estimating the initial 2016 abundance (another survey was conducted in early December 2015), Hostetler et al. questioned whether the estimated distribution of manatees between southwest and northwest Florida in early December 2015 was representative of the midwinter distribution and concluded it likely was (3)  29 . Therefore, instead of using the estimate for southwest Florida for December 2015, we multiplied the estimate for the entire west coast of Florida by an estimated proportion of those animals in the southwest, obtained from the synoptic surveys, including uncertainty reflected as temporal variability in synoptic proportions. We approximated this product and its uncertainty with a lognormal distribution with mean 8.16 and SD 0.132 on the log scale. Annual adult survival probability estimates were obtained from the Barker robust design mark-recapture model 52 . The data consisted of photo-documented sightings of live marked individuals under a traditional sampling framework at primary sampling locations during winter (Fig. 1), and identification of recovered carcasses. Substantial temporary emigration of individuals from study areas can bias estimates of survival rates for longlived species, particularly at the end of the time series. The Barker robust design models these types of emigration, reducing bias and increasing precision 30 .
Annual estimates of reproductive probability were obtained using a multistate robust design mark-recapture model 53 . The data for this model were photo-documented sightings of marked adult females and associated data indicative of calving. The reproductive state of individual female manatees can be uncertain even when the female is detected, resulting in biased estimates if one assumes they did not calf that year. This multistate model adjusts for reproductive state using hidden Markov processes, both increasing precision and reducing bias 54,55 . We used annual estimates from 1997 to 2015 for both survival and reproductive probabilities, using year as a random effect (method-of-moments variance component) on the sin-link coefficients 56 . Maximum likelihood shrinkage estimates, estimates of variance, and sampling covariances for these coefficients were obtained from program MARK.
The only direct estimates of Florida manatee survival rates for calves and subadults (s 1 -s 4 ) come from the Upper St. Johns River region and the years 1979/80-2000/01 26 . That study found no evidence of a difference between survival rates of subadults (s 3 and s 4 ) and adults (s p and s a ), but lower survival for each of the calf stages (s 1 and s 2 ). Following the CBM 11 , we constructed mortality ratios between adults and other age classes within the IPM, using as prior distributions the survival estimates (with associated uncertainty) of that analysis. The IPM assumes that those mortality ratios are the same in all years.
We used the estimate of the probability of reproduction for fourth-year females (γ 4 ) from Runge et al. 23 The CBM assumes no temporal variation in γ 4 ; we maintained that assumption. We diverged from the CBM, which assumes that γ p does not vary over time and is slightly less than the average γ a [t], by assuming that γ p [t] = γ a [t].
Carcass counts for southwest Florida in each year were obtained from FWC's manatee mortality database (http:// myfwc. com/ resea rch/ manat ee/ rescue-morta lity-respo nse/ morta lity-stati stics/). Carcasses were divided into coarse stage by carcass length (calves: 151-235 cm; subadults: 236-265 cm; and adults: 266 cm or more). Carcasses less than 151 cm in length (n = 1391) were defined as perinatal and were excluded from the analysis. Carcasses of unknown length (n = 79; mostly too decomposed or verified but not recovered) were also excluded from the analysis.

IPM.
The general methods of integrated population modeling have been well described 13,14 . Our IPM differs from most others in that we use previous estimates of abundance, adult survival probabilities, and female reproductive probabilities as prior distributions, an accepted alternative to directly integrating those analyses into the model [57][58][59] .
To provide a prior probability distribution on initial (1997) abundance, we ran several simulation-based hindcast models, some working backward from the 2011 abundance estimate and others from the 2016 abundance estimate (Supplementary Methods, online). We combined model runs across models and fit a lognormal distribution to all the 1997 abundances generated. Preliminary testing suggested that doubling the SD of the log from this estimate provided a better measure of our true uncertainty about initial abundance: mean 7.38 and SD 0.49 on the log scale.
In 1996 the manatee population in southwest Florida suffered at least two stochastic shocks: a cold winter and a red tide event classified as intense 11,60 . For this reason, we did not think that the stable stage distribution would necessarily be a good approximation of the true stage distribution in 1997. The latest CBM report provides estimates of the effects of these types of events on manatee mortality, by coarse stage 11 . We used these estimates to shift the initial stage distribution from a stable distribution in one of the simulation-based hindcast models (Supplementary Methods, online). We used the estimates of 1997 stage structure from that model, modeled with a Dirichlet distribution, and with variance quadrupled to allow for additional uncertainty (the population may not have reached the stable distribution at the beginning of 1996 either). We tested the sensitivity of abundance, stage structure, and age structure to uncertainty in the initial abundance and stage structure (Supplementary Methods, online).
A type of data we directly modeled in the IPM was recovered carcasses. We used a binomial model: where C s [t] is the number of carcasses recovered in year t of coarse stage s, M s [t] is the number dead in that year and stage, and r s [t] is the associated recovery probability. We tested several models for r and used binary inclusion factors 61,62 and parameter credible intervals to settle on additive effects of coarse stage and year, with year as a random effect.

Implementation.
We implemented the IPM using the R (version 3.5.1) package NIMBLE (version 0.8.0), an MCMC implementation using the BUGS language 63,64 . We ran each version of the model for three chains of 6 million iterations after a burn-in of 2 million each and thinned the results by 50. We tested for convergence using Gelman-Rubin statistic and visual examination of the chains 65 . All models presented converged successfully. www.nature.com/scientificreports/ To minimize mean absolute deviation, we used the medians of the posterior distributions as the point estimates for all parameters (except for binary inclusion factors, for which we used means) 66 . We present the 0.025 and 0.975 quantiles of the posterior distributions as the 95% Bayesian credible intervals for each parameter.
We used this IPM to estimate several parameters: total year-specific abundance for 1997-2016; realized yearspecific population growth rates for 1997-2015 ( [t] = N tot [t+1] N tot [t] ); the geometric mean population growth rate; the net population change in an average year and in 2013; annual year-specific age class structure for 2002-2016; age-class-specific and, in many cases, year-specific survival and reproductive probabilities for 1997-2015; and coarse stage-specific and year-specific carcass recovery probabilities for 1997-2015.