Livestock Science 251 (2021) 104629 Contents lists available at ScienceDirect Livestock Science journal homepage: www.elsevier.com/locate/livsci Low-coverage whole-genome sequencing in livestock species for individual traceability and parentage testing Joaquim Casellas a, *, Melani Martín de Hijas-Villalba a, Marta Vázquez-Gómez a, Samir Id-Lahoucine b a Departament de Ciència Animal i dels Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain b Animal and Veterinary Science Group, Scotland’s Rural College, Edinburgh EH9 3JG, United Kingdom H I G H L I G H T S • Traceability and paternity tests can adapt to low-coverage whole-genome sequencing data. • Testing performance depended on sequencing error rate and genotype frequencies. • Uncertainly had greater impact on false negatives than false positives. • 0.05× coverage sufficed to guarantee greater-than-99% success during testing. A R T I C L E I N F O A B S T R A C T Keywords: Procedures for genetic traceability of animal products and parentage testing mainly focus on microsatellites or Evidential statistics SNPs panels. Nevertheless, current availability of high-throughput sequencing technologies must be considered Paternity as an appealing alternative. This research focused on the evaluation of low-coverage whole-genome sequencing Sequencing for traceability and paternity testing purposes, within a context of evidential statistics. Analyses were performed Simulation Traceability on a simulation basis and assumed individuals with 30 100-Mb/100-cM chromosome pairs and ~1,000,000 polymorphic SNPs per chromosome. Ten independent populations were simulated under recombination and mutation with effective populations size 100 (generations 1–1000), 10,000 (generation 1001) and 25,000 (generation 1002), and this last generation was retained for analytical purposes. Appropriate both traceability and paternity tests were developed and evaluated on different high-throughput sequencing scenarios accounting for genome coverage depth (0.01×, 0.05×, 0.1× and 0.5×), length of base-pair reads (100, 1000 and 10,000 bp), and sequencing error rate (0%, 1% and 10%). Assuming true sequencing error rates and genotypic frequencies, 0.05× genome coverage depth guaranteed 100% sensitivity and specificity for traceability and paternity tests (n = 1000). Same results were obtained when sequencing error rate was arbitrarily set to 0, or the maximum value assumed during simulation (i.e., 1%). In a similar way, uncertainly about genotypic frecuencies did not impair sensitivity under 0.05× genome coverage, although it reduced specificity for paternity tests up to 85.2%. These results highlighted low-coverage whole-genome sequencing as a promising tool for the livestock and food in­ dustry with both technological and (maybe) economic advantages. 1. Introduction 2001), and is becoming more demanding by consumers and producers (Qian et al., 2020). On the other hand, parentage testing enables to Neutral genetic markers have been widely used for both traceability identify similar inheritance patterns between related individuals (Arana et al., 2002; Vázquez et al., 2004) and parentage testing (Heaton (Jamieson, 1965), and has a deep impact on breeding programs (Banos et al., 2014) in livestock populations. Traceability aims to maintain et al., 2001), where a moderate proportion of misidentified progeny can credible custody of identification for animals or animal products be anticipated (Geldermann et al., 1986; Visscher et al., 2002; Weller through various steps within the production and food chain (McKean, et al., 2004). Both approaches have relevant legal uses for animal * Corresponding author. E-mail address: joaquim.casellas@uab.cat (J. Casellas). https://doi.org/10.1016/j.livsci.2021.104629 Received 24 March 2021; Received in revised form 4 July 2021; Accepted 13 July 2021 Available online 17 July 2021 1871-1413/© 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). J. Casellas et al. Livestock Science 251 (2021) 104629 forensic determinations (Kanthaswamy, 2015) or pedigree certification 2.2. Sequencing and alignment simulation process regarding livestock breed societies. Genetic traceability and parentage testing rely on the fact that DNA is The number of reads per chromosome was defined as enormously variable among individuals despite the simple genetic C × (100×106) / L, where 100×106 was the assumed chromosome mechanisms ruled by Mendel’s laws of inheritance from parents to length in bp, C was the expected genome coverage, and L was the offspring. Moreover, DNA is present in every cell of the organism, does average read length in base-pairs. The length of each read was sampled not change during animal life, and is stable to different treatments of from a normal distribution with mean L and standard deviation L/10 to processed food (Dalvit et al., 2007). Current procedures for genetic account for variability on DNA sequencing products. Moreover, each traceability and parentage testing mainly focus on microsatellites or read was placed at random, both in the genome and chromosome phase. SNPs (Heaton et al., 2002), where standardized panels have already Following Fox et al. (2014) and Pfeiffer et al. (2018), an error rate be­ been defined to harmonize procedures worldwide (https://www.isag.us tween 10− 5 and 10− 2 was randomly assigned to each polymorphic SNP. /committees.asp, accessed March 18th, 2021). Nevertheless, current The same error rate applied to both alleles. Only the number of reads for advances in high-throughput technologies move towards partial or each allele was stored for further analyses. whole-genome sequencing procedures where closed SNP panels would be likely to have no future for further purposes. This requires additional 2.3. Evidential testing for single-individual traceability endeavors to elucidate the usefulness of sequencing data, mainly when low-coverage approaches are considered due to economic limitations. This research relied on evidential inference (Edwards, 1972) as a Although Zan et al. (2019) suggested that very low-coverage (<0.5×) way to compare two competing hypotheses (i.e., models). This approach sequencing data could be informative enough for inferring outbred relies on the likelihood function as the structure that contains all evi­ founder genotypes under an F2 design, little is known about their dence from the data relevant to the statistical model (Birnbaum, 1962), applicability in commercial populations of livestock. and compares hypotheses by calculating the ratio of their likelihood This research focused on three main objectives, (1) the development functions (Hacking, 1965). Within this context, an upper-than-1 likeli­ of both traceability and paternity tests for low-coverage sequencing data hood ratio favors the numerator model whereas a lower-than-1 likeli­ within the context of evidential inference (Bickel, 2012), (2) the vali­ hood ratio suggests the superiority of the denominator model, although dation of low-coverage sequencing for traceability and paternity testing a minimum likelihood ratio of 32 (or 1/32) is typically used in the in commercial livestock populations under full knowledge of population evidential literature (Blume, 2002; Royall, 1997), or even as high as (i.e., allele or genotype frequencies) and sequencing parameters (i.e., 1000, often used in genome-wide linkage studies (Morton, 1998). error rates), and (3) the evaluation of the impact of uncertainly about Traceability in the livestock industry can be defined as the ability to population and sequencing parameters on traceability and paternity identify animals or animals products through various steps within the tests for low-coverage sequencing data. food chain from the farm to the retailer (McKean, 2001). Within this context, the analysis of genetic polymorphisms must be viewed as a key 2. Materials and methods tool to verify the match between two independent samples. Take as a starting point a n × 2 matrix (S) to summarize sequence Animal Care and Use Committee approval was not obtained for this data from n polymorphic and biallelic sites of the genome. Once sorted study because analyses were performed on simulated data sets. Neither by chromosome and nucleotide within the chromosome, each column real animals nor biological tissues from alive animals were involved in stores the number of reads for alleles A and B, respectively. The analysis this research. of genetic traceability relied on two different samples (Sp and Sq), and two competing hypotheses, i.e., H0: samples belong to the same indi­ vidual (p = q), and H1: both samples belong to different individuals (p ∕= 2.1. Genome and population simulation process q). They can be tested through their likelihood ratio (Edwards, 1972) as follows, This research simulated an unspecific mammalian livestock popu­ lation. We took as a starting point a 100-Mb/100-cM chromosome with LR(H0, H1 | Sp, Sq) = p(Sp, Sq | H0) / p(Sp, Sq | H1), 5000,000 biallelic SNPs (one SNP each 20 base pairs and 2 × 10− 5 cM), and the whole genome consisted of 30 chromosome pairs. This gener­ where p(Sp, Sq | Hk) was the joint probability of obtaining data Sp and Sq ated a standard 3 Gb genome (Pérez-Encisoet al., 2015) with the same under hypothesis Hk. Under the H0 hypothesis, the likelihood must number of chromosomes as cattle and goat, and within the range of expand to other livestock species such as pig (19 pairs), sheep (27 pairs) and horse (32 pairs). The starting number of SNPs was assumed to guarantee more p(Sp, Sq | H0) = Πi=1,n p(sp,i | gp,i, εi) p(sq,i | gp,i, εi) p(gp,i) than 30,000,000 polymorphic SNPs at the end of the simulation process (see below), as reported by Daetwyler et al. (2014) in cattle. where sp,i was the ith row of Sp, gp,i was the genotype (i.e., AA, AB or BB) Populations started from a founder generation with 100 individuals of the pth individual for the ith polymorphic site, and εi was the that were heterozygous throughout the whole genome. They evolved sequencing error rate for the ith polymorphic site (we assume homo­ during 1000 non-overlapping generations under random mating and geneous error rates among alleles). Given that p and q were assumed to effective population size 100. Linkage disequilibrium between adjacent be the same individual and the parametric space accounted for three loci was generated based on Kosambi’s mapping function (Kosambi, genotypes (p(AA) + p(AB) + p(BB) = 1), the likelihood became 1944), and a mutation rate of 2.5 × 10− 3 per SNP was applied until generation 980 (Meuwissen et al., 2001), switching the allele state from p(Sp, Sq | H0) = Πi=1,n[Σα=AA,AB,BB p(sp,i | gp,i = α, εi) p(sq,i | gq,i = α, εi) p A to B, or vice versa. From generation 981 on, the mutation rate switched (gp,i = gq,i = α)] to 2.5 × 10− 8 (Hickey and Gorjanc, 2012). Only those populations with 1000,000 ± 10% (i.e., 900,000 to 1100,000) polymorphic (MAF > 0) Now, assume a reads for allele A and b reads for allele B in sp,i. The SNPs per chromosome in generation 1000 were retained for further following conditional probabilities can be straightforwardly calculated analyses. as binomial processes with trials, successes and success probability Populations expanded to 10,000 individuals in generation 1001 sequentially noted between parentheses, (1000 sires and 9000 dams), and 25,000 in generation 1002. A total of 10 independent populations were simulated. p(sp,i | gp,i = AA, εi) = Binomial(a + b, a, 1 – εi) 2 J. Casellas et al. Livestock Science 251 (2021) 104629 p(sp,i | gp,i = AB, εi) = Binomial(a + b, a, 0.5) × [Σγ=AA,AB,BB p(go,i = α | gs,i = γ) p(ss,i | gs,i = γ, εi) p(gs,1,i = γ)]}, p(sp,i | gp,i = BB, εi) = Binomial(a + b, a, εi), where p(go,i = α | gs,i = γ) can be obtained from Table 2. Finally, the probability of each genotype depends on its frequency in On the other hand, the alternative hypothesis could be defined on the the source population. following rationale, H1: only d was parent of o, whereas s was unrelated The same development can be applied to the alternative hypothesis to o and sampled from the same population. The likelihood expands to where p(So, Sd, Ss | H1) = Πi=1,n {p(so,i | go,i, εi) p(go,i | gd,i) p(Sp, Sq | H1) = Πi=1,n p(sp,i | gp,i, εi) p(gp,i) p(sq,i | gq,i, εi) p(gq,i), × p(sd,i | gd,i, εi) p(gd,i)} p(ss,i | gs,i, εi) p(gs,i), and p and q were assumed different and unrelated individuals from the and same population. Once accounted for all three possible genotypes, the previous expression expanded to p(So, Sd, Ss | H1) = Πi=1,n {Σα=AA,AB,BB p(so,i | go,i = α, εi) × [Σβ=AA,AB,BB p(go,i = α | gd,i = β) p(sd,i | gd,i = β, εi) p(gd,i = β)]} p(Sp, Sq | H1) = Πi=1,n {[Σα=AA,AB,BB p(sp,i | gp,i = α, εi) p(gp,i = α)] × Πi=1,n [Σγ=AA,AB,BB p(ss,i | gs,i = γ, εi) p(gs,i = γ)]. × [Σβ=AA,AB,BB p(sq,i | gq,i = β, εi) p(gq,i = β)]}. where p(go,i = α | gd,i = β) can be found in Table 2. As for previous hy­ pothesis, it was not mandatory to account for dam sequencing data if 2.4. Testing for parentage missing, Parentage testing relies on the use of biological markers to identify p(So, Ss | H1) = Πi=1,n p(so,i | go,i, εi) p(go,i) p(ss,i | gs,i, εi) p(gs,i), similar inheritance patterns between related individuals and traces back p(So, Ss | H1) = Πi=1,n Σα=AA,AB,BB p(so,i | go,i = α, εi) p(go,i = α) to the 1960s where blood typing was used as a regular part of some cattle × Σγ=AA,AB,BB p(ss,i | gs,i = γ, εi) p(gs,i = γ). breeding programs (Stormont, 1967). As seen with most domestic spe­ cies, the typical animal parentage case includes a dam, offspring, and 2.5. Uncertainly about population and sequencing parameters one or more alleged sires. The identity of the dam uses to be fairly certain, whereas the true sire must be identified from a set of m males. Single individual traceability and parentage testing were evaluated Our analytical approach will rely on this scenario, although it can be under different scenarios accounting for 0.01×, 0.05×, 0.1× and 0.5× straightforwardly generalized to test the other sex (i.e., dam). depth of genome coverage, with 100, 1000 and 10,000 base-pair reads. Paternity testing relied on data samples from the offspring (So), its Those read lengths were chosen to illustrate test performance under dam (Sd), and an alleged sire (Ss). The testing process started with the currently available sequencing platforms (Besser et al., 2018). definition of the null hypothesis such as H0,j: both s and d were parents of As noted above, both traceability and parentage tests relied on two o. Within this context, the joint likelihood of So, Sd and Ss was written as structural parameters, within-SNP sequencing error rate (εi) and geno­ p(So, Sd, Ss | H0) = Πi=1,n{p(so,i | go,i, εi) p(go,i | gd,i, gs,i) p(sd,i | gd,i, εi) p typing frequencies. The first mainly depends on the sequencing method (gd,i) and platform used (Fox et al., 2014) and uses to be estimated on an × p(ss,i | gs,i, εi) p(gs,i)}, where so,i was the ith row of So, go,i was the across-genome basis. Within this context, we compared test perfor­ genotype of the oth individual in the ith polymorphic site, and εi was the mances under three across-SNP homogeneous sequencing error rates: sequencing error rate for the ith polymorphic site (we assume homo­ 0%, 1% (the maximum sequencing error rate used for simulation of the geneous error rates among alleles). As for traceability tests, previous sequencing process), and 10% (i.e., ten times higher than the maximum likelihood expanded to account for biallelic genetic markers, sequencing error rate used for simulation of the sequencing process). On the other hand, genotypic frequencies could be approximated by p(So, Sd, Ss | H0) = Πi=1,n {Σα=AA,AB,BB p(so,i | go,i = α, εi) using sequence data generated for traceability and paternity tests. × [Σβ=AA,AB,BB Σγ=AA,AB,BB p(go,i = α | gd,i = β, gs,i = γ) p(sd,i | gd,i = β, Nevertheless, the number of sequenced animals could be small and εi) p(gd,i = β) contribute high uncertainly to estimated genotypic frequencies. To ac­ × p(ss,i | gs,i = γ, εi) p(gs,1 = γ)]}, count for this uncertaintly, the variance of the estimated A allele fre­ quency (π) can be calculated as (Cockerham, 1969) where p(so,i | go,i = α, εi), p(sd,i | gd,i = β, εi) and p(ss,i | gs,i = γ, εi) were V(π) = [π (1 - π)] / 2λ where λ was the number of sampled in­ binomial probabilities, p(gd,i = β) and p(gs,i = γ) were genotypic fre­ dividuals. We compared λ = 5, 10 and 100, and sampled the A allele quencies in the parental population, and p(go,i = α | gd,i = β, gs,i = γ) was frequency (π*) for each SNP from a truncated (0–1) normal distribution the conditional probability of the offspring’s genotype depending on with mean π and variance V(π). Genotypic frequencies were obtained parents’ genotype (Table 1). It is important to note that previous assuming Hardy-Weinberg equilibrium (Hardy, 1908). expression can also be applied when lacking of sequencing data from the dam as follows, p(So, Ss | H0) = Πi=1,n{p(so,i | go,i, εi) p(go,i | gs,i) p(ss,i | gs,i, εi) p(gs,i)}, p(So, Ss | H0) = Πi=1,n {Σα=AA,AB,BB p(so,i | go,i = α, εi) Table 2 Table 1 Conditional probability of the offspring’s genotype in a biallelic locus (alleles A Conditional probability of the offspring’s genotype in a biallelic locus (alleles A and B) when only one parent contributes to the paternity test. Each triad of and B) given the mother’s and the alleged father’s genotype. Each triad of numbers provides the probability for AA, AB and BB genotypes, respectively. numbers provides the probability for AA, AB and BB genotypes, respectively. Parent’s genotype Mother’s genotype Offspring’s genotype AA AB BB Father’s genotype AA AB BB AA p(A)1 0.5 p(A) 0 AA 1/0/0 0.5 / 0.5 / 0 0/1/0 AB 1-p(A) 0.5 p(A) AB 0.5 / 0.5 / 0 0.25 / 0.5 / 0.25 0 / 0.5 / 0.5 BB 0 0.5 [1 - p(A)] 1 - p(A) BB 0/1/0 0 / 0.5 / 0.5 0/0/1 1 p(A): allelic frequencies of A allele in parents’ generation. 3 J. Casellas et al. Livestock Science 251 (2021) 104629 3. Results 3.1. Simulated genomic data After 1000 non-overlapping generations, random mating and effec­ tive population size 100, we retained ten populations with 29,195,811 to 30,660,474 polymorphic SNPs. Allele frequencies widely distributed along with the parametric space, as shown in Fig. 1, and a remarkable percentage of SNPs had minimum allele frequency (MAF) below 0.05. Although this varied among chromosomes, between 36.7% and 51.9% of SNPs had MAF < 0.05. All these 10 simulated populations contributed equally to the subsequent analyses. After sequencing 10,000 individuals, the maximum number of reads per polymorphic SNP was 3 (0.01× genome coverage), 4 (0.05×), 5 Fig. 2. Average distribution of polymorphic SNPs depending on the number of (0.1×) and 7 (0.5×). Nevertheless, between 76.2% (0.5× genome reads when sequenced at 0.01X (black), 0.05X (blue), 0.1X (red) and 0.5X coverage) and 99.5% (0.01× genome coverage) of them had a single (white) genome coverage with 10,000 base-pair read length. The whiskers read, as shown in Fig. 2. The percentage of polymorphic SNPs with two extend to minimum and maximum estimates. (For interpretation of the refer­ reads increased with genome coverage, from 0.5% (0.01×) to 19.1% ences to color in this figure legend, the reader is referred to the web version of (0.5×), and a similar trend with smaller percentages was revealed for this article.) larger numbers of reads. Moreover, those percentages showed small variability across individuals, this uncertainly even reducing for smaller read length (Fig. 2). The same pattern was revealed when checking for shared SNPs among pairs of sequenced individuals. The longer the read length was, the wider the dispersion of the number of shared SNPs (Fig. 3). From the total of ~30,000,000 polymorphic SNPs, the average number of shared polymorphic SNPs decreased from 3355.6 ± 4.7 (100 base-pair read length) to 3093.2 ± 28.5 (10,000 base-pair read length). For SNPs with MAF ≥ 0.05, similar trends were observed, from 1748.4 ± 3.3 (100 base-pair read length) to 1586.9 ± 17.2 (10,000 base-pair read length). Within this context, subsequent results were reported based on the most uncertain (i.e., increased variability for the number of reads and shared SNPs) and less informative (i.e., reduced number of shared SNPs) scenario, this accounting for sequencing by 10,000 base- pair reads. Fig. 3. Shared SNPs among two unrelated individuals both sequenced at 0.01X 3.2. Traceability and parentage testing genome coverage with 100 (red dots), 1000 (blue dots) and 10,000 base-pair read length (black dots). The X-axis accounts for SNPs with non-zero mini­ As anticipated, the number of shared polymorphic SNPs among two mum allele frequency (MAF), whereas Y-axis accounts for SNPs with unrelated individuals quickly increased with genome coverage MAF>0.05. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) (Table 3). This generated a fast growth in terms of available information for traceability and paternity tests, as evidenced by the likelihood ratios provided in Fig. 4. Assuming true sequencing error rates and genotype Table 3 frequencies, 100% of traceability tests favored the true hypothesis when Mean ± SE of shared polymorphic SNPs among two unrelated individuals when genome coverage was 0.05× or deeper. The only exceptions were sequenced at different genome coverages with 10,000 base-pair read length. detected for 0.01X genome coverage, where 0.7% of false positives and Genome coverage Polymorphic SNPs MAF1>0.05 SNPs 0.04% of false negatives were reported (Fig. 4). The same pattern was 0.01X 3093.2 ± 28.5 1586.9 ± 17.2 revealed for paternity tests, they showing a 100% of true positive and 0.05X 77,007.2 ± 75.0 39,887.3 ± 46.4 true negatives under genome coverage 0.05X or deeper, and 1.1% 0.1X 290,845.5 ± 137 151,872.8 ± 86.5 (1.0%) of false positives and 0.5% (0.6%) of false negatives under 0.01X 0.5X 4965,993.9 ± 450.0 2589,401 ± 329.2 1 Minimum allele frequency. of genome coverage when the dam was known (unknown). In order to test for a more realistic scenario, different homogeneous sequencing error rates were evaluated. As shown in Fig. 5, 0.05X coverage sufficed to avoid false positives and negatives under both traceability and paternity tests when sequencing error rate was arbi­ trarily set to 0% or the maximum rate used during sequencing simula­ tion (i.e., 1%). The only assumption that generated wrong results under 0.05X coverage was when the sequencing error rate was unrealistically assumed 10 times higher than the upper bound during sequencing (i.e., 10%). In this case, 1.3% (traceability test), 24.5% (paternity test with known dam) and 22.8% (paternity test with an unknown dam) of false negatives were reported, whereas any test generated false positives. Fig. 1. Distribution of allele frequencies for the first chromosome of the first Higher genome coverage tested provided 100% of true positives and simulated population. 4 J. Casellas et al. Livestock Science 251 (2021) 104629 Fig. 4. Distribution of 1000 traceability tests (black dots) and paternity tests with known (blue dots) and unknown dam (red dots) under four different genome coverage, and assuming true SNP-specific sequencing error rate and true genotype fre­ quencies in parental generation. Traceability tests relied on the likelihood ratio between the null (H0: same individual) and the alternative hypothesis (H1: different individuals), and compared each individual against itself (X-axis) and against an unrelated individual (Y-axis). Paternity tests evaluated whether the alleged sire was the true sire (H0) or an unrelated male of the population (H1), and where applied on the true sire (X-axis) and on an unrelated male of the population (Y-axis). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) true negatives (results not shown), even under the assumption of 10% of require genotype-calling procedures to reach closed genotypes (Nielsen the sequencing error rate. et al., 2011). Nevertheless, genotype-calling approaches show little The other parameter accounting for uncertainly during traceability agreement when compared under low-coverage sequencing data (Liu and paternity testing was genotype frequencies. In this case, genotype et al., 2013; Vens et al., 2009; Yu and Sun, 2013), where heterozygous frequencies were assumed under Hardy-Weinberg equilibrium and genotypes cannot be adequately called with a single read (Brouard et al., calculated from allele frequency with uncertainly as sampled from 5, 10, 2017). Within this context, we omitted genotype-calling approaches in and 100 individuals. As shown in Fig. 6, the smaller the uncertaintly for our traceability and paternity tests and focused on genotype probabili­ allele frequency was, the larger the match with results was obtained ties within the context of appropriate likelihood functions. Although under true genotype frequencies. Nevertheless, 0.05× coverage sufficed these procedures were partially implemented in some genotype-calling to avoid both false positives and negatives in traceability tests, whatever approaches (Li et al., 2008, 2009; Martin et al., 2010), they summa­ the accuracy of allele frequencies. Paternity tests with known dam rized to the most probable genotype instead of keeping uncertainly for (unknown dam) revealed a similar pattern without false positives since further analyses. We kept uncertainly about genotypes along the whole 0.05× coverage, and 0.1% (0.1%), 0.6% (0.3%) and 10.2% (14.8%) of calculation of the likelihood ratio in order to avoid arbitrary decisions false negatives when the allele frequency was sampled with uncertainly when available information for each polymorphic site was very small in as calculated from 100, 10 and 5 individuals. Deeper genome coverage tested individuals (Fig. 2). provided 100% of true positives and true negatives at any uncertainly Our tests relied on the likelihood principle, a statistical proposition for allele frequency. that states that all the evidence in the data relevant to the statistical model is contained in the likelihood function (Birnbaum, 1962). Within 4. Discussion this context, a likelihood ratio must be viewed as an objective mea­ surement of the statistical evidence of one model against the other Current procedures for traceability and paternity testing rely on (Hacking, 1965), and establishes the foundations for the evidential SNPs where standardized panels have already been defined to harmo­ statistics (Edwards, 1972) in contrast with frequentist and Bayesian nize procedures worldwide (Heaton et al., 2002). Although their reli­ statistics. This inferential approach relies on two basic conditions that ability and statistical power fulfill the purpose for which they were are not completely fulfilled by frequentist and Bayesian inferences, ob­ created (Marshall et al., 1998), they depend on some dozens of a few jectivity (i.e., the strength of evidence does not vary from one researcher hundreds of SNP genotypes, too few to be reused for other purposes like to another) and interpretability (i.e., the strength of evidence has the genome-wide association analyses (Klein et al., 2005; Gilly et al., 2019) same practical interpretation for any sample size). The first condition or genomic evaluation (Meuwissen et al., 2001; Gorjanc et al., 2015, rules out Bayes factors that depend on subjective or default priors 2017). This is an important limitation because it drains the economic (Bickel, 2012), and the second rules out the frequentist p-value that capacity of food chain industries and breed societies and precludes forces the same type-I error percentage at any sample size (Bickel, additional investments in genomic techniques. The current explosion in 2011). By contrast, the likelihood ratio satisfies both of the necessary high-throughput sequencing technologies (Bansal and Boucher, 2019) conditions for a measure of the strength of statistical evidence. Within opens the door to more sustainable science where specificity and this context, the likelihoods used in our testing approaches had the same multiple-purpose data are not conflicting terms. Nevertheless, a first step mathematical structure than the likelihoods we could construct within a is required to verify that low-coverage whole-genome sequencing data frequentist scenario, as well as they are proportional to the joint pos­ can efficiently address both traceability and paternity tests in order to terior distributions with flat priors we could call in the Bayesian fulfill current standards at a similar economic cost. framework. The essential difference relies on the test itself and the as­ Theoretical approaches to test both traceability and paternity have sumptions carried out by the researcher. Within the context of evidential been widely developed in scientific literature on the basis of complete statistics, there are not additional assumptions apart from the statistical genotypes (Goffaux et al., 2005; Martin et al., 2010; Marshall et al., model itself and all the hypotheses have the same consideration during 1998), whereas high-throughput sequencing technologies provide a the analytical process. Indeed, paternity tests with panels of genetic variable number of random samples from each polymorphic site and markers were previously proposed by Marshal et al. (1998), and 5 J. Casellas et al. Livestock Science 251 (2021) 104629 Fig. 5. Distribution of 1000 traceability tests (upper panel), paternity test with Fig. 6. Distribution of 1000 traceability tests (upper panel), paternity test with known dam (mid panel) and paternity test with unknown dam (lower panel) known dam (mid panel) and paternity test with unknown dam (lower panel) under 0.05X genome coverage, 10,000 base-pair read length, and assuming true under 0.05X genome coverage, 10,000 base-pair read length, and assuming true genotype frequencies in parental generation. Tests assumed true sequencing sequencing error rates per SNPs. Tests assumed true genotyping frequencies error rates (black dots), null sequencing error rate (green dots), 1% sequencing (black dots), as well as genotyping frequencies under Hardy-Weinberg equi­ error rate (blue dots), and 10% sequencing error rate (red dots). Traceability librium after sampling the allele frequency from 5 (red dots), 10 (blue dots) and tests relied on the likelihood ratio between the null (H0: same individual) and 100 individuals (green dots). Traceability tests relied on the likelihood ratio the alternative hypothesis (H1: different individuals), and compared each in­ between the null (H0: same individual) and the alternative hypothesis (H1: dividual against itself (X-axis) and against an unrelated individual (Y-axis). different individuals), and compared each individual against itself (X-axis) and Paternity tests evaluated whether the alleged sire was the true sire (H0) or an against an unrelated individual (Y-axis). Paternity tests evaluated whether the unrelated male of the population (H1), and where applied on the true sire (X- alleged sire was the true sire (H0) or an unrelated male of the population (H1), axis) and on an unrelated male of the population (Y-axis). (For interpretation of and where applied on the true sire (X-axis) and on an unrelated male of the the references to color in this figure legend, the reader is referred to the web population (Y-axis). (For interpretation of the references to color in this figure version of this article.) legend, the reader is referred to the web version of this article.) evidential statistics have been growing attention in genetics and geno­ mild impact on the performance of both traceability and paternity tests, mics research (Strug et al., 2010; Strug, 2018). and it only impaired their results when an abnormally high sequencing The performance of both traceability and paternity tests was error rate was assumed (i.e., 10%). Indeed, results shown in Fig. 5 sug­ outstanding as evidenced in Figs. 4–6. Under the unrealistic assumption gested that the assumption of a null sequencing error rate provided the of known sequencing error rates and genotype frequencies, 0.01× most similar results to the ones obtained under true sequencing error genome coverage sufficed to guarantee ≥99% of true positives and true rates, simplifying both analytical models and subsequent calculations. negatives under traceability tests. In contrast, the minimum genome On the other hand, the impact of genotype frequencies was suggested as coverage for paternity tests must increase up to 0.05× genome coverage larger, where more accurate estimates were required to avoid false to reach the same rate of true positives and negatives. Nevertheless, our positives and negatives. current method works with low-coverage sequencing data and less false Statistical methodologies developed in this manuscript are ready to paternity assignments than previous methods found in the scientific use for both the food chain industry and breed societies. In fact, they literature (Snyder-Mackler et al., 2016; Whalen et al., 2019). The could also be useful for human studies. They do not need additional method design for very low sequencing coverage data from fecal-derived generalizations, as all required algorithms are detailed in the current DNA by Snyder-Mackler et al. (2016), which also performed peternity manuscript. It is important to highlight that 0.05X genome coverage tests wit known or unknown dam, was not available to assign paternity sufficed for traceability and paternity tests assuming null (or 1%) below 0.17×. On the other hand, results for paternity analyses by sequencing error rate and an accuracy for allele frequencies equal or Whalen et al. (2019) required greater coverage (0.4×) and larger higher to the ones obtained when sampling 10 individuals. This must be amount of genetic markers (50,000) to reach 100% sensibility. viewed as an outstanding result from technological, economic and sci­ In order to evaluate those procedures under more realistic scenarios, entific points of view. Moreover, the sequenceing data generated could different homogeneous error rates and accuracies for genotype fre­ have firther uses contributing more to sustainable science. The huge quencies were evaluated. In this case, the sequencing error rate had a amount of information available (even under very-low coverage) can be 6 J. Casellas et al. Livestock Science 251 (2021) 104629 exploited more in depth. Especially, with the structure of livestock Blume, J.D., 2002. Tutorial in biostatistics: likelihood methods for measuring statistical evidence. Stat. Med. 21, 2563–2599. species with dense family structures, large amounts of genomic data can Brouard, J.-.S., Boyle, B., Ibeagha-Awemu, E.M., Bissonnette, N., 2017. Low-depth accumulate across generations and years. This latter will open a new genotyping-by-sequencing (GBS) in a bovine population: strategies to maximize the window of animal breeding purposes, as the availability of whole selection of high quality gentypes and the accuracy of imputation. BMC Genet. 18, sequence for animal population may change the current animal breeding 32. Cockerham, C.C., 1969. Variance of gene frequencies. Evolution 23, 72–84. paradigm or even make a new revolution. Indeed, the exploration of Daetwyler, H.D., Capitan, A., Pausch, H., Stothard, P., van Binsbergen, R., Brøndum, R.F., sequence data at massive volume may allow to make animal breeding Liao, X., Djari, A., Rodriguez, S.C., Grohs, C., Esquerré, D., Bouchez, O., selection decisions more accurate by taking benefit of massive genomic Rossignol, M.-.N., Klopp, C., Rocha, D., Fritz, S., Eggen, A., Bowman, P.J., Coote, D., Chamberlain, A.J., Anderson, C., VanTassell, C.P., Hulsegge, I., Goddard, M.E., data (Knap, 2020). Thus, additional efforts to handle this new source of Guldbrandtsen, B., Lund, M.S., Veerkamp, R.F., Boichard, D., Fries, R., Hayes, B.J., partial genomic data may be of special relevance for the livestock in­ 2014. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and dustry (Knap, 2020). Evenmore, an additional investment to increase complex traits in cattle. Nat. Genet. 46, 858–867. Dalvit, C., De Marchi, M., Cassandro, M., 2007. Genetic traceability of livestock products: the sequencing coverage until 2×, which is still considered a review. Meat Sci. 77, 437–449. low-coverage, could allow to enhance animal breeding. Between the Edwards, A.W.F., 1972. Likelihood. Cambridge University Press, Cambridge, United possible options are the estimation of biological relatedness (Lipatov Kingdom. Fox, E.J., Reid-Bayliss, K.S., Emond, M.J., Loeb, L.A., 2014. Accuracy of next generation et al., 2015) and the imputation of the whole genome with high accuracy sequencing platforms. Next Gen. Seq. Appl. 1, 1000106. depending on the population size (Ros-Freixedes et al., 2020a, 2020b). Geldermann, H., Piper, U., Weber, W.E., 1986. Effect of misidentification on the This last step would be essential to implement whole-genome sequence estimation of breeding value and heritability in cattle. J. Anim. Sci. 63, 1759–1768. Gilly, A., Southam, L., Suveges, D., Kuchenbaecker, K., Moore, R., Melloni, G.E.M., data for genomic prediction and fine-mapping of causal variants. Hatzikotoulas, K., Farmaki, A.-.E., Ritchie, G., Schwartzentruber, J., Danecek, P., Kilian, B., Pollard, M.O., Ge, X., Tsafantakis, E., Dedoussis, G., Zeggini, E., 2019. 5. Conclusions Very low-depth whole-genome sequencing in complex trait association studies. Bioinformatics 35, 2555–2561. Goffaux, F., China, B., Dams, L., Clinquart, A., Daube, G., 2005. Development of a genetic Very low genome coverages in livestock species were enough to traceability test in pig based on single nucleotide polymorphism detection. Forensic guarantee ≥99% of true positives and true negatives for traceability Sci. Int. 1651, 239–247. Gorjanc, G., Cleveland, M.A., Houston, R.D., Hickey, J.M., 2015. Potential of genotyping- testing (from 0.01× coverage) and parentage testing (from 0.05× by-sequencing for genomic selection in livestock populations. Genet. Sel. Evol. 47, coverage). Even when 0.05× coverage sufficed for both tests, as genome 12. coverage increased, the percentage of reads per polymorphic SNPs and Gorjanc, G., Dumasy, J.-.F., Gonen, S., Gaynor, R.C., Antolin, R., Hickey, J., 2017. Potential of low-coverage genotyping-by-sequencing and imputation for cost- the certaintly of the estimate of its allele frequency increased, thus, effective genomic selection in biparental segregation populations. Crop Sci. 57, reducing the errors in the tests. Moreover, the length of the reads 1404–1420. affected the dispersion and number of shared SNPs among pairs of Hacking, I., 1965. Logic of Statistical Inference. Cambridge University Press, Cambridge, United Kingdom. sequenced individuals. Hardy, G.H., 1908. Mendelian proportions in a mixed population. Science 28, 49–50. Heaton, M.P., Harhay, G.P., Bennett, G.L., Stone, R.T., Grosse, W.M., Casas, E., Keele, J. Declaration of Competing Interest W., Smith, T.P., Chitko-McKown, C.G., Laegreid, W.W., 2002. Selection and use of SNP markers for animal identification and paternity analysis in US beef cattle. Mamm. Genome 13, 272–281. Authors declare no conflict of interest. Heaton, M.P., Leymaster, K.A., Kalbfleisch, T.S., Kijas, J.W., Clarke, S.M., McEwan, J., Maddox, J.F., Basnayake, V., Petrik, D.T., Simpson, B., Smith, T.P.L., Chitko- McKown, C.G., International Sheep Genomics Consortium, 2014. SNPs for parentage CRediT authorship contribution statement testing and traceability in globally diverse breeds of sheep. PLoS One 9, e94851. Hickey, J.M., Gorjanc, G., 2012. Simulated data for genomic selection and genome-wide Joaquim Casellas: Conceptualization, Data curation, Formal anal­ association studies using a combination of coalescent gene drop methods. G3 (Bethesda) 2, 425–427. ysis, Software, Writing – original draft. Melani Martín de Hijas-Vil­ Jamieson, A., 1965. The genetics of transferrins in cattle. Heredity (Edinb) 20, 419–441. lalba: Methodology, Writing – review & editing. Marta Vázquez- Kanthaswamy, S., 2015. Domestic animal forensic genetics – Biological evidence, genetic Gómez: Methodology, Writing – review & editing. Samir Id-Lahou­ markers, analytical approaches and challenges. Anim. Genet. 46, 473–484. Klein, R.J., Zeiss, C., Chew, E.Y., Tsai, J.Y., Sackler, R.S., Haynes, C., Henning, A.K., cine: Methodology, Writing – review & editing. SanGiovanni, J.P., Mane, S.M., Mayne, S.T., Bracken, M.B., Ferris, F.L., Ott, J., Barnstable, C., Hoh, J., 2005. Complement factor H polymorphism in age-related macular degeneration. Science 308, 385–389. Acknowledgments Knap, P.W., 2020. The scientific development that we need in the animal breeding industry. J. Anim. Breed. Genet. 137, 343–344. This research was partially supported by the grant CGL2016–80155- Kosambi, D., 1944. The estimation of map distances from recombination values. Ann. Eugen. 12, 172–175. R (Spain’s Ministerio de Economía y Competitividad, Madrid, Spain) and a Li, R., Yu, C., Li, Y., Lam, T.-.W., Yiu, S.-.M., Kristiansen, K., Wang, J., 2009. SOAP2: an fellowship granted to M. Martín de Hijas-Villalba (grant BES- improved ultrafast tool for short read alignment. Bioinformatics 25, 1966–1967. 2017–080596; Spain’s Ministerio de Economía y Competitividad, Madrid, Li, H., Ruan, J., Durbin, R.M., 2008. Mapping short DNA sequencing reads and calling Spain). variants using mapping quality scores. Genome Res. 18, 1851–1858. Lipatov, M., Sanjeev, K., Patro, R., Veeramah, K.R., 2015. Maximum likelihood estimation of biological relatedness from low coverage sequencing data. bioRxiv, References 023374. Liu, X., Han, S., Wang, Z., Gelernter, J., Yang, B.-.Z., 2013. Variant callers for next- generation sequencing data: a comparison study. PLoS One 8, e75619. Arana, A., Soret, B., Lasa, I., Alfonso, L., 2002. Meat traceability using DNA markers: Marshall, T.C., Slate, J., Kruuk, L.E.B., Pemberton, J.M., 1998. Statistical confidence for application to the beef industry. Meat Sci. 61, 367–373. likelihood-based paternity inference in natural populations. Mol. Ecol. 7, 639–655. Banos, G., Wiggans, G.R., Powell, R.L., 2001. Impact of paternity errors in cow Martin, E.R., Kinnamon, D.D., Schmidt, M.A., Powell, E.H., Zuchner, S., Morris, R.W., identification on genetic evaluations and international comparisons. J. Dairy Sci. 84, 2010. SeqEM: an adaptive genotyping-calling approach for next-generation 2523–2529. sequencing studies. Bioinformatics 26, 2803–2810. Bansal, V., Boucher, C., 2019. Sequencing technologies and analyses: where have we McKean, J.D., 2001. The importance of traceability for public health and consumer been and where are we going? Science 18, 37–41. protection. Rev. Sci. Tech. 20, 363–371. Besser, J., Carleton, H.A., Gerner-Smidt, P., Lindsey, R.L., Trees, E., 2018. Next- Meuwissen, T.H.E., Hayes, B.J., Goddard, M.E., 2001. Prediction of total genetic value generation sequencing technologies and their application to the study and control of using genome-wide dense marker maps. Genetics 157, 1819–1829. bacterial infections. Clin. Microbiol. Infect. 24, 335–341. Morton, N., 1998. Significance levels in complex inheritance. Am. J. Hum. Genet. 62, Bickel, D.R., 2011. A predictive approach to measuring the strength of statistical 690–697. evidence for single and multiple comparisons. Can. J. Stat. 39, 610–631. Nielsen, R., Paul, J.S., Albrechtsen, A., Song, Y.S., 2011. Genotype and SNP from next- Bickel, D.R., 2012. The strength of statistical evidence for composite hypotheses: generation sequencing data. Nat. Rev. Genet. 12, 443–451. inference to the best explanation. Stat. Sin. 22, 1147–1198. Pérez-Enciso, M., Rincón, J.C., Legarra, A., 2015. Sequence- vs. chip-assisted genomic Birnbaum, A., 1962. On the foundations of statistical inference. J. Am. Stat. Assoc. 57, selection: accurate biological information is advised. Genet. Sel. Evol. 47, 43. 269–306. 7 J. Casellas et al. Livestock Science 251 (2021) 104629 Pfeiffer, F., Gröber, C., Blank, M., Händler, K., Beyer, M., Schultze, J.L., Mayer, G., 2018. Strug, L.J., 2018. The evidential statistical paradigm in genetics. Genet. Epidemiol. 42, Systematic evaluation of error rates and causes in short samples in next-generation 590–607. sequencing. Sci. Rep. 8, 10950. Strug, L.J., Hodge, S.E., Chiang, T., Pal, D.K., Corey, P.N., Rohde, C., 2010. A pure Qian, J., Ruíz-García, L., Fan, B., Robla-Villalba, J.I., McCarthy, U., Zhang, B., Yu, Q., likelihood approach to the analysis of genetic association data: an alternative to Wu, W., 2020. Food traceability system for governmental, corporate, and consumer Bayesian and frequentist analysis. Europ. J. Hum. Genet. 18, 933–941. perspectives in the European Union and China: a comparative review. Trends Food Vázquez, J.F., Pérez, T., Ureña, F., Gudín, E., Albornoz, J., Domínguez, A., 2004. Sci. Technol. 99, 402–412. Practical application of DNA fingerprinting to trace beef. J. Food Prot. 67, 972–979. Ros-Freixedes, R., Whalen, A., Chen, C.-.Y., Gorjanc, G., Herring, W.O., Mileham, A.J., Vens, M., Schillert, A., König, I.R., Ziegler, A., 2009. Look who is calling: a comparison of Hickey, J.M., 2020a. Accuracy of whole-genome sequence imputation using hybrid genotype calling algorithms. BMC Proceed. 3, S59. peeling in large pedigreed livestock populations. Genet. Sel. Evol. 52, 17. Visscher, P.M., Woolliams, J.A., Smith, D., Williams, J.L., 2002. Estimation of pedigree Ros-Freixedes, R., Whalen, A., Gorjanc, G., Mileham, A.J., Hickey, J.M., 2020b. errors in the UK dairy population using microsatellite markers and the impact on Evaluation of sequencing strategies for whole-genome imputation with hybrid selection. J. Dairy Sci. 85, 2368–2375. peeling. Genet. Sel. Evol. 52, 18. Whalen, A., Gorjanc, G., Hickey, J.M., 2019. Parentage assignment with genotyping-by- Royall, R.M., 1997. Statistical Evidence: A Likelihood Paradigm. Chapman and Hall, sequencing data. J. Anim. Breed. Genet. 136, 102–112. London, United Kingdom. Weller, J.I., Feldmesser, E., Golik, M., Tager-Cohen, I., Domochovsky, R., Alus, O., Snyder-Mackler, N., Majoros, W.H., Yuan, M.L., Shaver, A.O., Gordon, J.B., Kopp, G.H., Ezra, E., Ron, M., 2004. Factors affecting incorrect paternity assignment in the Israeli Schlebusch, S.A., Wall, J.D., Alberts, S.C., Mukherjee, S., Zhou, X., Tung, J., 2016. Holstein population. J. Dairy Sci. 87, 2627–2640. Efficient genome-wide sequencing and low-coverage pedigree analysis from Yu, X., Sun, S., 2013. Comparing a few SNP calling algorithms using low-coverage noninvasively collected samples. Genetics 203, 699–714. sequencing data. BMC Bioinform. 14, 274. Stormont, C., 1967. Contribution of blood typing to dairy science progress. J. Dairy Sci. Zan, Y., Payen, T., Lillie, M., Honaker, C.F., Siegel, P.B., Carlborg, Ö., 2019. Genotyping 50, 253–260. by low-coverage whole-genome sequencing in intercross pedigrees from outbred founders: a cost-efficient approach. Genet. Sel. Evol. 51, 44. 8