Testing pipelines for genome-wide SNP calling from Genotyping-By-Sequencing (GBS) data for Pinus ponderosa

Abstract

Background Single Nucleotide Polymorphism (SNP) markers have rapidly gained popularity due to their abundance in most genomes and their amenability to high-throughput genotyping techniques. Reduced-representation restriction-enzyme-based sequencing methods (GBS or RADseq) have been demonstrated to be robust and cost-effective genotyping methods. While previous studies have shown that alignment of the short-read fragments to a genome sequence results in better SNP calling than de novo approaches, only a few tree species - and few conifers in particular - have an annotated sequence. While these could be used to align sequence fragments from related species, sequence divergence might result in SNPs being missed if they are in fragments that don't align properly. Producing a new annotated genome sequence for every conifer species before SNP analyses are conducted is still prohibitive, as many conifer genomes are huge (> 19 GB) and include a large proportion of repeat sequences, making assembly difficult. Here we compare four bioinformatics pipelines, two of which require a reference genome (TASSEL-GBS V2 and Stacks), two of which are de novo pipelines (UNEAK and Stacks). We used Illumina sequence data from 94 ponderosa pines, with loblolly pine as the reference genome. Results The number of SNPs called was much lower without a reference genome (62–196 thousand vs. 2.1–2.7 million SNPs). UNEAK was the fastest overall and identified more SNPs than Stacks de novo. Stacks with a reference genome produced the highest number of SNPs with lowest proportion of paralogs, while SNPs identified by TASSEL-GBS V2 exhibited the highest heterozygosity, minor allele frequency, and proportion of paralogs. More SNPs were uniquely identified by Stacks than TASSEL, though there was high overlap between methods. Conclusion The present case study provides a comprehensive comparison between four commonly-used SNP calling pipelines, and identifies the Stacks reference-based approach as the best overall for conifers (or other species with large repetitive genomes) that do not have a published reference genome for the same species. However, all four pipelines had distinct benefits and limitations, with Stacks for instance being less user-friendly than some of the other pipelines. In addition, researchers studying other conifer species using similar approaches should be prepared to analyze very large numbers of SNPs.

Haploid, diploid, and pooled exome capture recapitulate features of biology and paralogy in two non‐model tree species

“…There are several methods by which to detect such problematic sites, such as filtering by coverage (Dou et al, 2012 ), disomic models such as Hardy‐Weinberg proportions (Catchen et al, 2013 ; Chen et al, 2014 ; Hohenlohe et al, 2011 ), or gene annotation, though there are several shortcomings (see descriptions of these shortcomings in Table 1 of McKinney et al, 2017 ). When individual sequencing data is available for the same individuals or populations, such information can be used to isolate potentially paralogous sites from pool‐seq exome capture studies (e.g., Rellstab et al, 2019 ; Shu & Moran, 2020 ). However, a potentially cost‐saving alternative would be to sequence the haploid tissue of a single individual (if available).…”

“…This is particularly true for pool‐seq data sets relying on read counts for allele frequency estimation or population genetic inferences such as genotype‐environment associations. While individually sequenced data sets may be one path forward to identifying such problematic sites (as in Rellstab et al, 2019 ; Shu & Moran, 2020 ), the sequencing of sufficient quantities of DNA from haploid gametophyte tissue available for some plants, including conifers, seedless vascular plants, and bryophytes, offers an alternate path forward to balance sequencing cost and data reliability, particularly for organism using diverged and or highly fragmented reference genomes.…”

Haploid, diploid, and pooled exome capture recapitulate features of biology and paralogy in two non-model tree species

“…This is particularly true for pool-seq datasets relying on read counts for allele frequency estimation or population genetic inferences such as genotype-environment association (e.g., as implemented in baypass, Gautier et al 2015). While individually sequenced datasets may be one path forward to identifying such problematic sites (as in Rellstab et al 2019, Shu & Moran 2020, it will be cost-prohibitive for large projects with many populations. As such, the sequencing of sufficient quantities of DNA from haploid gametophyte tissue available for some plants, including conifers,…”

Evaluating the accuracy of variant calling methods using the frequency of parent‐offspring genotype mismatch

“…For instance, conifers often have exceptionally large genomes (20–40 Gbp; Neale et al, 2017) with histories of whole‐genome duplication (Zheng et al, 2015), gene family expansion (Scott et al, 2020; De La Torre et al, 2014), transposable element dynamics (Scott et al, 2020; Wang et al, 2020; Yi et al, 2018), and extensive repeat regions (Wegrzyn et al, 2014). These complexities present a major challenge for NGS data analysis and downstream hypothesis testing in conifers (Lind et al, 2022; Shu & Moran, 2020). Such challenges can be alleviated by quantifying the accuracy of SNP calling when using the above‐mentioned variant calling tools.…”