Y-SNPs as evolutionary markers

Y-SNPs are slowly mutating bi-allelic markers used to reconstruct a human phylogenetic tree, to predict ancestral origins and to study evolutionary migration patterns [5,7–9]. Based on the ISOGG YBrowse database (2019–2020), the 841 designed amplicons would target 13,812 known Y-SNPs (http://ybrowse.org/gb2/gbrowse/chrY). As reported by the ISOGG YBrowse database, they can be further divided into 5,927 Y-SNPs with a still unknown phylogenetic position and 7,885 haplogroup-specific Y-SNPs where 30 Y-SNPs have been identified as private SNPs. These haplogroup-specific Y-SNPs define 1,212 unique Y-subhaplogroups (covering 96% of the Minimal Y-tree) [10].

In practice, sequencing 130 males with our CSYseq panel and data analysis with Yleaf [42] successfully enabled the identification of 15,611 Y-SNPs. As reported by the ISOGG YBrowse database (2019–2020), they can be further divided into 6,597 Y-SNPs with a still unknown phylogenetic position and 9,014 evolutionary haplogroup-specific Y-SNPs where 32 Y-SNPs have been identified by ISOGG as private Y-SNPs. The haplogroup-specific Y-SNPs target 1,443 unique Y-subhaplogroups including all main haplogroups (from ‘A’ to ‘T’) divided across the entire phylogenetic tree (Table 1). The output of the panel covers 445 haplogroups (97%) of the 458 haplogroups included in the Minimal Y-tree. The 13 Y-subhaplogroups not covered by the panel are B1, C1b1a1a1a1a, I1a2a1a1d1a1a2b1a, K1a, K1b, K2a1a, K2b2, M3, N1a1a1a1a1a6, R1b1a1b1a1a2a1b1a, R1b1a1b1a1a2a2, R1b1a1b1a1a2b3b and S1a2. In total, 129 samples covered the 445 haplogroups and one sample targets only 403 Y-haplogroups. The latter sample was also observed to have the lowest output number of Y-SNPs (7,284) and the lowest chrY alignment percentage (13%). Even though this was a challenging sample, it is still able to target 88% of the haplogroups included in the Minimal Y-tree [10]. In Table 1, it can be observed that the CSYseq contains an equal subhaplogroup distribution per main haplogroup compared to the Minimal Y-tree [10]. A complete phylogenetic tree including all CSYseq typed Y-subhaplogroups can be found in S1 Table within the Supporting Information file.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. CSYseq Y-SNP and subhaplogroup coverage.

https://doi.org/10.1371/journal.pgen.1009758.t001

Targeted Y-SNPs contained between 10 and 5,218 reads per sample with an average of 717 reads. For the 65 non-related samples, the total number of reads per Y-SNP ranged from 10 till 339,189 with 70% between 10,000 and 100,000 reads (Fig 1A). On average, there are 12,281 Y-SNPs typed per sample and even the least extensive sample still contained 7,284 well-typed Y-SNPs (Fig 1B). The number of typed Y-SNPs per sample was significantly correlated with the total number of reads identified in the FASTQ files (p = 4.33×10⁻⁷) and the number of reads aligned against chrY (p = 8.85×10⁻⁴⁰) (Fig 1B). This was as expected, since the more reads the sample has in total (FASTQ) or aligned with chrY, the more reads it has per chrY amplicon containing the Y-polymorphisms of interest. Sample quality statistics revealed a slightly significant (at the margin of statistical significance) correlation between typed Y-SNPs with initial chrY concentrations measured before library preparation (p = 1.30×10⁻³) and their degradation index (DI, p = 1.27×10⁻²) (see section ‘CSYseq robustness’, S2 Fig). When MPS output is compared to the limited Y-SNP panel typed by the SBE SNaPshot PCR-CE technique used in most laboratories, a successfully deeper Y-SNP subhaplogroup was genotyped for 66% of the samples due to the massive number of typed Y-SNPs (Fig 1C). On average, four phylogenetic branches deeper were detracted in which a maximum of ten branches was observed: from ‘R1a1a’ (R-M198) with SNaPshot-CE to ‘R1a1a1b1a3a2b2b’ (R-AM00559) with MPS. In 32% of the samples, both techniques resulted in the same final derived Y-SNP, but for two samples with subhaplogroups ‘J-M92’ and ‘R-L2’, SNaPshot-CE surpassed MPS in typing one branch deeper. For these latter two cases, the CSYseq panel was able to sequence both final Y-SNPs, but the markers did not pass the selected sequencing criteria of at least 10 reads and the base calling percentage of 90%. This is sample specific and not a limitation of the panel. J-M92 was observed to be typed in 90 samples with an average depth of coverage of 43 reads. And R-L2 was typed in all the other samples with a high average depth of coverage of 1,095 reads.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. CSYseq targeted Y-SNPs.

A. The 15,611 genotyped Y-SNPs with number of paired-end reads for all samples (threshold = 10 reads). B. Correlation between typed Y-SNPs and the total reads per sample obtained from FASTQ files (▲) and chrY alignment (). C. Number of typed Y-SNP subhaplogroup branches through SNaPshot-CE (black) compared to MPS (grey) sorted by their main Y-haplogroup.

https://doi.org/10.1371/journal.pgen.1009758.g001

The female sample revealed output for 399 Y-SNPs, where 205 Y-SNPs (51%) have an unknown phylogenetic position. Detailed Y-SNP analysis revealed no unambiguous haplogroup determination. Although the depth of coverage of these Y-SNPs was above the threshold of 10 single-end reads, they only had a median of 31 reads. Additionally, for all 399 Y-SNPs, the average reads obtained for all male samples were 13 times higher compared to the female sample. This is valuable information for forensic genetics to set a valuable read threshold when DNA mixture analysis needs to be performed.

Y-STRs as patrilineage markers

Y-STRs provide a high degree of variability due to their fast mutating multi-allelic variations. This makes them highly interesting for population genetics. In theory, the 214 designed amplicons of our CSYseq panel target a total number of 228 Y-STR loci. In practice, after sequencing and Y-STR data analysis with FDSTools, 28 Y-STR loci were excluded from the panel due to not being sequenced (n = 20), low depth of coverage (n = 7) or strong chrX homology which could not be filtered out due to strong sequence similarities (n = 1). This was as expected and can be explained by the 90% primer design success rate of DesignStudio (see before). The excluded Y-STR loci with detailed information are listed in S2 Table. Through additional analysis of the high quality sequenced chrY reads with Tandem Repeat Finder (TRF) [43], two novel Y-STR loci were identified which are sequenced by the CSYseq panel. As no information about these specific Y-STRs is yet available in literature or within the ISOGG YBrowse database, they were named CSY1 and CSY2. Further, CSYseq analysis for the double sequenced male sample exposed equal data output and the female sample revealed no output. This indicates that our CSYseq panel is chrY-specific and possible output allele calls as a result of chrX homology were successfully filtered out using our in-house created ‘CSYseq.analYzer’ tool (see Materials and Methods).

In total, the CSYseq panel covers 202 well-targeted Y-STR loci. Table 2 provides detailed information concerning their repeat motif and discrimination capacity. HGVS nomenclature and Y-chromosome positions of these Y-markers can be found in S3 Table. The 202 Y-STRs from the CSYseq panel include 15 Y-STRs from the commercially available CE kits (PowerPlex Y23 and Yfiler Plus): DYS19, DYS389I/II, DYS390, DYS391, DYS392, DYS448, DYS456, DYS635, Y-GATA-H4, DYS533, DYS549, DYS570, DYS643 and DYS460. 17 Y-STRs targeted by the CSYseq are also present in the commercial kits developed for MPS (ForenSeq and PowerSeq): DYS19, DYS389I/II, DYS390, DYS391, DYS392, DYS448, DYS456, DYS460, DYS522, DYS533, DYS549, DYS570, DYS612, DYS635, DYS643 and Y-GATA-H4. As an internal control, 21 Y-STR loci sequenced by the CSYseq were compared to previously obtained PCR-CE results from our in-house YForGen kit (46 Y-STRs) and commercial Y-kits [44]. We observed that all Y-STR allele calls were in accordance with our previous results, which confirms that the results of MPS are reliable. A total of 188 Y-STRs are simple Y-STRs with one variable repeat motif, while 14 Y-STRs contain a more complex double repeat. For example, ‘DYS463’ exists of both AAAGG[n] and AAGGG[n] as variable repeat motifs which were easily discriminated using FDSTools. Furthermore, 156 Y-STR loci are single-copy (SC) Y-markers, whereas the other 46 are multi-copy (MC) Y-markers, including three Y-STRs with four loci (-abcd). For most MC Y-STRs, it remained difficult to discriminate the different loci due to sequence similarities of the flanking and repeat regions. The results of the MC Y-STR loci with indistinguishable genome alignment were grouped together for further analysis. Equal to CE analysis, if the exact sequence per locus remains unknown, we sort them from short to long Y-STR allele call. This makes it possible to still perform Y-STR mutation analysis using the principle of Parsimony: the least number of changes indicates the most likely event. DNA sequences of all included Y-STRs are publicly available in the ISOGG YBrowse database (http://ybrowse.org/gb2/gbrowse/chrY).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Detailed information concerning all 202 CSYseq targeted Y-STRs.

https://doi.org/10.1371/journal.pgen.1009758.t002

The CSYseq targets 57 di-, 38 tri-, 83 tetra-, 22 penta- and 2 hexanucleotide repeats. For autosomal DNA analysis, the two genuine allele calls of dinucleotide STRs can be difficult to interpret due to their stutter fragments. For Y-chromosomal STR analysis, this is different as it mostly results into one allele call due to its haploid nature. For the single-copy dinucleotide Y-STRs, stutter fragments and the true allele call can easily be identified using FDSTools. But the panel does include six multi-copy dinucleotide Y-STRs that cannot be distinguished by sequence variance in the flanking regions. These Y-STRs resulted in multiple stutter peaks. Therefore, a more complex separation of stutter alleles and genuine heterozygous alleles was necessary. The reported difficulties with these Y-STR stutters were taken into account within the CSYseq.analYser file. This file sorts the sequences according to their number of reads to additionally filter out all stutters. An example of the different allele and stutter output combinations to interpret multi-copy dinucleotide YCAII-ab within our sample is provided in S4 Table.

The average number of reads per Y-STR locus was 150 reads, ranging from 5 (DYS448) to 619 reads (TRF17200) (Fig 2A). Only 11 Y-STR loci had an average number of reads below 10, though eight of them were genotyped for the majority of the samples. These markers may have insufficient coverage for challenging forensic samples, but this needs to be confirmed by future research. However, they can still be interesting to include into the panel for genetic-genealogy purposes and mass-screening for forensic familial searching. 137 Y-STRs are well-typed in more than 90% of the samples, 54 Y-STRs in 60 to 90% and 11 Y-STRs in 30–60% of the samples. The number of Y-STRs typed per sample is visualized in Fig 2B. On average, there are 184 Y-STRs typed per sample and even the least extensive Y-haplotype still contained 115 well-typed Y-STRs. There was no significant correlation between the number of typed Y-STRs and the number of reads within their FASTQ file (Fig 2C). However, the number of typed Y-STRs was observed to correlate significantly with the number of reads aligned against chrY (p = 6.90×10⁻³). Sample quality statistics revealed a slightly significant (at the margin of statistical significance) correlation between typed Y-STRs and the initial chrY concentrations measured before library preparation (p = 2.80×10⁻²), but not with the DI (see section ‘CSYseq robustness’, S2 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. CSYseq targeted Y-STRs.

A. The average number of reads of the 202 targeted Y-STR loci of the CSYseq panel. B. The number of typed Y-STRs per sample. C. Correlation between the typed Y-STRs per sample and the total number of reads per sample obtained from FASTQ files (▲) and chrY alignment (×). D. A heatmap visualizing the allele ranges and frequencies per Y-STR.

https://doi.org/10.1371/journal.pgen.1009758.g002

All 202 Y-STR loci were investigated in detail using GenAlEx to determine the allele call frequencies with allele ranges (Fig 2D), discrimination capacity, and average repeat sizes (Table 2). The 14 Y-STR loci having multiple variable repeat units were divided into -M1 and -M2. The smallest variable repeat size contained only two repeats (DYS452-M2, DYS635-M2 and DYS19-M2), while the largest repeat number observed contained 30 repeats (DYS612). Detailed double repeat sequence variability with their allele call frequencies for the 14 variable complex and compound Y-STRs can be found in S5 Table within the Supporting Information file. Average discrimination capacity was 0.69 for complex and compound Y-STRs and 0.44 for simple Y-STRs. For 11 Y-STRs, no allele diversity was observed between the samples included in this study, wherefore consequently a discrimination capacity of zero was calculated.

Y-STR mutation analysis was conducted through Y-haplotype comparison between male relatives within the genealogical pairs and deep-rooting pedigrees. A detailed overview of the mutation statistics per Y-STR loci are listed in Table 3. The number of generations covered per Y-STR loci are on average 1,083 meioses and fluctuated between 218 and 1,279 meioses. This fluctuation can be explained by the fact that some Y-STR markers were not successfully typed in all samples. A total number of 910 Y-STR differences was observed over 214,859 allele transfers (Table 3). In total, 759 one-step, 98 two-step and 53 multi-step differences were observed. For 66 Y-STRs, no allele call differences within the sequenced genealogical pairs were observed. The mutation rates of the other 136 Y-STRs are listed with their 95% confidence interval (CI) in Table 3 and visualized in Fig 3A. An overall average mutation rate of 4.57×10⁻³ mpg (95% CI: 4.29×10⁻³–4.86×10⁻³) was observed for the CSYseq panel. When we exclude the Y-STRs without an observed mutation in our study, an average mutation rate of 6.64×10⁻³ mpg was obtained with a minimum of 4.15×10⁻⁴ mpg (DYS371-abcd) and a maximum of 4.13×10⁻² mpg (TRF14783). The mutating Y-STRs can be subdivided into 15 slow mutating Y-STRs (<10⁻³ mpg), 68 moderate mutating Y-STRs (≥10⁻³ to <5×10⁻³ mpg), 27 fast mutating Y-STRs (≥5×10⁻³ to <10⁻² mpg) and 26 rapid mutating Y-STRs (≥10⁻² mpg, Fig 3A, red line) [45]. The individual mutation rates of 101 Y-STRs were compared to literature [14,17,46–48]. In total, 95% of these Y-STRs were in accordance with literature, which means that a significantly different mutation rate was observed for only five Y-STRs (DYS390, DYS490, DYS525, DYS606 and DYS612). For the influencing molecular factors, a significant positive correlation between the individual Y-STR mutation rates with the average allele size (number of repeats) (p = 8.34×10⁻¹⁰) was observed (Fig 3B). Mutability rates had no significant difference between simple, compound or complex repeat Y-markers (Fig 3C). Further, significant differences were identified in the mutation rates between di-, tri-, tetra- and pentanucleotide Y-STRs, but no significant difference between tri- and tetranucleotide Y-STRs nor a linear correlation was observed (Fig 3D).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. CSYseq Y-STR mutation analysis.

A. Individual Y-STR mutation rates with their 95% CIs. The RM Y-STR treshold (10⁻² mpg, dashed line) and average mutation rate (4.57×10⁻³ mpg, black line) are indicated. B. Positive significant correlation between the Y-STR mutation rate and the average number of repeats. C. The average Y-STR mutation rates of the different repeat motif types (simple, compound and complex). D. Average mutation rates per length of the repeat unit (bp).

https://doi.org/10.1371/journal.pgen.1009758.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. CSYseq Y-STR mutation analysis.

https://doi.org/10.1371/journal.pgen.1009758.t003

Through detailed mutation analysis of complex and compound Y-STRs, it was observed that Y-STR differences occurred more frequently within the longest repeat sequence. For example, in DYS725-abcd, which has a compound repeat structure being GT[n]GTCT[n], the average number of repeats is respectively 20 and 4, and the observed number of mutations per motif is 32 and 17. Furthermore, we observed that three markers (TRF10691, DYS463 and DYS725) showed allele call differences in five genealogical pairs for both variable motifs at the same time. For three couples, these differences were found on DYS725 (NC_000024.10:g.24738202) e.g. one relative had GT[19]GTCT[4], while the other contained GT[20]GTCT[5] which reveals two independent one-step mutations in parallel. The other two genealogical pairs contained both a parallel mutation which would have remained hidden through CE as the two mutations resulted in the same allele call: 22 repeats for DYS463 (NC_000024.10:g.7775468) with AAAGG[7]AAGGG[15] ↔AAAGG[8]AAGGG[14] and 25 repeats for TRF10691 (NC_000024.10:g.15550131) with TG[21]N[10]TG[4] ↔TG[20]N[10]TG[5].

Through comparison analysis of the 202 Y-STR loci, it was possible to distinguish all non-related and related males, providing 130 unique Y-haplotypes. Using the Y-STR differences observed over the 136 mutating loci, the CSYseq succeeded in making a distinction between all paternally related males. On average, they were separated by 18 generations and discriminated by 13 Y-STR changes. A minimum number of four Y-STR differences was observed for a couple separated by 18 meioses, whereas a maximum of 22 Y-STR changes for two couples could be observed separated by 21 and 29 meioses. No significant correlation was observed between the number of generations and the number of mutations. This can be explained by the inclusion of fast and rapid mutating Y-STRs in the CSYseq panel and by the occurrence of back and parallel mutations which increases with the generational distance within genealogical pairs [23].

CSYseq robustness

A schematic overview of the MPS library quality and chrY data analysis steps is provided in S2A Fig. The TruSeq Custom Amplicon Low Input kit (Illumina, San Diego, CA, USA) recommends a DNA input of 10 ng and DNA concentration of 2.5 ng/μl [49]. The chrY DNA input concentration of all sequenced samples measured using PowerQuant qPCR was between 1.75 and 17.58 ng/μl (average 5.69 ng/μl) with a degradation index (DI) from 0.90 to 4.79 (average 1.92; S2B Fig). No significant correlation between chrY concentration with DI could be observed. Five samples did not fulfil the recommended input concentration from which two samples had a DI exceeding the manufacturer’s threshold of 2 [50]. The samples encountering the highest DI (4.79) and the lowest chrY concentration (1.75 ng/μl) are respectively indicated by the labels ‘d1’ and ‘c1’ throughout S2 Fig.

Library preparation quality control measured by the 2100 BioAnalyzer indicated a library peak size for all samples between 357 and 397 bp (average 380 bp) and a library concentration between 0.01 and 8.82 ng/μl (average 1.6 ng/μl; S2C Fig). BioAnalyzer library concentrations showed a significant correlation with the initial identified chrY concentrations (p = 4.09×10⁻³), but not with the DI. Additionally, normalized KAPA SYBR qPCR Ct values (S2D Fig) also revealed a significant correlation with the chrY concentrations (p = 1.27×10⁻⁶), but are only slightly significant with the DI (p = 0.006, R² = 0.062). As expected, normalized KAPA qPCR Ct values correlated significantly with the BioAnalyzer library concentrations (p = 1.24×10⁻¹⁰, R² = 0.293). The number of FASTQ file reads per library output after sequencing (S2E Fig) did not significantly correlate with both the initial chrY concentrations and the DI due to library normalization.

FASTQC software [51] flagged 63 samples with high per sequence base quality, meaning that, for both paired-end reads, the lower quartile of the first 150 bp did not have a FASTQC quality Phred score below 20. For all samples, the read position where FASTQC Phred scores went below 20 ranged from 85 to 278 bp (average 172 bp; S2F Fig). Again, only a slightly significant correlation could be observed with the initial chrY concentrations (p = 2.54×10⁻³), but not with the DI. Besides, the percentage of read alignment against GRCh37/hg19 reference genome and chrY also turned out to be only significant with the input chrY concentrations (p = 4.33×10⁻³ and 6.16×10⁻⁴) and not the DI. Remarkable was that both samples with the highest DI (4.79 and 4.69) showed a high and low quality. The FASTQC Phred score below 20, defining low per sequence base quality, are respectively at 258 bp and 104 bp and they have a chrY alignment of 61% and 36%. However, a high number of Y-markers (12,709 and 12,565 Y-SNPs; 186 and 192 Y-STRs) could still be sequenced (S2F and S2G Fig). In general, the number of typed Y-SNPs using Yleaf (S2G Fig) showed a slightly significant correlation with the input chrY concentrations (p = 1.30×10⁻³) and DI (p = 1.27×10⁻²). Yet the number of typed Y-STRs using FDSTools (S2H Fig) only turned out to be slightly significant with the input chrY concentrations (p = 2.80×10⁻²), but not the DI. Therefore, the success rate of the CSYseq panel was not clearly observed to be influenced by the initial chrY concentration or degradation index of a sample.

Further, we focus on sample d1 with the highest degradation index (4.79) and c1 with the lowest initial chrY concentration (1.75 ng/μl), both indicated throughout S2 Fig. We observed that they both encountered low concentrations after library preparation measured by the BioAnalyzer and KAPA qPCR (S2C and S2D Fig). Surprisingly, they both contained relatively high paired-end read outputs of respectively 1,065,926 and 939,462 reads. The paired-end read alignment against the GRCh37/hg19 reference genome differed strongly with respectively 85% and 25% and for chrY alignment this was respectively 61% and 16%. Consequently, d1 contained an overall higher FASTQC quality Phred score (from 258 bp below 20) compared to c1 (from 85 bp below 20; S2E and S2F Fig). As a result, the number of typed Y-SNPs was slightly higher for d1 (12,706 Y-SNPs; 81%) compared to c1 (11,139 Y-SNPs; 71%; S2G Fig). But, remarkably, with the high number of typed Y-SNPs, they both still resulted in well-typed deep Y-subhaplogrouping. The CSYseq kit even added eight branches in c1 from ‘I2a1b1’ with CE to ‘I2a1b1a2b1a1a1’. Additionally, a high number of typed Y-STRs for both d1 and c1 was still possible, respectively being 186 Y-STRs (92%, average = 245 reads per Y-STR) and 182 Y-STRs (90%, average = 51 reads per Y-STR; S2H Fig).

Y-SNPs as evolutionary markers

Y-SNPs enable the reconstruction of a well-preserved male phylogenetic tree. Neighboring populations represent a comparable evolutionary haplogroup distribution, while different continents can exhibit large differences. Y-SNPs are evolution markers which are more frequently present within specific geographical regions, for example ‘R-M269’ to West-Europe, ‘E1b1b’ to North-Africa and ‘Q’ to America [6,52–54]. This is of course without taking recent migration into account. The CSYseq panel successfully enabled the identification of 15,611 Y-SNPs, which is even more than the panel was estimated to target in theory due to primer homology. In total, 9,014 Y-SNPs are defined by the ISOGG YBrowse Database (2019–2020) as haplogroup-specific Y-SNPs targeting 1,443 evolutionary Y-subhaplogroups (Table 1). Since the Minimal Y-tree by Van Oven et al. is commonly used in Y-SNP genotyping, a large coverage is desired which facilitates Y-SNP genotyping [10]. With our CSYseq panel, 445 out of 458 subhaplogroups present in the Minimal Y-tree were genotyped, resulting in a successful coverage of 97%. In addition, 998 additional subhaplogroups were covered using the CSYseq panel. The number of phylogenetic Y-SNPs typed with the CSYseq is more than ten times higher than the current most extensive Y-SNP MPS kit on the market, namely the Ion AmpliSeq HID Y-SNP Research Panel v1 (Thermo Fisher Scientific) with 859 phylogenetic Y-SNPs [41]. Consequently, it analyzes 56% less Y-haplogroups (n = 640) than the CSYseq panel (n = 1,443). This large number of Y-SNPs spread over the euchromatic region of the Y-chromosome is interesting when we have to analyze challenging or degraded samples. Additionally, the CSYseq targets all main haplogroups (‘A’ to ‘T’) divided across the entire human Y-chromosome phylogenetic tree. This custom-made panel can therefore serve as a powerful tool to identify paternal evolutionary lineages and provide more information about males with any biogeographical background around the world. It is interesting to note that every main subhaplogroup present in our population (Belgium and the Netherlands) contained a high number of typed Y-SNPs in our sample (S3 Fig). Although sampling males with other biogeographical backgrounds still remains necessary, we are confident by the discrimination power of the CSYseq due to its high Y-SNP coverage across all main haplogroups of the human Y-chromosome phylogeny. In addition, about 6,597 Y-SNPs targeted by the CSYseq currently have no available information on their biogeographic paternal ancestry. As the list of Y-SNPs in the ISOGG index is constantly being updated with newly identified Y-SNPs, even more haplogroup-specific Y-SNPs will be identified using our CSYseq when more extensive geographic sampling will be performed in the future.

As a result of MPS, identification and mapping of Y-SNPs on the phylogenetic tree in population studies is facilitated as more Y-SNPs can be targeted and a large sample input (from a specific population) can be sequenced. Additional population data of newly discovered Y-SNPs still remains necessary to obtain as Y-SNP expansion gives rise to some complications. First, universal names for Y-SNPs are non-existing, meaning that equality and comparisons between studies and different phylogenetic trees remain extremely difficult to achieve [55]. For example, the Y-SNP R-L11 has also been described as R-S127 and R-PF6539. With the CSYseq, we present a universally exchangeable set of interesting Y-markers which can provide a basis for a uniform nomenclature between Y-SNP population studies worldwide. Second, many Y-SNPs are already mapped in a phylogenetic tree despite the lack of large-scale population data. Consequently, their position on the phylogenetic tree is deceitful and cannot be used as additional biogeographical background information [2,55,56]. For 183 Y-SNPs, a combination between ancestral and derivative alleles was observed within our samples indicating that these markers exhibit larger diversity within the population of the Low Countries (Belgium and the Netherlands). Herein, many unlisted, unknown Y-SNPs and Y-SNPs associated with multiple haplogroups are incorporated due to which they may not be allocated to a specific haplogroup yet be private to some families. The distribution of ancestral and derivative calls within a specific Y-SNP can serve as valuable information for population genetics since Y-SNPs with larger diversity have a higher discrimination power between individuals.

Y-STRs as patrilineage markers

Y-STRs are commonly used DNA polymorphisms to find distant or close relatives through patrilineage identification in interdisciplinary research fields. The current most extensive MPS kit for Y-STR haplotyping is the ForenSeq DNA Signature Prep kit (Illumina, 2015) targeting 24 Y-STRs [37]. As 24 Y-STRs is even smaller than the commercially available PCR amplification kit Yfiler Plus, including 27 Y-STR loci for fragment analysis by CE, the development of a more extensive MPS kit, sequencing more Y-STR markers, was definitely required. Our custom-made CSYseq panel successfully targets 202 Y-STR loci, where 15 Y-STR loci are present in today’s commercial PCR amplification Y-kits (PowerPlex Y23 and Yfiler Plus) and 17 Y-STR loci in commercial MPS kits (ForenSeq and PowerSeq). This is useful for chrY comparison with, for instance, the YHRD reference database where haplotypes from all over the world are gathered [57]. However, further inclusion of the other forensic commercially available Y-STR loci in the CSYseq will increase compatibility with the existing Y-haplotype kits. Next, the CSYseq output is highly chrY-specific as Y-STRs with primer homology on other chromosomes were successfully excluded by our in-house ‘CSYseq_analYser’ due to sequence or allele call differences. Only one Y-STR was excluded from the panel as it exhibited homology with chrX to such an extent that an extreme high number of reads was obtained for the homology allele calls. On average, the 202 Y-STRs within the CSYseq panel were well-typed in 90% of the samples, containing a high average depth of coverage from 150 reads. Unfortunately, for the majority of MC Y-STRs, FDSTools was still unable to distinguish between the different loci due to the extreme sequence similarity of the flanking regions for both loci. This causes non-distinguishable MC Y-STR loci to have a more complex separation of stutter alleles and genuine heterozygous alleles. However, the FDSTools software did succeed in genotyping more complex repeat structures. Moreover, the reported difficulties with these Y-STR stutters were taken into account within the CSYseq.analYser file. This file sorts the sequences according to their number of reads to additionally filter out all stutters. This is especially important for the dinucleotide Y-STRs, which result in multiple stutter peaks. Fortunately, no Y-STR loci had to be excluded due to unsuccessful genotyping caused by the complexity of the repeat structure, as was the case with previous STRaitRazor genotyping [58].

Y-STR discrimination power

The CSYseq includes 46 multi-copy Y-STRs and 26 RM Y-STR loci, which strongly increases the discrimination power of the panel. Our detailed population study revealed no variation for 11 Y-STRs, which was also observed before [59]. Their rather small allele size could explain the low variability. Moreover, our sample consisted exclusively of men from the Low Countries, which means that Y-marker variability in different populations cannot be excluded. Next, 25 Y-STRs from the CSYseq panel exhibited a high degree of variation (discrimination capacity >0.75) in our population. The most discriminating Y-STR in our panel is DYS712, which ranges from 11 to 22 repeats with a discrimination capacity of 0.85. This is presumably a result of the rapidly mutating nature of this locus (2.19×10⁻² mpg). Additionally, the compound and complex Y-STRs contain a higher discrimination capacity of 0.7 compared to 0.4 for the simple Y-STRs. Further, these markers provide the advantage to discriminate males containing equal total allele sizes through detailed repeat motif analysis. For instance, within the complex Y-STR DYS552, this resulted in twice the number of different allele calls typed with MPS (n = 12) compared to CE (n = 6). This provides a discrimination capacity of 0.8 for MPS compared to the 0.6 of CE. The inclusion of these highly discriminative markers to the CSYseq panel guarantees the robustness and reliability for kinship analysis or forensic familial searching. In general, the allele calls for the Y-STRs within the CSYseq ranged from 6 to 30 repeats. This wide range of allele sizes is valuable for the various demands of the panel. Y-STRs with larger allele sizes are known to exhibit larger mutation rates, which enlarges the panel’s discriminative power [14]. On the other hand, Y-STRs with short allele sizes are easier to target than long repeat stretches, since the latter are more susceptible to degradation by environmental factors and therefore often drop out of the profile in challenging samples [60]. However, this should be confirmed through additional sequencing analysis on degraded samples, such as forensically challenging or ancient DNA samples.

CSYseq robustness

Per sample, the CSYseq panel typed on average 12,281 Y-SNPs and 184 Y-STRs. This number was slightly dependent on the chrY concentration, but not with the initially identified degradation index (DI) (S2C and S2H Fig). For the five ‘challenging’ samples, with low chrY concentrations and high DI, over 10,800 Y-SNPs and between 127 and 182 Y-STRs were still well-typed. This reveals no clear influence of the initial concentration or DI on the success rate of MPS. Remarkably, despite the rather low chrY concentration for c1 and the high DI for d1, they both typed a high number of Y-SNPs (11,139 and 12,706) and Y-STRs (182 and 186). Therefore, it can be assumed that the input degradation and concentration requirement of the TruSeq Custom Amplicon Low Input kit is flexible and only the presence of DNA, and not the quantity, should be determined before library preparation. Although this statement needs to be confirmed by future studies with more challenging samples, this is in accordance to observations made by Poetsch et al. [62]. They concluded that only 0.00 ng/μl indicates the absence of DNA with the PowerQuant kit, which demonstrates that a concentration of >0.01 ng/μl could still provide a full Y-STR profile after multiplex PCR and CE. Forensic samples are often challenging due to the low amounts of DNA or high level of degradation. It can be concluded that the CSYseq Y-marker targeted resequencing is still successful with our lowest quantity DNA samples, making our CSYseq panel advantageous for challenging samples. Although high accuracy and sensitivity techniques such as MPS have already been observed as a solution compared to CE [63,64], additional analysis on degraded and low concentration samples is still recommended in order to have a more clear perspective on the precise CSYseq sensitivity.

CSYseq applications

As the CSYseq panel exclusively targets SNPs and STRs positioned on the Y-chromosome, the output data will be valuable for a wide range of genetic-genealogical applications in interdisciplinary research as it provides valuable paternal and biogeographical background information: family history, population genetics, evolutionary biology, forensic science and even medical diagnostics.

In-depth chrY genotyping helps to unravel family history by providing more detail in patrilineal relatedness between relatives. Genealogists could make their expanded chrY-profile public or available in a database hoping to find an unexpected patrilineal relation. For population genetics, the link between a surname and patrilineage provides the opportunity to detect signals of past or recent population stratification and migrations which are still undetectable within genomic analysis of the limited number of markers [65]. Additionally, the fixed set of Y-markers sequenced with the CSYseq could avoid problems with Y-SNP nomenclature and dataset differences as is currently observed in different population studies. For molecular biology, this panel will contribute to the knowledge concerning the Y-STRs mutation rate and the molecular mechanism of mutations, together with the general understanding of Y-STR evolution in the human genome. For evolutionary biology, the CSYseq enables deeper or equal subhaplogroup determination than CE in 98% of our samples, making Y-SNP haplogroup identification through the CSYseq panel definitely more convenient. Additionally, typing a large number of Y-STRs allows to study haplogroup phylogenetics in more detail and provide valuable information for tMRCA estimations and evolutionary dating to reconstruct phylogenetic trees in more detail. Furthermore, increasing Y-STR diversity by including sequence variation in the repeat region or in the flanking regions is beneficial for molecular, evolutionary and population genetic studies as this will increase their dataset resolution.

Y-chromosome knowledge gained with the CSYseq panel also provides applications in medical diagnostics concerning male infertility which affects one in five infertile couples and one in 20 men [66]. ChrY is essential for male fertility due to the presence of several spermatogenesis-related genes. High resolution mapping of hundreds of Y-SNPs and Y-STRs will allow to identify deletions or chromosome abnormalities, helping to identify the molecular mechanisms underlying male infertility. Next to infertility, studies have shown evidence that genetic variation within the NRY could also play a part in determining cardiovascular risks as well as immune and inflammatory responses in men [20]. Correlations between the subhaplogroup and an increased disease risk were already established. Haplogroup ‘I’ shows a correlation with coronary artery disease, while haplogroup ‘N’ with infertility [20,66]. Genotyping Y-SNPs could therefore serve as a prevention analysis to identify men with an increased risk. The availability of a molecular tool to type hundreds of Y-SNPs will therefore provide the opportunity to utilize the power of phylogenetic analysis which is currently not widely used in medical genetics to explore the potential chrY contribution to complex polygenetic traits.

For forensic science, this panel can resolve some complex paternity kinship questions or provide assistance with the identification of an unknown perpetrator [2,17]. The CSYseq panel enables to find both distant paternal relatives through approximately 81 slow mutating Y-STRs (<10⁻³ mpg) and distinguish closely related individuals through 26 RM Y-STRs (≥10⁻² mpg), which significantly increases the level of discrimination useful in forensic human identification processes [7]. This combination of Y-STRs is especially useful for familial searching, where the donor of an unknown trace has to be identified by searching for a male relative in a chrY database or through a large-scale voluntary DNA mass-screening [2,67]. Additionally, around 32 private Y-SNPs are typed which originated more recently, which makes it possible to link them to a specific population or even a single family [6]. If such a private SNP is found in a trace of the perpetrator, the number of suspects can be reduced significantly to that specific population or family. Thanks to the number of forensically interesting Y-chromosomal markers included in the CSYseq panel, a wide range of applications in the forensic field can be reached with this unique MPS panel. Through extensive Y-STR profiling, the CSYseq facilitates paternal kinship testing, disaster victim identification, cold case investigation and missing person identification [2].

MPS should become state-of-the-art in the near future to overcome the limitations of the traditionally used CE fragment analysis. CE analysis of about 15 Y-SNPs (one multiplex) and 46 Y-STRs (four multiplexes) takes multiple assays and approximately one day depending on the purification of the SNaPshot-PCR products. This does not provide detailed sequencing or subhaplogroup information as typically more than one Y-SNP multiplex has to be tested for deep subhaplogrouping. However, to implement MPS in routine DNA analysis, there are still some issues that need to be addressed concerning the nomenclature, minimum number of reads, sequencing errors, MPS strategy, data storage, minimal available MPS data and software adjustments towards new allele data [68]. To facilitate MPS in genetic research, there is also a need for specialized, expensive equipment and the improvement of the practical MPS work, as sometimes intensive laboratory effort is needed for library preparation and sequencing [69]. But, on the other hand, if this number of output needs to be obtained through CE-PCR, a much more intensive laboratory effort is needed. Furthermore, we are convinced that the low cost of less than €100 per sample (for the panel, library preparation and MiSeq run) of the CSYseq would help to apply this panel in population genetic studies worldwide.

In the end, with this study, we were able to successfully design the first extensive chrY sequencing kit. We tested the performance of the panel on samples with different concentrations and degradation levels. With the CSYseq, we offer a starting point for further investigation. To implement the kit in forensic chrY analysis, it will be necessary to further validate its repeatability and reproducibility in addition to determine its sensitivity in DNA mixture samples. Moreover, the inclusion of father-son pairs will be highly interesting to assess the discrimination power of the CSYseq panel in closely related males. This will provide additional mutability rate information for complex, compound, RM and MC Y-STRs.

Ethics statement

By means of written informed consents, permission for DNA analysis and scientific publication of the anonymized results was granted. Ethical approval has been allocated by the Ethical Commission of University Hospital Leuven (S55864, S59085).

The Y-chromosome is playing an increasingly important role in evolutionary biology and population genetics [6]. However, we currently lack a universal tool for sequencing these interesting Y-polymorphisms. Until now, genotyping mostly relied on fragment analysis for Y-STRs or a single-base extension assay for Y-SNPs, which is both based on capillary electrophoresis (CE). In order to fill this gap, we present the CSYseq, the first extensive chrY-specific targeted resequencing custom-made panel targeting both evolutionary Y-SNPs as familial Y-STRs. Therefore, the CSYseq is a unique tool to indicate both distant evolutionary lineages as close familial kinships, which provides biogeographical background information and paternal lineage identification.

Samples and DNA extraction

A total of 130 males with their residence in Belgium or the Netherlands (the Low Countries) were selected from previous studies investigating extra-pair paternity rates, haplogroup-specific Y-STR mutation rates, parallel Y-STR evolution and chrY-surname correlation [17,27,30,44,70,71]. These samples are further subdivided into five non-related males, 61 genealogical pairs with confirmed biological kinship and four extensive deep-rooted pedigrees enclosing three or four paternally related males. A total of 65 males within the study were confirmed to be unrelated. Judicial kinships between the relatives were certified per archival data and biological kinships were confirmed through 46 Y-STRs and 183 Y-SNPs genotyped using CE as described in detail within [17]. Relatives are separated by at least nine generations with an MRCA living between 1280 and 1900. Samples were collected via buccal swabs (Whatman OmniSwab, Sigma-Aldrich, USA).

For this study, DNA samples were re-extracted with the SwabSolution Kit (Promega, Madison, WI, USA). ChrY concentrations were quantified using an adapted protocol of the PowerQuant System kit (Promega, Madison, WI, USA). For 2 μl DNA-extract, 4 μl AmpSolution Reagent was added which was compensated by reducing the volume of Amplification Grade H₂O. Samples were quantified using the Applied Biosystems 7500 Real-Time PCR System and the HID Real-Time PCR Analysis Software v1.2 (ThermoFisher Scientific, Waltham, MA, USA). qPCR was performed under the following conditions: 2’ at 98°C followed by 39 cycles of 15” at 98°C and 35” at 62°C. Autosomal and Y-chromosomal DNA concentrations together with a degradation index (DI) per sample was obtained.

Custom panel, library preparation and sequencing

Regions of interest to target several pre-selected Y-STRs and Y-SNPs were carefully chosen based on literature. Y-STRs were selected on their sequence motif, Y-chromosome location and/or (if available) previously observed mutation rates and included all 46 Y-STRs of our in-house YForGen kit [14,17,47,59]. Y-SNPs were selected based on the Minimal reference Y-tree (www.phylotree.org/Y/tree/index.htm), containing 759 branch-defining Y-SNPs which targets 458 subhaplogroups distributed over the entire Y-SNP phylogenetic tree [10]. This resulted into a list of 865 defined chrY regions (cumulative 39,126 bp) containing 251 Y-STRs and 772 phylogenetic Y-SNPs. Selected regions ranged from 1 bp for Y-SNPs up to 950 bp for Y-STRs or a combination of Y-STRs and Y-SNPs. Primer pairs were designed using DesignStudio (Illumina, San Diego, CA, USA) with regard to the TruSeq Custom Amplicon Low Input kit with an amplicon length of 250 bp and avoiding SNPs in the primer positions (1000 Genomes as variant source). A final panel was obtained of 857 amplicons with their target length between 225 bp and 897 bp (average of 248 bp) encompassing 209,248 bp chrY sequence. Regions longer than 250 bp were covered by multiple amplicons with a combination of forward and reverse targets. In total, 228 Y-STRs and 757 Y-SNPs of the polymorphisms of interest could be genotyped with this panel.

The library preparation protocol was performed according to the guidelines of Illumina. With the TruSeq Custom Amplicon Low Input kit (Illumina, San Diego, CA, USA) a DNA input of 10 ng is recommended [49]. The custom designed primers (CAT, custom amplicon tube) were hybridized on the Veriti 96-well Thermal Cycler (ThermoFisher Scientific, Waltham, MA, USA). Unbound oligos were removed by magnetic sample purification beads (SPB) using four wash steps. PCR for primer extension and ligation was followed by library amplification while adding dual-index i7 and i5 adapters (TruSeq Custom Amplicon Index Kit, Illumina). To achieve the desired library yield and specificity, the optimal PCR cycle number for our oligo pool was 27 cycles (701–999 amplicons) according to the TruSeq kit guidelines. Libraries were again purified with SPB and washed three times to clean up other reaction components. Library quality was checked by DNA electrophoresis using the DNA 1000 kit (range 0.5–50 ng/μl) on the 2100 BioAnalyzer Instrument (Agilent, Santa Clara, CA, USA). An additional purification step was included for libraries with high primer or adapter dimers. In this step, purification beads were included in a ratio of 0.8:1 to filter out the shorter fragments (<200 bp). Libraries of approximately 370 bp and a sufficiently large peak height were normalized by qPCR DNA quantification using the KAPA SYBR Library Quantification kit for Illumina Platforms (KAPA Biosystems, Wilmington, MA, USA). P5 and P7 primers (Integrated DNA Technologies IDT, Leuven, Belgium) were diluted to 10 μM, libraries were diluted to a 1:10.000 ratio and analyzed in duplo with the Applied Biosystems 7500 Real-Time PCR System and the HID Real-Time PCR Analysis Software v1.2 (ThermoFisher Scientific, Waltham, MA, USA). A relative concentration of each sample was calculated using the comparative Ct method. For each sample, the highest Ct-value (the calibrator sample) was deducted from their Ct-value resulting in ΔCt. To determine the normalized volumes of each sample, all samples were divided by 2^-ΔCt and the calibrator volume was set at 10 μl to pool the libraries. The normalized pool of libraries was paired-end sequenced (2 × 300 bp), using the Miseq System and the MiSeq Reagent Kit v3 (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. In total 132 libraries were sequenced, including one female sample and one male sample twice as an internal control. Within the female sample, Y-markers can only be called by our CSYseq.analYser for Y-STRs or by Yleaf for Y-SNPs when the entire targeted region including primer positions is translocated. Sequencing and sample de-multiplexing using index-barcodes was done by Genomics Core (UZ Leuven, Belgium). The theoretical specification of 44–50 million reads per run (MiSeq Reagent Kit v3, Illumina) could be reached with 30 pM and a PhiX of 5–10%.

Data analysis

Sample quality control was executed using FastQC v0.11.8 software (Babraham Institute, Cambridge, UK) in which the per base sequencing quality was checked. Primers were trimmed using Galaxy tools and chrY alignment was done with SAMtools (SAM file) against the GRCh37/hg19 reference genome and visualized using the Integrative Genomics Viewer 2.8.0 (IGV, BAM/BAI file) [72]. Paired-end read alignment percentage per chromosome was calculated and compared to investigate primer homology.