Close-kin mark-recapture abundance estimation:
practical insights and lessons learned

Verena Trenkel, Grégory Charrier, Pascal Lorance, Mark Bravington

► To cite this version:

Verena Trenkel, Grégory Charrier, Pascal Lorance, Mark Bravington. Close-kin mark-recapture abundance estimation: practical insights and lessons learned. ICES Journal of Marine Science, 2022, 79 (2), pp.413-422. 10.1093/icesjms/fsac002. hal-03760486

HAL Id: hal-03760486

https://hal.inrae.fr/hal-03760486

Submitted on 25 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Trenkel, V. M., Charrier, G., Lorance, P., and Bravington, M. V. Close-kin mark-recapture abundance estimation: practical insights and lessons learned

Verena M. Trenkel ID<fnref n="1" />*, Grégory Charrier<fnref n="2" />, Pascal Lorance ID<fnref n="1" /> and Mark V. Bravington<fnref n="3" />

We present practical lessons learned from applying the recent close-kin mark-recapture (CKMR) abundance estimation method to thornback ray (Raja clavata). For CKMR, related individuals are identified from their genotypes and their number and pattern is used for abundance estimation. We genotyped over 7000 individuals collected in the Bay of Biscay using Single Nucleotide Polymorphism (SNP) markers finding 99 parent-offspring pairs. The estimated number of adult thornback rays in the central Bay of Biscay was around 135000 (CV 0.19) in 2013. In total, four lessons were drawn: (i) CKMR helps identifying metapopulation structure, which if ignored might affect abundance estimates and/or time trends. There was strong evidence for two distinct local populations of thornback ray with no demographic connectivity. (ii) Demographic sample composition can affect precision and needs to include a range of birth years, which turned out to be difficult for thornback ray. (iii) Reasonable age information for potential offspring is essential. (iv) The sex of potential parents is needed and might be identified from sex-related SNPs. Reliable abundance estimation by CKMR appears feasible for a wide range of species provided that: sampling adequately covers potential local population structure, has appropriate demographic composition, and the age of potential offspring is reasonably well-known.

**Keywords:** abundance estimation, Bay of Biscay, CKMR, metapopulation structure, parent-offspring pairs.

## Introduction

Reliable abundance estimates are essential for terrestrial wildlife management and the sustainable exploitation of marine living resources. Abundance estimates can be obtained using as input visual count and cull data (Trenkel et al., 2000), scientific trawling derived biomass indices (Hilborn and Walters, 1992; Marandel et al., 2016, 2019), or from fishery-derived information (Hilborn and Walters, 1992; Augustin et al., 2013). A fundamental problem in this context is the unknown detection probability underlying counts or trawl abundance indices, e.g. in aerial line transects (Trenkel et al., 1997), or when trawling (Krieger and Sigler, 1995; Trenkel and Skaug, 2005; Kotwicki et al., 2015).

The recently introduced close-kin mark-recapture (CKMR) method for absolute abundance estimation leads to a change in paradigm (Bravington et al., 2016b). Instead of surveying large areas of land or sea for counting individuals, a number of individuals are sampled, genotyped, and compared pairwise to look for specific types of kin, for example parent-offspring pairs (POPs). The rest is classical mark-recapture abundance estimation with the only difference that individuals are not physically marked, released, and recaptured, but rather mature individuals are “recaptured” either directly via sampling themselves and their offspring (Bravington et al., 2016b) or indirectly by sampling at least two offspring (Hillary et al., 2018). A key benefit is that there is no need to release an individual alive, as there is in classical mark-recapture. Thus,

sampling can be lethal, e.g. from fishery catches or cull data. Because an offspring's two parents were by definition alive at the offspring's birth (or conception), CKMR provides information about the abundance of mature individuals in the years of birth of the sampled offspring. The alternative view of recapturing offspring via their parents leads to the same result (Ruzzante et al., 2019).

Modern genomics provides the means to identify related individuals within a sample (e.g. Thompson, 2000; Anderson and Garza, 2006). Indeed, as each parent contributes exactly half of their DNA to their offspring in diploid species (i.e. for each locus, the descendant gets one copy from each of its parents), POPs have a distinct and unique genotype pattern which can be used to identify this kinship type (e.g. Bravington et al., 2016b). Full-Sibling Pairs (FSP), where both individuals have the same father and mother, on average share as much co-inherited (ibd; identical-by-descent) DNA as a POP, but with a different pattern that can be distinguished statistically, given enough loci. Half-siblings (HSPs) are less closely related; on average only 50% of their loci will have an ibd allele. Such second-order kin can still be distinguished reliably from weaker kin given reasonable numbers of loci (e.g. ~2000 Single Nucleotide Polymorphism (SNP) markers; Hillary et al. (2018)), but cannot be distinguished from grandparent-grandchild pairs (GGP) or full-thiatic pairs (FTPs) such as aunt-nephew, uncle-niece, etc.; i.e. where both parents of one animal are also grandparents of the other animal.

<footnote n="1">¹DECOD (Ecosystem Dynamics and Sustainability), IFREMER, INRAe, Institut-Agro - Agrocampus Ouest, rue de l'île d'Yeu, 44311 Nantes cedex 3, France</footnote>

<footnote n="2">²Univ. Brest, CNRS, IRD, Ifremer, LEMAR, Rue Dumont d'Urville, 29280 Plouzané, France</footnote>

<footnote n="3">³CSIRO, Castray Esplanade Battery Point, Tas 7004, Hobart, Australia</footnote>

<footnote>*Corresponding author. tel: +33 240374157; e-mail: verena.trenkel@ifremer.fr</footnote>
Rawding *et al.* (2014) and Bravington *et al.* (2016a) used versions of CKMR for estimating the abundance of adult chinook salmon and southern bluefin tuna, respectively, based on POPs. Davies *et al.* (2020) extended the southern bluefin tuna work to incorporate HSPs. Hillary *et al.* (2018) estimated abundance of white shark in eastern Australia and New Zealand using HSPs alone. Another application of CKMR POP-based abundance estimation is for brook trout populations (Ruzzante *et al.*, 2019).

Thornback ray (*Raja clavata*, Rajidae) is a medium-sized ray species, maturing at around 5 years (Serra-Pereira *et al.*, 2011). Based on one recapture after more than 16 years at liberty, maximum longevity might be 18–19 years (Bird *et al.*, 2020). For all Rajidae species, females spawn egg capsules that are attached to the substrate or buried (Clark, 1922; Maia *et al.*, 2015). Each adult female, independent of size, spawns 60–150 capsules per year (Holden, 1975; Serra-Pereira *et al.*, 2011), with no documented skipped spawning. Several capsules can be fertilized by the same father using stored sperm (Clark, 1922; Chevolot *et al.*, 2007). Juvenile thornback rays hatch after 4–6 months of incubation (Pawson and Ellis, 2005). This species is found in the Bay of Biscay from coastal and estuarine waters to the upper slope at around 300 m (unpublished survey data). The area encompassing the French shelf and slope of the Bay of Biscay and the Spanish Cantabrian Sea (ICES subdivisions 8a, b, c, and d) is currently managed as a single stock by the International Council for the Exploration of the Sea (ICES, 2020).

The original plan of this study was to sample juvenile and mature thornback ray individuals on the French shelf of the Bay of Biscay and then use the number of POPs and HSPs in the sample to estimate absolute abundance of adults. To this aim, a large number of SNPs were developed previously with HSPs in mind (Le Cam *et al.*, 2019; Marandel *et al.*, 2020). Sampling design is crucial for efficient abundance estimation with CKMR. Based on a rough estimate of mature population size and the assumption of a single population being present on the Bay of Biscay shelf, the initial plan was to sample roughly 3500 juvenile and 3500 adult thornback rays to efficiently achieve precise POP-based estimation, with information from HSPs as a bonus. Unfortunately, due to quota limits, the fishery did no longer land juveniles when the project started. Therefore, sampling ended up being concentrated primarily on the larger, and hence older individuals landed by the fishery, though some young individuals could be sampled from scientific surveys. In this context, age information became crucial for determining the year of birth of potential offspring sampled years later, and also for distinguishing HSP from GGP, which look the same genetically. Unfortunately, estimated ages for larger individuals were rather uncertain. Further, due to the particular biology of Rajidae with females storing the sperm of a given male, many FSP were found, as well as strong evidence for a large number of FTPs (e.g. aunt-nephew) which, like GGP, cannot be distinguished genetically from HSP. For these reasons, CKMR had to be concentrated on a subset of comparisons between individuals that could confidently be assumed to be potential POPs, and the use of HSPs for abundance estimation had to be abandoned (though we did still use HSPs or equivalently related pairs to investigate metapopulation structure). As a side effect of applying CKMR, using the geographic locations of related individuals can provide insight into metapopulation structure (Feutry *et al.*, 2020). Application of CKMR to thornback ray revealed unexpected metapopulation structure, which led to the need

of separate CKMR abundance estimation for local populations, creating further deviation from the planned sampling design. Despite these setbacks, in the end we were still able to produce solid abundance estimates and based on this experience, we summarize the lessons learned for the practical implementation of CKMR.

## Material and methods

### Sample collection

Tissue samples of over 7000 thornback rays were collected from the Bay of Biscay between 2011 and 2020 in landing ports, on board fishing vessels, and on scientific surveys. Sampling covered the main French fishing areas, primarily an offshore area in the central Bay of Biscay and the Gironde estuary (Supplementary Figure S1.1). Total length (lower cm) and sex were recorded for 2002 individuals (Supplementary Figure S2.5). For 5449 individuals, tail samples were provided by fishers and processed in the lab. To assign the sex of these individuals, an assignment method based on genetic markers located on the X chromosome was developed (Trenkel *et al.*, 2020). To estimate missing total lengths, allometric relationships between distance measurements on the tail and total length were derived (Supplementary Figure S2.1).

Growth curves were used to infer probability distributions of the year of egg fertilization (referred to as birth year) for each individual using measured or estimated total length. Age reading on tail vertebrae was attempted but turned out to be unreliable. The details are shown in Supplementary materials S2.1. The year of birth is more uncertain for larger individuals, i.e. close to asymptotic length (here, female $L_∞$ = 115 cm and male $L_∞$ = 105 cm). Therefore, for CKMR abundance estimation, only smaller individuals (< 75 cm) still in the linear growth phase were considered as potential offspring and all larger ones (≥ 75 cm) as potential parents. Sex-specific length-based maturity ogives were transformed to age-based maturities using the same growth curves (Supplementary materials S2.2).

### Genotyping and data filtering

All individuals were genotyped using an Infinium® XT iSelect-96 SNP-array with 9120 SNPs. The design of the SNP array and genotyping effort was performed by Labogena. SNPs and individuals were filtered in several steps to only retain reliably scored SNPs close to the Hardy-Weinberg equilibrium, without linkage disequilibrium, and with call frequency ≥ 98% as well as individuals with ≥ 98% of SNPs genotyped successfully. This resulted in a data set with 3668 SNPs for 6555 individuals (3412 females and 3143 males). The SNP development and filtering methods are described in Supplementary materials (S1.2 and S1.3).

### Kinship inference

POP s and HSPs (or other second-order kin) were found using standard statistical principles described in Thompson (2000). For POPs, the criterion we used was WPSEX, or Weighted PSeudo-EXclusion rate, an exclusion statistic designed to allow for null alleles, which is similar to the empirical exclusion rate when null alleles are negligible, which is the case in our study (see Supplementary materials S2.3.1). At an excluding locus, one individual scores as AA while the other scores as BB, which is an impossible combination for a true POP. Note that there should be no loci with null alleles (e.g. A0) for a true
POP (except for rare mutations), although with thousands of SNPs genotyped, even low genotyping error rates will lead to a small number of excluding loci in some true POPs. The number of excluding loci in Unrelated Pairs (UPs) is substantially larger (and its distribution can be calculated from allele frequencies); closely related non-POP pairs fall somewhere between UPs and POPs, and have non-zero expected WPSEX value even in the absence of error. In practice, the POPs are clearly distinguishable using the WPSEX statistics. Next, FSPs were identified using WPSEX as well as the number of identical genotypes (see Supplementary materials S2.3.2). All calculations were carried out using R and C code.

The cut off value for WPSEX to separate POPs from all other pairs was derived empirically. To evaluate the effect of the number of SNPs on the number of identified POPs, 500–3000 SNPs were randomly selected (10 replicate data sets). For a second group of replicate data sets, only SNPs with minor allele frequency (MAF) > 0.3 were used to evaluate whether such SNPs allowed using fewer SNPs. The same cut off value as for the full data set was then applied to these replicate data sets to identify POPs, which were compared to those obtained using all 3668 SNPs to calculate the number of false negative, i.e. missed POPs, and the number of false positive, i.e. falsely identified POPs due to not using enough SNPs.

For identifying HSPs and other second-order kin pairs for studying metapopulation structure, we calculated the log-likelihood ratio or “PLOD” described in Bravington et al. (2016b) and in the supplementary information of Hillary et al. (2018). For any single locus, the probability of the observed genotypes of a pair can be computed both under a “null hypothesis” that the pair is unrelated, and under the “alternative hypothesis” that the pair is a second-order kin, e.g. HSP. The log-likelihood for that pair and that locus is just the log-ratio of the two probabilities, and the pseudo log-likelihood (PLOD) for the pair is just the sum of the log-likelihoods across all loci. (The “pseudo” reflects that the loci will not be statistically independent unless the pair is really UP, because of linked inheritance between nearby loci in close relatives.) The expected PLOD can be predicted for any true kin-pair type, based on allele frequencies. It is negative for UPs, positive for HSPs, and largest for POPs and FSPs (for UPs, where the loci will not be linked because there is no shared inheritance, the variance and indeed the entire statistical distribution can also be predicted). The general idea is that the actual PLODs for the different kin-types of interest (principally HSPs and UPs) in the data will form clearly separated bumps centred on their predicted means, which can be used to assign kinship for individual pairs; see Bravington et al. (2016b) section 5 for further discussion of false negatives and false positives. This ability to compare observed and expected distributions of PLOD has proved extremely useful for checking whether genotyping and kin-finding has worked properly.

Geographic sampling positions of related pairs (POPs and other kin) were visually inspected to investigate metapopulation structure. Genetic differentiation between the two identified local populations was tested using the likelihood ratio G-statistic as recommended by Goudet et al. (1996) for unbalanced sampling. Calculations were carried out with the hirerfstat R package (Goudet, 2005).

## Abundance estimation

To estimate the abundance of mature thornback ray, we used a simple population dynamics model:

$$N_{s,t} = N_{s,0} e^{\lambda_{s,t} t}, \quad (1)$$

$$\lambda_{s,t} = \lambda_{s,t-1} + \varepsilon_t, \quad \varepsilon_t \sim N(0, \sigma^2), \quad (2)$$

where $N_{s,t}$ is the abundance of mature individuals of sex s in year t, $\lambda_{s,t}$ is the intrinsic population growth rate for this sex in year t, and $N_{s,0}$ is its abundance in the initial year $t = 0$. Considering mature females ($s = f$) and assuming equal fecundity on average, the probability female i with birth year $b_i$ is the mother of individual j (meaning i and j are a mother-offspring pair; MOP) depends only on the number of mature females $N_{f,b_i}$ in the birth year $b_j$ of individual j, given i was old enough to have been mature at that time. Further, since sampling was primarily lethal (port and onboard sampling), the sampling date $t_i$ of individual i needs to have been after the spawning season $s_j$ (egg-laying for thornback ray) in the birth year $b_j$ of individual j.

The general idea is that each comparison, say between animals i and j, is a Bernoulli yes/no random variable with very low expected value (on the order of the reciprocal of adult abundance). The Bernoulli distributions were approximated by Poisson distributions for computational efficiency. In other words, we need a length-based probability, given by

$$\begin{aligned} P_{\text{MOP}}^{[\text{len}]} (l_i, l_j, t_i, t_j) &\triangleq P(K_{ij} = \text{MOP} | l_i, l_j, t_i, t_j) \\ &= \sum_{a_i, a_j} P_{\text{MOP}}^{[\text{age}]} (a_i, t_i, b_j = t_j - a_j) \\ &\qquad \times P(a_i | l_i) \times P(a_j | l_j), \end{aligned} \quad (3)$$

where $l_i$ is the length of individual i and $a_i$ its age at capture date $t_i$ and the age-based probability is

$$P_{\text{MOP}}^{[\text{age}]} (a_i, t_i, b_j) = \begin{cases} \frac{P(\text{mat} = 1 | a_i, t_i, b_j)}{0.5 N_{f,b_j}} & t_i > b_j \\ 0 & \text{otherwise} \end{cases} \quad (4)$$

with $P(\text{mat} = 1 | a_i, t_i, b_j)$ the probability of i being mature at the birth of j.

We aggregate this probability across individuals with the same covariates (i.e. length-at-capture and capture-year) because the comparisons are effectively independent (see Bravington et al. (2016b) on “sparse sampling”), and everything stays Poisson. Let $m_{l,t}$ be the number of samples with covariates $l$ and $t$, so that the number of pairwise comparisons with $(l_i, l_j, t_i, t_j)$ is $m_{l,t_i} \times m_{l,t_j}$ for $l_i < l_j$. If $C_{\text{MOP}}(l_i, l_j, t_i, t_j)$ is the actual number of MOPs found amongst those comparisons, then let $\lambda_{l_i, l_j, t_i, t_j} \triangleq E[C_{\text{MOP}}(l_i, l_j, t_i, t_j)]$ so that

$$\mu_{l_i, l_j, t_i, t_j} = m_{l_i, t_i} \times m_{l_j, t_j} \times P_{\text{MOP}}^{[\text{len}]} (l_i, l_j, t_i, t_j). \quad (5)$$

The log-likelihood across all pairwise comparisons is just that of a set of independent Poisson distributions:

$$\begin{aligned} \Lambda_f &= \sum_{l_i, t_i, l_j, t_j} \log \left( e^{-\mu_{l_i, l_j, t_i, t_j}} \times \mu_{l_i, l_j, t_i, t_j}^{C_{\text{MOP}}} (l_i, l_j, t_i, t_j)! \right) / C_{\text{MOP}} (l_i, l_j, t_i, t_j)! \\ &= \sum_{l_i, t_i, l_j, t_j} \left\{ -\mu_{l_i, l_j, t_i, t_j} + C_{\text{MOP}} (l_i, l_j, t_i, t_j) \times \log (\mu_{l_i, l_j, t_i, t_j} + c) + c \right\}, \end{aligned} \quad (6)$$

where $C_{\text{MOP}}(l_i, l_j, t_i, t_j)!$ is constant for any given data set, hence set to c in the second line. The same assumptions and equations were used for mature males ($s = m$).
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Parent-offspring kinfinding statistics WPSEX for thornback ray. All pairs with WPSEX < 0.004 were considered POPs (left of vertical dashed line). FSP, 'other' related pairs includes HSPs and FTPs.</figcaption></figure>

Available studies on thornback ray and other Rajidae mostly reflect a 1:1 sex-ratio. Where deviation from this balanced ratio were observed, these were reported to be potentially biased by sex differences in spatial distribution, aggregative behaviour, or other factors, so more reflecting biased sampling than actual population parameters (e.g. Ellis and Shackley, 1995; Ebert, 2005; Frisk and Miller, 2006). Therefore, we assumed $N_{f,0} = N_{m,0} = N_0$ and common growth parameters $\lambda_{1...T}$.

Model parameters $N_0$, $\lambda_{1...T}$, and $\sigma^2$ were estimated maximizing $\Lambda$ as in Bravington *et al.* (2016b) and adding the log-prior on $\lambda$; $\sigma^2$ was estimated as usual in such random-effect frameworks by maximizing an automatic Laplace approximation. All calculations and analyses were carried out in R using the TMB R package for maximization (Kristensen *et al.*, 2016). TMB code is provided in Supplementary material S2.4 and online (Trenkel *et al.*, 2021).

# Results

## Kinship inference

Using the exclusion rate (WPSEX), 99 POPs were identified among the 6555 individuals; two offspring had both parents in the sample. The histogram of exclusion rates showed a clear gap between POPs and FSPs (Figure 1). In 58 of 99 POPs, the parent was female (Supplementary Table S3.1). Based on the exclusion rate value plus the number of identical genotypes, 431 FSPs were also identified. Restricting the comparisons to smaller potential offspring (< 75 cm) and larger potential parents (≥ 75 cm) reduced the number of POPs to 74 (57 offshore, 16 Gironde, and 1 Bay of Douarnenez, see below).

Subsampling SNPs showed that ~2000 SNPs were sufficient for robust identification of the number of POPs (Figure 2). As expected, using only SNPs with MAF > 0.3 decreased the number of missed (false negative) POPs more rapidly with increasing number of SNPs compared to using all SNPs (MAF > 0.1; Figure 2a), but at the same time slowed down the reduction of false positive POPs as the number of SNPs increased (Figure 2b).

The histogram of observed PLOD scores showed a clear bump with its mode at the value expected for HSPs and other second-order kin (PLOD = 112, Figure 3). The left-hand side of the bump (e.g. around PLOD = 50) clearly overlapped with weaker kin (third-order and higher), which is consistent with the experience of second-order kin in other CKMR applications (M. Bravington, pers. comm.). Nevertheless, it is clear that above some threshold of, say PLOD = 60, the great majority of the pairs is truly second-order. For quantitative use of HSPs in the CKMR model, it would be necessary to choose a PLOD threshold and carefully estimate its associated false-negative rate, but here we are just using second-order kin for qualitative insights on connectivity, so exact numbers are not crucial. Using a reasonably conservative cutoff PLOD = 60, we retained 3400 pairs of first- and second-order kin (from 3323 unique individuals, including 389 individuals from the Gironde and 2921 individuals from offshore) for studying metapopulation structure (family example in Supplementary Figure S3.1).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Number of false negative (a) and false positive (b) POPs as a function of the number of SNPs and the minimum MAF of selected SNPs for thornback ray in the Bay of Biscay. Results for 10 replicate data sets. For MAF > 0.3, only 1000 SNPs were available.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** PLOD scores for HSP, FSP, POP, and other kinship pairs for thornback ray in the Bay of Biscay. The vertical grey bar indicates the expected PLOD value for HSP, FTP, and GGP. The arrow on the x-axis indicates PLOD=60.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Sampling position of related individuals. (a) POPs. (b) Other pairs of related individuals including FSPs, HSPs, and FTPs (PLOD > 60, excluding POPs).</figcaption></figure>

Comparison of geographic sampling positions of related pairs revealed a signature of small-scale metapopulation structure (Figure 4). None of the 3400 close-kin pairs included both a Gironde and an offshore sample, compared to 245 Gironde-only and 3146 offshore-only pairs (further details in Supplementary Table S3.3). Since individuals were sampled primarily from landings, this indicates that two different local populations are exploited by fishers. From a population genetic point of view, the two local populations were significantly differentiated (G-statistic, *P*-value < 0.001).

## Abundance estimation

Although we identified 99 POPs overall, one quarter of those involved larger adults (> 75 cm) as the offspring. Unfortunately, uncertainty about the year of birth for large individuals was so large that they could not usefully be treated as potential offspring within a CKMR model. Their birth year lying far in the past, they would also have large leverage on the estimated abundance time trend. Accordingly, we restricted POP comparisons to pairs where the potential offspring was

< 75 cm at sampling, i.e. at sizes from which age can be estimated with adequate reliability. Further, given the identified metapopulation structure, abundance estimation had to be carried out by population. Thus, only POPs and samples from the Gironde estuary and the offshore area (Figure 1) were retained for abundance estimation. This left us with 0.15 × 10⁶ comparisons in the Gironde and 1.6 × 10⁶ comparisons offshore, among which we found 73 POPs (Supplementary Table S3.2). Among the 73 POPs, in the Gironde estuary, 12 POPs involved a mother and four a father, while for offshore, 31 POPs had a mother and 26 a father. Binomial tests confirmed that the proportion of POPs involving a mother was compatible with equal numbers of mature females and males in both local populations (Gironde *p* = 0.21; offshore *p* = 0.69; and combined *p* = 0.81).

Because of the unequivocal evidence for demographically separate populations, we fitted separate POP-based models to the Gironde samples on one hand, and the offshore samples on
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Estimated abundance of mature thornback ray in the Bay of Biscay. (a) Estimates by subarea (see map in Figure 4). (b) Total abundance treating subareas separately ("sum") or jointly in the estimation. Shaded areas are 95% confidence bands.</figcaption></figure>

the other hand (model parameter estimates in Supplementary Table S3.4).

For Gironde, the estimated adult abundance was around 19000 individuals (CV 41%) in 2013 (Figure 5a). In the offshore population, there were around 116000 individuals (CV 20%) in 2013. The estimated population growth rates were imprecise, because of the rather low number of POPs and the limited span of offspring cohorts, and the confidence intervals certainly include zero. There was no evidence of annual variations in growth rate (i.e. estimated $\sigma^2$ was zero in both populations), although statistical power to detect such variations would again be very low. The combined estimate of adults in the central Bay of Biscay was around 135000 (CV 0.19) in 2013 (Figure 5b). If metapopulation of the structure was ignored and total abundance was to be estimated treating all samples as coming from the same population, a similar overall abundance would be concluded for this case, though with a slightly different uncertain time trend (Figure 5b).

To study the potential effect of the unequal sex ratio in the sample on estimates, we compared abundance estimates obtained for resampled data with equal sex ratio of potential parents with results for a 1.33:1 sex ratio (females:males). For this, the combined data set was used to increase sample size. The resulting abundance estimates differed more between random data sets than between the two sex ratios (Supplementary Figure S3.3).

# Discussion

## Thornback ray in the Bay of Biscay

Kinship inference revealed 3400 clearly related individuals among the sampled thornback rays from the Bay of Biscay, forming 99 POPs, 431 FSPs, and at least 2870 second-order kin pairs (the true number of second-order pairs in the

dataset will be somewhat higher, because we used a conservative criterion for selection). Given that some of the second-order kin are actually GGP or FTPs rather than HSPs, the ratio of FSPs to HSPs might seem surprisingly high. In fact, though, it is entirely consistent with thornback ray biology. Within a given year, a female thornback ray is known to store sperm that she uses to fertilize several eggs when spawning on average every other day (Holden, 1975), so, even if she mates with more than one male that year, many within-year sibs will be FSPs rather than HSPs. Cross-year sibs will of course almost always be HSPs rather than FSPs so, provided that female adults live long enough to breed in several years, there will be more cross-year sibs overall than FSPs. However, the higher the adult mortality rate is, the fewer breeding occasions each female will have, and the lower the ratio of HSP:FSP will be in the population as a whole. Our existing data are too vague about individual ages to allow this to be made precise, but a ratio of somewhere around 7:1 HSP:FSP seems consistent with maximum adult lifespans of 18–19 years (Bird et al., 2020).

Our abundance results are based on the number (73) of usable POPs, though the large number of second-order kin are also potentially very useful for CKMR, providing information on adult mortality rates as well as on abundance (Hillary et al., 2018). However, the constitution of our samples (mostly large individuals), and the severe imprecision of length-based age estimates for those adults, have deterred us from trying to incorporate sib-pairs into our current model. Most of our pairwise comparisons are between two adults, and the birth-gap between those adults may be substantial, so it is to be expected on demographic grounds that the second-order kin will contain not just HSPs, but also a substantial number of GGPs and FTPs (since FSPs are reasonably common among adults). If the birth-gap between the two individuals in a pair is known with reasonable precision, then it is possible to restrict comparisons a priori to pairs separated by short gaps where GGPs and FTPs are unlikely or impossible (as in Farley et al., 2018;
Hillary *et al.*, 2018), but we cannot do that here. Despite the problems in using our second-order kin-pairs for quantitative abundance and mortality estimation, the sampling locations of a very large number of close-kin pairs (3000+) provided strong evidence for hitherto unknown metapopulation structure of thornback ray in the Bay of Biscay. In future work, it would be interesting to study more closely the fine-scale geographic distribution of first- and second-order kin within our sample, and to extend the study to other locations in the Bay of Biscay, e.g. the Cantabrian Sea along the Spanish north coast or the southern Brittany coast. We anticipate that further local populations will be found.

CKMR relies on identifying related individuals using genetic markers, with the number of SNPs not being too small. Bravington *et al.* (2016b) recommended to use several thousand SNPs. For thornback ray, around 2000 SNPs were found to be sufficient to reliably identify POPs. This is more SNPs than would be needed for finding POPs if all the non-POP pairs were unrelated or only weakly related (e.g. second-order kin), but the present dataset contains also a substantial number of FSPs. FSPs are statistically harder to separate from POPs, and our statistical criteria (WPSEX, which is really aimed at distinguishing POPs from weak kin, that being the situation we had expected *a priori*) may not be a fully efficient statistic for discriminating FSPs from POPs, so the number of SNPs required was fairly large. Thus, given data cleaning and filtering remove a substantial number of SNPs, we recommend to genotype individuals for more SNPs than thought to be needed. Of course, the method used for deriving and validating SNPs impacts the number of SNPs filtered out. For thornback ray, a Restriction Associated DNA Sequencing (RAD-Seq) protocol was used for identifying SNPs, followed by a SNP array genotyping method, with no intermediate validation step, which might explain why many SNPs were subsequently filtered out. From the set of 9120 SNPs on the array (including some duplicated loci for quality control), we ended up using 3668 for genotyping, which was adequate for identifying POPs, and would in fact have been adequate for second-order kin (using the thresholding/false-negative estimation approach in Bravington *et al.* (2016b)) if we had been using the latter directly in our models.

The sex ratio in the sample was biased towards females, presumably because the samples came primarily from commercial fisheries, which target larger individuals, and females grow larger. If, however, the sex ratio in the population was unequal and in addition the per capita chance of being sampled differed between sexes, biased abundance estimates would have been obtained unless the sexes were treated separately. However, published information suggests equal sex ratio for thornback ray at birth (Ellis and Shackley, 1995) and given the recent relatively low exploitation rate (Marandel *et al.*, 2019), the sex ratio should remain approximately equal for mature individuals even though some fishers might target larger individuals. Thus, we would not expect the bias to be large compared to the CV. More samples from male adults would resolve this in future, by allowing sex-specific models.

Before carrying out this study, insights into thornback ray population structure came from traditional tagging studies carried out around the British Isles (see review in Bird *et al.*, 2020) and a population genetic study covering a similar area (Chevolot *et al.*, 2006), while for the Bay of Biscay only a simulation study was available (Marandel *et al.*, 2018). Comparing release and recapture positions of tagged individuals sug-

gested 62.1% of movements were < 50 km after more than 50 d, with only 14.8% of movements > 100 km (Bird *et al.*, 2020). Release locations varied widely, so did the number of days at sea, but importantly, movements do not inform on local breeding and hence demographic connectivity in contrast to studying kinship co-locations. Using five microsatellite loci, Chevolot *et al.* (2006) detected weak but significant population differentiation in waters around southern England, but again were not able to infer contemporary connectivity patterns. The expected degree of demographic connectivity depends not only on the movement scenario, but also on the relative abundance of neighbouring populations (Marandel *et al.*, 2018), and of course on the assumption that individuals interbreed after arrival. A large number of comparisons between individuals sampled in the Gironde estuary and those sampled further offshore amounting to around eight times more comparisons than within the Gironde estuary did not return a single related pair. Thus, demographic connectivity between the two local populations, if it exists, must be small. It is, however, likely that in the past thornback ray abundance was higher in the Bay of Biscay (Marandel *et al.*, 2016, 2019), and therefore, individuals were more spread out as expected in the presence of an abundance-occupancy relationship; such a relationship has been observed for other Rajidae (Frisk *et al.*, 2011).

Subsequent application of CKMR to the two putative local populations revealed evidence for an increasing recent time trend of the number of mature thornback ray in the central Bay of Biscay, and most importantly provided the first absolute abundance estimates. Precision (especially of trend estimates) is currently limited, but could easily be improved by adding more samples. Previous attempts to estimate population abundance in the Bay of Biscay have been hampered by data uncertainty and strong effects of modelling assumptions on estimated abundance levels and time trends (Marandel *et al.*, 2016, 2019).

Application of CKMR comes with a number of challenges, from design, through sampling, through genetic analysis, through finding kin pairs, to modelling. Based on the experience gained by applying CKMR to thornback ray, we now summarize and discuss four main lessons learnt that should be of interest to scientists intending to apply the CKMR method for abundance estimation.

## Lessons learned

### Lesson 1: CKMR helps identifying metapopulation structure

The simple non-spatial model used here for abundance estimation (Equation (1)) assumes implicitly that all individuals (parents and offspring) are independently sampled from the same population (Bravington *et al.*, 2016b). Conn *et al.* (2020) using simulations studied the effect when this assumption is violated because of limited dispersal. Simulating incomplete mixing of individuals combined with spatially biased sampling, they found that simple CKMR abundance estimates can be substantially negatively biased; however, in such cases it was usually possible to detect statistically that the complete-mixing assumption was violated, provided that the samples had adequate geographical spread. For thornback ray in the Bay of Biscay, geographic position of related individuals indicated that mixing between two sampled local populations is currently probably rather limited (Figure 4). To account for this, separate population dynamics models were fitted for the two local populations. However, total abundance
estimates obtained ignoring metapopulation structure did not differ much from the sum of the two local estimates (Figure 5). Unbiased joint total estimates are expected if per capita sampling chance was similar for the two populations. By chance, this might have been the case for thornback ray.

## Lesson 2: demographic sample composition affects precision

For POP-based CKMR abundance estimation, both potential offspring and potential parents need to be sampled. For thornback ray, juveniles turned out to be difficult to sample as they are not much landed by fisheries and several dedicated surveys caught few individuals. The limited knowledge on metapopulation structure at the time of sampling meant we were targeting juveniles in coastal bays, while it appears now that there must also be juveniles further offshore. The compromise found was to consider as potential offspring all individuals < 75 cm for which length can most reliably be transformed into age estimates, and thus year of birth estimates (next section). However, this filtering reduced the number of POPs and potential parents used for abundance estimation and hence reduced estimation precision. A sample size effect is clearly visible as smaller confidence intervals for total abundance estimates obtained using the data for both areas jointly compared to the wider confidence intervals for the sum of local abundance estimates (Figure 5b). Further, for estimating the population time trend with POPs, offspring need of course to have been born in different years. Spread of birth years of sampled individuals is also essential when using the number of second-order kin for abundance estimation. Lastly, care needs to be taken that the sampled adults could actually have been mature at the time of birth of the sampled offspring, otherwise estimates will be biased (Waples and Feutry, 2021).

## Lesson 3: age information is essential

The age of individuals at sampling informs on their year of birth (egg-laying year for thornback ray), which is the reference year for the estimated number of mature individuals using CKMR. The older the offspring, the further in the past the reference year and the more consequential model misspecification, such as assumptions about growth rate variations, and age uncertainty are for CKMR abundance estimation. Further, for potential parents, age at sampling informs on their year of first maturity, which combined with information on maturity-at-age is used to calculate the expected number of POPs in the sample. Thus again, the further the reference year in the past, the more age uncertainty at sampling comes to play when determining whether an individual might have been mature or not. Overall, this means it is recommendable to restrict potential offspring to younger individuals at the time of sampling unless age is well-known for all sampled individuals.

For many species, age is, however, not easily accessible and length is used as a proxy. For thornback ray, we tried ageing individuals based on reading year rings in tail vertebrae but with no success. In such a case, length can be converted to age by applying a length-age curve, though with uncertainty, and potentially bias. Importantly, uncertainty becomes very large as individuals approach their asymptotic length. Thus, depending on the growth pattern of the species, length might not be a sufficiently reliable estimator of age except for the youngest age-classes. For applying CKMR to thornback ray, we only retained individuals still in the linear growth phase

as potential offspring to minimize age, and thus year of birth uncertainty.

Thus, the uncertainty associated with age (particularly of large individuals) prevented us from making full use of the very respectable total numbers of kin-pairs that we found, with 99 actual POPs reduced to 73 actual ones in the final model. Lastly, for studies with reliable age data, the use of HSPs provides not only abundance estimates but also total mortality estimates, given the HSPs are spread over several birth years.

## Lesson 4: sex markers are useful

The sex of sampled parents is needed to calculate the potential number of fathers and mothers in the sample. For a large proportion of thornback ray individuals, the sex was not recorded at sampling. The same issue might be faced when studying vulnerable species with no external sex attributes that need to be sampled non-lethally, or historic tissue samples for which not all individual related information was recorded. Fortunately, the RADSeq derived set of SNPs included SNPs on the sex chromosome that could be used for sex assignment (Trenkel et al., 2020). Further, many species have different reproductive dynamics for males and females, especially different fecundity-at-age schedules, so for those species, CKMR models should be constructed separately by sex, though perhaps sharing a few parameters. For this, however, sufficient samples and reliable age estimates are needed.

## Methodological comment

A fundamental issue of CKMR abundance estimation is that the denominators of kinship-probabilities involve “TRO” (Total Reproductive Output) rather than abundance per se as assumed in Equation (4) (section 3 in Bravington et al., 2016b). Thus, if fecundity varies among adults, POPs alone cannot distinguish between a large number of low-fecundity (small and young) adults and a smaller number of high-fecundity (large and old) ones that would produce the same TRO. And for HSPs alone, there is also a quadratic effect of systematic (within-lifespan, and/or between adults) variations in fecundity on the expected number of kin-pairs (section 3.10 in Bravington et al., 2016b), again meaning that numerical adult abundance *per se* cannot be separated from the extent of fecundity variation without further data sources. Both of these phenomena can, in fact, be dealt with inside a CKMR model (i.e. an unbiased and unambiguous abundance estimate is possible) provided that both POP and HSP data are used together, though the full details are too complicated to give here (see Bravington, 2017). However, if on biological grounds there is little reason to expect any such systematic variations in fecundity, as for thornback ray, then the “total reproductive output” can be expressed directly as numerical adult abundance, and there are no quadratic complications to worry about, so that either POPs alone or HSPs alone provide unambiguous abundance estimates. Nevertheless, we would always recommend using both POPs and HSPs if possible.

## Conclusion

The application of CKMR to thornback ray was successful, though a number of hurdles had to be overcome on the way, which at the same time provided valuable insights into metapopulation structure. To share this experience we have distilled the experience into four lessons, which we believe
should be of interest to scientists planning to apply CKMR abundance estimation. For many species it will be necessary to use both siblings and POPs, but for species where reproductive output does not change through adulthood (including thornback rays), POPs alone (or juvenile HSPs alone) may be sufficient.

## Funding

This study received funding from the French National Research Agency (project GenoPopTaille, contract ANR-14-CE02-0006-01) and the European Union's Horizon 2020 research and innovation programme under grant agreement number 773713 (PANDORA).

## Supplementary data

Supplementary material is available at the *ICESJMS* online version of the manuscript.

## Data availability

The RADSeq-derived sequences used for SNP development are available at https://doi.org/10.17882/70648 and the DNA sequences for the 9120 SNP array from https://doi.org/10.17882/70546. The genotypes of individuals used for identifying sex-related SNPs are provided at https://doi.org/10.17882/74390. The data and code used for CKMR abundance estimation are available at https://doi.org/10.17882/84524.

## Authors' contributions

MB, GC, PL, and VT designed the study, VT and PL analysed the data, and VT led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval.

## Acknowledgements

We thank E. Blanc, E. Stephan, APECS (Association Pour l'Étude et Conservation des Sélaciens) volunteers, and scientists for collecting tissue samples on the surveys carried out by INRAe in the Gironde estuary (Sturat survey), Ifremer on the Bay of Biscay shelf (Evhoe survey), and in the Bay of Douarnenez (RaiesJuves survey). We are grateful to the skippers of Fisher Golf and Le Battant, the crews of research vessels Albert Lucas and Thalia, as well as the fish auction markets in La Cotinière and Royan. We would also like to thank F. Marandel, A. Bidault, E. Borcier, and W. Handal for help with tissue preparation, S. Le Cam for contributing to the SNP identification process, V. Verrez-Bagnis for carrying out the BLAST analysis, and R. Hillary for providing R code to estimate age probabilities from length. The position data for commercial landings was provided by the Ifremer Système d'Information Halieutique. We thank R. Waples, the editor and an anonymous referee for insightful comments and suggestions.

## References

* Anderson, E. C., and Garza, J. C. 2006. The power of single-nucleotide polymorphisms for large-scale parentage inference. Genetics, 172: 2567-2582.

* Augustin, N. H., Trenkel, V. M., Wood, S. N., and Lorance, P. 2013. Space-time modelling of blue ling for fisheries stock management. Environmetrics, 24: 109-119.

* Bird, C., Burt, G. J., Hampton, N., Phillips, S. R. M., and Ellis, J. R. 2020. Fifty years of tagging skates (Rajidae): using mark-recapture data to evaluate stock units. Journal of the Marine Biological Association of the United Kingdom, 100: 121-131.

* Bravington, M. 2017. SBT kin-finding and genotyping update. ICES Document Technical report CCSBT-OMMP 8-12. Commission for the Conservation of Southern Bluefin Tuna.

* Bravington, M. V., Grewe, P. M., and Davies, C. R. 2016a. Absolute abundance of southern bluefin tuna estimated by close-kin mark-recapture. Nature Communications, 7: 13162.

* Bravington, M. V., Skaug, H. J., and Anderson, E. C. 2016b. Close-kin mark-recapture. Statistical Science, 31: 259-274.

* Chevolot, M., Ellis, J. R., Hoarau, G., Rijnsdorp, A. D., Stain, W. T., and Olsen, J. L. 2006. Population structure of the thornback ray (Raja clavata L.) in British waters. Journal of Sea Research, 56: 305-316.

* Chevolot, M., Ellis, J. R., Rijnsdorp, A. D., Stam, W. T., and Olsen, J. L. 2007. Multiple paternity analysis in the thornback ray *Raja clavata* L. Journal of Heredity, 98: 712-715.

* Clark, R. S. 1922. Rays and skates (Raïza) no. 1. Egg-capsules and young. Journal of the Marine Biological Association of the United Kingdom, 12: 578-643.

* Conn, P. B., Bravington, M. V., Baylis, S., and Ver Hoef, J. M. 2020. Robustness of close-kin mark-recapture estimators to dispersal limitation and spatially varying sampling probabilities. Ecology and Evolution, 10: 5558-5569.

* Davies, C. R., Bravington, M. V., Eveson, J. P., Lansdell, M., Aulich, J., and Grewe, P. M. 2020. Next-generation Close-kin Mark Recapture: using SNPs to identify half-sibling pairs in Southern Bluefin Tuna and estimate abundance, mortality and selectivity. Fisheries Research and Development Corporation. 52pp.

* Ebert, D. A. 2005. Reproductive biology of skates, Bathyraja (Ishiyama), along the eastern Bering Sea continental slope. Journal of Fish Biology, 66: 618-649.

* Ellis, J. R., and Shackley, S. E. 1995. Observations on egg-laying in thornback ray. Journal of Fish Biology, 46: 903-904.

* Farley, J., Eveson, P., Bravington, M., Aulich, J., and Grewe, P. 2018. Update on the SBT close-kin tissue sampling, processing, and kin finding. ICES Document CCSBT-ESC/1809/08. Australia's National Science Agency. 13pp.

* Feutry, P., Devloo-Delva, F., Tran Lu, Y. A., S., Gunasekera, R. M., Johnson, G., Pillans, R. D. et al. 2020. One panel to rule them all: DArT-cap genotyping for population structure, historical demography, and kinship analyses, and its application to a threatened shark. Molecular Ecology Resources, 20: 1470-1485.

* Frisk, M. G., Duplisea, D. E., and Trenkel, V. M. 2011. Exploring the abundance-occupancy relationships for the Georges Bank finfish and shellfish community from 1963-2006. Ecological Applications, 21: 227-240.

* Frisk, M. G., and Miller, T. J. 2006. Age, growth, and latitudinal patterns of two Rajidae species in the northwestern Atlantic: little skate (*Leucoraja erinacea*) and winter skate (*Leucoraja ocellata*). Canadian Journal of Fisheries and Aquatic Sciences, 63: 1078-1091.

* Goudet, J. 2005. HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Molecular Ecology Notes, 5: 184-186.

* Goudet, J., Raymond, M., deMeeus, T., and Rousset, F. 1996. Testing differentiation in diploid populations. Genetics, 144: 1933-1940.

* Hilborn, R., and Walters, C. J. 1992. Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty. Chapman and Hall, New York, NY. 570pp.

* Hillary, R. M., Bravington, M. V., Patterson, T. A., Grewe, P., Bradford, R., Feutry, P., Gunasekera, R. et al. 2018. Genetic relatedness reveals total population size of white sharks in eastern Australia and New Zealand. Scientific Reports, 8: 2661.
Holden, M. J. 1975. The fecundity of *Raja clavata* in British waters. ICES Journal of Marine Science, 36: 110-118.

ICES. 2020. Thornback ray (*Raja clavata*) in Subarea 8 (Bay of Biscay) (advice for 2019 and 2020). Report of the ICES Advisory Committee, 2020-rjc.27.8. CES Advisory Committee.

Kotwicki, S., Horne, J. K., Punt, A. E., and Ianelli, J. N. 2015. Factors affecting the availability of walleye pollock to acoustic and bottom trawl survey gear. ICES Journal of Marine Science, 72: 1425-1439.

Krieger, K. J., and Sigler, M. F. 1995. Catchability coefficient for rockfish estimated from trawl and submersible surveys. Fishery Bulletin, 94: 282-288.

Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., and Bell, B. M. 2016. TMB: Automatic Differentiation and Laplace Approximation. Journal of Statistical Software, 70: 1-21.

Le Cam, S., Bidault, A., Charrier, G., Cornette, F., Lamy, J.-B., Lapegue, S., Lorance, P. et al. SNPs for Thornback Ray *Raja clavata*, Validated from Genotyped Individuals. 2019. SEANOE. doi: 10.17882/70546.

Maia, C., Serra-Pereira, B., Erzini, K., and Figueiredo, I. 2015. How is the morphology of the ovi-ducal gland and of the resulting egg capsule associated with the egg laying habitats of Rajidae species? Environmental Biology of Fishes, 98: 2037-2048.

Marandel, F., Charrier, G., Lamy, J.-B., Le Cam, S., Lorance, P., and Trenkel, V. M. 2020. Estimating effective population size using RADseq: effects of SNP selection and sample size. Ecology and Evolution, 10: 1929-1937.

Marandel, F., Lorance, P., Andrello, M., Charrier, G., Le Cam, S., Lehuta, S., and Trenkel, V. M. 2018. Insights from genetic and demographic connectivity for the management of rays and skates. Canadian Journal of Fisheries and Aquatic Science, 75: 1291-1302.

Marandel, F., Lorance, P., and Trenkel, V. M. 2016. A Bayesian state-space model to estimate population biomass with catch and limited survey data: application to the thornback ray (*Raja clavata*) in the Bay of Biscay. Aquatic Living Resources, 29: 209.

Marandel, F., Lorance, P., and Trenkel, V. M. 2019. Determining long-term changes in a skate assemblage with aggregated landings and limited species data. Fisheries Management and Ecology, 26: 365-373.

Pawson, M. G., and Ellis, J. R. 2005. Stock identity of elasmobranchs in the Northeast Atlantic in relation to assessment and management. Journal of Northwest Atlantic Fishery Science, 37: 173-193.

Rawding, D. J., Sharpe, C. S., and Blankenship, S. M. 2014. Genetic-based estimates of adult chinook salmon spawner abundance from carcass surveys and juvenile out-migrant traps. Transactions of the American Fisheries Society, 143: 55-67.

Ruzzante, D. E., McCracken, G. R., Forland, B., MacMillan, J., Notte, D., Buhariwalla, C., Flemming, J. M. et al. 2019. Validation of close-kin mark-recapture (CKMR) methods for estimating population abundance. Methods in Ecology and Evolution, 10: 1445-1453.

Serra-Pereira, B., Figueiredo, I., and Gordo, L. S. 2011. Maturation, fecundity, and spawning strategy of the thornback ray, *Raja clavata*: do reproductive characteristics vary regionally? Marine Biology, 158: 2187-2197.

Thompson, E. A. 2000. Statistical inference from genetic data on pedigrees. NSF-CBMS Regional Conference Series in Probability and Statistics. 6: 1-169.

Trenkel, V. M., Boudry, P., Verrez-Bagnis, V., and Lorance, P. 2020. Methods for identifying and interpreting sex-linked SNP markers and carrying out sex assignment: application to thornback ray (*Raja clavata*). Molecular Ecology Resources, 20: 1610-1619.

Trenkel, V. M., Buckland, S. T., McLean, C., and Elston, D. A. 1997. Evaluation of aerial line transect methodology for estimating red deer (*Cervus elaphus*) abundance in Scotland. Journal of Environmental Management, 50: 39-50.

Trenkel, V. M., Elston, D. A., and Buckland, S. T. 2000. Fitting population dynamics models to count and cull data using sequential importance sampling. Journal of the American Statistical Association, 95: 363-374.

Trenkel, V. M., and Skaug, H. J. 2005. Disentangling the effects of capture efficiency and population abundance on catch data using random effects models. ICES Journal of Marine Science, 62: 1543-1555.

Trenkel, V., Charrier, G., and Lorance, P. 2021. Data for CKMR Abundance Estimation for Thornback Ray (*Raja clavata*) in the Central Bay of Biscay. SEANOE. doi: 10.17882/84524.

Waples, R. S., and Feutry, P. 2021. Close-kin methods to estimate census size and effective population size. Fish and Fisheries. Early access.

Handling Editor: Lorenz Hauser