Effective number of white shark (Carcharodon carcharias, Linnaeus) breeders is stable over four successive years in the population adjacent to eastern Australia and New Zealand

Abstract Population size is a central parameter for conservation; however, monitoring abundance is often problematic for threatened marine species. Despite substantial investment in research, many marine species remain data‐poor presenting barriers to the evaluation of conservation management outcomes and the modeling of future solutions. Such is the case for the white shark (Carcharodon carcharias), a highly mobile apex predator for whom recent and substantial population declines have been recorded in many globally distributed populations. Here, we estimate the effective number of breeders that successfully contribute offspring in one reproductive cycle (Nb) to provide a snapshot of recent reproductive effort in an east Australian–New Zealand population of white shark. Nb was estimated over four consecutive age cohorts (2010, 2011, 2012, and 2013) using two genetic estimators (linkage disequilibrium; LD and sibship assignment; SA) based on genetic data derived from two types of genetic markers (single nucleotide polymorphisms; SNPs and microsatellite loci). While estimates of Nb using different marker types produced comparable estimates, microsatellite loci were the least precise. The LD and SA estimates of Nb within cohorts using SNPs were comparable; for example, the 2013 age cohort Nb(SA) was 289 (95% CI 200–461) and Nb(LD) was 208.5 (95% CI 116.4–712.7). We show that over the time period studied, Nb was stable and ranged between 206.1 (SD ± 45.9) and 252.0 (SD ± 46.7) per year using a combined estimate of Nb(LD+SA) from SNP loci. In addition, a simulation approach showed that in this population the effective population size (Ne) per generation can be expected to be larger than Nb per reproductive cycle. This study demonstrates how breeding population size can be monitored over time to provide insight into the effectiveness of recovery and conservation measures for the white shark, where the methods described here may be applicable to other data‐poor species of conservation concern.


| INTRODUC TI ON
Assessing the size of natural populations is a key objective of monitoring programs which are vital for understanding the conservation status of species, the regulating effects of biotic and abiotic factors, and for the assessment of management efforts (Lindenmayer et al., 2020). However, for many marine populations there is a lack of consistent monitoring programs at appropriate spatial and temporal scales for conservation and policy needs (Papa et al., 2020). This presents a significant problem for chondrichthyans (sharks, skates, rays and chimaeras), where more than half of known species are characterized by insufficient data and one-quarter are estimated to be at risk of extinction . Within the elasmobranchs (sharks, skates, and rays), each contributes significantly to connect ecosystems and regulate marine food webs (Heupel et al., 2014).
However, habitat loss and continued pressures on mortality though bycatch and targeted fishing have resulted in many populations of elasmobranch being depleted at a rate that exceeds their natural recovery potential (Worm et al., 2013). Given the significant challenges facing elasmobranchs and the importance of their role in regulating marine ecosystems, improvements for monitoring changes in natural populations are critical.
Monitoring threatened elasmobranch species is particularly challenging for many reasons. In the case of the white shark, Carcharodon carcharias (Linnaeus, 1758), where monitoring is a both a social and conservation priority, efforts to evaluate long-term population trends have been hampered by issues including detectability [misidentification in photo-ID surveys (Burgess et al., 2014), lack of resightings in mark-recapture studies (Gore et al., 2016), effects of environment on heterogeneity of behavior (Jacoby et al., 2012)] and a lack of catch statistics (Roff et al., 2018). The need for alternate methods to index shark populations has therefore led to the increasing use of molecular markers to evaluate change and inform management (Blower et al., 2012;Bruce et al., 2018;Hillary et al., 2018). In this study, we focus on the concept of genetic effective population size (herein effective population size-Ne), which can be used to evaluate change in abundance from allele frequencies (Schwartz et al., 2007). When populations are small, genetic models predict that the evolutionary force of genetic drift (stochastic changes in allele frequencies) will predominate over other evolutionary forces such as natural selection, to reduce genetic diversity, population viability, and evolutionary potential (Frankham, 1996;Franklin, 1980). The extent to which a population is vulnerable to such effects is inversely related to the magnitude of Ne, where the effects of drift will occur more slowly in populations with larger effective sizes than those with smaller effective sizes (Wang, 2005). When a genetic sample contains only individuals from a single age cohort (a group of individuals having the same age-class), then the estimate of Ne corresponds to the effective number of breeders (Nb) which contributed offspring to that cohort Waples et al., 2013). For long-lived, iteroparous species, estimates of Nb are generally considered more useful for monitoring as they apply to a single breeding season and represent an accessible parameter for monitoring population trends at ecological timescales most relevant to conservation and management needs Schwartz et al., 2007;Waples & Do, 2008). Past research has confirmed the power and usefulness of Nb as a tool to monitor population trends (Antao et al., 2011;Luikart et al., 2020;Nunziata & Weisrock, 2018). For instance, quantifying changes in Nb over time provides high power to detect declines in Nb (Luikart et al., 2020) and has helped to identify factors relevant to shaping populations (i.e., management interventions, demographic parameters) with successful outcomes reported for populations of commercially important bony fishes. Examples include salmon (Bacles et al., 2018;Perrier et al., 2016), trout (Ruzzante et al., 2019;Whiteley et al., 2013;Wood et al., 2014), snapper (Jones et al., 2019), and tuna (Waples et al., 2018). In these examples, both Nb and Ne were used to investigate demographic (i.e., variance in reproductive success under commercial harvest conditions) and environmental (i.e., stream productivity, competition, habitat quality, year-of-theyoung development) effects on long-term population viability, with significant implications for management and conservation.
In this study, we trialed a sampling and genotyping protocol aimed at estimating Nb over time (four breeding seasons; 2010-2013) in a population of C. carcharias of conservation concern.
We focus on the east Australia-New Zealand population (EAP) of C. carcharias which, due to patterns of coastal residency and site fidelity (Bruce et al., 2019;Spaet et al., 2020), is genetically distinct from other identified populations in the North-Pacific, South-West Australia, Atlantic, South Africa, and Mediterranean (Andreotti et al., 2016;Blower et al., 2012;Gubili et al., 2010;O'Leary et al., 2015;Tanaka et al., 2011). The EAP has experienced a large (>90%) decline during the 20th century due to targeted fishing and mortalities associated with bather protection programs (Reid et al., 2011;Roff et al., 2018); however, recovery is now anticipated due to protection through international conventions and jurisdictional legislation [i.e., International Plan of Action for the Conservation and Management of Sharks (FAO, 2000) and the Environment Protection and Biodiversity Conservation (EPBC) Act of 1999(EPBC, 1999]. Previous efforts to detect population recovery using historical catch data (Roff et al., 2018) and genetic close-kin mark-recapture Hillary et al., 2018) found no significant evidence of population growth or recovery in the EAP. Updated bather protection programs along parts of east coast Australia (i.e., 'Shark Management Alert in Real Time' (SMART) drumlines, see Tate et al., 2019), aimed at minimizing unfavorable interactions with marine environment users, offer conservation, effective number of breeders, population genetics, effective population size, monitoring an opportunity for nonlethal tissue sampling and to determine the usefulness of this genetic monitoring method in the EAP.
Our specific objectives were to: (a) use two genetic methodologies to estimate Nb over time in the EAP [sibship assignment (SA) (Wang, 2009) and linkage disequilibrium (LD) (Hill, 1974, p. 197;Waples, 2006)]; (b) validate these results using two types of nuclear genetic markers (single nucleotide polymorphisms and microsatellites); (c) investigate Nb/N ratios using published estimates of the adult population size (Na); and (d) develop expectations for generational Ne in the EAP using life-history information and simulations. Our results for the EAP of C. carcahrias suggest that Nb has not changed significantly year-to-year and provides insight into the effectiveness of recovery and conservation measures following historical declines.

| Samples
To obtain genetic data to estimate Nb in the east coast Australia-New Zealand population (EAP) of C. carcharias, tissue samples (n = 247) were nonlethally collected during 2015 to 2018 from juvenile and subadult C. carcharias between Buckley Beach, Narrawallee (−35.29873, 150.48331), and Seven Mile Beach, Lennox Head (−28.76130, 153.62020) (Figure 1). Individuals were captured, restrained, tagged, and released as part of the New South Wales (Australia) Shark Management Strategy. Fin clips were collected for genetic purposes and fork length (FL) and total length (TL) measurements were taken from each individual. Since migration between populations can bias genetic estimates of both Ne and Nb (Macbeth et al., 2011), the population of origin for each individual was resolved through the inclusion of tissue samples of white sharks collected from other locations (Western Australia n = 3; South Australia n = 9; South Africa n = 20; total n including EAP samples = 279, see Table S1). All samples were used in the SNP discovery and genotyping pipeline.

| Cohort assignment
To group individuals into age cohorts, a year-of-birth was assigned to each sample using the year the individual was sampled minus the age of the individual in that given year. To estimate the age of individuals, we used the von Bertalanffy growth function (VBGF)

| SNP and microsatellite loci datasets
DNA was extracted from all samples (n = 279) using a standard salt precipitation procedure. For SNP data, the samples were genotyped by DArT P/L laboratory using DArTseq TM technology (Kilian et al., 2012). Sequencing steps followed Kilian et al., (2012) and were completed using an Illumina HiSeq 2500. Resulting sequences were processed using the proprietary DArT analytical software, DArTsoft14. DArTsoft14 uses technical sample replicates to optimize its algorithm parameters and ensure scoring consistency (see Georges et al., 2018). Postprocessing of SNPs was completed in R (R Core Team, 2018) using the R-Package radiator 0.0.5 (Gosselin, 2017) and custom R-scripts following current best practice (O'Leary et al., 2018;Shafer et al., 2017). A twostage postprocessing approach was applied to the SNP dataset to identify and remove 1) migrants and 2) outlier loci (candidate loci under selection, so as to retain only nonselective neutral loci).
SNP data representing all samples (East Australia, South Africa,  Table S3.1 (Supplementray Materials S3), and subsequently used for sample population assignment and initial outlier loci discovery. Strongly divergent individuals create strong mixture LD which downwardly bias estimates of Ne(LD) (Waples & England, 2011) and may contribute to upward bias in estimates using the SA method (Ackerman et al., 2017). To identify divergent individuals, we performed a discriminant analysis of principal components (DAPC) (Jombart et al., 2010) implemented in the R-package adegenet (Jombart et al., 2010). The optimal number of discriminant functions to retain was calculated using the function xvalDAPC using 80% of the data in the training set, and the number of PCs retained in the final DAPC was associated with the lowest mean squared error. As indicated in Figure S3.1, two samples collected from east Australia appeared distinct from other EAP samples (subsequently confirmed using tracking data from acoustic tagging, Spaet et al., 2020). These samples were removed from subsequent analysis. We also performed tests for putitative loci under selection which deviate from the assumptions necessary for estimating Ne (Waples & England, 2011). We used pcadapt (Luu et al., 2017) which identifies outlier loci in a multidimensional space (we used k = 3 principal components). We removed loci when the q-value (test statistic) was smaller than the false discovery rate ( = 0.05). In the second stage of SNP postprocessing, we used a dataset (herein Dataset-2) containing all SNP loci except those identified as selective outliers and including only samples representing genotypes of EAP origin only. We then filtered Dataset-2 using reproducibility greater than 98%, a minor-allelecount greater than three, coverage (minimum 5, maximum 25), retained only one SNP per locus and removed individuals missing greater than 20% of SNP loci. Loci were further removed where Hardy-Weinberg disequilibrium mid-p was less than 0.01 < 0. and if F IS was greater than or equal to +0.5 or less than or equal to −0.5 (see Table S3.2). Dataset-2 was then used to make estimates of Nb.

| Estimation of Nb
Two methods were used to estimate Nb from data derived from either SNP or microsatellite loci: (a) the linkage disequilibrium method (LD) (Hill, 1974;Waples, 2006) and (b) the sibship assignment method (SA) (Wang, 2009). Estimates are referred to as Nb(LD) and Nb(SA). Broadly, the LD method determines the size of the parental generation using a measure of the genetic association (or LD) in a given age cohort. In finite populations, random genetic drift leads to associations of alleles at different loci. The LD method uses the extent of nonrandom association between alleles at different loci to estimate genetic Ne and reflects the inbreeding Ne when loci are unlinked (Hill, 1981;Waples & Do, 2010). The formulation of the LD method uses the observed average disequilibrium between pairs of independent (i.e., nonlinked), neutral loci in a sample of individuals taken from a single, isolated, randomly mating population. Estimates of Nb(LD) are based on the theoretical relationship between r 2 and Ne as described in Hill (1981); Equation 2a from (Waples & Do, 2010) where r 2 is the mean squared correlation of allele frequencies at different gene loci adjusted for sampling error (i.e., the observed average disequilibrium) and S is the number of individuals sampled.
We implemented this method using the program NeEstimator v2.1 (Do et al., 2014). In contrast, the SA method uses the direct relationship between genetic relatedness and inbreeding Ne, such that any two individuals sampled randomly from a population with a small Ne will have a higher probability of sharing the same parent or parents (Wang, 2009). The SA method (Wang, 2009) determines the size of the parental generation by estimating the probability that dyad relationships are either full-or half-siblings in a sample from the same cohort, sharing two, one, or zero parents, respectively; Equation 10 from Wang (2009) where Q 1, Q 2 , and Q 3 are the paternal, maternal half-sibs, and fullsibs, respectively, N 1 and N 2 are the number of male and female parents, and is a measurement of the deviation from Hardy-Weinberg proportions in genotype frequencies (Wang, 2009). The SA method was implemented in the program COLONY (Wang, 2009).
Both Nb(LD) and Nb(SA) were estimated for the EAP across four year-of-birth cohorts (2010,2011,2012,2013) where sample size per cohort was greater than 25 individuals. Estimates of Nb were made using either SNP or microsatellite marker data. To estimate Nb(LD) with NeEstimator v2.1 (Do et al., 2014), a random mating model was specified, rare alleles which upwardly bias estimates were excluded using the criterion PCrit= Pcrit0.05 as recommended in Waples and Do (2010), and jackknife confidence intervals that accouns for pseudo-replication due to physical linkage and overlapping loci pairs were used Waples & Do, 2010

| Inference of Nb/Na ratios
To calculate Nb∕Na ratios, we used Na as described in Bruce et al., (2018), where Na is the number of adults in the population.
As C. carcharias has a low intrinsic capacity for population increase, low fecundity, and low lifetime variance in reproductive success (Bruce, 2008), the Na estimates from Bruce et al. (2018) are assumed to apply to the time period corresponding to our study; Na = 750 with an uncertainty range 470 to 1,030 . Our estimates of Nb(LD) and Nb(SA) were combined (herein Nb(LD+SA) to provide a single value of Nb with which to infer Nb/Na ratios. Nb (LD) and Nb(SA) were combined by taking the harmonic mean of the two values, weighted by the inverse of their variances as suggested in previous studies (see Waples & Do, 2010); see Appendix 1. In our study, the differences between the estimates from the LD and SA methods were not overly large, so using a combined estimate of Nb to determine the Nb/Na ratio would not change the conclusions described herein. Furthermore, in our study the SA estimate generally has a lower variance and provided a less bias and more precise estimate of Ne and therefore contributed more (around 2/3 rd ) to the final combined estimate. We also note that when two methods with approximately comparable performance provide an estimate of Ne, then the variance of the combined estimate will be smaller than for either estimate alone (Waples et al., 2016).

| Expectations for Ne
To develop expectations for generational Ne in the EAP of C. carcharias, we use a simulation-based approach. This route was taken as the assumptions of single-sample genetic estimators of Ne, including LD and SA methods used herein, dictate that data used to make estimates represent a random sample of a population across an entire generation (Hare et al., 2011). Since the white shark is long-lived and samples in this study were mostly juvenile or subadults, we instead characterize the expected Nb/Ne ratio using simulations based on published methods and parameterized using the life-history of white shark. This indirectly allowed the inference of an expected Nb/Ne ratio to permit a better understanding of inbreeding and implied fitness of the population. We use both deterministic and forward-time population simulations following methods described in , to determine Ne and Nb. First, we implemented the discrete-time, deterministic hybrid Felsenstein-Hill method for calculating Ne in iteroparous species . The model was implemented in the software AgeNe , herein Nb(ageNe), and parameterized using life-history information from white sharks in the EAP (Supplementray Materials S7). Furthermore, since the Felsenstein-Hill method assumes the probability of reproduction is not affected by events in previous time periods, we also use forward-time population simulations implemented in simuPOP (Peng & Kimmel, 2005), to create a single, isolated, randomly mating population to further characterize the Nb/Ne ratio under two intermittent-breeding scenarios as in . Each simulation was parameterized using outputs from AgeNe, including total population size and stable age distribution in the population, given the specified vital rates and a specified number of offspring produced per cycle that survived to age 1 (N1), here N1 = 1,000.
Each individual was represented by 100 microsatellite-like loci, each having 10 possible allelic states, no mutation, and data were tracked for 50 reproductive cycles after a burn-in period of 50 cycles. We forced a number of females to skip either zero, one, or two cycles of breeding (a proportion of females) hypothesized in this species (Domeier & Nasby-Lucas, 2013;Mollet & Cailliet, 2002). Intermittent or skipped breeding occurs when sexually mature adults skip breeding opportunities (Last & Stevens, ;Shaw & Levin, 2011), in this case likely due to the costs of reproduction or prolonged gestation period in females (Bruce, 2008). We directly calculated mean (k) and variance (Vk) lifetime reproductive success, and Ne and Nb directly from simulation demographic data (not genetic data) for each reproduc-

| Cohort assignment
Using the relationship between TL and age, we found that one in- As low sample sizes can bias estimates of Nb using the methods of this study, only age cohorts containing greater than 25 samples were used (cohorts 2010, 2011, 2012 and 2013).

| SNP and microsatellite loci data
The DArTsoft14  (loci Cca1419, Cca1072, CcSA2). These markers were removed from further analysis (LD method only). One locus (CcSA5) was not polymorphic (see Table S6.1) and was also excluded. Per individual, 97% had no missing loci while the remaining 3% of samples had three or less missing loci.

| Estimates of Nb
Using SNP data, Nb estimates per year-of-birth cohort were similar between the LD and SA methods and had overlapping 95% confidence intervals ( Table 1). Estimates of Nb(SA) were not sensitive to changes in model parameters such as the sibship prior, inbreeding settings, error rate, and polygamy settings ( Comparing between the SA and LD method using data from microsatellite loci, estimates of Nb(LD MSAT ) were higher than the equiv-   Table S7.2).

| D ISCUSS I ON
Monitoring Ne can inform management decisions in populations of conservation concern, where Nb is analogous to Ne except that it represents the effective number of breeders per year rather than per generation. Using data from SNP and microsatellite loci and two sin- Indeed, as past studies have shown, monitoring Nb for as few as five consecutive reproductive cycles could be used to detect change in Nb (declines), even in species with long generation intervals (Antao et al.,2011;Wang, 2005;Luikart et al., 2020) such as the white shark, with implications for both N and Ne in some circumstances. For example, if Nb were to decline significantly for multiple reproductive cycles, then both Ne and Nc may be affected (Luikart et al., 2020).
Alternatively, if Nb/N ratios are observed to be stable over many generations, then it has been suggested that N (or Na) may be inferred from Nb, which would be of use to population monitoring and evaluation of conservation and management actions (Luikart et al., 2020). To this end, we recommend using Nb to track year-to-year changes in the effective number of breeders as a timely assessment of population status over time to provide insights into the effects of current management actions and co-occurrences such as environmental changes. As in this study, future tissue samples for Nb monitoring could be obtained as part of existing bather protection programs (i.e., SMART drumlines; see Tate et al., 2019).
In this study, we used two genetic marker types (SNPs and microsatellites) and two single-sample genetic estimators of effective  (Table 1). This has been noted in previous studies (Ackerman et al., 2017;Wang, 2009) which have demonstrated that false sibships (type I errors) occur with a higher frequency compared to false nonsibships (type II errors) when either genetic information or true sibship within a sample is insufficient (i.e., few loci, low polymorphism, small sample size relative to total population size, low inclusion of siblings). Nb estimated using SNPs differed between methods such that Nb(LD) was lower compared to Nb(SA), although differences were not significant having overlapping CIs. Nb(LD) estimated using SNPs showed those cohorts with larger numbers of samples (i.e., 2013) provided more precise estimates, a result expected given genetic methods for estimating contemporary effective size depend on signals that are proportional to 1/Ne (Waples et al., , 2018. Monitoring studies are often focused on the number of individuals in a population; however, the relationship between effective size and population size (i.e., Ne/Na, Nb/N) is important to understand since genetic drift results in the loss of neutral genetic variation at a rate rate inversely proportional to Ne per generation, not N (Wright, 1931) and can therefore be useful for examining how different ecological factors influence genetic variation (Nunney, 1996). In this study, the ratio of Nb/Na was approximately 1/3 of mature adults for a single reproductive cycle, where mature adults represent different adult age classes. This is somewhat comparable to ratios inferred for other Carcharhinidae, including C. plumbeus (sandbar shark) in Delaware Bay, North Atlantic, which ranges between 0.50 (95% CI 0.45) and 0.63(95% CI 0.57) (Portnoy et al., 2009). Since a ratio of Nb/Na applies to a single reproductive cycle; when ratios are close to 1, we can infer that the majority of the adult population contribute to the next generation and that the offspring number per adult approaches the standard scenario of binomial distribution (Hedgecock, 1994). In contrast, when ratios are < 1, we can infer there is some deviation from the ideal (Hare et al., 2011). A number of factors will affect this relationship (Ne/Na, Nb/Na) including fluctuations in population size and several important life-history factors that change variance in reproductive success (e.g. mating system, generation time, sex difference including sex ratio, survival, recruitment age). In one case, Nb may be expected to be reduced relative to Na if females with high fecundity skip reproductive cycles after giving birth, resulting in different females breeding in different cycles (Waples & Anato, 2014). This should decrease both lifetime Vk and Nb, while increasing Ne. The ratio reported herein appears to be consistent with expectations for the breeding behavior of C. carcharias, suspected to undergo intermittent breeding (Bruce, 2008). Observations suggest the gestation period of C. carcharias females may approach 18 months from fertilization to parturition (Bruce, 2008;Mollet et al., 2000), resulting in the unavailability of a portion of adult females to produce offspring each cycle. However, we caveat that Nb/Na ratios determined in this study used estimates of Na from best available information from close-kin-mark-recapture estimates for the EAP in 2017 (see Bruce et al., 2018). Thus we assumed temporal stability of N over time; an assumption which would be violated if N has increased (or decreased) over the time period to which our Nb estimates apply (2010-2013).
Since neutral genetic variation is lost at a rate of 1/2Ne per generation (Wright, 1931), even numerically large populations can be at genetic risk if Ne is small (Waples et al., 2018). Although important, due to sampling restriction (i.e., difficulty sampling across a generation as required by estimators) and uncertainty of breeding histories, we could not estimate Ne directly nor did we consider the linear relationship between Nb and Ne which requires either true or estimated Nb/Ne to be quantified (Waples et al., 2013). Blower et al., (2012)   . This result is important as it suggests the study population in the EAP at least exceeds the inbreeding avoidance goal (Nb100) (Frankham et al., 2014). Ne > 100 (previously the Ne-50 rule) describes the short term goal required to avoid inbreeding, which results in excess homozygosity for deleterious and recessive alleles, leading to inbreeding depression and reduced fitness (Frankham et al., 2014). However, in relation to the long-term viable population benchmark, Ne 1,000 (previously Ne > 500) (Frankham et al., 2014) our results are less certain, and we caveat that this rule (Ne > 1000) refers to the loss of additive genetic variation that may negatively effect adaptation in response to changes in selective regimes, not inbreeding effective size as estimated herein. We suggest any genetic effects of a recently and significantly reduced population size in the EAP, such as a decline in Ne or loss of heterozygosity, may not be fully realized until adequate benchmark studies can be completed (i.e., historical or ancient DNA). However, genetic bottlenecks in white sharks have been recorded elsewhere (O'Leary et al., 2015). Given this, together with the lack of evidence from other studies to date of an expected recovery (except see Department of Primary Industries, 2019), our results emphasize the importance of continued monitoring, improved protections, and interventions to reduce mortality.
Indeed, the vulnerability of chondrichthyan fishes to exploitation has been comprehensively documented (Hutchings et al., 2012) and relative to other marine fish, the intrinsic capacity for population increase and rebound potential in white shark is low (Cortés, 2002) (i.e., long-lived, late age to maturity, high juvenile survival). In addition, shark species often travel large distances and use different habitats throughout their lives (Fujioka & Halpin, 2014), where they may be vulnerable to environmental changes (density, food availability, climate, illegal fishing). Regrettably, mortalities continue to occur in the EAP driven by action taken to mitigate human-shark in- along the east coast of Australia, which has been lethal for sharks despite catch-and-release programs (Roff et al., 2018). The recent modeling of the recovery of the North West Atlantic white shark population provides a useful principal in this regard; "every fish counts" (Bowlby & Gibson, 2020, p.9).

| CON CLUS ION
We

ACK N OWLED G M ENTS
Project funding and primary support for sample collection were pro- useful comments, and we extend our sincere thanks to the reviewers of this manuscript.

WO R K ED E X A M PLE
This is a worked example using data from the cohort 2013. This weighted mean will give the lowest variance of any weighted mean of the values. As with nearly all Ne calculations the harmonic mean must be used as the real quantities of interest are proportional to 1/Ne.

The Variances
Unfortunately, neither COLONY or NeEstimator (Do et al., 2014;Jones & Wang, 2010) provides the raw variance figures required; however, these can be approximated by working backwards from the provided confidence intervals.

COLONY
Colony generates 95% confidence intervals using the following for-

mula [cite]:
where V * is the variance of the estimate of 1/2N e . Knowing the upper and lower bounds of this confidence interval, we can estimate V * as, or, where L * and L are the lower confidence bounds in terms of 1/2N e and N e respectively. An identical argument follows for the upper bounds.
However, we desire the variance of N e , which we can approximate using a first order Taylor

N E E S TI M ATO R
Similar to COLONY Ne Estimator does not provide raw variances, and we need to work in terms of the confidence intervals for r 2 − drift which we will call here r*. The drift term is ~ 1/S, where S is the sample size. The 95% confidence interval for r* is explicitly normal in the case of the jackknife confidence interval . Thus, V*, the variance of r * is approximated by where U and L are the upper and lower bounds provided for r * .
However, again, we wish to have the variance in terms of N e . Using the same approach as for COLONY, we will approximate this using

WE I G HTE D M E A N A N D FI N A L VA R I A N CE
Now we follow the formula for the weighted harmonic mean, where the weights, w i , sum to 1. In this case, we need to normalize the inverse variances to sum to 1. , r * can be found in Table 1 in Waples and Do (2008), here we use the simpler original form (Equation 1 in Jones et al., 2016 and others) as our f (X) for this estimate, that is Working as before, However, we still need to obtain U * and L * the upper and lower confidence bounds for r*. This can be achieved by inverting the function of r * and S from Table 1 in Waples and Do (2008). In this particular case, we cannot use the simple approximation. In our case, the inverted function is N e = 1 3(r 2 − drift) = 1 3(r * ) .