Grandparent inference from genetic data: The potential for parentage-based tagging programs to identify offspring of hatchery strays

Fisheries managers routinely use hatcheries to increase angling opportunity. Many hatcheries operate as segregated programs where hatchery-origin fish are not intended to spawn with natural-origin conspecifics in order to prevent potential negative effects on the natural-origin population. Currently available techniques to monitor the frequency with which hatchery-origin strays successfully spawn in the wild rely on either genetic differentiation between the hatchery- and natural-origin fish or extensive sampling of fish on the spawning grounds. We present a method to infer grandparent-grandchild trios using only genotypes from two putative grandparents and one putative grandchild. We developed estimators of false positive and false negative error rates and showed that genetic panels containing 500 - 700 single nucleotide polymorphisms or 200 - 300 microhaplotypes are expected to allow application of this technique for monitoring segregated hatchery programs. We discuss the ease with which this technique can be implemented by pre-existing parentage-based tagging programs and provide an R package that applies the method.


Introduction 31
Fisheries managers have long used hatcheries to increase angling opportunity and to 32 compensate for anthropogenic impacts that have decreased population sizes (Waples et al. 2007). 33 In some situations, it has been observed that hatchery-origin fish have lower fitness in the wild 34 relative to natural-origin conspecifics, potentially due to selection in the hatchery environment or 35 the use of hatchery strains that are not locally adapted (Ford 2002;Miller et al. 2004;Araki et al. 36 2007aAraki et al. 36 , 2007bAraki et al. 36 , 2008Christie et al. 2014). In order to prevent negative effects of hatchery-and 37 natural-origin fish interbreeding, it has been recommended that hatcheries operate as either 38 In order to evaluate the efficacy of segregated hatchery programs, a method of 47 monitoring gene flow from hatchery-origin fish to nearby natural-origin populations is needed. 48 This has been previously estimated in several ways. Some studies estimate the proportion of fish 49 on the spawning grounds that are of hatchery origin (pHOS) through observation of marks and/or 50 tags present on hatchery fish (Dauer et al. 2009;Tattam and Ruzycki 2020). This metric is 51 important as it demonstrates the potential for competition between hatchery-and natural-origin 52 fish on the spawning grounds and suggests that hatchery introgression may be occurring. 53 7 were excluded from consideration. The value of m was chosen so that the probability of 123 excluding a true grandparent-grandchild trio was less than 0.0001 for a trio with no missing 124 genotypes. With a constant value of m, the probability of rejecting a trio with one or more 125 missing genotypes is smaller, as a locus with a missing genotype cannot be considered a 126 Mendelian incompatibility. The probability of rejecting a true trio given a value of m was 127 calculated by representing the number of MIs as a Markov Chain, following the method detailed 128 by Anderson (2012) but extended to the relationship of grandparent-grandchild trio. 129 Assigning grandparent-grandchild relationships can be treated as a hypothesis test by 130 comparing the calculated log-likelihood ratio (LLR) to a critical value, c (SanCristobal and 131 Chevalet 1997;Anderson and Garza 2006). If the LLR was greater than or equal to c, the trio 132 was considered related; otherwise, the trio was considered unrelated. The value of c can be 133 chosen to achieve a desired balance of false negative and false positive error rates. 134 135 False positive error rates 136 Per-comparison false positive error rates (probability of assigning a relationship to a trio 137 that is not a grandparent-grandchild trio) are dependent upon the true relationship for a trio. We 138 have implemented methods to assess false positive error rates for unrelated trios (typically the 139 most important) as well as 13 other types of relationships. In these relationships, the two putative 140 grandparents were considered unrelated to each other, but the putative grandparents had different 141 relationships to the putative grandchild. The relationships considered represented combinations 142 of true grandparents, individuals unrelated to the putative grandchild, great-aunts, half-great-143 aunts, and first cousins of the putative grandchild's true grandparent. 144 Estimating the per-comparison false positive error rate for a given value of c has been 145 8 demonstrated for parent-offspring and sibling relationships using importance sampling 146 for trios of a given relationship (such as unrelated) and record how many of them fit the criteria 159 to be considered related. Importance sampling can be thought of as focusing the simulation on 160 producing mostly genotypes that do assign and then correcting for this modification. To 161 implement importance sampling, genotypes were simulated from the distribution of genotypes in 162 true grandparent-grandchild trios. Missing genotypes were accounted for utilizing the forward-163 backward algorithm described below with the state of the Markov chain representing the number 164 of missing genotypes in one individual. If the simulated trio had fewer than m MIs and the 165 calculated LLR was greater than or equal to c, the observation was recorded as a false positive 166 with the appropriate importance sampling weight (Owen 2013). 167 168 Stratified sampling 169 Similar to importance sampling, the goal of stratified sampling was to focus simulation 170 effort on categories (strata) that produce false positives. We stratified the distribution of trio 171 genotypes by the number of observed MIs. This was a natural choice because the algorithm we 172 use explicitly filters possibilities based on the number of observed MIs. Therefore, we can 173 eliminate simulating genotypes for most strata. Sampling effort can then be focused on trios with 174 m or fewer MIs. This method requires calculating the probability that a trio of given relationship 175 has a given number of observed MIs and the ability to simulate genotypes for a trio given a 176 relationship and number of MIs. Genotypes for trios were simulated in each stratum and the false 177 positive rates were recorded. Utilizing the probabilities that a trio has each number of MIs (i.e., 178 the size of the strata), the overall false positive rate was then calculated. 179 To calculate the probability a trio has a given number of observed MIs, we represented 180 the observation of MIs and missing genotypes as a Markov chain and utilized the forward step of 181 the forward-backward algorithm. We extended the approach described by Anderson (2012) to 182 account for the common practice of only analyzing samples given a maximum number of 183 missing genotypes and to fit the target relationships. The maximum number of missing 184 genotypes allowed in the current analyses was 10% of loci. Let be the state, describing the 185 number of observed MIs and missing genotypes, after observing locus i. Prior to observing any 186 loci, 0 = (0,0,0,0), representing the number of observed incompatibilities and number of 187 missing genotypes for the three individuals. Let be a vector indicating whether an 188 incompatibility or any missing genotypes are observed at locus i, in the same order as . 189 Assuming HWE, known allele frequencies, known probability of a genotype being missing, and 190 observation of missing genotypes is independent across loci and individuals, the probabilities of 191 each possibility for can be calculated according to standard probability arguments for a given 192 true relationship. The probability of being in state after a given locus can then be calculated as 193 P(s i+1 = ) = P(s i ) ∑ P(a i+1 )I{s i a + a i+1 = }. (1) This can be evaluated recursively to obtain the probabilities of each final state. The probability a 194 trio has a given number of observed MIs can then be calculated given a maximum number of 195 missing genotypes. However, memory constraints can make saving all the probabilities at each 196 step impractical even with moderate numbers of loci. The probability of observing an individual 197 with more than the allowed number of missing genotypes can be obtained through the same 198 algorithm, but with and now only representing missing genotypes in one individual. Because 199 we assumed that missing genotypes are independent between individuals, the probability of all 200 three samples having a valid number of missing genotypes is then straightforward to calculate 201 and only the probabilities of states with m or fewer MIs and with all three individuals having an 202 allowable number of missing genotypes need to be saved. 203 To utilize stratified sampling, genotypes need to be simulated for trios with a specified 204 number of MIs. The backwards step of the forward-backward algorithm fills this need. Given L 205 loci, a value of is chosen given the number of MIs by sampling a categorical distribution with 206 probabilities proportional to the probability of each that has the specified number of MIs (and 207 allowable number of missing genotypes). Next, for each locus and iterating backwards, a value 208 for is chosen by sampling a categorical distribution with 209 (2) Once is chosen, genotypes are sampled using the genotype frequencies for the true 210 A key question for designing experiments implementing this method is, how many loci 228 need to be genotyped to obtain reliable assignments? To help answer this question, we estimated 229 error rates for panels of different size. We simulated panels containing 100, 300, 500, 700, and 230 900 biallelic single nucleotide polymorphisms (SNPs) and panels containing 100, 200, 300, and 231 400 triallelic microhaplotypes. The SNPs and microhaplotypes were simulated with expected 232 heterozygosity of 0.22 (allele frequencies of 0.125 and 0.872) and 0.42 (allele frequencies of 233 0.13, 0.13, and 0.74), respectively. This choice reflects the mean expected heterozygosities for 234 SNPs and microhaplotypes in a comparison of both marker types for relationship inference in 235 rockfish (Baetscher et al. 2018). The probability of a missing genotype at a locus was set at 3% 236 in these analyses. It is common practice to remove any samples with more than a threshold 237 number of missing genotypes, and so we restricted all simulated genotypes to have 90% or more 238 genotypes present. False negative and false positive, given a true relationship of unrelated, error 239 rates were estimated for a range of c values and the relationship between error rates was 240 compared between panels. Integer values of c were chosen starting at 0 and increasing until the 241 estimated false negative rate was above 0.05. Estimates of false negative rates and estimates 242 from the importance sampling method were derived from 10,000 and 1,000,000 Monte Carlo 243 iterations, respectively. Estimates from the stratified sampling routine were derived from 244 1,000,000 iterations for each stratum (number of observed Mendelian incompatibilities) less than 245 or equal to m. 246 To compare performance of the two methods for estimating false positive error rates, the 247 estimated error rates for unrelated trios were compared between the methods with all four 248 simulated microhaplotype panels. Estimates for other true relationships were compared using the 249 300 locus microhaplotype panel. Finally, to examine the importance of modelling the presence of 250 missing genotypes, we compared importance sampling estimates using the 300 locus 251 microhaplotype panel with the probability of a genotype being missing equal to 3% and 0%. 252 The scripts used to perform these simulations and their outputs are available at 253 https://github.com/delomast/gpError2021. analysis. For each hatchery program, we consider the effect on the desired per-comparison false 266 positive rate of using no data, phenotypic sex, spawn day, phenotypic sex and spawn day, or 267 cross records. 268 In the example analysis, we assume that natural-origin juveniles are sampled and that 269 exact age is unknown but is constrained to one, two, or three years old. The effect of this 270 assumption is that three potential years of parents must be considered for each juvenile. necessarily harmful as the purpose is to identify fish with recent hatchery-origin ancestry. 284 Additionally, the vast majority of trios considered are likely unrelated, and so the false-positive 285 rate for unrelated trios will often be the most impactful. 286

287
Results 288 Estimated error rates declined with increasing panel size, and the microhaplotype panels 289 showed lower error rates than SNP panels of similar size (Figure 2). For example, false positive 290 error rates (true relationship of unrelated) below 10 -10 were achieved at false negative error rates 291 below 0.05 with SNP panels containing 700 or more loci and microhaplotype panels containing 292 300 or more loci. False positive rates for related trios decreased with decreasing relatedness 293 between individuals in the trio (Supplemental Figure 1). 294 False positive rates estimated by importance sampling and stratified sampling were 295 practically identical when estimates were above approximately 10 -8 ( Figure 3). As the false 296 positive rate decreased below this, the importance sampling method estimated false positive rates 297 higher than the stratified sampling method. In these cases, the stratified sampling method 298 estimated false positive rates of 0 within one or more strata (i.e., no false positives were sampled 299 out of the 1,000,000 iterations). Similar results were obtained for false positive rates estimated 300 for trios with relationships other than unrelated (Figure 4) except for two trio types (true 301 grandparent and unrelated; true grandparent and cousin of grandparent) that had noticeably 302 different estimates between the two methods and did not have a low false positive rate. 303 Incorporating a 3% missing genotype rate in the error rate estimation had a moderate 304 15 effect on the results ( Figure 5). For example, at a false negative rate of 0.05, the false positive 305 (true relationship of unrelated) error rates were approximately (derived from linear interpolation 306 with neighboring points) 1.3 • 10 -11 and 2.1 • 10 -12 when missing genotypes had rates of 3% and 307 0%, respectively. 308 The example hatchery programs considered show that analysis of smaller programs (e.g., 309 USB) require per-comparison false positive rates on the order of 10 -6 -10 -8 depending on what 310 data, if any, is available to reduce the number of comparisons (Table 1). Larger programs, similar 311 to DNFH, require rates on the order of 10 -7 -10 -10 if the number of comparisons can be reduced 312 by one of the data sources considered. Error rates estimated for the simulated panels ( Figure 2) 313 suggest that these rates can be achieved with panels containing 500 -700 SNPs or 200 -300 314 microhaplotypes. If no data is available to reduce the number of pairs of grandparents that must 315 be considered, then hatchery programs similar in size to DNFH (1,700 fish spawned / year) will 316 require error rates on the order of 10 -11 , which is achievable with panels containing 700 -900 317 SNPs or 300 -400 microhaplotypes. 318 319 Discussion 320 The methods developed here facilitate direct monitoring of introgression between 321 hatchery-and natural-origin populations without relying on genetic differentiation. This fills an 322 unmet need for monitoring the segregated hatchery programs upon which many fisheries rely. 323 The method identifies fish with an unsampled, hatchery-origin parent using genetic samples from 324 the hatchery broodstock and estimates error rates given allele frequencies for a population. that PBT will continue to grow in usage, and the method described here will become feasible for 336 a larger number of hatchery programs. 337 The simulated microhaplotype panels demonstrated that error rates low enough for 338 relatively large analyses of segregated hatchery programs are achievable with panels containing 339 300 -400 microhaplotypes. The size/power required from the genetic marker panel will depend 340 on the data available to limit the number of comparisons. As has been previously demonstrated, 341 the availability of accurate hatchery cross-records greatly reduces the power required for 342 grandparent inference (Letcher and King 2001), but we note that simply having sex information 343 associated with a broodstock individual's genotype can have a sizeable effect (Table 1) We expect this method to be most applicable for assessing segregated hatchery programs. 357 For some closely related trios, false positive error rates were relatively high with the 300 358 microhaplotype panel (Supplemental Figure 1). For all analyses, the impact of these error rates is 359 mediated by the infrequency of comparisons involving closely related individuals. When 360 analyzing segregated programs, false positives from closely related trios may not negatively 361 impact conclusions as they indicate recent hatchery ancestry. However, when analyzing 362 integrated hatchery programs, distinguishing trios with different relationships can be important. 363 Application of this method could then either be infeasible or require a more powerful genetic 364 panel. Additionally, the number of relationships for which we provide estimators covers trios 365 with a range of relatedness but is not exhaustive. In some specific cases, other relationships may 366 be present at impactful frequencies. For example, if generations overlap substantially then the 367 impact of trios containing aunts, half-aunts, and first-cousins may need to be considered. 368 The method developed here addressed shortcomings of previously developed methods for 369 grandparent inference with moderately sized genetic panels. Previous methods required either 370 that all four grandparents were sampled (Letcher and King 2001) or used an exclusionary method 371 that inherently had lower power than likelihood-based methods (Christie et al. 2011). 372 Additionally, the exclusionary method did not provide a formal treatment of genotyping error. methods showed that when the true relationship was unrelated, they estimated essentially the 389 same false positive rate when that rate was above approximately 10 -8 . This suggests the 390 importance sampling routine performed well for unrelated trios. At lower false positive rates, the 391 stratified sampling estimates were lower (and in some cases 0) than those from importance 392 sampling. In these cases, one or more strata had estimates of 0, indicating that the 1,000,000 393 samples taken were not sufficient to observe one or more false positives. This demonstrates one 394 of the strengths of importance sampling compared to stratified sampling -some situations will 395 have false positive rates small enough they cannot be efficiently estimated by this stratified 396 19 sampling routine. This is further emphasized by the greater computational effort devoted to 397 stratified sampling in this study (1,000,000 iterations in each stratum vs. 1,000,000 iterations 398 total). 399 For false positive error rates under relationships other than unrelated, the importance 400 sampling and stratified sampling estimates were again largely the same. In a few cases where one 401 potential grandparent was the true grandparent and the other was not, the estimates were 402 noticeably different and the stratified sampling method achieved at least 150 false positive 403 observations in each strata. This suggests that for these relationships the importance sampling 404 method may have performed suboptimally. Similar observations were made for importance 405 sampling estimates of false positive error rates in parentage inference, where performance of a 406 given importance distribution varied depending on the true relationship for which an error rate 407 was estimated (Anderson and Garza 2006). One strategy would be to design an alternative 408 importance distribution, but for closely related trios, error rates for most panels will likely be 409 high enough that they are amenable to stratified sampling. 410 The current method assumes loci are in linkage equilibrium. If loci are physically linked 411 (but still in equilibrium) the methods described here can be applied. Physical linkage of two loci 412 will result in grandchildren being more likely to inherit alleles at both loci from one of the two 413 grandparents in a trio. The estimated false positive (when true relationship is unrelated) error 414 rates are not affected because unrelated individuals are not impacted by physical linkage. The 415 Monte Carlo method we implement for estimating false negative error rates simulates trios 416 assuming no linkage, but given the alleles present at two loci in the grandparents, the probability 417 of a given combination of alleles being inherited by the grandchild is identical regardless of 418 physical linkage (Supplemental file 1). Therefore, as long as the set of alleles present in the 419 20 grandparents (ignoring which grandparent is assigned which alleles) is simulated accurately, the 420 frequency of combinations of grandparent alleles and grandchild alleles will be simulated 421 accurately by assuming the loci are unlinked. The calculated LLR does not rely on assignment of 422 genotypes to specific grandparents, and so the distribution of LLR, and thus the false negative 423 error rate, is also simulated accurately by assuming the loci are unlinked. 424 One drawback of the method as implemented is that computational efficiency decreases 425 with increasing numbers of alleles per locus. This is partly due to the flexibility of the 426 genotyping error model and the need to marginalize over all possible true genotypes. In the 427 current implementation, computation is sped up by precomputing likelihood values for observed 428 trio genotypes at each locus. This works well when the number of alleles per locus is small, but 429 as the number of alleles increases, this can become impractically slow. As such, the current 430 implementation will be most suitable to panels containing biallelic SNPs and microhaplotypes, 431 which in our experience typically have three to five alleles. For panels of highly variable loci, a 432 different implementation of this method, particularly with a more streamlined genotyping error 433 model, would be necessary.   relationship labels indicate the relationship of the two putative grandparents to the putative 587 grandchild. Where the label is only one relationship, both putative grandparents have the same 588 relationship to the putative grandchild. True is true grandparent, GAunt is great-aunt, HGAunt is 589 half great-aunt, GpCous is first-cousin to the putative grandchild's true grandparent, Unrel is 590 unrelated. The inset gives a magnified view of the region containing error rates close to 1. 591 Figure 5. Error rates estimated by importance sampling for the simulated 300 locus 592 microhaplotype panel with two rates of missing genotypes: 3% (Missing) and 0% (No missing). 593 False positive error rate is the rate for unrelated trios. 594 Each locus is treated separately and each genotype is treated as two observations (one for each allele in the 4 presumed diploid). For a given locus, the error model first relies upon a "per allele error rate", , which 5 represent the probability that, when observing one allele, you observe any allele other than the true one.

6
The default value is 0.005, or 0.5%. This value is used along with a measure of similarity between alleles to 7 calculate the probability of observing each allele given the true allele. For a locus with I alleles, given the 8 true allele i, the probability of correctly observing i is 1 − and the probability of observing allele j where where s ij is the similarity between alleles i and j and s ii = 0. The measure of similarity between two alleles 11 of a microhaplotype used here was the reciprocal of the number of base pair differences between them.

12
These probabilities are then used to calculate the probability of observing each genotype given a true according to basic probability arguments, where P (B|A) is the probability of observing B when the true allele 15 is A. The probability of observing each possible genotype given a true heterozygous genotype is calculated 16 similarly, but also incorporates a probability of allelic dropout for each allele. Each allele is given a probability 17 of dropping out, d i for allele i, and in the current study a probability of 0.005 was used for all alleles.

18
It is assumed that allelic dropout events are disjoint, and so the probability of no dropout occurring is The probability of observing each genotype given no dropout can be calculated using the 20 observation probabilities for each allele calculated in equation 1. The probability of observing a given 21 genotype if a dropout has occurred is taken to be the probability of observing that genotype given a true 22 homozygous genotype for the remaining allele. The total probability of observing a genotype is the sum of 23 1 the probabilities of all three cases (no dropout, allele 1 dropout, allele 2 dropout). The probability of case 1 is 0.25, because the loci are in linkage equilibrium, and so B 1 is equally likely to 32 be on the same chromosome as any of the alleles at locus 1. Similarly, the probability of cases 2 and 3 are 33 0.25 and 0.5, respectively.

34
Under case 1, the probability the unsampled parent inherits A 1 and B 1 is 0.5(1 − r), where r is the 35 probability of a recombination event and the probability the grandchild then inherits A 1 and B 1 is 0.5(1 − r).

36
So, the probability the grandchild inherits A 1 and B 1 under case 1 is 0.25(1 − r) 2 .

37
Under case 2, the probability the unsampled parent inherits A 1 and B 1 is 0.5r and the probability the 38 grandchild then inherits A 1 and B 1 is 0.5(1 − r). So, the probability the grandchild inherits A 1 and B 1 under 39 case 2 is 0.25r(1 − r).

40
Under case 3, the probability the unsampled parent inherits A 1 and B 1 is 0.5 · 0.5 = 0.25 and the 41 probability the grandchild then inherits A 1 and B 1 is 0.5r. So, the probability the grandchild inherits A 1 42 and B 1 under case 3 is 0.125r.