Estimates of regional annual abundance and population growth rates of white sharks off central California

Determining population trends is critical for evaluating management actions and prioritizing species protections. In this study, we used empirical data to produce an estimate of the population trend for sub-adult and adult white sharks in central California. We used the unique dorsal fin morphology to build a mark-recapture data set in a modified Jolly-Seber model (POPAN formulation) to estimate annual abundance and then investigate population growth rates using parametric bootstrapping methods for sub-adult and adult sharks (males and females). For all demographic groups combined, we found equivocal evidence for a positive regional population growth ( λ = 1.07 (95% CI = 0.91 to 1.23)). However, sex-and size-specific population growth rate estimates provided some evidence of population increases for reproductively mature males ( λ = 1.06 (95% CI = 0.99 to 1.13)) and females ( λ = 1.06 (95% CI = 0.95 to 1.17)). For sub-adult male and female white sharks, point estimates of λ were positive but uncertainty prevents strong inference ( λ = 1.07 (95% CI = 0.85 to 1.29)) and ( λ = 1.08 (95% CI = 0.88 to 1.28)), respectively. Our findings of a potential increase in reproductive-aged white sharks in central California may be a result of regional fluxes in density or attributed in part to current protection efforts and subsequent increase in abundance of pinnipeds as well as reduced gill-net fisheries mortality of juveniles. A trend estimate for the entire northeastern Pacific will require obtaining similar data across known aggregation areas along the west coast of North America.


Introduction
Determining whether the number of individuals in a population is increasing, decreasing or stable over time is a central goal for assessing the status of wildlife populations (Williams et al., 2002). However, population trends can be difficult to estimate because long-term monitoring is required, especially for long-lived species. In marine systems, where observations beneath the sea surface are challenging, assessing non-target species remains problematic (Hilborn and Walters, 1992;Pope et al., 2000). This is even more so for large, highly mobile species with extensive home ranges and relatively low abundance. As a result, ~20% of marine vertebrates globally remain in the 'Data Deficient' assessment category under the International Union for Conservation of Nature (IUCN, 2019).
While catch data is traditionally used in marine fisheries to assess populations, the proliferation of high-quality and low-cost photographic and video technologies, along with pattern recognition techniques, has increased opportunities to assess more cryptic species through photo identification mark recapture (Drechsler et al., 2015). Simultaneously, the increasing accessibility of genetic (Miller et al., 2005) and electronic biologging data (Dudgeon et al., 2015) is helping define the relevant geographic scope of populations targeted for assessment using imagebased mark recapture.
Elasmobranchs are a taxa of conservation concern due to intrinsically low reproductive rates and vulnerability to direct and indirect fishing mortality, which has decimated many populations (Myers et al., 2007;Worm et al., 2013) and numerous species lack robust population assessments due to insufficient data (Dulvy et al., 2008;Queiroz et al., 2019). However, many shark and ray species possess morphological (Holmberg et al., 2009;Tienhoven et al., 2007), behavioral (Barnett et al., 2011;Kneebone et al., 2012) and life history traits (Cortés, 2000) that are particularly suited for assessment using photo identification mark-recapture.
The northeast Pacific population of sub-adult and adult individuals (hereafter called white sharks) primarily ranges from aggregations in Guadalupe Island, Mexico (Domeier and Nasby-Lucas, 2008) and central California, USA out to the north Pacific subtropical gyre ('White Shark Café') (Jorgensen et al., 2010;Weng et al., 2007). Additional aggregation areas may exist but to date none have been described in the north Pacific. In central California, a predictable migratory pattern reflects high fidelity for both sexes to the region and white shark individuals can be identified by a unique pattern of notches on the trailing edge of the dorsal fin that does not change significantly through time . Taken together, these attributes have led to multiple markrecapture studies to estimate essential population parameters (Chapple et al., , 2016Kanive et al., 2015Kanive et al., , 2019).
An initial photo ID based estimate of abundance for white sharks off central California from a three-year study estimated 219 individuals (95% credible intervals, 130-275) . That study provided a proof of concept for photographic mark-recapture as a viable population assessment tool for white sharks, however no baseline value pre-dated this estimate so no temporal trend could be inferred. White shark populations have a low resilience to mortality (Chapple and Botsford, 2013) and have been protected since 1994 in California waters (Heneman and Glazer, 1996). In addition, pinniped populations in North America have increased (Boeuf and Laws, 1994;Laake et al., 2018) following the Marine Mammal Protection Act (1972) potentially influencing white shark population growth as pinnipeds are a preferred prey source while white sharks are coastal. Furthermore, policy measures such as banning the use of gillnets within 4.83 km (three miles) of the California coast (California proposition 132, 1990) may be benefitting white sharks by reducing incidental mortality of juveniles (Benson et al., 2018;Lowe et al., 2012). Therefore, we hypothesize that sub-adult and adult white sharks may have experienced a recent positive population growth.
To test this, we estimated the abundance and population growth rates of a regional sub-population of white sharks in central California (hereafter simplified as 'population') using 8 years of data from all 3 known aggregation sites. We tested competing models within a Jolly-Seber mark-recapture framework, based on size and sex-specific demographic assumptions derived from previous work (Chapple et al., 2016;Kanive et al., 2015Kanive et al., , 2019 to estimate annual abundance and determine whether there have been changes in the population size over time. This was done for the overall population, and across demographic groups (based on sex and maturity) within the sampled population.

Data collection
From 2011 to 2018, we sampled at known white shark aggregation sites (Southeast Farallon Island, Año Nuevo Island, and Tomales Point) along the central coast of California, USA ( Fig. 1) during periods of peak coastal residency (September-February; Klimley, 1985;Jorgensen et al., 2010). We pooled data from the three sites given their proximity (<150 km; total study area =~1500 km 2 ) and the species' mobility.
White sharks were attracted to a research boat using a seal decoy and a localized scent of marine mammal (Kanive et al., 2015). Individual sharks were identified from photographs and video images of the unique patterns on the trailing edge of their dorsal fin. The dorsal fins of white sharks have unique trailing edge markings shown to be stable for >22 years , establishing these morphological differences as a reliable means of identifying individuals. This methodology has been used in previous population studies of white sharks (Chapple et al., 2016Towner et al., 2013;Kanive et al., 2015Kanive et al., , 2019. To estimate abundance, we used seasonal sighting records of all sampled individuals to generate mark-recapture encounter histories. We determined sex by unambiguous confirmation of presence (male) or absence (female) of claspers using underwater video. In this study, because sex assignment was imperfect, some sharks remained unsexed throughout the study and were thus classified as "unknown sex". To deal with this uncertainty, we used the program LOLASURVIV (Nichols et al., 2004) to estimate the sex ratio of the sample population. Unsexed individual sharks were then randomly assigned to sex based on the estimate of the sex ratio from the LOLASURVIV analysis. Multiple iterations of random assignment in the dataset were run against the model suite to check for sensitivity of the resulting annual abundance estimates. Data used to estimate sex ratio extended from 2006 to 2018 (https://zenodo. org/record/4022914#.X2IxTGdKjOR).
Total body length was estimated as the shark swam near a vessel of known size and maturity was assigned based on published values (adult males ≥380 cm, sub-adult males < 380 cm, adult females ≥440 cm, and sub-adult females < 440 cm based on Pratt (1996) and Francis (1996), respectively). These length-based classifications resulted in four groups: adult males and females, sub-adult males and females.

Jolly-Seber abundance estimation
To estimate annual abundance of the four groups, we constructed Jolly-Seber (Jolly, 1965;Seber, 1965) open mark-recapture models using the POPAN (Schwarz and Arnason, 1996) formulation in Program MARK (White and Burnham, 1999). We implemented MARK using the RMark package (Laake, 2013) in the R software environment (R Core Team, 2018). The POPAN parameterization includes: φ i = probability of apparent survival from year i to i + 1, p i = probability of capture (observation) on occasion i, N super = super-population consisting of all sharks that were alive and available for detection at some point during the sampling periods of the study, β i = probability that a shark from the super-population, N super , would enter the population between occasion i and i + 1.
Jolly-Seber models use information from the number of unmarked individuals observed on each sampling occasion and assume that unmarked animals have the same capture probabilities as marked animals, which allows for the estimation of abundance (Schwarz and Arnason, 1996). We analyzed data from 2011 through 2018 when sharks at all three sampling sites were consistently sampled. Assumptions for this model include: (1) sharks retain their fin pattern throughout the study and are recorded and matched properly, (2) sampling is instantaneous (or realistically conducted over a short time), (3) survival probabilities are the same for marked and unmarked sharks between each sampling occasion (homogeneous survival), (4) catchability is the same for marked and unmarked sharks (homogeneous catchability), (5) the study area remains constant throughout the study period, and (6) any temporary emigration is random (Kendall et al., 1997) based on previous work with this dataset that assessed random versus Markovian temporary emigration (Chapple et al., 2016;Kanive et al., 2015). Multiple combinations of model structures for φ, p, and β were tested, and the super-population size (N super ) was estimated for each group in each model. Following evidence that φ increases significantly with total length (Kanive et al., 2019), we allowed φ to vary by maturity (mat) where adult males and females were pooled and sub-adult males and females were pooled or allowed to vary by maturity and time via an interaction (mat*t). For capture probability, p, we assessed the effects of sex, effort, time (t), group (g), as well as the following combinations: sex*effort, sex*t, g*t. We also defined group effects where p for sub-adult males and females was constrained to be the same but p for adults was allowed to vary by sex and year via an interaction [((subs(.)ads(sex))*t] or held constant across years [(subs(.)ads (sex)]. An alternative model structure for p was also considered in which p for sub-adult males and females was allowed to vary by sex while the rate for adult males and females was constrained to be the same but either (a) allowed to interact with year [((subs(sex)ads(.))*t] or (b) was held constant across years [(subs(sex)ads(.)]. For models where p, β, and φ were all allowed to vary over time, not all parameters in the model were estimable (Schwarz and Arnason, 1996). Therefore, for models where p and β were both time variant, the first and second capture probabilities were constrained to be equal for sharks of each sex (p 1 = p 2 ) and other capture probabilities were allowed to vary. To eliminate inestimable parameters in models where p and φ were both time variant, the last two capture probabilities were constrained to be equal for sharks of each sex (in this case p 7 = p 8 ) and other capture probabilities to vary. Probability of entry, β, was allowed to differ by group and was either held constant over time (g), allowed to vary independently over time in each group (g*t), or to vary by maturity class while either being held constant over time (mat) or allowed to vary independently over time in each maturity class (mat*t). The initial β, (β 0 ), was derived from subtraction (β 0 = 1 − β 1 − β 2 − … − β 7 ) representing the estimated proportion of the superpopulation present before sampling occurs. Finally, N super was estimated for each group (g) for every model. The models were fitted using the logit-link function for φ and p, the log-link function for N super , and the multinomial-link for β to constrain the set of parameters to have a sum that was ≤ 1.00 (Schwarz and Arnason, 1996;Williams et al., 2002). We considered models that contained all-possible combinations of parameterizations for each of the parameters, which yielded 88 different competing models. Model selection was based on Akaike's Information Criterion for small sample sizes (AICc) (Burnham and Anderson, 2002). Given issues of model selection uncertainty, we also used model averaging to obtain estimates for parameters of interest (Burnham and Anderson, 2002).
We derived the number of new animals (B i ) that joined the superpopulation during each year and the size of the annual population sizes (N t ) using standard methods for POPAN results (Schwarz and Arnason, 1996). To estimate total abundance for each year, we summed the group abundances for and quantified uncertainty using the delta method (Powell, 2007).

Population growth rate calculation
To obtain estimates of population growth rates (λ), we used parametric bootstrapping (Buckland and Garthwaite, 1991;Efron and Tibshirani, 1986) using the derived abundance estimates and associated variance from the top model and model averaged results for each of the P.E. Kanive et al. Biological Conservation 257 (2021) 109104 four groups of white sharks. Specifically, we generated 5000 random samples of annual abundance (N t ) for each group, for each year. Population growth rates (λ) were then calculated using the 5000 randomly generated estimates of annual abundance for each group from the top model. Specifically, from 2011 to 2018, we estimated λ by dividing subsequent abundance estimates (i.e. N 2012 /N 2011 , N 2013 /N 2012 , …, N 2018 /N 2017 ) yielding seven inter-annual estimated λ's representing the growth rate from one time step (year) to the next. For each of the 5000 sets of 7 estimated inter-annual λ's, we estimated a geometric mean λ that accounts for temporal variation in rates of annual population growth (Cooch and White, 2019;Lewontin and Cohen, 1969). These geometric means yielded approximately normal distributions and allowed us to calculate 95% confidence intervals for each group's estimated geometric mean.

Data summary
During the sampling seasons from 2011 to 2018, we logged 2587 h attracting White sharks to research vessels at Southeast Farallon Island (478 h), Año Nuevo Island (755 h) and Tomales Point (1354 h). A total of 419 individual White sharks (241 male; 147 female; 31 unknown) ranging from 200 to 550 cm total length (mean = 390 cm, SE = 3.71) were identified (https://zenodo.org/record/4022914#.X2IxTGdKjOR). Unique fin identifications were logged from 1533 images. The observed sex ratio of the sample population was 1.6 males per female. The number of unique and resights per unit effort is summarized in Table 1.
Of the sample population of 419 white sharks, 388 (241 males, 147 females) were confirmed to be either male or female at some point in the study period and 31 sharks remained unsexed. From the LOLASURVIV model (Supplementary materials S1), the proportion of males in the sample population was estimated to be 0.60 (SE = 0.02). Therefore, for this sample population to equate to 0.60 male, we randomly assigned 10 of 31 unsexed sharks to male and 21 as female resulting in 251 males (144 adults, 107 sub-adults) and 168 females (71 adults, 97 sub-adults) (Fig. 2).

Model results
Our best-supported model among a suite of 88 models had an AICc score with a difference of 2.63 units from the next best model and carried 54% of the AICc weight. The second and third best-supported models each carried 14% of the AICc weight (Table 2). Given the resulting model-selection uncertainty, we report the parameter estimates from the top model alongside model-averaged results obtained from all models that received any model weight (Table 2) (Burnham and Anderson, 2002). In the best-supported model, φ differed by maturity class but not by year and was estimated as 0.70 (SE = 0.03) for subadults and 0.85 (SE = 0.02) for adults. In the best-supported model, p was allowed to vary by sex in the mature class, but held constant by sex for sub-adults, and allowed to vary by year [((subs(.)ads(sex)) * time]. Time-varying estimates of p ranged from 0.21 (SE = 0.08) to 0.66 (SE = 0.08) for sub-adult sharks (both sexes), from 0.25 (SE = 0.05) to 0.87 (SE = 0.13) for adult males, and 0.17 (SE = 0.06) to 0.44 (SE = 0.12) for adult females. Estimated annual probability of entry into the population, β i was best estimated as being different by maturity class but constant over time was 0.09 (SE = 0.01) for adults and 0.10 (SE = 0.01) for subadults. The probability of being present at the start of the study, β 0 , was 0.34 (SE = 0.06) for adults and 0.28 (SE = 0.08) for sub-adults. The super-population (N super ) for each group based on results from the bestsupported model was 161.15 (SE = 13.56) for sub-adult males, 146.04 (SE = 12.57) for sub-adult females, 177.57 (SE = 7.99) for adult males, and 114.01 (SE = 12.50) for adult females.  (Fig. 3).

Annual abundance estimates
Seasonal total abundance estimates summed over all groups of subadult and adult male and female White sharks in our study area provided some evidence that overall abundance increased over the study (Fig. 3) Results were consistent for three versions of the dataset each with different random sex assignments for the 31 individuals of unknown sex. Across demographic groups, resulting estimates of annual abundance were within the 95% CI's of the presented results. Given the relatively small number of unknowns (7.4%), the modest influence of the assignment of unknowns is perhaps not surprising.

Estimates of population growth rate (λ)
Estimates of overall geometric mean λ and associated quantiles from the top model and from model averaging were calculated from empirical bootstrap distributions as 1.05 (95% C.I. = 0.91-1.16) and 1.07 (95% C. I. = 0.91-1.23) respectively. Of the 5000 bootstrap estimates generated for λ for the overall shark population, 83.4% of values from the bestsupported and 83.5% of values from model averaging were >1.00. Taken as a whole, the bootstrap results provide evidence that the population's stochastic growth rate, as measured by geometric mean λ was adequate for population growth over the study period. However, the uncertainty associated with the estimates does not rule out other possibilities.
Group-specific estimates of geometric mean λ from the top model provide strong evidence that adult males and females had a stochastic growth rate >1.00 (98.8% of estimates of bootstrapped values were >1.00 for each sex) indicating an increase in abundance over the study period. In contrast, results for sub-adults from the top model were more equivocal (68.5% and 68.6% of values for sub-adult males and females, respectively, were >1.00) (Fig. 4). Annual geometric mean estimates for each group from top model were calculated as 1.07 (95% C.I. = 1.01 to 1.13) for adult males, 1.07 (95% C.I. = 1.00 to 1.15) for adult females, 1.04 (95% C.I. = 0.90 to 1.17) for sub-adult males, and 1.04 (95% C.I. = 0.90 to 1.17) for sub-adult females. Model-averaged results were similar to those from the top model but slightly more uncertain (96.1% and 87.0% of estimates of bootstrapped values were >1.00 for adult males and females, respectively). For subadults, results from model-averaging provided slightly more evidence than results from the top model that the stochastic growth rate was >1.00; however, results from both had notable levels of uncertainty (71.4% and 78.3% of estimates of bootstrapped values were >1.00 for sub-adult males and females, respectively). Annual geometric mean estimates based on results from model averaging were 1.06 (95% C.I. = 0.99 to 1.13) for adult males, 1.06 (95% C.I. = 0.95 to 1.17) for adult females, 1.07 (95% C.I. = 0.85 to 1.29) for sub-adult males, and 1.08 (95% C.I. = 0.88 to 1.28) for sub-adult females. Consequently, we found considerable evidence for positive population growth in adult males and females during the 8-year study period and more equivocal evidence for positive population growth among sub-adult males and females, although a stable population growth rate can be conservatively inferred for all groups.

Discussion
Using visual mark-recapture techniques based on dorsal fin morphology, we estimated the regional abundance over eight years for sub-adult and adult white sharks that seasonally aggregate off central California. Importantly, our results are indicative of a relatively small group of white sharks (<300 individuals) along with evidence that the adult demographic appears to be stable or slightly increasing (see Fig. 3); however, the possibility of decline cannot be definitively ruled out. Our analysis of separate demographic classes suggests a potential modest increase in the abundance of adults (male and female) yet an equivocal estimate for sub-adults (also both male and females) ( Table 3). These findings represent the first empirical estimate of population trend of white sharks in central California.
Importantly, the abundance estimates and inferences of this study are specific to the region of central California encompassing the three primary aggregation sites sampled (Fig. 1). However, electronic tagging data have shown sharks tagged within our study area move to, and return from, coastal areas both to the north and south of the study area

Table 2
Model suite with Akaike weights >0.01 where K is the number of parameters estimated in the model, AICc is the Akaike Information Criterion score corrected for small sample size, ΔAICc is the difference in the AICc values from the top model, and w i are the calculated Akaike weights. Parameters are defined as follows: φ (apparent survival), p (capture probability), β (entry probability), and N (super-population). Notation for parameters is as follows: (mat) -time invariant, adults and sub-adults allowed to vary, (g) -time invariant, each demographic (group) allowed to vary, (sex*t) -time variant, sexes allowed to vary, (g*t) -time variant, each demographic (group) allowed to vary, (mat*t) -time variant, adults and sub-adults allowed to vary, [((subs(.)ads(sex))*t] -time variant, sub-adults grouped, adults vary by sex.  Jorgensen et al., 2010). Additionally, opportunistic dorsal fin photographs of white sharks taken outside the study area including Hawaii, Mexico and throughout California, have been matched to white sharks identified in central California. The hypothesis that the central California white shark population may comprise a substantial proportion of 'transient' individuals  has been robustly ruled out (Kanive et al., 2019) based on empirical data testing (Pradel et al., 1997). Whether this potential regional increase of reproductive-age white sharks in central California is representative of an overall increase in the total northeastern Pacific white shark population (Jorgensen et al., 2010) or a result of regional fluxes in density remains to be determined. A trend estimate for the entire northeastern Pacific will require integration of data across their known aggregation areas along the west coast of North America. Although our results contain uncertainty as outlined above, there are two key protections that may have contributed to the success of white sharks in central California. First, white sharks were first protected in California waters in 1994 and federally in 2004 (CITES, appendix ii), and bans on coastal gillnetting in 1990 (California Proposition 132) likely reduced early incidental mortality of juveniles. Secondly, pinnipeds, their main prey, were protected under the Marine Mammal Protection Act in 1972 and have greatly increased in abundance since (Boeuf and Laws, 1994;Laake et al., 2018;Lowry, 2014). Taken together, these protections may influence the prospect for population increase of this important marine predator.
Our findings provide evidence that mature individuals may be increasing in this region, while the probability of an increase in subadults remains less clear. Ancillary information from increased sightings of YOY and juvenile sharks (Lowe et al., 2012) and an increase in sea otter mortality due to white shark bites (Tinker et al., 2016) likely inflicted by both adults and juveniles (Moxley et al., 2019), provide support to our findings. An increase in incidental fisheries landings of juvenile white sharks between the mid 1990's and 2008, following the 1994 California near shore gill net ban has frequently been cited as evidence of population increase (Lowe et al., 2012). However, this potential abundance index, which depends on consistent reporting by fishers, fell to zero for many years following the 2012 petitions for this segment of the white shark population listing under the Endangered Species Act. Public reports of juvenile white sharks sightings in California have also increased substantially in recent years, but the number of sightings has been highly variable in magnitude and location year-toyear, and has been linked to climatic changes in optimal temperature habitat across their latitudinal range (Tanaka et al., 2021;White et al., 2019). These observations suggest a possible increase in the abundance of white sharks in the juvenile size class, smaller than the sizes sampled in this study area. However, distributional changes, climate variation, effort consistency, social media, and fisheries reporting accuracy are confounding factors that need to be resolved before firm conclusions can be made.
Among adults and sub-adults, length-based apparent survival rate estimates show that adults have a substantially higher probability of surviving and reusing sites in the study area (Kanive et al., 2019). The 'apparent' term means that mortality is confounded with potential permanent emigration, where the animal is alive and leaves the study area for the remainder of the study period. As sub-adults recruit to prime foraging grounds occupied by adults, intraspecific competition could force out smaller, less competitive individuals (Kanive et al., 2019). It is therefore possible these sub-adult sharks are actually alive but dispersed outside of the prime hunting grounds of our study area. Future research should examine the coastal movement of sub-adult white sharks tagged at the three aggregation sites described in this study to determine if they remain alive and outside of the study area. Resolution of the movement of smaller and less competitive sub-adults could reveal whether density dependent recruitment is a factor at prime adult aggregation sites (Kanive et al., 2019). This is the first empirical assessment of regional population trend for white sharks in the northeastern Pacific. A similar regional increase in white shark numbers in the US Economic Exclusive Zone of the northwest Atlantic showed a recovery towards early historical numbers (Curtis et al., 2014). While there are no historic baselines of white shark abundance that preclude management policies implemented in 1994 in the Pacific northeast, the results of this study suggest recent favorable management in California (including marine mammal protection, and white shark bycatch mitigation) may be contributing to a pathway for the positive population growth for this important marine top predator. Given the vulnerable life history traits of white sharks, insuring their future will take a continued effort to monitor populations and identify potential future threats such as climate change and unregulated highseas fisheries.

CRediT authorship contribution statement
Paul Kanive, Salvador Jorgensen, Taylor Chapple, and Jay Rotella conceived the ideas and designed methodology; Paul Kanive, Salvador Jorgensen, Taylor Chapple, Scot Anderson, Timothy White and Barbara Block collected the data; Paul Kanive and Jay Rotella analyzed the data; Paul Kanive led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Declaration of competing interest
The authors declare no conflict of interests.  Table 3 Estimated seasonal geometric mean of λ with associated 95% confidence intervals for each group of white shark as well as the overall estimate. permission to salvage marine mammal parts from the National Marine Fisheries Service.

Data availability statement
Data will be made available directly from online repository, https://zenodo.org/record/4022914#.X2IxTGdKjOR.