Estimating Pacific walrus abundance and survival with multievent mark−recapture models

: Arctic marine ecosystems are undergoing rapid physical and biological change associated with climate warming and loss of sea ice. Sea ice loss will impact many species through altered spatial and temporal availability of resources. In the Bering and Chukchi Seas, the Pacific walrus Odobenus rosmarus divergens is one species that could be impacted by rapid environmental change, and thus, population assessments are needed to monitor changes in the status of this ecologically and culturally important marine mammal. We conducted a 5 yr genetic mark−recap-ture study to estimate demographic parameters for the Pacific walrus. We developed a Bayesian multievent mark−recapture model to estimate walrus survival and abundance while accounting for age misclassification. We estimated the probability of juvenile annual survival as 0.63 (95% credible interval [CrI]: 0.39−0.87) and adult female annual survival as 0.90 (95% CrI: 0.74−1.00). We estimated total abundance as 257193 (95% CrI: 171138−366366). We provide the first estimate of total Pacific walrus abundance since an aerial survey in 2006, which generated a substantially less precise total population size estimate (129 000; 95% CI: 55000−507000). The emerging ecosystem state in the northern Bering and Chukchi Seas will likely result in a decline in Pacific walrus abundance, but there is substantial uncertainty regarding the magnitude of the anticipated decline. Our demographic estimates provide critical information to evaluate future population trends of this subsistence resource vital to communities that border the Bering and Chukchi Seas in the USA and Russia.


INTRODUCTION
In the Pacific Arctic marine environment, sea ice loss, altered currents, and increased sea surface temperatures have precipitated a major shift in ecosystem state with direct impacts on biological communi-ties (Baker et al. 2020, Huntington et al. 2020). In the northern Bering and Chukchi Seas, pelagic−benthic coupling combined with other physical properties of the Pacific Arctic has historically generated persistent areas of elevated benthic biomass (Grebmeier et al. 1989(Grebmeier et al. , 2015. Currently, these benthic hotspots re -present a predictable food resource for upper trophic level predators with locally high abundances of seabirds and marine mammals (Kuletz et al. 2015). However, sea ice loss will likely reduce access to these benthic hotspots for marine mammal species that depend on ice as a platform (e.g. bearded seal Erignathus barbatus; Breed et al. 2018). Furthermore, the emerging physical state of the ecosystem (Baker et al. 2020) could cascade to alter the spatial distribution and seasonal availability of these hotspots in the future (Grebmeier et al. 2015). Thus, the pending ecosystem shift in the Pacific Arctic environment will have indirect and direct consequences on space use for many species. Species also could respond to an altered Pacific Arctic marine environment at a population level via reduced reproductive rates, decreased survival, and reduced overall abundance.
The Pacific walrus Odobenus rosmarus divergens is a sea ice-associated pinniped that inhabits continental shelf waters in the Bering and Chukchi Seas. Pacific walruses (hereafter, walruses) primarily feed on benthic invertebrates, which are acquired via short dives to the ocean floor (Fay 1982, Sheffield & Grebmeier 2009). Walruses haul out to rest on sea ice in singlets or groups between foraging bouts when sea ice is available (Beatty et al. 2016), but they also rest on land at coastal haulouts. Walruses breed in the winter in the Bering Sea, and as sea ice melts in spring, adult females, juveniles, and some adult males follow the retreating sea ice north to summer in the productive offshore waters of the Chukchi Sea (Fay 1982) (Fig. 1). Most adult males remain in the Bering Sea through summer, resting at coastal haulouts between foraging bouts (Fay 1982). As sea ice continues to retreat north into the Arctic basin, walruses in the Chukchi Sea move to coastal haulouts in the USA and Russian Federation (Russia) in the late summer and autumn (Udevitz et al. 2013). When sea ice re-forms in late autumn and early winter, walruses migrate back to the Bering Sea (Fay 1982).
Sea ice loss will alter the spatial and temporal distribution of available foraging habitat for walruses (Beatty et al. 2016, MacCracken et al. 2017. Sea ice loss is predicted to result in increased use of coastal haulouts, and individuals that travel from coastal haulouts to access productive offshore foraging areas could face increased energetic costs . Alternatively, walruses that feed near coastal haulouts may forage in nearshore areas that are less productive than offshore areas. Walruses may exhibit seasonal reductions in body mass in response to a changed environment, which could reduce reproductive rates, survival, and abundance ). Consequently, estimates of current population abundance and vital rates could provide valuable information to quantify the impact that emerging and future ecosystem changes may have on this species.
Walrus abundance is difficult to estimate because the species is heterogeneously distributed over a vast and remote area (Speckman et al. 2011). Aerial surveys to estimate walrus abundance were conducted in 1975 (Estes & Gol'tsev 1984), 1976 (Krogman et al. 1979), 1980(Johnson et al. 1982), 1985 (Gilbert 1989), and 1990 (Gilbert et al. 1992). Estimates from these surveys ranged from 140 000 in 1975 (Estes & Gol' tsev 1984) to 328 000 in 1976 (Krogman et al. 1979). However, none of the original estimates accounted for walruses that were not available to be counted because they were underwater during the survey (Fay et al. 1997). In 2006, a survey conducted jointly by the USA and Russia corrected for availability and accounted for uncertainty, yielding a population size estimate of 129 000 walruses (Speckman et al. 2011). However, the high spatial and temporal variability in the distribution of walruses generated an estimate with relatively low precision (95% CI: 55 000−507 000; Speckman et al. 2011), which warranted a new approach to estimate walrus abundance.
In 2013, we initiated a 5 yr genetic mark−recapture study to estimate walrus abundance and annual survival. Specifically, our objectives were to (1) estimate walrus abundance and age-specific survival and (2) evaluate the efficacy of a shorter (i.e. 3 yr) genetic mark−recapture study on this species to reduce costs associated with field work. We developed a Bayesian multievent mark−recapture model to estimate survival and (re)capture probabilities and used a ratio estimator to estimate walrus abundance. We specifically used a multievent model because these models are generalizations of multistate models that can estimate age-dependent recapture probabilities and survival while accounting for the natural aging process and age (i.e. state) misclassification. We then conducted a simulation exercise where we compared the bias and precision of a 5 yr study to a hypothetical 3 yr study. We consulted with the Eskimo Walrus Commission before we initiated the study to minimize any effects to subsistence hunting or other community activities in the region. Alaska Native walrus hunters participated in research cruises as team members and shared traditional ecological knowledge of walrus behavior and sea ice dynamics, which was applied to the field methodology at great benefit to the study. We also provided annual updates on the project to the Eskimo Walrus Commission to inform them of progress and for feedback.

Sample collection, laboratory methods, and sample matching
We conducted a research cruise in June of each year from 2013 to 2017 to collect skin biopsy samples from live walruses hauled out on sea ice in the Bering and Chukchi Seas. We focused sampling on groups of adult females and juveniles, as these classes are demographically important segments of the population because this species is polygynous (Fay 1982). Sea ice dynamics in the Chukchi Sea in June and international restrictions precluded us from a sam-pling approach that could be consistently implemented among years. In 2013, 2014, and 2016, all sampling was conducted in the US exclusive economic zone (EEZ). In 2015 and 2017, sampling was conducted in both the Russian and US EEZs. We attempted to completely sample the portion of the marginal ice zone that was accessible on each cruise. In addition, we used aerial reconnaissance in both Russian and US EEZs to locate walrus aggregations to increase efficiency. We collected biopsy samples in pulses in space and time because walruses are heterogeneously distributed throughout their range (Beatty et al. 2020). We typically encountered large aggregations comprised of many groups ranging in  Beatty et al. 2020). We biopsied as many individuals as possible in each group and within each aggregation until no animals remained on the ice to biopsy. Biopsy sampling methods, laboratory methods, and the genotype matching algorithm were detailed by Beatty et al. (2020). Briefly, we used crossbows and darts outfitted with removable biopsy tips to sample animals hauled-out on sea ice. In the field, we assigned sampled animals to 1 of 7 age categories (calf, 1, 2, 3, 4−5, or ≥ 6 yr old, and unknown) based on the width and/or depth of the snout compared to the length of the tusks (Fay 1982, Citta et al. 2014. All procedures were approved by the US Fish and Wildlife Service Region 7 Animal Care and Use Committee (permits 2011Committee (permits -02, 2017 and the Alaska Department of Fish and Game Animal Care and Use Committee (permits 2013Committee (permits -20, 2014Committee (permits -03, 2015, and carried out under Marine Mammal Protection Act permits (MA801652-7 and MA039386-2).
We generated multilocus genotypes with 116 single-nucleotide polymorphisms (SNPs) (Cook et al. 2020, GenBank accession numbers MH399912− MH400027) using TaqMan real-time PCR. We assigned sex to each sample with a sex-linked SNP locus (Table S1 in the Supplement at www.int-res. com/articles/suppl/m697p167_supp.pdf). Briefly, we quantified DNA using a commercial assay kit (Thermo Fischer Scientific) and fluorometer (Perkin -Elmer). We diluted samples to 50 ng μl −1 , and reextracted samples with lower DNA concentrations. We processed the samples using prepared OpenAr-ray® TaqMan SNP genotyping assays on a real-time PCR system (QuantStudio 12K Flex; Thermo Fischer Scientific) using standard software and default running specifications (QuantStudio 12K Flex Software v1.2.2). Each locus was independently genotyped with the following conditions: pre-PCR hold at 93°C for 10 min at 100% ramp rate; 50 PCR cycles of 95°C for 45 s at an 84% ramp rate, 94°C for 13 s at a 100% ramp rate, 53.5°C for 2 min 14 s at 44% ramp rate; followed by post-PCR hold at 25°C for 2 min at 100% ramp rate. Resulting genotype cluster plots for all runs were combined and scored with commercial software (TaqMan Genotyper Software v1.3; Thermo Fischer Scientific), and samples with call rates below 80% or excessive homozygote calls were re-run. We verified the sex-linked SNP with a 3-step process that compared sex assignments from the SNP to a previously published marker that sequenced fragments of the zinc finger intron (Text S1; Fischbach et al. 2008). We identified individuals using a maximum likelihood approach that allowed for genotyping error and missing data with a custom script in R 4.1.2 (Kalinowski et al. 2006, Wang 2006, Sethi et al. 2016, R Core Team 2020).

Data and model development
Our study was designed to model (re)capture events among annual research cruises to estimate probability of recapture, annual survival, and abundance. However, our approach resulted in multiple biopsy samples from individuals within the same cruise, which provided an opportunity to evaluate the repeatability of the aging method. We first condensed animals assigned to 1 of the 4 juvenile age classes (1, 2, 3, or 4−5 yr old) into 1 juvenile age class because of low recapture rates (see Section 3). Thus, we retained 3 age classes for all data analysis: calf (young of the year), juvenile (1−5 yr), and adult (≥6 yr). We then evaluated the consistency of age class assignments for all animals with ≥ 2 biopsy samples collected during the same research cruise. We examined all pairwise comparisons of age class (e.g. an animal with 2 biopsy samples had 1 comparison, an animal with 3 biopsy samples had 3 comparisons, an animal with 4 biopsy samples had 6 comparisons, etc.), and observed a high discordance rate (0.18) in age class assignments. Individuals were often assigned discordant ages within cruises and assigned inconsistent ages among cruises. Consequently, we required an approach to estimate age-specific survival and recapture probabilities while accounting for error in the age assignment process.
Multievent mark−recapture models represent a framework to estimate survival and abundance while accounting for state classification error (Arnason 1972, Hestbeck et al. 1991, Pradel 2005. Multievent models are a generalization of multistate models, and multistate models are a generalization of the traditional Cormack-Jolly-Seber (CJS) model (Cormack 1964, Jolly 1965, Seber 1965, Arnason 1972, Pradel 2005. The traditional CJS model estimates timedependent recapture and survival probabilities for 1 stratum (e.g. age class), and multistate models generalize this model to multiple strata or states (Lebreton et al. 1999). Multievent models further generalize multistate models by assuming that observed data represent 'events' rather than true states, and these events provide information on the latent, true state of an animal (Pradel 2005). Nevertheless, one common feature of CJS, multistate, and multievent models is that the observation process is typically conditioned on the first capture occasion, which means these models estimate the probability of recapture rather than the probability of capture.
We conceptualized multievent models as a type of hidden Markov model where observed data for an animal (i.e. the event) is dependent on its true state (Royle 2008, Gimenez et al. 2012, Mc Clintock et al. 2020. We specifically implemented the model with a discrete ap proach in a Bayesian framework (Gim enez et al. 2012). We defined 8 states for our model, i.e. the following 4 states for each sex: calf (1), juvenile (2), adult (3), and dead (4). The calf age class comprised exclusively young of the year. The juvenile age class comprised animals ages 1 through 5 yr, and the adult age class was comprised of animals ≥ 6 yr old. At first capture, individuals were in 1 of 6 living states be cause we did not sample dead animals. We defined a sampling occasion as lasting for the duration of 1 research cruise; thus, we had 5 sampling occasions (T = 5) for this study, corresponding to the 5 years in which we conducted cruises (2013−2017). We defined an enumeration as the process of collecting a biopsy from an animal and assigning an age class to that animal. Thus, a walrus could be enumerated multiple times during a sampling occasion. We defined recaptured animals as walruses that were enumerated at least once on ≥ 2 sampling occasions. Consequently, a recaptured animal on sampling occasion t was enumerated at least once on sampling occasion t and enumerated at least once on a previous sampling occasion. Observed events (i.e. data) were integers that represented a combination of the number of times an animal was enumerated during a sampling occasion and age assignments that could arise for an animal given it was enumerated that number of times during a sampling occasion (Table 1).
We used a Bayesian framework, and our model included 4 matrices for each sex, with each matrix conditioned on the true state of the animal. The first matrix was a state transition matrix that represented the annual survival and aging process. The second matrix was either the capture matrix or the recapture matrix. The third matrix was a conditional enumeration matrix that represented the probability that an animal was enumerated exactly once, exactly twice, or at least thrice on a sampling occasion given that it was (re)captured on that occasion. The final matrix accounted for error in the age (state) assignment process and represented age (state) assignment probabilities given the number of times the animal was enumerated on that sampling occasion. The state transition, recapture, and age (state) assignment matrices are typical components of a multievent model (e.g. Gimenez et al. 2012). Our model included an ad ditional matrix (i.e. the conditional enumeration matrix) to model the probability that animals were enumerated multiple times within a sampling occasion.
We specifically describe the 4 matrices for females below, and we indicate differences between male and female matrices when applicable. We assumed all parameters were constant for the 5 yr study, but we use t as a subscript to provide clear interpretation of parameters when appropriate. The female state transition matrix (Φ F ) estimated the probability of transitioning from state (calf, juvenile, adult female, dead) on occasion t −1 (rows) to state on occasion t (columns):  Table 1. Event definitions used for the Pacific walrus multievent mark− recapture analysis. The number of times an animal was enumerated represented the number of times the animal was biopsied and assigned an age class during a sampling occasion. Only 3 animals were enumerated more than 3 times, and we dropped the fourth and fifth enumerations for these animals because enumerations for these 3 walruses were concordant. Age assignment sequence was irrelevant where φ C , φ J , and φ F represent annual calf, juvenile, and adult female survival, respectively, and ψ represents the transition probability from juvenile to adult (Table S2). We set ψ as a function of juvenile survival based on a formula in Caswell (2001): (1) We set λ = 1 because we assumed a stationary population (see Section 2.3), and K = 5 represented the number of years in the juvenile age class (1−5 yr). We specifically estimated ψ as a derived parameter with this approach rather than estimating it as a free parameter because we had a sparse dataset (see Section 3). We conditioned on sampling occasion of first capture, which means that an animal initially captured as a calf can only be recaptured in subsequent sampling occasions as a juvenile or adult. Thus, the first column in the state transition matrix contained all zeros. The male state transition matrix (Φ M ) was identical to the female state transition matrix (Φ F ), except we replaced φ F with φ M . We allowed survival and other parameters (see below) to vary between adult males and adult females because the majority of adult males summer in the Bering Sea. However, some adult males may spend occasional summers in the Chukchi Sea. Thus, most adult males likely ex hibit permanent emigration, whereas other adult males may exhibit temporary emigration. Nevertheless, adult male survival is likely confounded with emigration out of the study area. Our small sample size for adult males and limited number of recaptures (see Section 3) precluded us from accounting for variation in adult male survival from temporary or permanent emigration.
The second matrix was either the annual capture matrix P 0 or the annual recapture matrix (P F ). The annual capture matrix applied to an animal only on its first occasion of capture: where rows represent true states (calf, juvenile, adult, dead), the first column represents not captured on sampling occasion t, the second column represents calf captured on occasion t, the third column represents juvenile captured on occasion t, and the fourth column represents an adult female captured on occasion t. The annual recapture matrix applied to sampling occasions after the first capture occasion and estimated the probability that a female walrus was recaptured (i.e. enumerated at least once on occasion t given it was enumerated at least once on a previous occasion) on occasion t: where rows are defined as in P 0 , and columns are: not recaptured on occasion t, calf recaptured on occasion t, juvenile recaptured on occasion t, and adult female recaptured on occasion t. Thus, ρ J represents the probability that a juvenile walrus was recaptured on occasion t, and ρ F represents the probability that an adult female was recaptured on occasion t (Table S2). The annual recapture matrix for males (P M ) was identical except it contained a different parameter for adult males (ρ M ). We made a conventional assumption associated with CJS models of no recapture heterogeneity within modeled groups (Williams et al. 2001). Specifically, we assumed that biopsied juveniles had the same recapture probability, biopsied adult females had the same recapture probability, and biopsied adult males had the same recapture probability.
The third matrix was the conditional enumeration matrix (Θ), which represented the probabilities that an animal that was (re)captured on occasion t was enumerated a certain number of times on occasion t. We set the maximum number of enumerations within a sampling occasion (i.e. research cruise) to 3, which was based on our observed data. The overwhelming majority of animals were enumerated once on an occasion, and successively smaller proportions were enumerated exactly twice, exactly thrice, and more than thrice on an occasion. Only 2 animals were enumerated 4 times within an occasion and only 1 animal was enumerated 5 times. Age classifications for these 3 animals were concordant, so we simply removed the fourth and fifth enumerations for data analysis. We did not extend our matrices to 4 and 5 enumerations because age classifications for these 3 animals were concordant, and extending our model to 4 or 5 enumerations would introduce matrices that had effective zeros in many columns. Thus, we estimated probabilities that a walrus was enumerated exactly once, exactly twice, or at least thrice on sampling occasion t given it was (re)captured on occasion t with: where rows represent: not (re)captured on occasion t, calf captured on occasion t, juvenile (re)captured on occasion t, and adult female (re)captured on occasion t. Columns represent: (1) not (re)captured, (2) calf enumerated once on occasion t given it was captured on occasion t (θ C1 ), (3) calf enumerated twice on occasion t given it was captured on occasion t (θ C2 ), (4) calf enumerated at least thrice on occasion t given it was captured on occasion t (θ C3 ), (5) juvenile enumerated once on occasion t given it was (re)captured on occasion t (θ J1 ), (6) juvenile enumerated twice on occasion t given it was (re)captured on occasion t (θ J2 ), (7) juvenile enumerated at least thrice on occasion t given it was (re)captured on occasion t (θ J3 ), (8) adult female enumerated once on occasion t given it was (re)captured on occasion t (θ F1 ), (9) adult female enumerated twice on occasion t given it was (re)captured on occasion t (θ F2 ), and (10) adult female enumerated at least thrice on occasion t given it was (re)captured on occasion t (θ F3 ) ( Table S2). The conditional enumeration matrix for males (Θ M ) was identical except θ F1 , θ F2 and θ F3 were replaced with θ M1 , θ M2 and θ M3 , respectively, to allow conditional enumeration probabilities to vary between adult males and adult females. We assumed that conditional enumeration probabilities (Θ M , Θ F ) were constant for the 5 yr study. We also assumed that all biopsied calves had the same conditional enumeration probabilities (θ C = θ C1 , θ C2 , θ C3 ), all biopsied juveniles had the same conditional enumeration probabilities (θ J = θ J1 , θ J2 , θ J3 ), all biopsied adult females had the same conditional enumeration probabilities (θ F = θ F1 , θ F2 , θ F3 ), and all biopsied adult males had the same conditional enumeration probabilities (θ M = θ M1 , θ M2 , θ M3 ). We had confidence that we met these assumptions for 2 reasons. First, we allocated similar biopsy sampling effort into all walrus groups/aggregations encountered during field work. Second, the overwhelming majority of repeat enumerations for a single animal occurred within 24 h (see Section 3), suggesting that conditional enumeration of walruses was a function of sampling effort at the level of the aggregation, and did not require population mixing or repeated passes through the entire population. The matrix products P F Θ and P M Θ represented the joint probability that a female or male walrus, respectively, was recaptured and enumerated exactly once, exactly twice, or at least thrice on occasion t (Fig. S1). Consequently, our analysis assumed that all biopsied juveniles had the same set of joint recapture/enumeration probabilities on each occasion (ρ J θ J1 , ρ J θ J2 , ρ J θ J3 ), all biopsied adult females had the same set of joint recapture/enumeration probabilities on each occasion (ρ F θ F1 , ρ F θ F2 , ρ F θ F3 ), and all biopsied adult males had the same set of joint recapture/ enumeration probabilities on each occasion (ρ M θ M1 , ρ M θ M2 , ρ M θ M3 ).
In multievent mark−recapture models, the assignment matrix estimates error in the state assignment process where rows represent true states, and columns represent events (i.e. observed data). A simplified age assignment matrix can be written with rows that represent true states (not enumerated, calf enumerated once, juvenile enumerated once, and adult enumerated once) and columns that represent not assigned, assigned calf, assigned juvenile, and assigned adult with a row-stochastic matrix: Thus, δ C|C represents the probability that a calf was correctly assigned calf, δ J|C represents the probability that a calf was incorrectly assigned juvenile, and δ A|C represents the probability that a calf was incorrectly assigned adult (δ C|C + δ J|C + δ A|C = 1). The above assignment matrix would be applicable to a study where animals were enumerated once. However, we required an expansion of the above assignment matrix to account for multiple enumerations of individuals within a sampling occasion. Thus, our fourth matrix was an expanded age assignment matrix that defined observed events (i.e. data) as the different possible combinations of unordered age classes that an animal could be assigned given the number of times it was enumerated on sampling occasion t. Our assignment matrix (Δ) was a 10 × 20 matrix where rows represent the number of enumerations for each true state on a sampling occasion, and columns represent possible age classifications (Fig. 2). All rows in Δ sum to 1.0. For example, cell (3, 7) in Δ gives the probability (2 δ C|C δ J|C ) that a calf enumerated twice on an occasion was classified once as a calf and once as a juvenile, which could occur in either order. We removed walruses that were enumerated once and assigned unknown age class. However, we retained animals that were enumerated twice and assigned a valid age class for one enumeration and an unknown age class for the other enumeration. For these ani- J|A A|A mals, we removed the unknown age enumeration. Thus, our age assignment matrix did not include parameters for animals assigned to an unknown age class.

Abundance
We used a Horvitz-Thompson ratio estimator (Horvitz & Thompson 1952) to estimate mean juvenile abundance as a derived quantity (N J ), which assumed that recapture probability equaled capture probability for all juvenile walruses: ( 2) where I i,J,t is an indicator (0/1) that denotes whether the latent, true age (z i,t , see Section 2.4) of animal i was a juvenile on sampling occasion t, ρ J is recapture probability for juveniles, n t is the total number of animals captured on sampling occasion t, and T = 5 is the number of sampling occasions. We estimated adult female abundance as a derived quantity (N F ) similarly by re placing I i,J,t with I i,F,t , and ρ J with ρ F . Consequently, we assumed that recapture probability equaled capture probability for all adult female walruses. In a typical multistate model, the indicator I i,J,t is conditioned on information provided to model as data (i.e. it is known without error) whereas in our implementation the indicator I i,J,t is conditioned on random variable z i,t (see Section 2.4).
We could not estimate calf abundance (N C ) with a Horvitz-Thompson estimator because we could not estimate calf recapture probability (i.e. calves always transitioned from calf to juvenile or dead after the first capture occasion). As an alternative approach, we used age structure data (Jay et al. 2020) to estimate the proportion of the calf−juvenile−adult female walrus population that was comprised of juveniles and adult females, which we then used to estimate calf, juvenile, and adult female abundance (N CJF ) as a derived parameter. Our approach (see below) assumed a stable age structure and a constant population size for the 5 yr period of our study. Age structure data were collected on the same research cruises where biopsy samples were collected in 2013−2017. Observer teams sampled groups of walruses hauled out on sea ice and categorized all walruses in a group into 1 of 9 age classes (calf, 1, 2, 3, 4−5, 6−9, 10−15, >15 yr of age, and unknown) and assigned sex to all animals ≥ 6 yr (Fay 1982, Citta et al. 2014. We removed all groups that had at least 1 walrus assigned adult enumerated exactly once, (9) adult enumerated exactly twice, and (10) adult enumerated at least thrice. Columns represent the age assignments listed in Table 1 to an unknown age class and removed all adult males from the dataset. To estimate the proportion of the calf−juvenile−adult female population comprised of juveniles and adult females, we reduced the age structure data to 2 age classes: calf and non-calf. The non-calf age class included juveniles and adult females, whereas the calf age class contained exclusively calves. Thus, our age structure data were reduced to a 5 × 2 matrix, L, where rows corresponded to 5 yr of age structure data (2013−2017) and columns represented the number of animals in each of the 2 age classes (calf, non-calf) for each year. Thus, l t,1 and l t,2 are the number of animals observed in the calf and non-calf age categories in year t of the age structure surveys, respectively. Furthermore, H t = l t,1 + l t,2 is the total number of animals aged in year t, which is the number of trials and provided as data.
We used a beta-binomial likelihood function to estimate the proportion of the population that was in the non-calf age class (Citta et al. 2014): where α and β are parameters to be estimated. H represents a vector of length 5 comprised of the number of trials for each year age structure data were collected. The beta-binomial pro bability (ξ) is α/(α + β), which is the expected proportion of non-calves based on all 5 yr of data. We used this approach because we were specifically interested in obtaining a mean estimate of the beta-binomial probability ξ for animals in the non-calf age category for all 5 yr of age structure data while accounting for overdispersion among years. To estimate walrus abundance in the Chukchi Sea, we divided the total number of juveniles and adult females by ξ: Thus, the estimated number of calves was derived as Although we occasionally sampled adult males, we did not estimate adult male abundance with the Horvitz-Thompson estimator because most adult male walruses summer in the Bering Sea. Thus, the adult male recapture rate (ρ M ) may be confounded with emigration. Although skewed sex ratios have been proposed for walruses (e.g. 1:3 to 1:5), these represented sex ratios for breeding adults, and male walruses do not reach breeding age until approximately 15 yr (Fay 1982). Consequently, we assumed a 1:1 sex ratio for walruses ≥ 6 yr old to estimate total walrus abundance: N = N CJF + N F. Total abundance (N ) represents the estimated number of walruses that reside in the Bering and Chukchi Seas. In addition, we estimated minimum abundance (N min ) as the 20% quantile of the posterior distribution for total abundance, which is a value that agencies are required to report to execute provisions of the Marine Mammal Protection Act (Wade 1994).

Model description and implementation
We jointly fit the multievent model and derived abundance estimates in a Bayesian framework in Just Another Gibbs Sampler (JAGS) v.4.3.0 (Plummer 2017) via jagsUI (Kellner 2021). For the multievent likelihood function, we defined x as a vector of length n (the total number of walruses enumerated at least once) that identified the first capture occasion for each walrus, and we considered sex as known based on the genetic marker. We modeled true state for individual i on the occasion of first capture with a categorical distribution: where π is a vector of length equal to the number of true states (calf, juvenile, adult, dead) containing 4 probabilities that sum to 1. Thus, π represents the probabilities that an animal was in states calf, juvenile, adult, or dead, respectively, on first capture. We did not sample dead animals, and walruses could not enter the study on first capture in the dead state. Thus, we set the fourth probability to 0 and used a Dirichlet prior for π with parameters (1, 1, 1, 0). In preliminary runs of the model, informative priors on π did not meaningfully alter the posterior distribution of π.
We modeled true state for individual i for all subsequent occasions conditioned on state at t − 1. The row of the transition matrix (Φ F or Φ M ) that corresponded to the state of individual i at t − 1 provided arguments for the categorical distribution: We introduce the observation matrix (Ο) for convenience, which is the matrix product of the (re)capture matrix (P 0 , P F , or P M ), the conditional enumeration matrix (Θ F or Θ M ), and the age assignment matrix (Δ). Thus, Ο is a 4 × 20 matrix (Fig. S2) that represents the probabilities of observing any possible combination of (re)capture, conditional enumeration, and age assignments on occasion t conditioned on true state. The observation matrix allows us to link the true state of the animal to our data. We modeled the observations for the first capture occasion for individual i with a categorical distribution with arguments: where О 0 = P 0 ΘΔ. We modeled subsequent occasions as: where О = PΘΔ, P was P F or P M as appropriate.
For the derived estimate of Chukchi Sea walrus abundance, we fit the beta-binomial likelihood function using its component parts (Citta et al. 2014) because JAGS does not have a function for the betabinomial distribution: where ω t is the estimated proportion of calves−juveniles−adult females in the not-calf age category in year t.
We used vague uniform priors between 0 and 1 for all annual survival rates and all recapture probabilities. We used vague Dirichlet(1,1,1) priors for all conditional enumeration probabilities (θ C , θ J , θ F , θ M ). We defined prior distributions for age assignment probabilities based on uniform distributions between 0.5 and 1.0 for the correct classifications (δ C|C , δ J|J , δ A|A ), and conditional Dirichlet distributions for the 2 incorrect classifications for each age class. For example, we assigned δ C|C ~ Uniform(0.5, 1.0) and δ * J|C , δ * A|C ~ Dirichlet(1,1) and δ J|C = δ * J|C (1-δ C|C ) and δ A|C = δ * A|C (1-δ C|C ). Thus, priors for δ J|C and δ A|C were uniformly distributed in bivariate space, and we assumed we were more likely than not to correctly classify an individual. We used vague gamma priors for α and β with shape = 0.001 and rate = 0.001. We ran 19 parallel chains with 125 500 iterations for each chain and discarded the first 20 000 iterations as adaptation and burn-in. We thinned chains at a rate of 500 to obtain a total of 4009 draws from the posterior distribution. To assess model convergence, we inspected trace plots and evaluated parameters with the potential scale reduction factor (Brooks & Gelman 1998) with values <1.05 interpreted as indicating convergence. To assess model fit, we conducted posterior predictive checks on the number of recaptures in each age class (conditioned on true state, z i ) and the number of animals in the non-calf age class from the age structure surveys for each year. We considered Bayesian p-values 0.10 ≤ p ≤ 0.90 to indicate adequate fit. We ran the model on the USGS supercomputer Yeti (Falgout et al. 2021).

Study duration evaluation
Research cruises to collect biopsy samples required extensive investments of human and financial capital. Thus, we evaluated the efficacy of a 3 yr mark− recapture study compared to a 5 yr study with a series of simulations to inform whether useful walrus abundance estimates could be obtained with a shorter-duration study. We specifically identified these 2 study durations because a study period > 5 yr would require intractable investments of financial and human capital. We simulated 5 yr datasets (n = 50) and 3 yr datasets (n = 50) using a data-generating process that mirrored the model outlined above. We calculated 'true' abundance of juveniles and adult females based on the mean number of first captures + recaptures for all 5 years in the empirical dataset and the posterior mean values for juvenile and adult female recapture probability (Text S2). For the remaining parameters, we set 'true' parameter values to the posterior means from the above model fit to observed data. We calculated the expected number of first captures for subsequent years based on the number of animals in the population that had not been captured previously (Text S2). We generated simulated data with custom scripts in R 4.1.2 (R Core Team 2020).
We fit the above model to all simulated datasets and calculated standard error and coefficient of variation for each parameter in the 5 and 3 yr simulations to evaluate precision. We also calculated absolute and relative bias for each parameter in the 5 and 3 yr simulations. Markov chain Monte Carlo details were identical to those used in the analysis of the real data. We ran all simulations to evaluate bias and precision on the USGS supercomputer Yeti (Falgout et al. 2021).

Sample collection, laboratory methods, and sample matching
We collected tissue samples from 8303 unique individuals. We removed 7 individuals that did not have accompanying field data, and 81 animals that did not have sex assigned from the genetic marker. We removed 50 individuals that were biopsied once and assigned an unknown age class. We had 10 animals that were biopsied twice within an occasion and assigned a valid age class for one sample and an unknown age class for the other sample. We removed the unknown classification from the dataset for these 10 walruses and retained the valid age class assignment for analysis. Thus, our final dataset included 8165 unique walruses with 5890 females and 2275 males with 82 recaptures (Fig. 1). Furthermore, 669 individuals were enumerated within an occasion at least twice, and 52 of those 669 walruses were enumerated within an occasion at least thrice. We enumerated 492 male walruses that had at least one valid age class assignment as adult, and we recaptured only 3 of these. Approximately 83% of enumerations for animals enumerated exactly twice or at least thrice occurred within 12 h of one another, and approximately 87% occurred within 24 h of one another.

Model results
The model converged with all potential scale reduction factors for fitted parameters ≤ 1.05, and all Bayesian p-values indicated adequate fit. We estimated the probability that a juvenile was recaptured (ρ J ) as 0.008 (95% credible interval [CrI]: 0.004−0.015). We estimated the probability that an adult female was recaptured (ρ F ) as 0.011 (95% CrI: 0.007−0.016) and the probability that an adult male was recaptured (ρ M ) as 0.010 (95% CrI: 0.001−0.042). The probability that a juvenile transitioned to an adult was estimated as 0.068 (95% CrI: 0.014−0.147). We observed a consistent pattern in conditional enumeration probabilities with the highest values observed for animals enumerated once and the lowest values for those enumerated at least thrice as expected (Table 2). Estimates for parameters in the age assignment matrix (Δ) also exhibited consistent patterns with probabilities > 0.80 for correct assignments (δ C|C , δ J|J , δ A|A ) and values < 0.11 for all incorrect assignments (Table 3).
We observed heterogeneity in survival among age classes and sexes. Juvenile and adult female survival were estimated as 0.626 (95% CrI: 0.388− 0.865) and 0.904 (95% CrI: 0.744−0.996), respectively (Fig. 3). The 95% CrI for calf survival was 0.094−0.958, indicating that this parameter was not estimable given available data. The 95% CrI for adult male survival was 0.117−0.964. Thus, we do not report the posterior mean of adult male survival because our estimate was imprecise. Furthermore, our study was not designed to sample adult males, and we recognize that this survival estimate is likely conflated with emigration from the study domain. We estimated the mean abun-dance of juveniles and adult females as 163 565 (95% CrI: 104 455 − 250 520) and the mean abundance of calves, juveniles, and adult females as 187 912 (95% CrI: 118 590 − 287 380) (Fig. 4). We estimated total walrus abundance as 257 193 animals with a 95% CrI of 171 138 − 366 366 and coefficient of variation of 0.19 (Fig. 4). In addition, we estimated the minimum abundance (N min ) as 214 008 (Fig. S3).

Study duration evaluation
Standard errors for ρ J , ρ F , φ C , φ J , and φ F ranged from < 0.01 to 0.15, while coefficients of variation  Table 3. Estimates for parameters from the age assignment matrix (Δ) to account for age misclassification for Pacific walrus calves (C), juveniles (J), and adults (A). The first subscript refers to assigned age, and the second subscript refers to true age ranged from 0.07 to 0.34 for the 5 yr simulations (Table S3). Standard error and coefficient of variation for juvenile and adult female abundance were 28 085 and 0.21, respectively. Absolute bias for the simulated 5 yr studies ranged from −0.05 to < 0.01 while relative bias ranged from −0.06 to 0.24 for ρ J , ρ F , φ C , φ J , and φ F (Table S3). Generally, parameters that estimated recapture probabilities (ρ J , ρ F , ρ M ) exhibited high coefficients of variation and relative bias because true values for these parameters were low (≤ 0.02).
The 3 yr simulations did not perform as well overall as 5 yr simulations. Standard errors for ρ J , ρ F , φ J , φ F , and juvenile and adult female abundance were higher in the 3 yr simulations compared to the 5 yr simulations (Table S4). Coefficients of variation for ρ J , ρ F , φ C , φ J , φ F , and juvenile and adult female abundance were 1.2 to 3.5 times higher in the 3 yr simulations compared to the 5 yr simulations. Absolute values of both bias metrics for ρ J , ρ F , φ C , φ J , φ F , and juvenile and adult female abundance were 5.1 to 18.7 times higher for the 3 yr simulations compared to the 5 yr simulations (Table S4).

DISCUSSION
We estimated age-specific walrus survival and abundance with multievent mark−recapture models that accounted for age misclassification. We specifically accounted for age (i.e. state) misclassification because it can lead to biased parameter estimates in multistate models, which could lead to biased estimates of population size (Kendall 2009). We generated the first total abundance estimate for the walrus population since the most recent aerial survey was conducted in 2006 (Speckman et al. 2011). The contemporary abundance estimate (257 193, 95% CrI: 171 138−366 366) has a substantially higher point estimate than the 2006 aerial survey-based estimate (129 000, 95% CI, 55 000−507 000). However, these 2 estimates are statistically indistinguishable, as our credibility intervals are completely contained within the confidence intervals of the 2006 survey. Speckman et al. (2011) acknowledged their abundance estimate was biased low because only half of the available habitat was surveyed. Therefore, we cannot address population trends because of the known bias in the previous estimate and the difference in methodology between the historical aerial surveys and our genetic mark−recapture estimate.
Although we could not estimate annual calf survival with available data, we achieved reasonably precise estimates of juvenile and adult female survival. Our posterior mean estimate of juvenile survival (0.63, 95% CrI: 0.39−0.87) was lower than previous estimates. Taylor et al. (2018) developed an integrated population model to examine Pacific walrus demography for the same population, and re ported juvenile and reproductive adult female survival estimates. Their most parsimonious model held juvenile survival constant over time (1975−2015). Our juvenile survival estimate (0.63) was lower than that reported by Taylor et al. (2018), and our CrI (0.39−0.87) did not encompass their posterior mean for juvenile survival (0.90, 95% CrI 0.83−0.97; Taylor et al. 2018). Similarly, our posterior mean estimate of adult female survival (0.90) was lower than that estimated by Taylor et al. (2018), but our CrI (0.74−1.00) included the posterior mean estimates for reproductive adult survival from the model of Taylor et al. (2018) (0.99, 95% CrI: 0.98−1.00). Notably, Taylor et al. (2018) estimated separate survival rates for reproductive adult females and senescent females, whereas we estimated a single survival rate for all adult fe males. Thus, we attribute some of the differences in adult female survival estimates between Taylor et al. (2018) and our study to the fact that we estimated survival for a single adult female age class. In addition, we attribute at least part of the apparent difference in juvenile survival rate estimates (i.e. the large difference between point estimates) to the relatively low precision for this parameter in our model compared to the analysis by Taylor et al. (2018). We estimated survival from a 5 yr mark−recapture study with a low recapture rate, whereas Taylor et al. (2018) estimated survival from a population model that included age structure data from 9 ship-based surveys spanning 35 yr.
Permanent emigration is confounded with mortality in our model as is typical of many mark−recapture models (Lebreton et al. 1992). For example, some 4− 5 yr old males may spend summers in the Bering Sea, which would represent permanent emigration from the Chukchi Sea sampling area and result in negatively biased juvenile survival. To acknowledge that we cannot distinguish permanent emigration from death, our survival estimates are best described as apparent survival estimates (Lebreton et al. 1992, Williams et al. 2001, Gilroy et al. 2012. Furthermore, movement between Russian and US EEZs could be interpreted as a form of temporary emigration, which may have impacted our survival and recapture estimates. Random movement between the 2 EEZs would generate unbiased estimates of survival (not withstanding bias from permanent emigration) whereas Markovian movement between the 2 EEZs would generate biased estimates of survival (Kendall et al. 19 97). A multievent model could be developed that included additional unobservable states for each age class in the Russian EEZ in 2014 and 2016 and observable states for each class in the Russian EEZ in 2015 and 2017. The model would require high recapture probabilities and large sample sizes to obtain precise estimates because of the complex structure (Bailey et al. 2010), which was infeasible for this sparse walrus dataset. Nevertheless, we conclude that our survival estimates are biased low because of permanent and/or temporary emigration.
Results from the 5 yr simulations suggest that our abundance estimate for juveniles and adult females may be biased low (Table S3). However, we assumed a 1:1 sex ratio for adult females to adult males to derive a total abundance estimate, which likely introduced positive bias for adult male abundance. Pacific walruses are one panmictic population (Sonsthagen et al. 2012, Beatty et al. 2020; therefore, sampling in the Russian EEZ on 2 of 5 sampling occasions for the study should not bias our abundance estimate per se. Thus, our total abundance estimate is likely broadly applicable to the entire Pacific walrus population. We developed a new multievent model formulation for studies that enumerate (i.e. observe) animals mul-tiple times during a sampling occasion. Although we considered a robust design to account for temporary emigration, our study differs from the robust design in several important ways (Pollock 1982, Kendall et al. 1997, Kendall & Bjorkland 2001. The robust design requires 2 types of distinct sampling occasions: primary periods and secondary periods (Pollock 1982) . In the robust design, the population is open to births, deaths, emigration, and immigration between primary periods and closed or open between secondary periods (Pollock 1982, Kendall & Bjorkland 2001. In contrast, our model had 1 sampling occasion type (i.e. years) between which the population was open, and our study design did not have distinct occasions that could be used to delineate secondary periods. Additionally, the robust design assumes that secondary sampling occasions represent independent and random samples of the population, but we could not meet this assumption because we could not make repeated passes through the entire population in a single cruise. Furthermore, we wanted to develop a model that could leverage all of the information on age misclassification that is contained within each enumeration (i.e. observation) regardless of the time difference between enumerations. In our study, we observed multiple enumerations on the same animal while sampling a group because of the challenges in monitoring which animals had already been biopsied in a group.
Age misclassification in our study likely arose from multiple sources. First, field personnel may have incorrectly identified which animals were biopsied, and this is evident in the misclassification probabilities for adults and calves (Table 3), which are easily distinguished in the field by size. Second, field personnel may have recorded data incorrectly in field notebooks or in data storage software. Although we ensured notebook records coincided with electronic records for all animals with discordant ages, this did not preclude errors when records were entered into notebooks. Third, age misclassification also might arise from natural variation in snout depth/width and tusk length. Sample sizes used to develop the walrus age classification method were relatively small, and standard deviations for morphometric parameters imply a substantial amount of overlap among some age classes (Fay 1982, Citta et al. 2014). Finally, walruses were often sleeping or facing the opposite direction of field personnel when biopsied. In such instances, age assignments must be made quickly, and errors should be expected. However, we developed a framework to deal with age misclassification errors, regardless of their origin.
The Pacific walrus population likely underwent a large population decline that began no later than 1981 and moderated in the 1990s (Taylor & Udevitz 2015, Taylor et al. 2018). The population decline was likely the result of high harvests in the 1980s and a food-limited population that had reached carrying capacity sometime in the 1970s or early 1980s (MacCracken et al. 2017). As of 2015, the population was likely stable, and was approximately 42% (95% CrI: 23−64%) of the 1981 population size (Taylor et al. 2018).
Pacific walrus population assessments are challenging due to the species' broad geographic range, diving behavior, and sea ice habitats. Our genetic mark−recapture study provided a novel opportunity to generate a relatively precise abundance estimate for the walrus population. The 5 yr field effort required a significant investment of financial resources and human capital, beginning with collection of biopsy samples from adult males in 2011 for SNP marker discovery and development (Cook et al. 2020). Planning and execution of research cruises re quired extensive project management and administration from 2012 through 2018. After the field effort was completed, samples were genotyped and genetic analyses were conducted to evaluate study design (Beatty et al. 2020) before the current analysis could be developed. Unfortunately, our study duration simulations indicated that a 3 yr study instead of a 5 yr study would likely lead to unacceptably imprecise abundance estimates. Our genetic mark−recapture dataset could be used in a combined model with data from age structure surveys, which could im prove precision of survival and abundance estimates. An integrated population model that included all available data, including past aerial surveys, age structure, mark−recapture, and harvest (sensu Taylor et al. 2018) could improve our understanding of the historical fluctuations of the walrus population. An integrated model could also provide insight into how the species may respond to a rapidly changing Arctic.
Walruses are a vital cultural and subsistence resource to communities that border the Bering and Chukchi Seas in the USA and Russia. Thus, estimates of walrus survival and abundance can provide information to these communities on the status of this species in the region. The walrus population likely will be under increased stress in the future due to sea ice loss and its related effects (e.g. increased shipping activity, altered spatial distribution and availability of benthic biomass), which will almost certainly result in a population decline (Grebmeier et al. 2015, Mac-Cracken et al. 2017. Sea ice facilitates access to productive offshore foraging areas for walruses in the Chukchi Sea in summer , and decreasing sea ice extent during the summer months will either reduce the amount of accessible foraging habitat, alter the energy budgets of individuals that travel from coastal haul-outs to offshore foraging areas , or result in some combination of the two. Furthermore, benthic biomass could be impacted throughout the range of the walrus because of declining sea ice extent, ocean acidification, and an altered ecosystem state (Huntington et al. 2020). The severity and magnitude of the projected walrus population decline is uncertain (MacCracken et al. 2017); however, our abundance estimate provides a benchmark to assess trends in walrus population abundance in the future.