An evaluation of sequencing coverage and genotyping strategies to assess neutral and adaptive diversity

Abstract

Whole genome sequences (WGS) greatly increase our ability to precisely infer population genetic parameters, demographic processes, and selection signatures. However, WGS may still be not affordable for a representative number of individuals/populations. In this context, our goal was to assess the efficiency of several SNP genotyping strategies by testing their ability to accurately estimate parameters describing neutral diversity and to detect signatures of selection. We analysed 110 WGS at 12× coverage for four different species, i.e., sheep, goats and their wild counterparts. From these data we generated 946 data sets corresponding to random panels of 1K to 5M variants, commercial SNP chips and exome capture, for sample sizes of five to 48 individuals. We also extracted low-coverage genome resequencing of 1×, 2× and 5× by randomly subsampling reads from the 12× resequencing data. Globally, 5K to 10K random variants were enough for an accurate estimation of genome diversity. Conversely, commercial panels and exome capture displayed strong ascertainment biases. Besides the characterization of neutral diversity, the detection of the signature of selection and the accurate estimation of linkage disequilibrium (LD) required high-density panels of at least 1M variants. Finally, genotype likelihoods increased the quality of variant calling from low coverage resequencing but proportions of incorrect genotypes remained substantial, especially for heterozygote sites. Whole genome resequencing coverage of at least 5× appeared to be necessary for accurate assessment of genomic variations. These results have implications for studies seeking to deploy low-density SNP collections or genome scans across genetically diverse populations/species showing similar genetic characteristics and patterns of LD decay for a wide variety of purposes.

Keywords: SNP chip; depth of coverage; genotyping-by-sequencing; mammals; population genomics; whole genome sequencing.

Figure 1. Flow-chart describing sampling random and non-random panels of variants and individuals.

Whole genome sequences are denoted by WGS.

Figure 2. Nucleotide diversity (π) in sheep calculated from WGS data and from random and non-random panels of variants.

Nucleotide diversity (π) was estimated for each replicate of the different numbers of variants of the random panels and for each non-random panel. Sample sizes varied for each estimate from 10 to 30 individuals. Random panels are denoted by their number of variants (from 1K to 5M) and non-random panels by: 50K.Chip (Illumina® ovine 50K SNP Beadchip), HD.Chip (Illumina® ovine HD Beadchip) exome (exome capture simulation), WGS (all variants extracted from whole genome sequences). For each panel of variants the sample sizes are from left to right: 10 (red), 20 (green) and 30 (yellow) individuals.

Figure 3. Site frequency spectra (SFS) in sheep inferred from WGS data and from random and non-random panels of variants.

Site frequency spectra were estimated using different random and non-random panels for 30 sheep. Random panels are denoted by their number of variants (from 1K to 5M) and non-random panels by: 50K.Beadchip (Illumina® ovine 50K SNP Beadchip), HD.Beadchip (Illumina® ovine HD Beadchip) exome (exome capture simulation), WGS (all variants extracted from whole genome sequences). Pearson correlation coefficients with the WGS inferences are shown for each panel.

Figure 4. Nucleotide diversity (π) estimated in two Ovis groups with random and commonly used panels of variants.

Plot of Nucleotide diversity (π) estimated with a random set of 10K variants sampled in sheep data (10K), and with Illumina® ovine 50K SNP Beadchip (50K.Chip), Illumina® ovine HD Beadchip (HD.Chip), and variants extracted from whole genome sequences (WGS).

Figure 5. Estimates of individual inbreeding coefficient (F) and observed heterozygosity (Ho) from different panels of variants compared to WGS data estimates in sheep.

Plot of individual inbreeding coefficient (F; top) and observed Heterozygosity (Ho; bottom) estimated with variants extracted from whole genome sequences (WGS) versus inferences with Illumina® ovine 50K SNP Beadchip (50K.Chip), Illumina® ovine HD Beadchip (HD.Chip), and 1 set of 10K variants defined in Moroccan sheep (random 10K). The red lines represent the relationship for which the estimates of the different panels are identical to the ones of WGS inferences.

Figure 6. Fixation index (Fst) between Moroccan goats and Bezoar ibex for different panels of variants and different samples of individuals.

The fixation index Fst (Weir & Cockerham, 1984) was estimated for each random panel for the 5 independent replicates, and for each non-random dataset for each sample size. Random panels are denoted by their number of variants (from 1K to 5M) and non-random panels by: 50K.Chip (Illumina® caprine 50K SNP Beadchip), WGS (all variants extracted from whole genome sequences). For each panel of variants the sample sizes are from left to right: 18 (red), 33 (green) and 48 (yellow) individuals.

Figure 7. XP-CLR scores calculated along the 20M-40M bp segment on chromosome 10 in a horned-polled Moroccan sheep comparison for different sets of variants.

The two peaks of XP-CLR scores showed in the WGS data plot are located respectively in the two genes NBEA (chr 10: 26,007,917 - 26,592,574) and MAB21L1 (chr 10: 26,231,353 - 26,232,432) and in the RXFP2 gene (chr 10: 29,454,677 - 29,502,617 bp). The horizontal dashed line represents a XP-CLR score of 15 to represent a scale among the different plots.

Figure 8. Efficiency and accuracy of different genotyping strategies

For each purpose, the different strategies are rated according to the accuracy of the estimates taking as a reference the WGS 12x depth inferences. Grey dots indicate that the genotyping approach allow detecting some selection signatures but could miss some further signals detected by high density panels and WGS (12x depth) data. MD chip = 50K SNP BeadChip (caprine and ovine); HD chip = 600K SNP Ovine BeadChip. Low and medium re-sequencing coverages are represented by: (i) classical variant calling and filtering denoted by 1x, 2x and 5x and (ii) variant discovery based on genotype likelihoods denoted by 1xGL, 2xGL and 5xGL.