Insights into the genetic architecture of haematological traits from deep phenotyping and whole-genome sequencing for two Mediterranean isolated populations

Haematological traits are linked to cardiovascular, metabolic, infectious and immune disorders, as well as cancer. Here, we examine the role of genetic variation in shaping haematological traits in two isolated Mediterranean populations. Using whole-genome sequencing data at 22× depth for 1457 individuals from Crete (MANOLIS) and 1617 from the Pomak villages in Greece, we carry out a genome-wide association scan for haematological traits using linear mixed models. We discover novel associations (p < 5 × 10–9) of five rare non-coding variants with alleles conferring effects of 1.44–2.63 units of standard deviation on red and white blood cell count, platelet and red cell distribution width. Moreover, 10.0% of individuals in the Pomak population and 6.8% in MANOLIS carry a pathogenic mutation in the Haemoglobin Subunit Beta (HBB) gene. The mutational spectrum is highly diverse (10 different mutations). The most frequent mutation in MANOLIS is the common Mediterranean variant IVS-I-110 (G>A) (rs35004220). In the Pomak population, c.364C>A (“HbO-Arab”, rs33946267) is most frequent (4.4% allele frequency). We demonstrate effects on haematological and other traits, including bilirubin, cholesterol, and, in MANOLIS, height and gestation age. We find less severe effects on red blood cell traits for HbS, HbO, and IVS-I-6 (T>C) compared to other b+ mutations. Overall, we uncover allelic diversity of HBB in Greek isolated populations and find an important role for additional rare variants outside of HBB.


Results
Novel genetic associations with haematological traits. Following quality control, data were available for 1,457 individuals for MANOLIS and 1,617 individuals for the Pomak population. Five rare, previously unreported loci were associated with haematological traits in either cohort after adjusting for multiple testing (p < 5 × 10 -9 ) ( Table 1).
The G-allele (minor allele frequency (MAF) = 0.009) of rs1320751535 at 15q26 was associated with an increase in red blood cell count by 1.52 units of standard deviation error (SE = 0.23, p = 6.2 × 10 -10 ) in MANOLIS. This single nucleotide variant (SNV) is extremely rare in reference samples. It is not observed in the 1000 Genomes Project data and is carried by 2 out of 125,568 individuals in TopMed 18 . Of note, the credible set in this region, a group of variants that is likely to contain the causal one, included only two markers, neither of which has been reported in previous GWAS for blood traits. In the Pomak population, we identified four novel loci, including an association at 2p11.2 with increased white blood cell count (WBC) (beta = 1.72, SE = 0.29, p value = 4.1 × 10 -9 ). The locus contains 16 highly associated variants in perfect linkage disequilibrium that span a region of 3 Mb (Supplementary Material, Supplementary Table 2). One of the seven variants in the credible set was associated with WBC at nominal significance in a large cosmopolitan GWAS 1 Table 3). There was also a novel association with platelet distribution width (PDW) at 20q13. 32 Table 2). Variants in the credible set, including rs201343203 intronic to FSD1L and rs186868542 intronic to TMEM38B, overlap with enhancer sites that are active in blood. In a previous study, an association with platelet distribution width has been reported for another variant in SVEP1: missense mutation rs6175193 1 , but not with RDW. There was another novel association on chromosome 9. The rare G-allele of a variant at 9q31.3, rs145221983, had a large effect on RDW with a beta of − 2.62 (SE = 0.39, p value = 3.9 × 10 -11 ) in the Pomak population. It is intronic to lysophosphatidic acid receptor 1 (LPAR1). The encoded membrane protein belongs to a group known as endothelial differentiation gene receptors which mediate platelet aggregation. Common variants in LPAR1 have previously been linked to haematological traits 1 .

(Supplementary Material, Supplementary
We replicated several previously-reported associations 1 at genome-wide significance, mostly with platelet traits. Variant rs11553699 at 12q24.31 was associated with mean platelet volume (beta = 0.39, SE = 0.06, p value = 4.9 × 10 -11 ) and platelet distribution width (beta = 0.39, SE = 0.06, p value = 2.3 × 10 -10 ) in MANOLIS, and rs1354034 at 3p14.3 and rs342293 at 7q22.3 with large platelet distribution ratio (beta = 0.23, SE = 0.037, p value = 7.0 × 10 -10 and beta = 0.22, SE = 0.04, p value = 1.4 × 10 -9 , respectively), mean platelet volume (beta = 0.24, SE = 0.04, p value = 1.4 × 10 -10 and beta = 0.23, SE = 0.04, p value = 3.8 × 10 -10 , respectively) and platelet distribution width (beta = 0.23, SE = 0.04, p value = 2.1 × 10 -9 and beta = 0.22, SE = 0.04, p value = 3.9 × 10 -9 , respectively) in the Pomak population (Supplementary Table 4). HBB mutations. Across the entire genome, a region on chromosome 11 displayed the strongest association with all measured red blood cell traits, except haematocrit (HCT) in Pomak and mean corpuscular haemoglobin concentration (MCHC) in MANOLIS (Supplementary Material, Supplementary Fig. 1). Conditional analyses demonstrated that the chr11 peak in the Pomak population could be explained by the three most frequent mutations in the HBB gene, c.364C>A, IVS-II-745 and IVS-I-6 ( Fig. 1). The same was true for MANOLIS, again with three independent signals from HBB mutations, IVS-I-110, CD8/9+G and CD39C>T. We identified all mutations in the HBB gene that have previously been classified as pathogenic, as described in the methods section. In both groups, we observed high proportions of carriers. A total of 99 individuals in MANOLIS (6.8% of the individuals) carried a pathogenic HBB mutation and 162 individuals in the Pomak population (10.0% of the individuals). Different mutation spectra were observed in the two populations (Table 2). There were ten muta- Table 1. Association results of the lead SNPs of novel genome-wide significant loci with haematological traits. Betas are reported in units of standard deviation of the traits. The lead SNP for the locus on chromosome 2 is represented by 16 variants in perfect linkage disequilibrium covering an area of 3 Mb.*A0 is the reference allele and A1 is the effect allele for which the allele frequency (AF) is reported in the current sample and in TopMed.  www.nature.com/scientificreports/ the sickle cell HbS mutation c.20A>T (rs334). CD8/9+G in MANOLIS was located on a haplotype ranging from 0 to 15,000,000.
We tested whether the effects on red cell traits differ between specific HBB mutations. Across the mutational spectrum, there was a clear separation between carriers and non-carriers (Fig. 2), and each mutation was individually associated with red cell traits, such as RDW levels ( Table 2). A case-only analysis demonstrated differences between the structural variants and thalassemia mutations (Supplementary Material, Supplementary Table 5). Moreover, in both Pomak and MANOLIS, IVS-I-6 had significantly less severe effects on red cell distribution width, mean corpuscular volume and mean corpuscular haemoglobin compared to the most common b+ mutation in each group.

Discussion
We present a first detailed characterisation of the mutational spectrum of haemoglobin in two isolated Mediterranean populations, HELIC MANOLIS (Crete) and the HELIC Pomak cohort. A large proportion of the variation in the haematological traits can be explained by HBB mutations in these populations. We provide the first effect estimates for two large non-clinically ascertained samples. We also demonstrate an important role for rare noncoding variation. Finally, we replicate associations of common variants with smaller effects that have previously been reported for cosmopolitan populations with European ancestry. We discover novel associations of variants at 15q26, 2p11.2, 20q13.32, and 9q31.3 with red and white blood cell count, platelet and red cell distribution width, respectively. All of these variants are rare and located outside of coding regions, making it difficult to understand the exact mechanisms through which these variants affect blood traits. Several genes in the 9q31.3 region have been associated with abnormal lipid profiles and coronary artery disease, including SVEP1 21 and ATP Binding Cassette Subfamily A Member 1 gene (ABCA1) 22,23 . Links between haematological traits and lipid profiles have been previously demonstrated 24 . Another variant at the locus with likely regulatory function, rs201343203, is intronic to Fibronectin Type III and SPRY Domain Containing 1 Like (FSD1L), a gene that has been previously linked to red cell distribution width 25 .
For the PDW-associated variants at 20q13.32, results from chromatin interaction experiments 26 further implicate APCDD1L as a likely target gene. An important paralog of this gene (APCDD1) is a negative regulator of the Wnt signaling pathway that is involved in the regulation of platelet function 27 . Chromatin immunoprecipitation sequencing experiments have also shown that the Histone-lysin N-methyltransferase SETDB1 protein binds in the lead SNP in this region. Furthermore, female heterozygous mutant mice have abnormal peripheral blood lymphocytes data 28 . Further investigation is required to fully elucidate the underlying mechanisms of these novel associations.
In each of the two populations studied here, moderate to large effects are observed, with alleles conferring effects of 1.44-2.63 units of standard deviation. This contrasts with other studies. For example, a study based on whole-genome sequencing of 3,781 individuals from a cosmopolitan European population did not discover any novel associations with blood traits 29 . We hypothesize that the observed enrichment of novel haematological association of rare variants is a consequence of population history 30 . Genetic drift due to the founder event in these isolated groups may have resulted in increased allele frequencies for the associated variants, which provides better statistical power for discovery of rare variants with large impact 16 . In fact, most of the lead variants at the novel loci have risen in frequency compared to large reference populations. For example, the rare allele of rs1320751535 at chromosome 15 is only carried by two individuals in TopMed (MAF = 0.000016) but has a frequency of almost 1% in MANOLIS. Limitations of genetic association studies include the possibility of false positive associations. Relative to some array-based genotyping efforts, our sample size was smaller. To fully confirm these novel associations, replication studies would be warranted. However, the low frequency of these variants in other populations and the need for large-scale sequencing to detect them has prevented replication testing using data from other available studies.
Previous reports of high levels of mutational diversity were limited to cosmopolitan populations. We found six different pathogenic HBB mutations in MANOLIS and six in the Pomak group, two of which were seen in both populations. We compared the mutational spectrum observed in the two Greek isolated populations to published data from 3,796 individuals, who were referred to genetic counselling at the NTC from all over Greece 31 . The two most common variants seen in the Greek NTC samples, IVS-I-110G>A (42.1% of carriers) and CD39C>T (18.8% of carriers), were both relatively common in MANOLIS (44% and 13% of carriers, respectively)). However, some  Table 3. Differences in haematological, cardiometabolic and anthropometrics traits between carriers and non-carriers of HBB mutations. width RDW-SD in Pomak, RDW-a in Manolis. We grouped HBB mutations into those that either reduce (b +) or abolish (b0) expression of beta-globin (see Table 2). We also separated c.364C>A (HbO-Arab) in Pomak and sickle cell HbS in MANOLIS. Associations significant after Bonferroni correction (p value < 0.0013) were bolded. www.nature.com/scientificreports/ www.nature.com/scientificreports/ of the mutations from the NTC sample, such as IVS-I-1 (G>A) (12.8% of carriers in NTC) and IVS-II-745C>G (6.3%), were not present in MANOLIS. Conversely, the second most frequent variant in MANOLIS, CD8/9+G, was rare in the NTC sample (0.1% of carriers). It should be noted, however, that recruitment for NTC was for symptomatic cases and relatives which may affect the mutation spectrum. Therefore, the carrier frequencies may not be comparable. There were marked differences in frequency patterns between the Pomak population and NTC sample with the top three mutations from NTC sample, IVS-I-110G>A, CD39C>T, and IVS-I-1, either not present or carried by only one individual in the Pomak population. The most common mutation in Pomak, c.364C>A (HbO), was not observed in the NTC sample. It has been previously postulated that this variant originated in the Pomak population 19 . We provide detailed additional information to characterise the effects of this variant based on data of heterozygous and homozygous carriers. Firstly, we confirm a high allele frequency of 4.4% in the Pomak sample. In line with previous reports, we observe decreased levels of mean corpuscular volume but increases in mean corpuscular haemoglobin concentration 20,32 . This can be explained by the strong positive charge of the HbO molecule which results in their accumulation below the inner surface of the negatively charged erythrocyte membrane, leading to more space that can be filled with haemoglobin as well as denser, more spherical cells 33 . However, in discordance with previous research 20,32 , we do not find statistically significant differences between heterozygotes and homozygotes for these traits. Moreover, we observe novel effects of HbO on white blood cell traits and platelets, including increased platelet distribution width, mean platelet volume, large platelet distribution ratio, as well as white blood cell, lymphocyte and neutrophil count.
The isolated Cretan population has low levels of cardiometabolic complications despite exposure to risk factors such as obesity 15 . In line with previous reports [34][35][36] , we found decreased levels of total and LDL-cholesterol in carriers of thalassemia variants in MANOLIS. This has previously been linked to a decreased risk of atherosclerotic cardiovascular disease [37][38][39][40][41][42] . In MANOLIS we also found an association between the sickle cell HbS mutation with increased height. However, it cannot be ruled out that carrier status tags a specific ancestral group with taller stature. There was also a novel association of b+ mutations with increased risk of being born pre-term.
We show that c.364C>A (HbO-Arab), IVS-II-745 (C>G), CD8/9 + G and IVS-I-110 (G>A) are located on a very large haplotype each, extending over up to 15 Mb, a range that has previously only been reported for the major histocompatibility complex (MHC) region. This pattern likely represents a trace of the natural history of haemoglobin disorders. Positive natural selection can lead to increased linkage disequilibrium 43 . There is strong evidence that heterozygous carrier status of certain haemoglobin mutations provides a protective effect against malaria 9 and Greece is one of the regions with a history of long-standing endemic malaria 10 .
In conclusion, whole-genome sequencing enabled a detailed characterisation of the spectrum of mutations, providing important insights into the allelic architecture of medically-relevant haematological traits. This can provide important guidance for mutation screening in these regions. Future research should extend this work to other populations with a high prevalence of haemoglobin disorders.

Materials and methods
Samples. The HELIC cohorts (https:// www. helmh oltz-muenc hen. de/ itg/ proje cts-and-cohor ts/ helic/ index. html) have previously been described in detail [15][16][17] . Briefly, MANOLIS includes individuals from the mountainous Mylopotamos villages on Crete. Individuals for the Pomak cohort were recruited at the Pomak villages, a set of mountainous villages in the North of Greece. Genetic isolatedness has been demonstrated for both cohorts 16 . A wide range of phenotypic information was collected including anthropometric and biometric measurements, biochemical and haematological blood measures, medical history, demographic, socioeconomic and lifestyle information. All participants provided written informed consent. Ethical approval was obtained from the Harokopio University Bioethics Committee. All methods were performed in accordance with the relevant guidelines and regulations.
Phenotype data. The distribution of each haematological trait was assessed. For those not sufficiently approximating a normal distribution, log-or rank-based inverse normal transformation was applied to phenotype measures and outliers were excluded (Supplementary Material, Supplementary Table 1). Values were adjusted for sex, age and squared age if any of these were significantly (p < 0.05) associated with the trait in a linear regression analysis. Standardised residuals from these analyses were used as outcome for the genome-wide association analysis. For the characterisation of pathogenic mutations in haemoglobin genes, we used un-transformed values as outcomes in the regression analyses to retain interpretability in original units. These analyses were conducted in R v3.4 44 . Sequencing. Whole-genome sequencing was carried out for 1,482 samples from MANOLIS and 1,642 samples from Pomak using Illumina's HiSeqX platform with a target depth of 30x. Processing followed the GATK best practice guideline and has been described in detail elsewhere 45 . Alignment was carried out using BWA mem 0.7.8 using hg38 as the reference. Picard was used to mark duplicates. HaplotypeCaller v.3.5 was used to call genotypes. VQSR was used to for variant quality control (QC) using a tranche threshold of 99.4% for single nucleotide polymorphisms (SNP). For indels, we used the recommended threshold of 1%. We also filtered out 14% of variants with call rates < 99%.
We excluded 25 samples from MANOLIS that failed one or more of the following checks: four samples failed sex checks, eight had low concordance with previous genotyping efforts 17  www.nature.com/scientificreports/ Genome-wide association analyses. The association of genetic variants with each haematological trait was evaluated using a linear mixed model implemented in GEMMA 46 . This approach accounts for relatives in the sample as well as any population substructure. GEMMA was used to estimate the genetic relatedness matrix after filtering for minor allele frequency (MAF) < 0.05, missingness < 1% and linkage disequilibrium (LD)-based pruning. We considered associations of variants with minor allele count of at least 10. To determine the multipletesting burden, we estimated the effective number of traits by carrying out a principal component analysis for the correlation matrix of traits. The first 10 principal components explained 99% of the variants, therefore the effective number of traits was estimated to be 10. We also accounted for number of variants which resulted in an adjusted p-value threshold of 5 × 10 -9 . Measures of linkage disequilibrium, D' and r 2 , were calculated using plink 47 .
For each of the previously unreported variants significantly associated with haematological traits, we considered all variants in a ± 500 Kb distance. In order to identify potentially causal variants, we excluded SNPs with a likelihood of being causal of less than 1:100, by comparing the likelihood of each SNP from the association analysis with the one of the most strongly associated SNP 48 . The remaining variants at each locus, henceforth called credible set, were annotated used FUMA 49 , Ensembl including VEP 50 , HaploReg 51 and Open Targets Genetics 52 to characterise their putative functional impact.

HBB mutations.
We identified all mutations in the Haemoglobin Subunit Beta (HBB) gene that are classified as pathogenic with a review status of at least one star in the ClinVar data base 53 . These were referenced against the HbVar data bank 12 . The HELIC sequence data were queried for these 83 variants.
In addition to single variant association analyses as described above, we also carried out conditional analyses where we included either the most strongly associated variant or the pathogenic HBB variants as covariates in the model.
We used burden testing to evaluate the combined effect of variants in HBB and linked regulatory elements on blood traits. We followed the approach outlined in 45 . Briefly, we used an extended SKAT-O model 54 to account for relatedness or population structure as implemented in MONSTER 55 . Boundaries for HBB were extracted from GENCODE v25. We applied eleven different conditions: regions of interest (coding regions only, coding and regulatory regions and regulatory regions only), variant filters (inclusion criteria based on severity of predicted consequence) and weighting schemes.
Ethics approval and consent to participate. All participants provided written informed consent. Ethical approval was obtained from the Harokopio University Bioethics Committee.