Whole genome sequencing of Turkish genomes reveals functional private alleles and impact of genetic interactions with Europe, Asia and Africa

Whole genome analysis and variant discovery

We recruited 16 healthy volunteers from across Turkey (Figure 1A). The individuals were included in the study irrespective of their mother-tongue/ethnicity; we refer to them collectively as “Turkish”.

We isolated DNA from blood samples obtained from each individual, and generated whole genome shotgun (WGS) sequence data at high coverage ranging between 32 to 48X, using the Illumina platform (Methods, Table 1, Additional file 1: Table S1). Using BWA [23] and GATK [24], we identified an average of 3,254,782 SNPs per individual, and a non-redundant total of 8,161,894 SNPs, of which 651,936 (7.99%) were found to be novel when compared to dbSNP v132 and the 1000 Genomes Project (Figure 1B). We observed that 3,583 of the novel SNPs found were identified as non-synonymous mutations in protein coding sequence (Figure 1C). We also indentified a total of 1,765,584 short (<50 bp) indels. Additionally, using VariationHunter [25] we identified 1,751 large (50 bp to 100 Kbp) deletions per individual on average (Table 1).

To estimate the false discovery and false negative rates (FDR and FNR) of our SNP calls, we analyzed the same samples using Affymetrix 6.0 SNP arrays (Methods) to genotype for common SNPs. We then compared the microarray results with the SNP calls we generated from WGS. Assuming the SNP arrays as the gold standard, we estimated the FDR to be at 0.174% and FNR as 0.209%.

Common genetic variation in contemporary Turkey

Next we studied the 16 genomes with respect to major patterns of population diversity, structure, and migration, comparing these profiles with those in the 1000 Genomes Project [2]. We found that genome-wide nucleotide diversity (π) in Turkey is comparable to that observed in Europeans, while lower than in Africans or admixed Native American populations, and higher than in East Asians (Additional file 2: Figure S1). Even though we expect lower values of π in Turkish populations as compared to African populations, it was surprising to find π in European populations equal to or slightly higher than that observed in Turkey. We expected that Turkish genomes might exhibit significantly higher nucleotide diversity, given Turkey’s location at the crossroads of out-of-Africa migrations, as well as more recent population movements [9]. Still, this result should be taken with caution given differences in SNP calling procedures and power between the Turkish and 1000 Genomes Project datasets, which could potentially affect estimated diversity levels.

To obtain an overview of population relationships, we conducted principal component analyses using 16 individuals randomly selected from each one of the 14 populations within the 1000 Genomes Project, combined with the profiles from Turkey (Methods). The first two principal components of the combined dataset replicated the global genetic structure identified in earlier studies (e.g. [26]). Within this picture, the genetic variation observed among Turkish population clusters with variation observed in European populations (Figure 2A). The genome profiles from Turkey also overlapped with those of Native American populations. This is interesting, because the genetic affinity of Native American populations to both East Asia and Europe has traditionally been attributed to their Asian origins and subsequent admixture with Europeans [27]. However, this picture has recently been reinterpreted as evidence of shared Ancient North Eurasian ancestry in Native Americans and Europeans [28]. This result implies that significant Ancient North Eurasian ancestry might also be found in Turkish genetic profiles; this requires further study.

We then estimated the maximum likelihood population tree with migration using the Treemix software [29] (Methods). The Treemix program models populations as having ancestry from multiple parental populations, it can then calculate migration weights [29]. These weights are correlated with the fraction of alleles contributed by a parental population; however, they are not an unbiased predictor of this fraction.

In the Treemix analysis, Turkish samples clustered together with South Europe/Mediterranean populations: Iberians from Spain and Tuscans from Italy (Figure 2B). Within the worldwide population tree, two of the four predicted migration events involved Turkey. The strongest predicted migration event represents admixture from the root of the East Asian branch into Turkey, which could be reflecting Central/South Asian population migration [10]. The second migration pattern represents admixture from the root of the Eurasian branch (close to the African branches) to the common node between Turkey and Italy. The algorithm predicts a parallel branch from Africa into Spain. These patterns plausibly reflect South Mediterranean admixture into North Mediterranean, as observed for other populations in the Mediterranean basin [30].

The weight for the migration event predicted to originate from the branch ancestral to all Eurasians (presumably Middle East and North Africa), to the Turkey-Tuscan clade, was only 0.048. In comparison, the weight for the migration event predicted to originate from the branch ancestral to East Asia (presumably Central Asia) into current-day Turkey was 0.217. Although this implies a major population event from the East to West Asia, we note that these weights are not direct estimates of the migration rates. First, the original contributing populations to the ancestral population in Turkey are not known. For instance, we do not know the exact genetic relationship between current-day East Asian populations and the Turkic speakers from Central Asia who migrated into Anatolia about 1,000 years before present. In fact, Hodoglugil and Mahley, using HGDP genotyping data, predict that South Asian contribution to Turkey's population was significantly higher than East/Central Asian contributions [13], suggesting that the genetic variation of medieval Central Asian populations may be more closely related to South Asian populations, or that there was continued low level migration from South Asia into Anatolia. Another possibility is Ancient North Eurasian genetic contribution to both the historical Anatolian and East Asian populations [28], which might have been interpreted as migration in this dataset. Second, Pickrell and Pritchard [29] also note that in their simulations, the weights underestimate relatively high admixture proportions. Data from more closely related populations coupled with extensive population genetic simulation may eventually allow determining the relative contributions of migration events that shaped population variation in Turkey.

Finally, we investigated possible population structure within Turkey. Principal component analyses did not reveal any subclustering among the 16 individuals, and we found no evidence for longitudinal or latitudinal divergence (Figure 2C). Supporting this lack of structure, we found no correlation between genetic distances among subjects and geographic distances among subject locations (Mantel test p > 0.10). One explanation for this observed pattern is the well-documented high-levels of recent population mobility within this geographic region [31, 32]. However, this pattern might not be exclusive to Turkey; for example, a neighbor joining tree of 98 Tuscan individuals in the 1000 Genomes Project similarly displays a star-like phylogeny (Additional file 3: Figure S2), i.e., we observe long terminal branches coalescing at about the same time, indicating lack of structure within the population.

Genetic variants exhibiting unusual frequency in the Turkish samples

The assessment of whole genome and exome sequencing data at the population level has shown that due to recent population growth, human populations carry an excess of rare genetic variants [17], many of which can be functional [33]. Indeed, we identified 5,523 synonymous or non-synonymous rare –i.e., variants that were found in a single Turkish chromosome but not in any other population - and 7 synonymous or non-synonymous private –found in more than 3 Turkish chromosomes, but not in any other population- single nucleotide variants (SNVs) in the Turkish population (Additional file 4: Table S2). As expected, rare SNVs are significantly more likely to be non-synonymous than SNPs observed in multiple populations (Figure 3). We further scrutinized 3 non-synonymous private alleles. Interestingly, one of these genes, CCDC82, was shown to evolve under positive selection in humans, and diverged from Neandertals [34]. Indeed, the non-synonymous private SNP in the Turkish population observed at chr19:50832152, is a T- > C mutation and homologous to Denisovan haplotype at the orthologous site.

We then explored potentially functional variants at higher or lower frequency in the Turkish population relative to the closely related European populations (Figure 4). For this, we compiled the frequencies of published GWAS SNPs among samples resequenced by the 1000 Genomes Project [2]. We then compared the frequencies of these SNPs in continental populations with what we observed in the Turkish population (Figure 4A). As expected, the frequency distribution of SNPs in the Turkish population is, by and large, similar to that observed among European populations. To identify potential outliers to this expectation, we compared the frequency distributions of the European and Turkish populations and identified 7 SNPs displaying the highest frequency differences between these populations (Figure 4B, 0.1st percentile, corresponding to > ~0.35 absolute difference in frequencies).

Among these 7 SNPs, 2 exist at lower frequencies in Turkey compared to Europe. We found that both are related to pigmentation [36], and hair color [37]. Pigmentation is a genetic and variable trait in humans [38], with lighter skin color being associated with better vitamin D absorption, but higher incidence of skin cancer [39]. As such, human pigmentation correlates with latitude and the low frequency of this SNP in Turkish population as compared to Europe may reflect positive selection for these alleles in Northern Europe and/or selection against the alleles in Southern Europe and Turkey.

To our surprise, 2 of the 5 SNPs that show higher frequency in Turkish population as compared to Europeans are associated with lower total cholesterol counts. This is particularly interesting, given that the Turkish population has indeed been shown to have lower cholesterol than Western European populations, including total cholesterol, and high and low density lipoprotein cholesterol measures, while carrying relatively high triglyceride levels [40]. Although cholesterol levels are strongly shaped by diet and lifestyle, they are also under the influence of genetic factors [41]. The frequency distribution of one of these SNPs (rs7570971), associated with total cholesterol counts [42], is fixed in Asian continent and almost fixed among African populations, decreased to slightly above 90% in Turkish populations, further decreasing to 49% among Europeans. Therefore, similar to what is observed for pigmentation related SNPs, rs7570971 is best explained by clinal adaptation or drift patterns. More striking is the frequency distribution of rs1169288, another SNP associated with total cholesterol count [42], which reaches almost 70% in Turkish populations, but is at lower than 10% frequency in Africa and remains a minor allele in other Eurasian populations (except Japan, where it reaches 54%) (Figure 4C). Intriguingly, however, for both cholesterol level related SNPs, the variant at high frequency in Turkey (A for rs7570971 and C for rs1169288) is reported to elevate total cholesterol levels [42]. Hence, the contribution of these variants to phenotypic differences between populations needs further study.

Deletion polymorphisms

We used the mrFAST aligner [43] together with the VariationHunter [44] algorithm to discover deletions in the genomes of the samples resequenced in this study. We merged individual sample-based callsets, then used the Genome STRiP tool [45] to both genotype and in silico validate our discovery callset (Methods). In summary, we identified 1,751 deletions per individual, with a non-redundant total of 3,292 deletions. To assess novelty, we compared against the 1000 Genomes Project datasets (Methods), and found 199 novel deletions per sample on average (494 total). We did not find any whole-gene deletions, but 21/494 deletions were predicted to delete coding exons where 10/21 were predicted in more than one chromosome (allele frequency >3%; Additional file 5: Table S3).

We also analyzed the genomic structure at the 17q21 locus, a region with increased plasticity [46–49], where eight different haplotypes are characterized [46]. Previous studies have shown that all eight haplotypes show themselves as an inversion polymorphism that exist as two main variants: a direction-orientation haplotype, H1, prevalent in most human populations and an inverted haplotype, H2 which predominantly occurs in European populations [46–48] and presents itself as an 970-kbp inversion [49]. This complex region has been shown to be associated with increased fecundity and an increase in global recombination and it was shown that the H2 haplotype is enriched in Europeans [46, 49]. We have found an allele frequency of 31.25% for 17q21.31 H1/H2 inversion polymorphism in the samples we analyzed, which shows an increase when compared to European populations at 25% allele frequency [48]. One of the samples, sequenced at 34X coverage, was also homozygous for H2 inversion allele (Additional file 1: Table S1 and Additional file 6: Figure S3).

Ethics statement and sample collection

Institutional review board permission was obtained from INAREK (Committee on Ethical Conduct in Studies Involving Human Subjects at the Boğaziçi University). We collected blood samples from 16 volunteers after receiving signed informed consent forms from each of the individuals. The individuals that participated in this study were selected to represent different geographical regions of Turkey.

DNA extraction from human blood sample – NaCl extraction

Approximately 10 ml peripheral blood was collected from each subject into a tube containing K3EDTA to prevent coagulation. Each blood sample was then transferred into a 50 ml centrifuge tube and 30 ml of ice cold red blood cell (RBC) lysis buffer was added. The contents were mixed thoroughly and the mixture was kept at 4°C for at least 20 minutes to lyse the cell membranes. Centrifugation was performed at 5000 rpm at 4°C for 10 minutes. The supernatant containing the RBC debris was removed and the pellet containing the leukocyte nuclei was washed with 3 ml RBC lysis buffer to remove the cell debris. The pellet was suspended in 10 ml of cold RBC lysis buffer by vortexing. After centrifugation at 5000 rpm at 4°C for 10 minutes, the supernatant was discarded, the pellet was cleaned with 3 ml RBC lysis buffer and centrifuged again. 3 ml nuclei lysis buffer was added and the pellet was dissolved by vortexing. In order to digest nuclear proteins, 30 ul Proteinase K (20 mg/ml) and 50 ul of 10% SDS were added, mixed gently and incubated at 37°C overnight or at 56°C for 3 hours. After the incubation step, 10 ml 2.5 M NaCl solution was added to the mixture, shaken thoroughly and centrifuged at 5000 rpm at 20°C for 30 minutes. The supernatant was transferred into a new 50 ml centrifuge tube and approximately 30 ml absolute ethanol was added. The tube was inverted gently until the DNA threads become visible. DNA was fished out with the aid of a micropipette and was left to dry for 2 hours. It was dissolved in 300 ul TE and stored at -20°C. The purity and quality of the extracted DNA were measured at A260/A280, and were found to be within the acceptable range (1.8 - 2.0).

Sequencing

Whole genome shotgun sequence data were generated using the Illumina HiSeq2000 platform (paired-end 101 bp reads clipped to 90 bp) at BGI-Shenzhen. The sequencing reads that contained Illumina adapter sequences or high number of low-quality bases were removed. We finally obtained sequence data coverage between 32- to 48-fold per sample.

SNP and indel discovery

We aligned the reads to the human reference genome (NCBI GRCh37) using the BWA aligner [23] with default options and paired-end mode (sampe). For easier comparison and mitochondrial DNA correspondence, we used the GRCh37 version that is used by the 1000 Genomes Project. Next we used the GATK tool [24] to realign indel-containing reads to the reference genome. We then used GATK UnifiedGenotyper in multi-sample mode to generate SNP and indel callsets separately by pooling the data from 16 genomes. We applied the Variant Quality Score Recalibration filter using the GATK resource bundle version 1.2 to help minimize false positives. We then removed any SNP and indel calls that overlap with segmental duplications to help further reduce false positives. We used the experimental validation data (described below) to approximately calculate the sensitivity and specificity of these call sets.

Validation with SNP microarrays

We hybridized the DNA extracted from 15 of the samples with the Affymetrix 6.0 SNP arrays to test for sensitivity and specificity of the call sets. We compared the Affymetrix 6.0 SNP array results with SNPs we obtained from WGS. Assuming that microarray results provide the gold standard, we calculated false discovery rate to be 0.174% on average and false negative rate to be 0.209%.

Deletion polymorphism discovery and validation

Using the BAM files generated in the previous step, we filtered read pairs that map concordantly to the reference genome within 4 standard deviations around the mean span size. In order to discover all possible mapping locations of the remaining discordant read pairs and unmapped reads in GRCh37, we used the mrFAST aligner [43] in paired-end mode allowing edit distance of at most 4 for both reads. Next, we used the VariationHunter algorithm [44] to call putative deletions (>50 bp) less than 100 Kbp. To calculate the non-redundant set of deletions, we pairwise merged any deletion calls that overlapped at least 80% reciprocally. We then used the GenomeSTRiP tool [45] to genotype the deletions discovered by VariationHunter. We also used the GenomeSTRiP results as a surrogate for in silico validation, and we filtered out any calls that were not successfully genotyped using split-read or read-depth information. We reasoned that any true variation discovered with VariationHunter using read-pair information will be supported by read-depth and split-read information used by GenomeSTRiP to genotype the variant in the genomes of multiple individuals. Although it is still possible that GenomeSTRiP may have false negatives in genotyping, therefore incorrectly invalidating a true variant detected by VariationHunter, this effect is acceptable to obtain a more reliable (lower FDR) call set.

As a second in silico validation strategy, we counted the number of heterozygous SNPs detected by GATK within the novel deletion intervals. We applied this validation to only the novel deletion calls when compared to the 1000 Genomes Project. We reasoned that, if a deletion is real, then it couldn’t contain heterozygous SNPs within, as at least one of the copies is deleted. The only exception to this rule is if the deletion is found within duplication. However, such deletions are too difficult to detect to start with, and filtered out by both VariationHunter and GenomeSTRiP. We found that >85% of our novel deletion calls contain no heterozygous SNPs; where only 10% contains exactly one heterozygous SNP.

We determined novelty in our call set by filtering those deletions that overlap (>50% reciprocal) with any of the 1000 Genomes Project releases [1, 2], or calls discovered within the 1000 Genomes Project but not genotyped. We used the liftOver tool to convert the NCBI Build36 coordinates to GRCh37 for the 1000 Genomes Project Pilot Project results [1].

PCR validation of the 17q21.31 haplotypes

Genotypes of each individual were confirmed by PCR by using the primer specific for Tau deletion region using primers 5’GGAAGACGTTCTCACTGATCTG, and 3’AAGAGTCTGGCTTCAGTCTCTC and subsequent capillary based sequence analysis [49, 50].

Variant annotation

We used ANNOVAR (version 2012Oct23) to annotate the detected variants. For gene based and filter based annotation, we used the April 2012 version of the annotation database (hg19-1000g2012apr_all, hg19_ALL.sites.2012_04.txt), and the dbSNP version 132. For comparison against the 1000 Genomes Project, we used the data released by the consortium as Phase I of the project [2].

Population genetic analyses

To compare Turkish populations with worldwide populations, we used the phase 1 version of the 1000 Genomes Project dataset [2]. Average nucleotide diversity was calculated in each dataset separately, using biallelic autosomal SNPs passing the respective quality thresholds. The sum of π values across all SNPs was divided by sum of ungapped chromosome lengths (GRCh37). Next, to compare populations in the two datasets directly, we combined the datasets. For simplicity, as well as to avoid technical biases, we restricted our analysis to biallelic autosomal SNPs polymorphic in both datasets, totaling 7,134,695. We pruned this dataset to remove SNPs in high linkage using the PLINK software, using the recommended parameters (window size = 50, shift = 5, r^2 = 0.5) [51], which yielded 1,090,112 SNPs. To further limit ascertainment bias, we randomly chose 16 individual subsets (the same number as in the Turkish sample) from each of the 14 populations represented in 1000 Genomes Project, and required a SNP to be polymorphic in each of these 14 population subsets (except for the IBS population, which originally had 14 samples). This step further reduced the number of SNPs to 208,816. Genotypes were coded as follows: 1 for homozygous reference, 2 for heterozygous, and 3 for homozygous alternative. Using this dataset we first conducted principal component (PC) analyses using the ‘prcomp’ function in the R environment, after scaling the matrix to unit variance. To further resolve relationships among populations, we used a recently developed method, Treemix, that allows constructing maximum likelihood phylogenetic trees while allowing migration from edges to nodes [29], providing an efficient alternative to STRUCTURE [52] and similar programs for studying population structure. We removed populations with known large admixture (Native American and African American populations) to focus on admixture events in Turkey. We ran the program with Yorubans as root, grouping SNPs in bins of 100, and allowing for four migrations. We further performed principal component analyses and constructed a neighbor joining tree (using the ‘nj’ function in the R ‘ape’ package) on the 16 individuals from Turkey. We additionally calculated the matrix of genetic distances between individuals and compared these distances with geographic distances calculated from longitude and latitude of subject locations, using the Mantel test.

GWAS analysis

To compare the frequencies of SNPs associated with phenotypes and disease between Turkish population and other world populations, we calculated the frequencies of SNPs listed in the GWAS Catalog [53] among world populations using data from 1000 Genomes Project [2]. Then, we have conducted a pairwise comparison of frequencies of these SNPs observed in Turkish population with populations included in 1000 Genomes Project (Additional file 7: Figure S4). Our initial results, concordant with our population genetics analysis, indicated that the frequency distribution of GWAS SNPs in European populations is closest to that observed in the Turkish populations. As such, we have calculated the distribution of absolute frequency differences between Turkish and European populations and identified the outliers (0.1 percentile) (Figure 4).

Availability of supporting data

Sequence reads are deposited to the SRA read archive (SRP021510). The VCF file that lists all genomic variation characterized in this study is available at: http://turkiyegenomprojesi.boun.edu.tr/tgp_vcf/TGP.integrated_callset.vcf.gz.