2.1 SNP Discovery, quality assessment and selection
After processing raw sequencing data from a DArTseq GBS library of 650 individuals (ten individuals were excluded during library preparation), a dataset of 24,683 RAD-tags containing 31,262 SNP markers (31K SNP dataset), was returned. An average of 1.37 SNPs ± 0.6 SD were present in each RAD-tag, with a maximum of six SNPs observed in a single RAD-tag (19,288 tags with one SNP, 4,300 tags with two SNPs, 915 tags with three SNPs, 155 tags with four SNP, 22 tags with five SNP, three tags with six SNPs). For this raw 31K SNP dataset, the average genotype call rate was 0.86 ± 0.14 SD and the average MAF was 0.11 ± 0.15 SD.
An average of ~ 2.5 million reads were allocated to each individual in library preparation and sequencing; however, after removal of low quality sequences, monomorphic loci, and SilicoDArT markers (presence/absence variants), an average of 459,987 ± 88,493 SD reads were associated with each individual (366,609,683 total reads over 650 individuals). The average read depth over all non-missing genotype calls was 17.0 ± 18.3 SD, however, 27.0% of non-missing genotype calls had 5 or less reads associated.
To ensure only the highest quality markers were available for bates probe selection and assay development, a series of SNP filtering thresholds were implemented to remove individual genotype calls with low confidence, remove less informative markers and remove erroneous data (Fig. 1). Filtering for minimum read depth (< 5) removed 20.4% of the non-missing genotype calls (4,155,990), increasing the overall missing-ness of the 31K dataset from 14–33%. Subsequently filtering for minimum minor allele frequency of ≥ 0.02 removed 16,391 SNPs. A further 4,646 SNPs were removed due to having a minimum call rate equal to or below 0.5; no SNPs were discarded due to having a repeatability score of less than 0.9 due to pre-filtering on repeatability before data was provided. For SNPs that were derived from identical clones (100% rad-tag sequence similarity), the SNP with the highest MAF was retained from each clone resulting in the removal of 1,500 SNPs. Similarly, redundancy clustering of clone sequences at 95% sequence similarity removed a further 572 SNPs derived from highly similar clones. Of the markers retained, tests of conformity to Mendelian inheritance patterns and HWE were conducted. A total of 364 SNPs were removed due to MI errors (> 9%) and 247 SNPs were removed due to HWE deviations. Finally, 3,101 genotype calls identified as MI errors across the remaining SNPs were silenced (Fig. 1).
In total 23,720 SNPs were removed from the 31K dataset, leaving 7,542 high quality SNPs for further ranking and final SNP selection. After filtering, the average MAF increased from 0.11 ± 0.15 to 0.22 ± 0.15. While the average per genotype call rate decreased from 0.86 ± 0.14 to 0.84 ± 0.14, the average read depth per genotype call increased from 17.0 ± 18.2 SD to 30.7 ± 27.9 SD. Furthermore, the rate of average Mendelian inheritance errors detected across SNPs decreased from 3.2–1.4% after filtering.
2.1.1 SNP number and GRM analysis
To determine the number of markers required for accurate genomic relationship calculations and selection of the density of the DArTcap panel, a number of marker subsets were modelled. An increasingly consistent GRM estimate was achieved by including progressively more random markers from 100 to 1,000 markers (Fig. 2). With 1,000 markers a correlation of 0.95 was consistently achieved when compared to the full 7,452 marker set (Fig. 2). Increasing the number of markers from 1,000 to 4,000 further increased the correlation between marker sets, with the correlation exceeding the desired minimum cut off for future assay applications of 0.98 at 4,000 markers (Fig. 2).
Since GRM correlations indicated an optimal assay size of 4,000 SNPs, a QC score was assigned to the set of 7,452 high quality DArTseq markers allowing the prioritisation of 4,236 SNPs for DArTcap probe systhesis. For this selected subset of markers, the average call rate was 80.2% ± 14.7%, MAF was 0.35 ± 0.28 and the average read depth per non-missing genotype call was 20.4 ± 15.6. Furthermore, genomic relationship values calculated with the 4,236 markers showed high correlation (r2 = 0.987) to the full 7,452 marker panel (Fig. 3a). A number of pairwise relationships were estimated to be negative. These values are a result of the distinct underlying population structure between individuals from East Australina Coast populations and those from Northern Territory populations [29].
2.1.2 Linkage disequilibrium
Limited evidence of substantial linkage disequilibrium was observed with the DArTseq datasets. A total of 4 marker pairs were found to be linked (r2 > 0.80) in the 7,452 SNP dataset when considering all available individual genotypes, while only two of these marker pairs persisted within the selected 4,236 SNPs (r2 values of 0.96 and 0.85). When assessing genotypes of the Northern Territory individuals only (Tiwi Island, Joseph Bonaparte Gulf, Gulf of Carpentaria) from the 7,452 SNP dataset, there were three marker pairs with r2 values greater than 0.8. Likewise, for East Australian Coast individuals only (Bramston Beach, Townsville, Etty Bay), there were 13 r2 values greater than 0.80, however, only one of these was also observed in the Northern Territory population. Using only Northern Territory individuals from the 4,236 SNP dataset, there was one marker pair identified to be in LD (r2 > 0.80). Likewise, when considering genotypes of East Australian Coast individuals, there were 12 r2 values greater than 0.80, with one pair shared with those observed in the Northern Territory population.
2.2 DArTcap assay validation
A number of samples resubmitted for genotyping did not pass sample digestion QC, and were excluded from sequencing (Table 1). This was most prominently observed in the initial set of 2nd generation farm stock (44 of 90 submitted), and was identified to be due to DNA degradation during storage. Across the remaining 462 individuals, 21 samples from across various populations failed QC and were excluded (Table 1). Processing of raw sequencing data obtained from a DArTcap target-capture library including a total of 485 individuals produced a dataset of 15,880 RAD-tags containing a total of 26,751 SNP markers (raw DARTcap dataset). An average of 1.7 SNPs ± 1.0 SD were present in each RAD-tag, with a maximum of 18 SNPs observed in a single RAD-tag (9,110 tags with 1 SNP, 4,142 tags with 2 SNPs, 1,731 tags with 3 SNPs, and 897 with > 3 SNPs). For the raw DArTcap dataset, the average individual call rate was 72.0% ± 6.8% SD, average genotype call rate was 72.0% ± 25.17% SD, and average MAF was 0.11 ± 0.13 SD. Approximately 600,000–700,000 reads were allocated to each individual in library preparation and sequencing, and after removal of low quality sequences, an average of 670,610 ± 185,993 reads were associated with each individual (328,598,913 total). The average read depth per genotype call was 34.8 ± 61.8 SD, however 32.0% of non-missing genotype calls had 5 or less reads associated.
Within the raw DArTcap dataset, 4,194 (99.0%) of the selected 4,236 panel were observed, with only 52 bates probes failing to capture the desired rad-tag sequences. For the 4,194 DArTcap dataset, the average individual call rate was 93.0% ± 3.5% SD, and average MAF was 0.23 ± 0.15 SD. An average of 312,343 ± 83,235 reads were associated with each individual (153,048,102 total). An average of 46.8% ± 1.6% of the total DArTcap sequencing effort was successful in obtaining target sequences per individual. The average read depth per non-missing genotype call was 80.0 ± 88.7 SD, while 5.78% of non-missing genotype calls had 5 or less reads associated.
Table 1
Individuals genotyped with DArTseq and DArTcap
Population | Region | DArTSeq | DArTCap |
# submitted | # passing library preparation | # submitted | # passing library preparation |
Townsville | East Coast, Australia | 22 | 22 | 10 | 10 |
Etty Bay | East Coast, Australia | 50 | 50 | 15 | 7 |
Bramston Beach | East Coast, Australia | 60 | 60 | 12 | 9 |
Gulf of Carpentaria | Northern Territory, Australia | 42 | 35 | 14 | 14 |
Tiwi Island | Northern Territory, Australia | 56 | 56 | 10 | 10 |
Joseph Bonaparte Gulf | Northern Territory, Australia | 34 | 34 | 13 | 13 |
Nickol Bay | Western Australia, Australia | - | - | 19 | 18 |
1st Generation | Farm Stock | 165 | 162 | 87 | 86 |
2nd Generation Set 1 | Farm Stock | 231 | 231 | 90 | 46 |
2nd Generation Set 2 | Farm Stock | - | - | 282 | 272 |
Total | 660 | 650 | 552 | 485 |
2.3 DArTcap assay utilisation
2.3.1 Parentage assignment
The results of parentage assignment success using Colony are presented in Table 2. At a conservative genotyping error rate of 10%, parentage assignment across the three sets of markers (7,452 DArTseq, 4,236 DArTseq, 4,194 DArTcap) were essentially identical, where only one individual (from 72 known parent-progeny relationships) had a single parent incorrectly left unassigned (false exclusion) when using the 4,194 DArTcap marker set (1.39% assignment error). The occurrence of false exclusions increased to 4.2% and 8.3% for the DArTcap marker set when using genotyping error rates of 5% and 1% respectively. The rate of false exclusions (2.3%) remained consistent irrespective of the error rate used in Colony for both the 7,452 and 4,236 DArTseq marker sets; however, the individual false exclusions did vary between either exclusion of the mother or father for some progeny.
Across all analyses, progeny containing higher missing data (> 10%) accounted for the majority of false parent exclusions. Furthermore, irrespective of the markers used in the analysis, there was a single individual where at least one of the two true parents could not be assigned consistently. For this individual, across analyses the unassigned parent was not consistent with either both parents unassigned, only the mother unassigned, or only the father unassigned. There were no untrue parent assignments (false assignments) observed across any of the assignment analyses, with all putative parents correctly excluded where the true parent/s were absent from the dataset tested.
Table 2. Success rate of parentage assignment analysis using three SNP marker sets (7,452 DArTseq, 4,236 DArTseq, 4,194 DArTcap) and three genotyping error rates that range from strict to conservative (0.01, 0.05 and 0.1) in Colony [30]. There were no untrue parent assignments for any dataset at any error rate. |
SNP dataset | Estimate genotyping error rate |
0.01 | 0.05 | 0.10 |
7,452 DArTseq | 97.2% | 97.2% | 97.2% |
4,236 DArTseq | 97.2% | 97.2% | 97.2% |
4,194 DArTcap | 91.8% | 95.8% | 98.6% |
2.3.2 Genomic relationship calculations
Genomic relationship values calculated with the 4,194 DArTcap markers showed high correlation (r2 = 0.98) to the 4,236 DArTseq marker panel (Fig. 3c), importantly indicating high concordance between the genotypes obtained across 195 samples through the two GBS approaches. As with the comparision of DArTseq 7,452 and 4,236 marker panels (see above) a number of pairwise relationships were less than 0 (no relationship). This is likely to be due to the structured nature of the wild samples included in the dataset [29, 31, 32].
To further explore the utility of the 4,194 DArTcap SNP panel, genomic relationship values were calculated for an additional independent set of G2 samples (n = 272) that were produced from novel Northern Territory sourced broodstock (i.e. not genotyped with either DArTseq or DArTcap). Utilising the GRM values of progeny alone, in the absence of parental genotypes, it is possible to clearly obtain sib-ship information (block structure in Fig. 4), including delineating full-sib and half-sib relationships. Furthermore, when compared to routine pedigree based relationship matrixes, it is possible, by assessing the values off the diagonal of the heat-map (within and outside blocks), to identify cryptic relatedness between otherwise traditionally unrelated individuals (Fig. 4).
2.3.3 Population segregation
The ability of the 4,194 DArTcap SNP array to segregate samples obtained from different populations across the Australian distribution of P. monodon was assessed by comparing these individuals through discriminant analysis of principle components (DAPC). By comparing the first two principle components (with PC1 explaining 52.8% and PC2 explaining 44.3% of the variation), three distinct clusters of individuals are clearly identified (Fig. 5). Specifically, samples obtained from within East Australia Coast, Northern Territory and Western Australia form their own clusters. The G1 broodstock, G1 progeny and G2 progeny also cluster with their source population (Northern Territory) showing no clear separation as a result of breeding practices (Supplementary Fig. 2). The associated regional structure is consistent with previously identified population structure [22].