Material and Methods

Ethics Statement
This project was excluded from institutional-review-board review by the Office of Human Subjects Research Protections, National Institutes of Health (OHSRP ID# 17-NHGRI-00282).

Sequence Data
We retrieved whole-genome-sequence data from 2,504 individuals in the 1000 Genomes Project,25 320 individuals in the African Genome Variation Project,26 and 108 individuals from Qatar.27 As previously detailed by the 1000 Genomes Project, a three-stage approach was used to establish phased haplotypes: (1) given array-based genotype data and family information, SHAPEIT2 was used to estimate a scaffold of long-range phased haplotypes, (2) a combination of Beagle and SHAPEIT2 was used to jointly analyze the long-range phased haplotypes and di-allelic variants, and (3) MVNcall was used to place multi-allelic and structural variants onto the haplotype scaffold.25 The average switch error rate was 0.56%, and the average distance between phasing errors was 1,062 kb.25 For comparison, the length of the β-globin locus, including the locus control region, is 67 kb. All data were accessed with VCFtools version 0.1.14.28 We assessed linkage disequilibrium in the 1000 Genomes Project data by using the --hap-r2 function for phased haplotypes.

Measuring Multi-locus Association
We calculated Cramér’s V for the association between two nominal variables. In our analyses, the number of levels equaled the number of haplotypes. V2 is the square of the mean canonical correlation and is equivalent to Pearson’s r2 if at least one variable has only two levels. We implemented a bias-corrected version of Cramér’s V given by V˜=φ˜2/m˜, φ˜2=max(0,φ2−((r−1)(c−1)/n−1)), φ2=χ2/n, m˜=min(r˜−1,c˜−1), r˜=r−((r−1)2/n−1), c˜=c−((c−1)2/n−1), where n is the sample size and χ2 is the chi-square statistic without a continuity correction for a contingency table with r rows and c columns.29

African Ancestry
We used YFitter to call Y chromosome haplogroups30 and HaploFind to call mitochondrial DNA haplogroups.31 We used projection analysis in ADMIXTURE version 1.332 with a reference panel of 21 global ancestries33 to analyze autosomal ancestry. To determine standard errors for the proportions of ancestral components for each individual, we reran ADMIXTURE with the addition of 200 bootstrap replicates. Accounting for both within and between individual variances, we calculated the proportions for average ancestry by using inverse variance weights. We then calculated 95% confidence intervals for each ancestry and individual, zeroed out any average proportions for which the 95% confidence intervals included 0, and renormalized the remaining averages to sum to 1.

Balancing Selection
Let the genotype frequencies of the sickle homozygote, heterozygote, and wild-type homozygote be p2, 2pq, and q2, respectively. Let the corresponding relative fitnesses be 0, 1+s, and 1, respectively.2 Then, at equilibrium, s=(p/1−2p). For each of the five continental African samples in the 1000 Genomes Project Phase 3 release version 5a, we estimated the effective population size Ne on the basis of the heterozygosities of all single-nucleotide polymorphisms (i.e., diallelic, triallelic, and quadrallelic) by assuming a mutation rate of 0.97×10−8 mutations per site per generation.34 We then took the harmonic mean of the five Ne estimates. Assuming one initial copy of the mutant allele, we simulated 1,000 generations under a combination of random genetic drift and balancing selection. We repeated this process 1,000 times.

Phylogenetic Network Analysis
Phylogenetic trees are based on models of evolution that assume a purely bifurcating process. This assumption is violated by several evolutionary processes, including recombination, gene conversion, and recurrent mutation. A single phylogenetic tree is unable to represent incompatible signals across sites. In contrast, a phylogenetic network is a more general form of graph that relaxes the assumption of bifurcation. Split decomposition is an inferential method that computes a set of incompatible splits from a given distance matrix. Trivial splits separate a set of taxa into two sets, one set containing a single taxon and the other set containing all other taxa, thereby defining terminal branches. Non-trivial splits separate a set of taxa into two sets, both containing at least two taxa, thereby defining internal branches. A split network depicts incompatible splits as parallel branches. We used the Neighbor-Net method implemented in SplitsTree version 4.13.1 to perform split decomposition analysis of haplotypes.35

Inferring the Ancestral Recombination Graph
We used ARGweaver to infer the ancestral recombination graph.36 ARGweaver is based on the standard coalescent model and is sensitive to balancing selection, such that regions under balancing selection have older times to the most recent common ancestor than comparable neutral regions. To account for uncertainty in both the mutation and recombination rates, we used a grid approach. We tested mutation rates of 0.72×10−8, 0.97×10−8, and 1.44×10−8 mutations per generation per site.34 We tested recombination rates of 5×10−9, 1×10−8, 1.5×10−8 (the default value for human data), 2×10−8, and 2×10−7 recombinations per generation per site.34 We ran the sampler for 3,000 iterations and saved the ancestral recombination graph from every tenth iteration. We used the functions heidel.diag and geweke.diag in the coda library of R, version 3.2.3, to assess convergence diagnostics on the basis of the posterior distribution of the number of recombination events.37 To convert generations into years, we assumed a generation interval of 28 years.38, 39