The genetic dissection of fetal haemoglobin persistence in sickle cell disease in Nigeria

Abstract

The clinical severity of sickle cell disease (SCD) is strongly influenced by the level of fetal haemoglobin (HbF) persistent in each patient. Three major HbF loci (BCL11A, HBS1L-MYB, and Xmn1-HBG2) have been reported, but a considerable hidden heritability remains. We conducted a genome-wide association study for HbF levels in 1006 Nigerian patients with SCD (HbSS/HbSβ0), followed by a replication and meta-analysis exercise in four independent SCD cohorts (3,582 patients). To dissect association signals at the major loci, we performed stepwise conditional and haplotype association analyses and included public functional annotation datasets. Association signals were detected for BCL11A (lead SNP rs6706648, β = -0.39, P = 4.96 × 10-34) and HBS1L-MYB (lead SNP rs61028892, β = 0.73, P = 1.18 × 10-9), whereas the variant allele for Xmn1-HBG2 was found to be very rare. In addition, we detected three putative new trait-associated regions. Genetically, dissecting the two major loci BCL11A and HBS1L-MYB, we defined trait-increasing haplotypes (P < 0.0001) containing so far unidentified causal variants. At BCL11A, in addition to a haplotype harbouring the putative functional variant rs1427407-'T', we identified a second haplotype, tagged by the rs7565301-'A' allele, where a yet-to-be-discovered causal DNA variant may reside. Similarly, at HBS1L-MYB, one HbF-increasing haplotype contains the likely functional small indel rs66650371, and a second tagged by rs61028892-'C' is likely to harbour a presently unknown functional allele. Together, variants at BCL11A and HBS1L-MYB SNPs explained 24.1% of the trait variance. Our findings provide a path for further investigation of the causes of variable fetal haemoglobin persistence in sickle cell disease.

Keywords: fetal haemoglobin; genome-wide association study; haplotype analysis; heritability; sickle cell disease.

Graphical abstract

Figure 1

(A) Manhattan plot of -log₁₀(P-values). The upper line indicates the threshold for genome-wide significance (P < 5 × 10⁻⁸), the lower line for suggestive association (P < 1 × 10⁻⁶). (B) Quantile-quantile plot of -log10(P-values) for all SNPs tested. Population structure is implied to be effectively controlled since there is no deviation from the null hypothesis line except at the tail end. The λ_GC value (lambda genomic control inflation factor) of 1.004 indicates that population stratification has been adjusted effectively by the model and relationship matrix. Alt-text. (A). This shows strong genetic association signals for BCL11A (P = 4.96 × 10⁻³⁴) and HBS1L-MYB (P = 1.18 × 10⁻⁹) on chromosomes 2 and 6 respectively. In addition, we identified three novel loci, SLC28A3 on chromosome 9 (P = 2.52 × 10⁻⁸), TICRR/KIF7 on chromosome 15 (P = 3.34 × 10⁻⁸), and PIEZO2 on chromosome 18 (P = 8.04 × 10⁻⁸). (B). The quantile-quantile plot aligns with the null hypothesis line, deviating only at the tail end, representing the associated loci.

Figure 2

Regional association plots for BCL11A (A) and HBS1L-MYB loci (B). The dashed line shows the threshold for genome-wide significance (p < 5 × 10⁻⁸). The black rectangular blocks show key erythroid regulatory elements, detected as DNase I hypersensitive sites for BCL11A [29] and LDB1 binding sites for HBS1L-MYB [33]. Alt-text. (A). The genome-wide significant BCL11A SNPs overlap with the erythroid enhancer elements and DNase I hypersensitive sites (DHSs), denoted as +55 kb, +58 kb, and + 62 kb sites relative to the transcription start sites with their proxy SNPs: rs7606173, rs6706648, and rs1427407 respectively. (B). Similar to the BCL11A SNPs, the HBS1L-MYB intergenic region contains proxy SNPs that fall within the erythroid regulatory elements with rs61028892 and rs66650371 tagging +84 kb DHS.

Figure 3

Major haplotypes at the BCL11A locus are defined by critical markers rs1427407, rs7565301, and rs7606173. A: shows the relationship of these tagging SNPs with other variants populating these haplotypes, i.e. associated alleles (top) and LD relationships (bottom); the effect alleles ‘T’, ‘A’, and ‘G’ for rs1427407, rs7565301, and rs7606173 are highlighted in red . LD plot: Each diamond shows the value of D′, and the LD indicates white (r² = 0), shades of grey (0 < r² < 1), and black (r² = 1). B: shows haplotype effects on HbF levels. Alt-text. (A). Two haplotypes (haplotypes 1 and 3) with strong HbF-increasing effects, as defined by rs1427407, rs7565301, and rs7606173, were detected at the BCL11A locus. Haplotype 1: TGG carries HbF-boosting alleles at rs1427407 (‘T’) and rs7606173 (‘G’) while haplotype 3 carries rs7565301 variant (‘A’ allele) and rs7606173 (‘G’). These two haplotypes showed significant association with HbF compared with haplotype 4, carrying low-HbF variants (P-values < 0.0001). (B). BCL11A haplotypes 1 and 3 showed the highest levels of fetal haemoglobin.

Figure 4

Major haplotypes at the HBS1L-MYB locus are defined by critical markers rs61028892 and rs9399137. A: shows the relationship of these tagging SNPs with other variants populating these haplotypes, i.e. associated alleles (top) and LD relationships (bottom). rs66650371: D and I represent 3-bp deletion and insertion, respectively. The effect allele ‘C’ for rs61028892 and rs9399137 is highlighted in red. LD plot: Each diamond shows the value of D′, and the LD indicates white (r² = 0), shades of grey (0 < r² < 1), and black (r² = 1). B: shows haplotype effects on HbF levels. Alt-text (A). Two major HbF-boosting haplotypes were identified at the HBS1L-MYB locus as defined by rs9399137 and rs61028892: namely, haplotypes 1 (carrying the high-HbF C-allele at rs9399137) and 2 (carrying the high-HbF C-allele at rs61028892). (B). This demonstrates levels of fetal haemoglobin according to the haplotypes identified; haplotypes 1 and 2 show significantly (P < 0.0001) higher levels of HbF compared to the low-HbF-carrying haplotype 3 (GT).

Figure 5

(A) Manhattan plot showing meta-analysis of 3582 SCD patients. The y-axis indicates the −log(P-value) for each variant in the meta-analysis, and the x-axis shows the chromosomal position. (B) Quantile-quantile plot of the observed vs expected P-values with the genomic control inflation factor (λ) of 1.002. Alt-text. (A). This shows genetic association signals for BCL11A on chromosome 2 (P = 1.83 × 10⁻¹⁰⁹), HBS1L-MYB on chromosome 6 (P = 1.77 × 10⁻²⁷), SLC28A3 on chromosome 9 (P = 2.73 × 10^–8) and OR52S1P on chromosome 11 (P = 2.03 × 10⁻¹⁰). (B). This is a quantile-quantile plot of -log10 (P-values) for all SNPs tested, showing the distributions of the observed P-values from the genome-wide association study against the null hypothesis.