It has previously been found that calling tetraploid allele dosages from sequence-based genotyping data is challenging from low-sequence coverage studies (Li et al., 2014), and it has been suggested that a coverage of 60–80× is appropriate to secure concordance between sequence-based genotyping calls and KASP analysis (Uitdewilligen et al., 2013). However, the recent finding that elite autotetraploid potato has an observed average ploidy of 3.2 across all loci (Sun et al., 2022) complicates this problem even further. An even higher resolution is demanded for distinguishing between all the possible heterozygous genotypes of both tetraploid and triploid loci, which can be expected to occur at high frequencies: AAAB, AABB, ABBB, AAB, and ABB. When plotting the observed frequency ratio of heterozygous SNPs across all MASPOT genotypes, peaks corresponding to 0.25, 0.33, 0.5, 0.67, and 0.75 are apparent ( Figure 1A ), supporting the finding that many loci deviate from a tetraploid state. Furthermore, peaks are observed at 0.125, 0.375, 0.625, and 0.875, which seem to suggest octoploid heterozygous loci, as well as peaks at 0.167 and 0.833 would suggest hexaploid heterozygous loci. We speculate that these may stem from the mapping of duplicated gene loci in the genotyped cultivars to a single locus in the reference. In such a case, the observed octoploid distribution is expected for the mapping of two tetraploid loci, and the observed hexaploid distribution for the mapping of two triploid loci. Following this notion, another two loci scenario (tetraploid at one locus and triploid at another) would give rise to peaks in the histogram at 0.143, 0.286, 0.428, 0.571, 0.714, and 0.857 and, indeed, peaks can be observed at these positions. A diploid locus and a triploid locus mapping to one locus would give rise to peaks at 0.2, 0.4, 0.6, and 0.8, which is also observed. Overall, this is in agreement with the finding by Sun et al. (2022) that allele copy numbers of 1–4 were observed across loci in the Otava genome. On the other hand, the read coverage of SNPs included in the analysis was in the range of 5–60, with a mean of 19.2 and a median of 14.6 ( Supplementary Figure S1 ). Therefore, we cannot exclude that the observed structure in the data stems from binomial sampling at low sampling depth which is also expected to give rise to structure in the data. However, the frequency ratios of the subset of markers with the highest coverage (between 20 and 60) also display the distinct 0.25, 0.33, 0.5, 0.67, and 0.75 peaks representing the cooccurrence of loci of varying ploidy in the potato genome ( Figure 1B ), supporting that the observed non-tetraploid ratios are not caused by sampling errors. Regardless, it is apparent from Figures 1A, B that, with the read coverage in our study, the main distributions, representing triploid and tetraploid loci, are overlapping and hence, there is insufficient resolution to reliably determine the per-loci ploidy state using the observed data directly for most loci. Therefore, we sought to simulate how the ability to distinguish between the “worst case scenario”—ABB and ABBB is developing with sequence coverage by simulating a random sampling of N depth from a distribution of 0.33 A/0.67 B and 0.25 A/0.75 B, respectively. In Figure 1C , the 95th and 5th percentiles of the simulated distribution of SNP frequency ratios are shown for tetraploid and triploid loci, respectively. It is obvious, that even in this ideal case in the absence of sequencing errors, extreme sequence depth in excess of ~320X is needed to distinguish between a heterozygous triploid and a heterozygous tetraploid state with 95% confidence. This underscores that it is not practically feasible to distinguish between tetraploid and triploid loci with sequence-based genotyping methods, at least with the current cost of sequencing technologies. Similarly, Grandke et al. (2016) found that heterozygous states could not successfully be called in hexaploid chrysanthemum, because allele signals did not segregate into distinguishable clusters—an assumption of the applied genotype caller SuperMASSA. The poor resolution of heterozygous signals in polyploids, an issue of escalating complexity with increased ploidy (de Bem Oliveira et al., 2019), indicates that the continuous genotypes are currently the best option for genotyping polyploids. This adds to the relevance of investigating whether calling discrete genotypes affects GP and GWAS performance in potatoes compared to continuous genotype values.

Single-trait GP of the five traits using the true marker dosage (continuous allele frequency SNPs) and forced tetraploid markers performed almost identically ( Figure 2 , Supplementary Figure S2 , Table 1 ), only chipping quality (with the highest amount of missing phenotype data [30%]) (p = 1.3*10⁻⁸) and senescence (p = 0.018) being modeled to improved mean correlation coefficients using the tetraploid compared to continuous SNPs. However, the numerical increase of the mean correlation is only a diminutive rise from 0.534 to 0.543 for chipping quality and from 0.343 to 0.347 for senescence. This could be a result of a very small trait-dependent response to the genotyping method, similar to previous reports on diploidization in potato (Gemenet et al., 2020). However, the difference in mean correlation for chipping quality could also be attributable to the effectively smaller training population for this trait, due to the large amounts of missing phenotypic data. Overall, based on these results genomic prediction is largely insensitive to calling genotypes as continuous or using a tetraploid model.

It has already been shown for loci with relatively low read coverage that it is not possible to distinguish between all three heterozygous states of an autotetraploid (Uitdewilligen et al., 2013; Li et al., 2014; Byrne et al., 2020), and our simulation shows that it is practically unattainable to confidently discriminate between heterozygous triploid and tetraploid loci with coverage below <320×. Irrespective of this, forcing tetraploid SNP calls in our case with loci of read depth greater than 5× and smaller than 60× did not affect model performance compared to using the observed allele frequencies—effectively validating previous genomic prediction model results based on tetraploid SNP calls in potato. Notably, while high coverage is needed to confidently call heterozygous tetraploid genotypes and discriminate ploidy, using read depths >60× to determine allele frequencies has been reported to result in changes in allele frequency variances (Ashraf et al., 2016), possibly influenced by signal from abundant repetitive elements, why using loci with extremely high read coverage is generally not recommended or employed. Results from autotetraploid blueberry further indicates that for accurate predictions using continuous genotypes, it is sufficient to use low-mid read depths of 6–12× compared to 60× (de Bem Oliveira et al., 2020).

We previously established that prioritizing continuous SNPs according to type (i.e., nsSNP, sSNP, ncSNP, or all combined) also did not affect GP model performance for most of the traits analyzed here (except yield) (Aalborg et al., 2024). Likely, the reasons for the invariability of the prediction performances observed in these two cases are similar and both related to the LD decay of the potato genome and the distance between the anonymous marker and the causal variant. With an estimated LD block size of 0.6–1.5 Mb and 2–5 Mb in introgressed genome regions (Vos et al., 2017), it can be extrapolated that, for an LD block size of 1 Mb, all markers in a 500 kb window on either side of a causal variant can be used as a reasonably reliable marker, that is, with a redundant SNP set of ~31k genome-wide markers there is likely to be found correctly called variants in LD with causal trait variants when using a tetraploid model. Hence, the linkage effect is likely sufficiently large to predominate over any erroneous signals of miscalled allele dosages. While we did not examine the effect of read depth on prediction accuracy in potato, it is plausible that a similar effect could support the use of low read depth data. As genotype quality is proportional to read depth, the effective pool of accurately sampled SNPs is lowered with decreased sequencing depth but is likely still sufficiently large to include reliable LD block markers. Previous results have shown that GP models are resilient to marker reductions in a trait-dependent manner in potato, but that across several traits only 1–10k markers of read depth between 5–60x were required to reach optimal model performance (Aalborg et al., 2024).

An elaborate description of the major QTLs identified is presented in (Aalborg et al., 2024). The GWAS analyses with tetraploid SNPs reproduced the major QTLs found with continuous genotypes ( Figure 3 ; Supplementary Files S1, S2 ). Differences constituted possible minor QTLs, major QTL significance, and therefore also the number of significant SNPs constituting the major QTL signals. The latter is observed for dry matter content, where continuous SNPs generate a major QTL signal on chromosome 10 from 48 to 58 Mb constituted of 13 significant markers, whereas the chromosome 10 signal is reduced to four SNPs in the 53–56 Mb range when using tetraploid SNPs. For major trait QTLs, the tetraploid SNPs hence hold sufficient true information about the trait genetic variation to accurately map those QTLs, albeit with slightly reduced significance, observed as variation in signal intensity in terms of both range and amplitude. However, for the identification of candidate minor QTL effects, a differential outcome is observed, underscoring the need for validation of such QTLs. The higher sensitivity of GWAS performance compared to GP is in agreement with what we previously found when analyzing the effect of marker type (Aalborg et al., 2024).

This indicates that, while the overall amount of additive genetic variance explained is similar across the two genotyping methods, as it is captured to the same level by the genomic prediction model, continuous genotypes contribute alternative information when resolved to chromosomal signal location. We speculate that this may vary across traits for two reasons. First, this could reflect uneven dispersion of the gene copy number variations, a likely result of progressive phase deletions during successive meiotic recombination events, effectively generating a locus-dependent response to the genotyping method. In regions where the local LD-decay is low due to hot spots for recombination events, the local concentration of true markers becomes site-dependent using the tetraploid model. Secondly, our inability to reliably call the missing phase for any locus and, hence, the potential haploinsufficiency effects of a triploid versus tetraploid site becomes cryptic. If there indeed is an effect, it is not specially specified as a potential feature in the statistical models and, hence, we have no marker for the specific difference of triploid/tetraploid and therefore it is not fully included in the analysis with either genotyping method.

The resilience of population/family structure-corrected GP models to marker reduction (Sverrisdóttir et al., 2018; Selga et al., 2021; Aalborg et al., 2024), their indifference to marker type filtration (Aalborg et al., 2024), and, as seen here, their general tolerance to allele dosage parametrization (Endelman et al., 2018; de Bem Oliveira et al., 2019; Matias et al., 2019; Gemenet et al., 2020; Yadav et al., 2024) all suggest that efforts such as SNP pruning, high-depth sequencing (de Bem Oliveira et al., 2020), and so forth can neither notably improve model performance in potato for complex, quantitative traits such as yield nor detection of single major QTL traits such as senescence and length/width ratio beyond a plateau level, which is reached with relatively shallow analyses. Rather, future efforts in increasing model performance should be directed toward increasing the size of phenotype/genotype datasets and capturing the non-additive genetic variation of traits, for example, dominance, epistasis [though this has generated mixed results in previous studies (Su et al., 2012; Endelman et al., 2018; Wilson et al., 2021)], pleiotropy, as well as gene–environment interactions. Such studies are further complicated by the fact that a full tetraploid model is not appropriate for modeling these effects across all loci. To capture the effects of haploinsufficiency of phase deletions these should also be specified in models for continuous markers capable of capturing ploidy variations.

In practice, as documented in this study, using a tetraploid model for calling SNPs in potato had largely no effect on genomic predictions compared to using observed allele frequencies and only resulted in very minor variations in GWAS results. This constitutes a relief for existing genomics-assisted breeding programs since existing GP models that are used need not be remade. Particularly, because determining gene allele dosage from sequence data directly is practically infeasible. Nonetheless, a general conservative statistical modeling approach is generally advocated where models are kept as simple as possible unless an explicit argument for more complicated models is presented and the differential effects documented. This dictates that tetraploid SNP calls should not be used, but rather the observed allele frequencies computed directly from sequencing data as continuous genotypes. As results from (de Bem Oliveira et al., 2020) indicates that low-mid sequencing depth is sufficient for accurate prediction with continuous genotypes, this does not necessarily entail inflated genotyping costs (Ferrão et al., 2021). In extension, pseudo-diploidization is also discouraged based on results from (Endelman et al., 2018), showing deflated prediction accuracies with this parametrization. Using continuous genotypes also entails a reduction in computational time, as the genotype calling step is circumvented by calculating allele frequencies directly from sequencing data (de Bem Oliveira et al., 2019).

Genotyping of the MASPOT panel was performed by genotyping-by-sequencing (GBS) as described in (Sverrisdóttir et al., 2017). GBS libraries were prepared following a protocol adapted from (Elshire et al., 2011). Illumina 5′ and 3′ adaptors for sequencing were designed for a 96-multiplexing system. Leaf tissue DNA was extracted and digested with ApeKI, the fragments were ligated to adaptors and pooled in 96-plex libraries, purified, and amplified by PCR. The libraries were then sequenced on a HiSeq 2000 (Illumina, San Diego, USA) with single-read sequencing (100 bp). Each 96-plex library was sequenced on three channels of a flow cell.

Custom bash scripts ( Supplementary File S3 ) were used to simulate a random sampling of N depth at 2,500 loci from parental distributions of ABB or ABBB, respectively. N was varied from 3 to1,000. The SNP ratio average, the 95th and 5th percentiles was calculated from either case and plotted using Microsoft Excel. The crossing of the 5th percentile of the ABBB case with the 95th percentile of the ABB case was taken as 95% confidence threshold.

The sequenced reads were processed as described in (Sverrisdóttir et al., 2017). The reads were demultiplexed, trimmed, and mapped onto the reference genome of double monoploid S. tuberosum Group Phureja DMv4.03 (Sharma et al., 2013). Tetraploid SNPs were called using the Genome Analysis Toolkits (McKenna et al., 2010) UnifiedGenotyper tool with ploidy set to 4, a minimum phred-scaled confidence threshold of 50 for variant calling and of 20 for variant omission (and filtered with LowQual for below the calling threshold), cf (Sverrisdóttir et al., 2017). SNPs were filtered to a minimum root mean squared quality of 30 and reduced to only biallelic variants. In addition to the tetraploid variant call, observed allele frequency estimated directly from sequencing data was also used as genotyping method, cf. (Ashraf et al., 2016) to accommodate the gene copy number variation across loci in the potato genome (Sun et al., 2022). The allele frequencies were determined as the ratio between allele counts of the alternative allele and the total allele count, giving a continuous allele ratio between 0 and 1, but not constricted to tetraploid dosages ( Equation 1 ).

This provided two sets of genotypes for the same SNPs, that is, continuous genotypes and tetraploid genotypes. The five possible tetraploid genotypes [AAAA, AAAB, AABB, ABBB, or BBBB] were recoded to [0, 1, 2, 3, 4] allele dosages. Minor alleles frequency (MAF) was calculated from read coverage, and SNPs were filtered to an MAF of 1% (i.e., mean allele frequency<0.99 and >0.01), missing rate of maximum 50%, and finally filtered to read coverage >5 at all positions, introducing a missing genotype for low-quality positions in individual samples. Based on the continuous genotypes, the SNPs were further filtered to SNPs with read coverage between 5 and 60 and samples were filtered to clones with <70% missing data. A final SNP set of ~31k positions of the filtered SNPs spanning the entire genome was used for further analysis. The SNP set originates from (Aalborg et al., 2024) and consists of 31,032 SNPs that are approximately evenly composed of non-synonymous SNPs, synonymous SNPs, and non-coding SNPs ( Supplementary Figure S3 , Supplementary File S4 ). Both the continuous and tetraploid genotypes were reduced to that SNP set. SNPs that were monomorphic in either genotype matrix were filtered from both sets, yielding two final genotype matrices of 31,007 markers. Finally, the tetraploid genotypes were recoded from [0, 1, 2, 3, 4] to [0, 0.25, 0.50, 0.75, 1] to match the continuous genotype matrix format.