Ancestry-Dependent Enrichment of Deleterious Homozygotes in Runs of Homozygosity.

Runs of homozygosity (ROH) are important genomic features that manifest when an individual inherits two haplotypes that are identical by descent. Their length distributions are informative about population history, and their genomic locations are useful for mapping recessive loci contributing to both Mendelian and complex disease risk. We have previously shown that ROH, and especially long ROH that are likely the result of recent parental relatedness, are enriched for homozygous deleterious coding variation in a worldwide sample of outbred individuals. However, the distribution of ROH in admixed populations and their relationship to deleterious homozygous genotypes is understudied. Here we analyze whole-genome sequencing data from 1,441 unrelated individuals from self-identified African American, Puerto Rican, and Mexican American populations. These populations are three-way admixed between European, African, and Native American ancestries and provide an opportunity to study the distribution of deleterious alleles partitioned by local ancestry and ROH. We re-capitulate previous findings that long ROH are enriched for deleterious variation genome-wide. We then partition by local ancestry and show that deleterious homozygotes arise at a higher rate when ROH overlap African ancestry segments than when they overlap European or Native American ancestry segments of the genome. These results suggest that, while ROH on any haplotype background are associated with an inflation of deleterious homozygous variation, African haplotype backgrounds may play a particularly important role in the genetic architecture of complex diseases for admixed individuals, highlighting the need for further study of these populations.

Introduction

Runs of homozygosity (ROH) are long stretches of identical-by-descent (IBD) haplotypes that manifest in individual genomes as the result of recent parental relatedness. Originally conceived to improve the accuracy of homozygosity mapping of recessive Mendelian diseases, ROH have formed the foundation of studies investigating the contribution of recessive deleterious variants to the genetic risk for complex diseases and to the determination of complex traits. Moreover, they have provided unique insights into the demographic and sociocultural processes that have shaped genomic variation patterns in contemporary worldwide human populations,2, 3, 4, 5, 6, 7, 8 ancient hominins,9, 10, 11, 12 non-human primates,13, 14 woolly mammoths, livestock,16, 17, 18, 19, 20, 21 birds,22, 23 felines, and canids.25, 26, 27, 28, 29, 30, 31 Recent population bottlenecks, cultural preferences for endogamy or consanguineous marriage, and natural selection can create increased rates of ROH in individual genomes, substantially increasing overall homozygosity in such populations.

Several studies of the distribution of ROH in ostensibly outbred human populations have shown that ROH are common and range in size from tens of kilobases to several megabases in length.2, 3, 4, 5 Furthermore, total length and prevalence of ROH are correlated with distance from Africa,3, 4, 5 with more and longer ROH manifesting in individuals from populations a longer distance away. These patterns likely reflect increased IBD among haplotypes as a result of the serial bottlenecking process that humans experienced as they migrated out of Africa.

The prevalence of ROH in individual genomes has also been an important factor for understanding the genetic basis of complex phenotypes.32, 33, 34 High levels of ROH have been associated with heart disease,35, 36 cancer,37, 38, 39 blood pressure,40, 41 LDL cholesterol, various mental disorders,42, 43, 44, 45 human height,46, 47 and increased susceptibility to infectious diseases. Indeed, these results are consistent with the idea that many rare alleles of small effect may be the cause of increased risk for complex diseases,49, 50, 51 especially if these mutations are recessive.

We have previously shown that ROH, especially long ROH, are enriched for deleterious homozygous variation.52, 53 Whereas an overall increase in homozygotes is expected with increasing genomic ROH, we have shown that the rate at which deleterious homozygotes accumulate outpaces the rate at which benign homozygotes accumulate52, 53 in long ROH (ROH on the order of several megabases). This is a consequence of young (long) haplotypes containing low-frequency variants getting paired IBD. As low-frequency variants are more likely to be deleterious than common variants, the processes that create very long ROH can also generate unusually high numbers of deleterious homozygotes within these regions.

Although a few studies describing the worldwide distribution of ROH patterns have included a small number of admixed populations,3, 4, 5 the number of individuals per admixed population has been fairly small. Even as the number of admixed individuals continues to grow in the United States, they are still relatively understudied, which translates to disparities in our understanding of population-specific genetic factors that may influence complex phenotypes. Indeed, admixed populations have unique features compared to other populations, in that genomes from these populations are recent combinations of two or more ancestral populations.

This ancestral mosaicism has been exploited to make inferences about the natural history of human populations56, 57, 58, 59, 60, 61, 62, 63 and to search for ancestral haplotypes that influence complex phenotypes.64, 65, 66, 67, 68 Here we add to the body of work on admixed populations by examining the relationship between ROH, local ancestry, and the accumulation of deleterious alleles. We use 1,441 recently published whole-genome sequences distributed roughly equally across three admixed populations in the Americas: African American (n = 475), Mexican American (n = 483), and Puerto Rican (n = 483). Each of these populations is three-way admixed, with distinct contributions from European, Native American, and African ancestral populations.

Among the ancestral populations that contributed haplotypes to these admixed populations, it has been shown that the distribution of deleterious heterozygotes and deleterious homozygotes changes with distance from Africa.70, 71, 72, 73 With this in mind, we propose that accumulation of deleterious homozygotes via increased genomic ROH may also differ within admixed populations based on differing ancestral haplotypes. Indeed, with high deleterious heterozygosity, we propose that African ancestral haplotypes may be most susceptible to large increases in deleterious homozygotes when subjected to harsh bottlenecks or inbreeding, as these low-frequency deleterious alleles will be paired into homozygotes as a result of increased genomic ROH.

Material and Methods

Sample Selection and Quality Control

We used 1,441 whole-genome sequences (dbGaP accession numbers phs000920 and phs000921) from three different admixed populations: African American (n = 475), Mexican American (n = 483), and Puerto Rican (n = 483). These data are an unrelated (up to third-degree relative) set that were previously published by Mak et al., who previously identified and removed third-degree (and closer) relatives and conducted all QC. These genomes all had mean genome coverage >30× with >95% of genome covered at >10× and were called with GATK HaplotypeCaller. Site-level QC was conducted via GATK Variant Quality Score Recalibration, filtering at the 99.8% tranche. Individual genotypes were filtered if they did not have a minimum read depth of 10 and genotype quality of 20. Full details are available in Mak et al.

Calling Local Ancestry

We used 90 African (YRI) individuals and 90 European (CEU) individuals for ancestry references (genotypes obtained from the Axiom Genotype Dataset, see Web Resources) and SNPs with less than 95% call rate were removed. For Native American reference genotypes, we used 71 Native American individuals previously genotyped on the Axiom Genome-Wide LAT 1 array. These samples are unrelated and unadmixed individuals including 14 Zapotec, 2 Mixe, and 11 Mixtec from the southern Mexican state of Oaxaca and 44 Nahua individuals from Central Mexico. Although these individuals are unlikely to exactly match the Native components of all the individuals in our sample, they act as a reasonable proxy for inferring those components, just as our YRI and CEU reference populations act as a reasonable proxy for inferring the African and European components, respectively.

We then subset our 1,441 whole-genome sequences corresponding to sites found on the Axiom Genome-Wide LAT 1 array, leaving 765,321 markers. We then merge these data with our European (CEU), African (YRI), and Native American (NAM) reference panels, which overlapped at 434,145 markers. After filtering multi-allelic SNPs and SNPs with >10% missing data, we obtained a final merged dataset of 428,644 markers. We phased this combined dataset using SHAPEIT2 and called local ancestry tracts jointly with RFMix under a three-way admixture model based on the African, European, and Native American reference genotypes described above.

Calling Runs of Homozygosity

We called runs of homozygosity using the program GARLIC v.1.1.4, which implements the ROH calling pipeline of Pemberton et al. for each population separately on the full whole-genome call set, filtering only monomorphic sites. For the 475 African American (AA) individuals, this left 39,517,679 segregating sites; for the 483 Puerto Rican (PR) individuals, this left 31,961,900 segregating sites; and for the 483 Mexican American (MX) individuals, this left 30,744,389 segregating sites. Instead of asserting a single constant genotyping error rate (as in Pemberton et al.), we used genotype quality scores provided with the WGS data to give GARLIC a per-genotype estimation of error. Using GARLIC’s rule of thumb parameter estimation, we chose analysis window sizes of 290 SNPs, 250 SNPs, and 210 SNPs and overlap fractions of 0.3688, 0.3553, and 0.3528 for the AA, PR, and MX populations, respectively. GARLIC chose LOD score cutoffs of −47.5169, −70.1977, and −60.9221 for the AA, PR, and MX populations, respectively. Using a three-component Gaussian mixture model, GARLIC determined three size classes: small class A, medium class B, and long class C ROH. Class A/B and class B/C size boundaries were inferred as 38,389 bps and 142,925 bps for AA; as 50,618 bps and 230,079 bps for PR; and 46,979 bps and 217,054 bps for MX.

Computing Ancestry Enrichment in ROH

To determine whether the ROH covering a gene region is overrepresented for a particular ancestry, we first compute, for each gene region, the quantities and , which represent the mean proportion of ancestry i in ROH at gene region R and the “number” of ROH in each gene region, respectively. Note that if an ROH only covers part of a gene region, then only that fraction is counted, thus N is continuous and not a whole number. We also compute the mean proportion of ancestry i in the population, A. If we consider the fraction of ancestry type i in ROH () as a random sample from the distribution of ancestry in the population , then we can model the ancestry-specific ROH sampling process with a beta distribution. This is conceptually similar to a binomial sampling process, where sampling ancestry i in an ROH is considered a “success” but in continuous space. Here we wish to compute the probability of sampling ROH regions of ancestry i (or more) given that the population admixture fraction of ancestry i is A and that we have N ROH total. We can do this by computing , where is the regularized incomplete beta function.

Calling Deleterious Alleles

Using the Whole Genome Sequencing Annotation (WGSA) pipeline to generate annotation data, we extracted PolyPhen 2, SIFT, Provean, and GERP scores for deleteriousness, as well as high-confidence ancestral allele states (from Enredo-Pecan-Ortheus alignments) and synonymous annotations and for all mutations in coding regions (WGSA pre-computed annotations available online, see Web Resources).

PolyPhen 2 generates three deleteriousness categories: Probably Damaging, Possibly Damaging, and Benign. If a mutation has more than one PolyPhen2 classification (e.g., Benign and Probably Damaging), it is reassigned to have only the most damaging category of the group. All mutations that have a PolyPhen 2 prediction or that are synonymous are then pooled into two separate categories: “damaging” and “benign.” All Probably Damaging or Possibly Damaging mutations are pooled into the “damaging” category, and all Benign and synonymous mutations are pooled into the “benign” category.

SIFT generates two deleteriousness categories, Intolerant and Tolerant, which we relabel “damaging” and “benign.” If a mutation has more than one SIFT classification, it is reassigned to have only the most damaging category of the group.

Provean generates two deleteriousness categories, Deleterious and Neutral, which we relabel “damaging” and “benign.” If a mutation has more than one Provean classification, it is reassigned to have only the most damaging category of the group.

GERP generates a numerical score at a given locus where a higher score indicates more deleteriousness for a derived allele at that locus. Here we focus on derived alleles that are very likely to be deleterious and combine all derived mutations at sites with GERP ≥ 6 into the category “damaging.” We form our “benign” category with all derived mutations with GERP ≤ 2.

Defining Gene Sets

We sought to define three sets of genes for further analysis based on the probability of intolerance to loss of function (pLI) predicted as part of the gnomAD project (Web Resources). This score ranges from 0 to 1, with high scores suggesting an intolerance to inactivation and low scores suggesting a tolerance for inactivation. The distribution of these scores is bimodal, with most genes having a pLI near 0 or 1. Of the 18,451 autosomal genes with a pLI score, we create a “low-pLI” category consisting of 13,128 genes with a and a “high-pLI” category consisting of 3,241 genes with a . We finally create an “all” category consisting of all 18,451 autosomal genes reported as part of the gnomAD project.

Computing Minor Allele Frequencies

In order to determine minor allele frequency (MAF) category, we use frequencies computed from all TOPMed Freeze 3 whole-genome sequencing datasets (dbGaP accession numbers phs000920, phs000921, phs001062, phs001032, phs000997, phs000993, phs001189, phs001211, phs001040, phs001024, phs000974, phs000956, phs000951, phs000946, phs000988, phs000964, phs000972, phs000954, and phs001143) forming a total sample size of n = 18,581. We then categorize variants in the dataset analyzed here as common (MAF ≥ 0.05) and rare (MAF < 0.05) based on these “global” allele frequencies.

Simulations

We perform simulations to examine how demographic history affects the concentration of deleterious homozygotes in ROH. We use the forward simulation program SLiM 387, 88 to simulate deleterious mutations within a complex demography in conjunction with the coalescent simulator msprime to simulate neutral mutations conditional on the forward simulation genealogy. This allows us to efficiently simulate very large genomic regions, which is a requirement for analyzing the distribution of long ROH that typically extend several megabases. We complete 500 replicates of the following simulations.

We simulate a three-population demographic history after Gravel et al. in SLiM 3, introducing recessive mildly deleterious alleles with selection coefficients drawn from . We simulate a 100 Mbps region, where deleterious alleles are allowed to occur in designated “coding regions.” These regions are defined based on the hg19 exon coordinates of all CCDS genes in the first 100 Mbps of human chromosome 1. Similarly, we simulate a variable recombination rate based on the HapMap phase II inferred map. We allow a mutation rate based on the Gravel et al. inferred mutation rate of 2.36 × 10−8, setting the deleterious mutation rate at one-tenth of this value. At the end of the forward simulation, a list of segregating deleterious mutations and their genomic locations is output along with the full tree sequence88, 89 of the entire simulation. Neutral mutations are then added with msprime.

To simulate neutral mutations conditional on the forward simulation history, we load our population tree sequence with the pyslim package, recapitate to ensure all lineages fully coalesce, and then lay down neutral mutations at a rate of 90% of the Gravel-inferred rate (so that the neutral plus deleterious mutation rate equals the inferred rate). Finally, we sample 500 diploid individuals from each population in the simulation for analysis.

Simulation code is available online (see Web Resources).

Results

Admixture

Using the subset of sites from our whole-genome sequencing data that intersected with our African, European, and Native American reference panels, we called 3-way local ancestry tracts in all 1,441 samples (see Material and Methods). We also estimated global ancestry proportions by summing the length of all haplotypes inferred to be from a given ancestry and dividing by the total genome length. Figure 1 summarizes the global ancestry proportions for all individuals from each population on a ternary plot. The admixture proportions largely accord with previous results in these populations, with Puerto Ricans having mostly African and European ancestry, Mexican Americans having mostly European and Native American ancestry, and African Americans having mostly African and European ancestry to the near exclusion of any Native American ancestry. However, although African Americans are frequently treated as a 2-way admixed population between European and African sources, we show that several AA individuals have non-trivial proportions of Native American ancestry. This suggests that, in general, a 2-way admixture model may not be uniformly appropriate for studying admixture patterns among self-identified African American individuals.

Figure 1

A Ternary Plot of Global Ancestry Proportions

Each point represents a single individual, with their global ancestry proportions shown on each of the three axes (European, EUR; African, AFR; and Native American, NAM). Individuals are colored based on their reported ethnicity, with African Americans (AA) colored gray, Puerto Ricans (PR) colored purple, and Mexican Americans (MX) colored green.

Runs of Homozygosity

We followed the ROH calling pipeline of Pemberton et al. as implemented in the software GARLIC to call ROH from the full whole-genome sequencing data (see Material and Methods). This method identifies three classes of ROH based on the length distribution in each population. We refer to these size classes as short, medium, and long. These classes roughly correspond to ROH formed of IBD haplotypes from different time periods from the population history. Short ROH are tens of kilobases in length and likely reflect the homozygosity of old haplotypes; medium ROH are hundreds of kilobases in length and likely reflect background relatedness in the population; and long ROH are hundreds of kilobases to several megabases in length and are likely the result of recent parental relatedness. Total length of ROH in the genome is correlated with distance from Africa.2, 4 In the case of our admixed populations, we therefore expect the total length of ROH to be correlated with increased European and Native American admixture fraction. Figure 2A illustrates this pattern, with AA individuals having lowest total ROH, PR individuals having intermediate total ROH, and MX individuals having the highest total ROH (all pairwise Mann-Whitney U tests p < 2.2 × 10−16). Indeed, if we do multiple regression of total ROH coverage (in Mbps) onto total European and total Native American coverage (in Mbps), we find a significant positive association with both ancestry backgrounds in all three populations (Table S9). Breaking down ROH by size class, we find that the total length of short ROH is similar but still significantly higher in PR than in MX individuals (p < 2.2 × 10−16; Figure 2B), but the total length of both medium ROH (p < 2.2 × 10−16; Figure 2C) and total long ROH (p < 2.2 × 10−16; Figure 2D) is highest on average in MX individuals.

Figure 2

The Distribution of Summed ROH Lengths across Size Classes

Shown are (A) all ROH, (B) short ROH, (C) medium ROH, and (D) long ROH. AA, African American; PR, Puerto Rican; MX, Mexican American.

As it has been previously noted that ROH do not occur uniformly across the genome,4, 5 we also examined the proportion of ROH coverage of each of 18,451 coding genes from the gnomAD project across all individuals in each population to discover whether certain genes or sets of genes were enriched for ROH coverage. For each gene region (exons plus introns), we compute the fraction of basepairs that are covered by ROH in each individual and take the mean of this fraction across individuals. Next, we look at the top 0.1% of genes with the highest overall ROH coverage across individuals in each population (Table S8). This corresponds to genes with greater than 0.661, 0.891, and 0.971 ROH coverage across individuals in the African American, Puerto Rican, and Mexican populations, respectively. Although none of these gene sets were enriched for any gene ontology terms, four gene regions were found in all populations: CCDC189, PDCD7, PHKG2, and TMEM139.

We also examine whether certain gene sets may have more enrichment for ROH than others. In particular we create two gene sets based on the gnomAD project’s predicted intolerance to loss of function (pLI) measurement (see Material and Methods). The high-pLI gene set consists of 3,241 genes predicted to be most intolerant to loss of function in humans, and the low-pLI gene set consists of 13,128 genes predicted to be least intolerant to loss of function in humans.

Table 1 lists the means and ranges for ROH coverage across individuals for both high-pLI and low-pLI gene sets. Although the ranges tend to span most of the [0,1] interval, we do observe a small but significant difference in the mean ROH coverage between high-pLI and low-pLI gene sets (as tested by a two-sided Mann-Whitney U test) across all populations, with high-pLI genes having slightly more ROH on average. This may be a result of high-pLI genes experiencing stronger background selection, as high-pLI genes are intolerant to loss of function in humans and mutations in these genes may therefore be more deleterious on average. This, in turn, may contribute a non-trivial amount of homozygosity to the patterns of ROH we observe.

Table 1

Range and Mean ROH Coverage of High-pLI and Low-pLI Gene Sets by Population

Population	High-pLI Genes		Low-pLI Genes		Difference of Means (p value)
	Range	Mean	Range	Mean
AA	[0.013,0.699]	0.195	[0,0.818]	0.181	∗∗∗<2×10^−16
PR	[0.023,0.914]	0.346	[0,0.974]	0.329	∗∗∗1.196×10^−9
MX	[0.019,0.977]	0.428	[0,0.992]	0.414	∗∗∗1.586×10^−5

p value for difference of means computed by two-sided Mann-Whitney U test. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

We also tested whether ROH in certain gene regions are overrepresented with one ancestry background relative to the distribution of ancestries at that gene region population-wide. We compute the probability of observing as much or more of each ancestry among the set of ROH at a gene region for all populations (see Material and Methods) for each 18,451 gene regions from the gnomAD project. Significance was determined via Bonferroni correction, and we find numerous genes in each population enriched for various ancestries (Tables S10, S11–S14, S15, S16, and S17). Each population had at least one gene enriched for each ancestry, except African Americans, where we found no genes enriched for Native American ancestry (though the proportion of Native American ancestry in this population is low, ∼2%, so power may be limited).

We conduct a gene ontology (GO) enrichment analysis using PantherDB, as some of the enrichment lists were large. We find among genes enriched for African ancestry in Mexican Americans significant enrichment of GO terms related to nucleosome assembly (FDR = 2.13 × 10−5), cellular response to unfolded protein (FDR = 6.95 × 10−3), and cellular response to heat (FDR = 1.22 × 10−2). Among genes enriched for Native American ancestry in Mexican Americans we find significant enrichment of GO terms related to spindle assembly (FDR = 1.27 × 10−2) and detection of chemical stimulus involved in sensory perception (FDR = 3.59 × 10−3). Finally, we also find among genes enriched for African ancestry in Puerto Ricans significant enrichment of GO terms related to cytokine-mediated signaling pathways (FDR = 2.27 × 10−2).

Deleterious Alleles

We used multiple approaches to predict the deleteriousness of all sites in the genome (see Material and Methods), but focus on missense mutations classified as Probably Damaging, Possibly Damaging, or Benign using PolyPhen 2. As in Szpiech et al., we combine the Probably Damaging and Possibly Damaging mutations into a single “damaging” class, and we combine all Benign mutations with synonymous mutations into a single “benign” class. For individual i across all sites, we denote by and the total number of sites with alternate alleles classified as damaging or benign, respectively. In Figure 3A we plot the distribution of deleterious heterozygotes per individual, , split by population. Consistent with previous work,70, 71, 72, 73 we see an increased number of deleterious heterozygotes in populations with more African ancestry, with AA individuals having the most and MX individuals having the fewest (patterns replicate with other deleterious categories, see Figures S5–S10). Conversely, we would expect an increase of deleterious homozygotes per individual in populations with more non-African ancestry. Indeed, in Figure 3B we plot the distribution of deleterious homozygotes per individual, , split by population and observe AA individuals with the fewest and MX individuals having the most (these patterns also replicate with other deleterious categories, see Figures S5–S10). Figure 3C plots the total number of deleterious alleles per individual (). Contrary to other work, we find a total deleterious load highest on average in AA individuals. Although this pattern replicates across several other deleterious calling methods (Figures S5–S9), when using GERP scores (as in Henn et al.), the pattern reverses (Figure S10) and is consistent with Henn et al.

Figure 3

The Distribution of Deleterious Alleles across Populations

The number of (A) deleterious heterozygotes, (B) deleterious homozygotes, and (C) total deleterious alleles per individual using PolyPhen2 classifications. AA, African American; PR, Puerto Rican; MX, Mexican American.

Deleterious Alleles across Local Ancestry

We next investigate whether there are any differences in deleterious load by local ancestry. Although our local ancestry calls provide us with phased local ancestry inferences, we were limited to a small subset of sites for our reference populations. Since the vast majority of our deleterious alleles come from our unphased whole-genome data, we do not have phase information for the deleterious alleles and cannot assign a specific ancestral haplotype in regions of discordant ancestry. Therefore, we calculate total load based on six different ancestry backgrounds. AFR, EUR, and NAM ancestry regions represent regions that are homozygous for African, European, and Native American ancestries, respectively, and AFEU, EUNA, and AFNA ancestry regions represent regions that are called heterozygous for African/European, European/Native American, and African/Native American ancestries, respectively. We then calculate for each population the number of deleterious alleles per basepair for each ancestry background.

Table 2 shows the number of deleterious alleles per basepair for each population and each ancestry background using PolyPhen 2 deleterious calls (results were qualitatively similar across all other deleterious call sets). We perform two types of tests for independence in order to determine whether there are significant differences in the number of deleterious alleles per basepair. First, we test for independence of the count of deleterious alleles on an ancestry background and the count of basepairs covered by that ancestry across populations. We find that neither African ancestry nor European ancestry have statistical differences in the number of deleterious alleles per MB across populations. Further, while NAM, EUAF, and AFNA exhibit statistically differences across populations, it appears to be driven by one of the two populations (AA, MX, and PR, respectively). Next, we test for independence of these counts across ancestries within each population. Here we find that all populations have statistically significant differences in the distribution of deleterious alleles across ancestry backgrounds (AA p < 2.2 × 10−16; MX p < 2.2 × 10−16; PR p < 2.2 × 10−16), with NAM ancestry having the lowest rate in AA and PR individuals and EUR having the lowest rate in MX individuals. However, we note that the overall differences were very small (a difference of <0.1 deleterious alleles per Mbp).

Table 2

The Number of Deleterious Alleles per Megabase Partitioned by Population and Local Ancestry Background

	AFR (p=0.160)	EUR(p=0.452)	NAM∗∗∗ (p=3.314×10−7)	EUAF∗∗ (p1.131×10−3)	EUNA (p=0.123)	AFNA∗∗ (p=4.392×10−3)
AA∗∗∗ (p <2×10^−16)	0.335 (1.642×10^6)	0.284 (1.009×10^5)	0.237 (8.648×10^2)	0.311 (7.943×10^5)	0.280 (2.491×10^4)	0.315 (8.364×10^4)
PR∗∗∗ (p <2×10^−16)	0.337 (1.603×10^5)	0.282 (1.064×10^6)	0.275(5.395×10^4)	0.313 (7.517×10^5)	0.286 (4.912×10^5)	0.308 (1.700×10^5)
MX∗∗∗ (p <2×10^−16)	0.341 (7.651×10^3)	0.282 (4.585×10^5)	0.286 (8.275×10^5)	0.317 (1.154×10^5)	0.287 (1.142×10^6)	0.314 (1.393×10^5)

Total number of megabases, summed across all individuals, in parentheses. A significant difference (Pearson’s chi-square test, p value in parentheses) across populations for a given ancestry background is denoted at the beginning of a column. A significant difference across ancestry backgrounds for a given population (Pearson’s chi-square test, p value in parentheses) is denoted at the beginning of a row. Population codes: AA, African American; PR, Puerto Rican; MX, Mexican American. Local ancestry codes: AFR, homozygous African; EUR, homozygous European; NAM, homozygous Native American; EUAF, heterozygous European/African; EUNA, heterozygous European/Native American; AFNA, heterozygous African/Native American. ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001.

Deleterious Alleles in ROH

Next, we turn to examining the distribution of deleterious homozygotes within ROH. It was previously reported52, 53 that there is a higher proportion of deleterious homozygotes per unit increase of ROH than expected from the proportion of benign homozygotes. Naturally, as the total amount of genomic ROH increases, we expect more homozygotes to fall within ROH. However, Szpiech et al. and Pemberton and Szpiech found that the rate of increase of the proportion of deleterious homozygotes was greater than for benign homozygotes. This effect was strongest for long ROH, which are likely the result of recent parental relatedness.

For each individual i and for each ROH class (A, short ROH; B, medium ROH; C, long ROH; R, all ROH; and N, outside ROH), we define the number of damaging or benign sites with alternate alleles as and , respectively. Thus, we calculate the proportion of damaging homozygotes in ROH class j asand the proportion of benign homozygotes in ROH asrespectively. We also compute, for each individual i and each class j, the fraction of the genome covered in ROH as

We plot the proportions of ROH homozygotes versus genomic fraction of ROH in Figure 4, which is analogous to Figure 4 from Szpiech et al. In order to determine whether there is a statistically significant difference in the accumulation of deleterious homozygotes versus benign homozygotes, we construct a linear regression model (as in Szpiech et al. and Pemberton and Szpiech), , where is a vector of length 2,882 containing the proportions of both damaging and benign homozygotes in ROH class j for all individuals, is a vector of genomic class j ROH proportions, and D is an indicator variable taking a value of 1 when the response represents damaging homozygotes and 0 for benign homozygotes. In this framework, a statistically significant suggests an overall higher proportion of damaging homozygotes in ROH compared to benign homozygotes, e.g., means that an extra 10% of genome-wide deleterious homozygotes fall in ROH compared to the distribution of benign homozygotes. A statistically significant suggests a difference in the rate of accumulation per unit increase of ROH, e.g., means that for a 10% increase in genomic ROH, 10% more deleterious homozygotes fall in ROH compared to benign homozygotes. Inferred coefficients for the four regressions corresponding to each are given in Table S1.

Figure 4

Deleterious and Benign Homozygotes in ROH Classes

The proportion of damaging (red) and benign (blue) homozygotes falling in ROH of different size classes: (A) all ROH, (B) short ROH, (C) medium ROH, and (D) long ROH. Data shown is across all populations. Gray line plots Y = X.

Figure 4A plots these proportions versus total ROH for all ROH classes combined. In agreement with Szpiech et al., we find that there is an overall greater proportion of damaging homozygotes in ROH compared to benign homozygotes (, p < 2 × 10−16), but in contrast the overall rate of accumulation is not different (, p = 0.0671). When we partition ROH by size class, the distribution of homozygotes in short ROH (Figure 4B) also differs from Szpiech et al. Whereas previously there were no statistically significant differences in or , here we find a significant positive (p < 2 × 10−16) and a statistically significant negative (p < 1.10 × 10−8), suggesting that ROH comprised of old haplotypes accumulate deleterious homozygotes at a slower rate that benign homozygotes. As we expect short ROH to be comprised of old haplotypes that have been segregating for a long time, it is reasonable to think that only haplotypes with relatively few deleterious alleles remain segregating in the population. Our results for medium (Figure 4C) and long ROH (Figure 4D) are consistent with previous work;52, 53 in particular we find that the difference in rates of gain of deleterious versus benign homozygotes is greatest in long ROH (; p < 2 × 10−16).

We also consider whether we can detect a difference in concentration of deleterious homozygotes in our high-pLI and low-pLI gene sets. For this analysis we only consider predicted deleterious homozygotes, and we wish to compare the genome-wide proportion of these genotypes between high-pLI and low-pLI genes. To do this we construct the following linear regression, , where and are as above and is an indicator variable taking a value of 1 or 0 if the response comes from the high-pLI gene set or the low-pLI gene set, respectively (Table S4). Here represents the difference in rate of accumulation of deleterious homozygotes in high-pLI genes versus low-pLI genes. We find a significant difference in the accumulation of deleterious homozygotes in high-pLI genes versus low-pLI genes for total ROH (, ) and short ROH (, ), although not for long ROH (, ) or medium (, ). In this analysis we compare damaging alleles across two gene sets (instead of comparing damaging to non-damaging), where we might expect mutations in loss-of-function intolerant genes (high-pLI) to be more deleterious compared to mutations in loss-of-function tolerant genes (low-pLI). In this case, the effect size may be much smaller, and by restricting our high-pLI gene set to such a small number of genes we may lack power to detect it. However, in aggregate these results suggest that a higher proportion of genome-wide deleterious homozygotes fall within high-pLI genes versus low-pLI genes.

Deleterious Alleles in ROH Partitioned by Local Ancestry

Now we turn to analyzing the distribution of deleterious homozygotes in ROH comprised of only one particular ancestral haplotype. As shown in Figure 3A and in other work,70, 71, 72, 73 populations with more African ancestry tend to have high numbers of deleterious heterozygotes genome-wide. This contrasts with populations that have more European and Native American ancestry, which tend to have more genome-wide deleterious homozygotes (Figure 3B) as a result of the serial bottlenecks they experienced since migrating out of Africa.

We have already shown (Figure 4) that as total genomic ROH increases the proportion of deleterious homozygotes falling in ROH increases faster than the proportion of benign homozygotes, but here we want to know whether the ancestral background of the IBD haplotypes matters. We propose that haplotypes sourced from ancestral populations with high deleterious heterozygosity have highest rates of accumulation of deleterious homozygotes when paired IBD to generate ROH.

To test this proposition, we first partition ROH based on the ancestral background of the underlying IBD haplotypes. Then we compute for each individual (i) the fraction of all deleterious (d) and benign (b) homozygotes across the genome that fall into each ROH class (j) as:andwhere and are the number of deleterious and benign homozygotes, respectively, in individual i in ROH class j on ancestral haplotype background . Similarly, and are the genome-wide fraction of deleterious and benign homozygotes, respectively, in individual i in ROH class j that fall on haplotype background A. Finally, we fit a linear model similar as above,, in order to test for differences in the rate of accumulation of deleterious homozygotes compared to benign homozygotes as a function of , the genomic fraction of ROH on ancestral background A. The results are plotted in Figure 5 for total ROH (; Figures 5A–5C) and for long ROH (; Figures 5D–5F), and the regression coefficients are also summarized in Table S2.

Figure 5

Deleterious and Benign Homozygotes in ROH Classes Separated by Ancestry

The proportion of damaging (red) and benign (blue) homozygotes falling in ROH comprised of different ancestral haplotypes and size classes: (A) all NAM ROH, (B) all EUR ROH, (C) all AFR ROH, (D) long NAM ROH, (E) long EUR ROH, and (F) long AFR ROH. EUR, European; AFR, African, and NAM, Native American. Gray line plots Y = X.

For total ROH, we find significant differences in the rate of accumulation of deleterious homozygotes on all ancestry backgrounds (Figures 5A–5C). Furthermore, consistent with our expectations, we find that ROH on African ancestral haplotypes have the highest rate difference (, p < 2 × 10−16; Figure 5C), whereas ROH on European ancestral haplotypes have an intermediate rate difference (, p < 2 × 10−16; Figure 5B) and ROH on Native American ancestral haplotypes have the lowest rate difference (, p < 2 × 10−16; Figure 5A). This pattern is repeated when we consider only long ROH comprised of young haplotypes (Figures 5D–5F) and also when we analyze smaller ROH (albeit with weaker effects; Figure S1).

We also perform a variation of this analysis to compare the rate of gain of deleterious homozygotes in high-pLI versus low-pLI genes in ROH across different ancestral backgrounds. We fit the regression , which is similar to above except that is an indicator variable taking a value of 1 or 0 if the response comes from the high-pLI gene set or the low-pLI gene set, respectively (Table S5).

For all ROH combined, we find a significantly higher rate of gain of deleterious homozygotes in high-pLI genes versus low-pLI genes on Native American haplotypes (,) but not for European (,) or African (,) haplotypes. Considering only long ROH, there is a significant difference for Native American (,) and European (,), but again not for African (,). Since we have restricted our dataset by gene set, ROH class, and ancestral background, we may lack power to detect small effect sizes in this African case. Alternatively, there may be more complicated dynamics relating deleteriousness to demography and inbreeding.

We next directly compare the rate of increase of deleterious homozygotes across different ancestral haplotype backgrounds. To do this we compute the following regression, , where is a vector representing the proportion of damaging homozygotes in ROH class j on each local ancestry background across all individuals. represents the genome-wide fraction ROH class j falling on each local ancestry background across all individuals, and is an indicator variable which takes the value 1 if the associated response is on ancestral background and takes the value 0 otherwise. Here we analyze each ROH class: all, long, medium, and short.

We plot the results for “all” and “long” in Figure 6 (“medium” and “short” in Figure S2) and summarize the inferred regression coefficients for all classes in Table S3. We focus on the regression coefficients and , which represent the difference in rate of gain of deleterious homozygotes in ROH on European or Native American haplotypes compared to African haplotypes, respectively. Graphically, in Figures 6 and S2, a significant corresponds to a significant difference in the slope of the orange and blue line, and a significant corresponds to a significant difference in the slope of the orange and red line. Since we expect that the rate of gain of deleterious homozygotes to be lowest in ROH on European and Native American haplotypes compared to ROH on African ones, we expect significant negative values for both and .

Figure 6

Deleterious Homozygotes in ROH Classes Compared across Ancestry

A direct comparison of the proportion of damaging homozygotes falling in ROH comprised of different ancestral haplotypes for (A) all ROH and (B) long ROH. EUR, European, colored blue; AFR, African, colored orange; and NAM, Native American, colored red. Gray line plots Y = X.

Consistent with our expectations, when analyzing all ROH (Figure 6A) we find a significant negative () and (), indicating that the gain rate of damaging homozygotes in ROH on African ancestral haplotypes outpaces that of ROH on the other ancestral haplotypes. This pattern continues when considering only long ROH (, ; , ; Figure 6B) and smaller ROH (Table S3 and Figure S2).

We repeat a similar analysis to compare the rate of gain of deleterious homozygotes in high-pLI genes directly across ancestry backgrounds. In this case, although African ancestral backgrounds do not show a significant difference in the accumulation of deleterious homozygotes between high- and low-pLI genes, they show a clearly higher rate of gain in high-pLI genes compared to European and Native American ancestral backgrounds (Table S6).

To check the robustness of these results, we reran these analyses using several other deleterious classification methods including SIFT,82, 93 Provean, and GERP. Since GERP scores sites and not mutations, we restricted the GERP analysis to loci where the ancestral and derived states were inferred to high confidence. As this ancestral polarization results in discarding a large number of loci with ambiguous ancestral allele state, we also reran these analyses for PolyPhen 2, SIFT,82, 93 and Provean restricted only to loci for which we have ancestral/derived state information. Figure S3 plots the inferred for each of these analyses for each ROH size class and demonstrates qualitatively similar patterns as shown above.

We further re-analyzed a subset of the ROH and deleteriousness calls from Pemberton and Szpiech, which contains data on six admixed populations from the 1000 Genomes Project and used CADD scores as a deleteriousness prediction (Supplemental Material and Methods). After extracting the data relating to the admixed individuals from Pemberton and Szpiech and calling local ancestries, we again find qualitatively similar patterns as above (Figure S4).

Since Pemberton and Szpiech showed that these enrichment patterns appear to be driven by an abundance of homozygotes in ROH comprised of low-frequency alleles, we re-analyzed our data using categories of minor allele frequency (MAF) instead of deleteriousness (see Material and Methods for how we determined MAF category). Using these allele frequencies, we categorize each polymorphic locus in a gene region (exons plus introns) into one of two categories: common (MAF ≥ 0.05) and rare (MAF < 0.05). We then fit the same models as above, except that instead of comparing the proportion of deleterious alternate allele homozygotes to benign homozygotes as a function of ROH coverage, we compare the number of minor allele homozygotes in the rare class to the common class.

We summarize the results of these analyses for each ancestral background, each ROH size class, and each low-frequency class in Figure 7. We find that ROH on African haplotype backgrounds gain more low-frequency minor allele homozygotes per unit increase of ROH (and especially long class C ROH) compared to common minor allele homozygotes. Since low-frequency alleles are enriched for deleterious variants relative to high-frequency alleles, this result accords with our previous analyses.

Figure 7

Enrichment of Low-Frequency Variants across ROH Sizes

The difference in rate of gain of low-frequency minor allele homozygotes (MAF < 0.05) compared to common minor allele homozygotes (MAF ≥ 0.05; from regression analysis). ROH size classes: A, short; B, medium; C, long; R, all sizes. EUR, European, colored blue; AFR, African, colored orange; and NAM, Native American, colored red. Error bars represent standard error of the regression coefficient.

Simulating Deleterious Alleles in ROH

We have proposed that autozygosity of haplotypes with recent ancestry from high-heterozygosity source populations concentrate deleterious homozygotes at a higher rate per unit increase of ROH coverage (Figure 6). We wish to test via simulations whether these differences in ancestral demographic history can account for this pattern. To this end, we simulate recessive deleterious alleles in a complex three population demographic history, corresponding roughly to African, European, and Asian human populations (see Material and Methods). Although our other analyses considered haplotypes from African, European, and Native American ancestral populations, this three-population demographic model has been well studied and is readily available. As this three-population model contains a high-heterozygosity source population with two population splits undergoing multiple bottlenecks, we feel this will provide a set of simulated data with a qualitatively similar demographic history.

For each of 500 simulation replicates, we sample 500 diploid individuals from each population, call ROH, and then compute the proportion of genome-wide deleterious homozygotes falling within each ROH class. We then compute a regression, similar to the previous section where we analyzed the differences between deleterious homozygotes in ROH on different ancestral backgrounds. We compute, , where is a vector representing the proportion of damaging homozygotes in ROH class j in each population across all individuals. represents the genome-wide fraction ROH class j in each population across all individuals, and is an indicator variable which takes the value 1 if the associated individual is from population and takes the value 0 otherwise. Here AFR corresponds to the simulated African population, EUR corresponds to the simulated European population, and ASN corresponds to the simulated Asian population. We analyze each ROH class: all, long, medium, and short, and within each class we combine our regression coefficients across replicates with inverse-variance weighted meta-analysis.

In this formulation, the regression terms and represent the difference in rate of gain of deleterious homozygotes in ROH on European or Asian haplotypes compared to African haplotypes, respectively. For example, a would represent a scenario where an increase of 1% ROH genome-wide in the simulated European population concentrated 1% more genome-wide deleterious homozygotes in those regions compared to the simulated African population. Similarly, a would represent a scenario where an increase of 1% ROH genome-wide in the simulated Asian population concentrated 1% less genome-wide deleterious homozygotes in those regions compared to the simulated African population. Since we hypothesize that the simulated African population will have the highest rate of gain of deleterious homozygotes as a function of genomic ROH coverage, we expect both of these terms to be negative. Indeed, this is what we find across all ROH classes (Table S7). Considering all ROH together, we find (p < 2 × 10−16) and (p < 2 × 10−16), and when analyzing only long ROH we find (p < 2 × 10−16) and (p < 2 × 10−16).

Discussion

The distribution of runs of homozygosity in individual genomes has provided insights into evolutionary, population, and medical genetics. By examining their genomic location and prevalence in a population, we can learn about the history and adaptation of natural populations,2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 96, 97 and we can make discoveries about the genetic basis of complex phenotypes.32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48 Given the importance of demographic history and socio-cultural practices in the generation of ROH in individual genomes, and their relationship to complex phenotypes including many genetic diseases, it naturally follows to study the distribution of deleterious alleles and their relationship to ROH.

Previous work has described the effect of demographic history on the distribution of deleterious alleles,31, 70, 71, 72, 73, 98 including a few specifically investigating their relationship with runs of homozygosity.17, 29, 31, 52, 53, 99, 100 However, little work has been done on the relationship between deleterious alleles and ROH in admixed populations (although see Mooney et al.). Since there is evidence of very recent bottlenecks (which generate ROH) within admixed populations living in the Americas,63, 100 the relationship between ROH and the accumulation of deleterious homozygotes may provide valuable insights into the genetic basis of complex phenotypes in these individuals.

Here we analyzed 1,441 individuals across three admixed populations: African American, Puerto Rican, and Mexican American. We found that, consistent with other studies, the proportion of deleterious homozygotes found in ROH increases faster than the proportion of benign homozygotes as a function of total genomic ROH (Figure 4 and Table S1). We also found that the genome-wide proportion of deleterious homozygotes in ROH on African ancestral haplotypes increased faster per unit ROH than on ether European or Native American ancestral haplotypes (Figures 5, 6, and Tables S2 and S3). These patterns are also consistent with population-specific worldwide patterns of deleterious homozygotes in ROH, where three of the five African populations analyzed had among the highest rates of enrichment in long ROH.

To explain this observation, we propose that ancestral haplotypes from populations with high deleterious heterozygosity would exhibit even greater increases of deleterious homozygotes per unit ROH. We reason that, under random mating, the larger number of low-frequency deleterious alleles in the population would largely segregate as heterozygotes, whereas, when a harsh bottleneck or consanguinity occurs, these mutations get paired IBD as homozygotes, concentrating more deleterious homozygotes within ROH. Indeed, via simulation of a realistic human demographic history, we found that the rate of gain of deleterious homozygotes was significantly higher in high heterozygosity source populations compared to others (Table S7).

The idea that population bottlenecks and inbreeding can concentrate more deleterious homozygotes on haplotype backgrounds from a high heterozygosity founder population has also been proposed as a reason for the deterioration of the wolf population on Isle Royale, MI, USA. This population, numbering around 50 at its height, was founded by two to three animals from a large and genetically diverse source population on mainland Minnesota. The extreme bottleneck and inbreeding have manifested numerous conspicuous phenotypes among these wolves, and several extremely long ROH have been identified in its members. This can be contrasted with the historically small wolf populations in Ethiopia, which have successfully avoided the pitfalls of inbreeding depression. Robinson et al. further demonstrate through simulations that although historically small populations tend to have a higher burden of deleterious alleles, there are fewer strongly deleterious alleles segregating compared to large populations. Thus, in the event of a population size crash or inbreeding, smaller populations have reduced risk of severe fitness consequences compared to large populations.

This suggests that ROH on haplotypes from high-heterozygosity populations (e.g., African populations) may generate more homozygotes of strong deleterious alleles compared to other haplotype backgrounds. In the context of human health, this may mean that ROH on those haplotype backgrounds are relevant for understanding the genetic basis of various diseases.

Whereas ROH on any haplotype background are associated with an increased rate of deleterious homozygotes, we show that ROH on African haplotypes tend to have a larger share of the genome-wide deleterious homozygotes. Indeed, this accords with recent work that has independently associated increased ROH and increased African ancestry with reduced lung function. This suggests that these ROH on African haplotypes may play a particularly important role in the genetic architecture of complex phenotypes in admixed individuals, especially for populations with African ancestry that have undergone very harsh bottlenecks in the recent past.

Declaration of Interests

The authors declare no competing interests.