Fine-scale population structure and demographic history of British Pakistanis

Abstract

Previous genetic and public health research in the Pakistani population has focused on the role of consanguinity in increasing recessive disease risk, but little is known about its recent population history or the effects of endogamy. Here, we investigate fine-scale population structure, history and consanguinity patterns using genotype chip data from 2,200 British Pakistanis. We reveal strong recent population structure driven by the biraderi social stratification system. We find that all subgroups have had low recent effective population sizes (N_e), with some showing a decrease 15‒20 generations ago that has resulted in extensive identity-by-descent sharing and homozygosity, increasing the risk of recessive disorders. Our results from two orthogonal methods (one using machine learning and the other coalescent-based) suggest that the detailed reporting of parental relatedness for mothers in the cohort under-represents the true levels of consanguinity. These results demonstrate the impact of cultural practices on population structure and genomic diversity in Pakistanis, and have important implications for medical genetic studies.

Fig. 1. Population structure within British Pakistanis from Bradford.

a, b Principal components analysis of 2200 unrelated Pakistani mothers, with the self-reported subgroups with >20 samples indicated in different colours. Plots show PC1 versus PC2 (a) and PC2 versus PC3 (b). Proportion of overall variation explained by each PC is noted in brackets on the axis label. c ADMIXTURE plot (K = 4) illustrating different ancestral components making up the various subgroups, with the largest self-reported subgroups indicated. We have indicated the two distinct subgroups amongst the Rajput.

Fig. 2. Fine-scale population structure inferred amongst 1520 Bradford Pakistani mothers using fineSTRUCTURE.

a The tree illustrates the results of hierarchical clustering of the co-ancestry matrix using patterns of haplotype sharing from ChromoPainter. Each ‘leaf’ on the tree contains multiple individuals. The coloured bars represent the composition of each of the major clusters indicated by the thick black vertical lines. The length of each coloured bar is proportional to the number of individuals in that cluster, with the proportion of each colour representing the fraction of individuals from each self-reported subgroup that make up the cluster. The labels on the edges are the posterior assignment probabilities from fineSTRUCTURE. b The barplot illustrates the fraction of individuals from each self-reported subgroup that fall into the indicated clusters defined by fineSTRUCTURE. The numbers on the barplot indicate the clusters that were used to define a homogeneous subset of individuals from that self-reported subgroup for subsequent demographic analyses.

Fig. 3. Divergence times and historical effective population size changes of Bradford Pakistani subgroups.

Note that these are the homogeneous subgroups defined using the fineSTRUCTURE clusters. a Divergence times estimated using NeON (n = 850 independent individuals). Within each panel, the points show the estimated divergence time between the group indicated on the left on the y-axis and those indicated in the legend. The horizontal lines indicate 95% confidence intervals of the estimated divergence time. b Change in effective population size (N_e) through time estimated with IBDNe. The coloured lines indicate the mean estimate and the grey shading indicates 95% confidence intervals of the N_e estimates.

Fig. 4. Effect of endogamy and consanguinity on IBD sharing and ROH patterns in BiB Pakistani subgroups.

Note that these are the homogeneous subgroups defined using the fineSTRUCTURE clusters. a IBD scores were calculated as the average total length of IBD segments between 5 and 30 cM shared between individuals from the indicated group (n = 964 independent samples), standardised by the value for the 1000 Genomes Finns (FIN). The error bars indicate the IBD score standard error. YRI Yoruba, CEU Western Europeans. b, c Stacked bar plots showing inferred parental relatedness b in the Pakistani mothers broken down by self-reported parental relatedness or c in all white British versus all Pakistani mothers, and in Pakistani mothers broken down by subgroup. d Boxplots showing the fraction of the genome in regions of homozygosity (F_ROH) broken down by subgroup and inferred degree of parental relatedness (n = 850 independent samples). In each boxplot the centre is equal to the median, the upper and lower bounds of the box correspond to 25th and 75th percentiles and the whiskers represent 1.5 times the inter-quartile range (IQR) from the bounds of the box. Outliers are represented as points. Lines indicated expected F_ROH for individuals with parents who are first cousins (red), first cousins once removed (grey) and second cousins (brown). The two-way ANOVA tests whether F_ROH differs between the subgroups.

Fig. 5. Observed ROH and IBD footprints compared to the expectation under a coalescent model.

The footprint is the average fraction of the genome covered by segments of a given length interval. The top lines represent the ROH footprint and the bottom lines the IBD footprint. Points are plotted at the beginning of each 1 cM interval. The N_e profile used to compute the expected footprints was determined by IBDNe using a the Pathan from fineSTRUCTURE’s Cluster 8 or b Jatt/Choudhry from Cluster 10. In a, for the IBD segments, the pink line is beneath the green one. We used kinship values calculated using either the self-reported relationships or the inferred estimates of consanguinity for the relevant group. The observed footprints (black lines) were determined using filtered IBD and ROH calls from IBDseq and bcftools/roh, respectively. Supplementary Fig. 16 shows equivalent plots using the upper and lower bounds of the 95% confidence interval for the IBDNe estimate as the N_e trajectory, and using observed footprints from different ROH or IBD calling/filtering strategies.

Fig. 6. Estimated fraction of couples from intra-biraderi, inter-biraderi and first cousin unions who would be at risk of having a child with a recessive developmental disorder.

These estimates are based on simulations that use P/LP variants ascertained from exome-sequence data in Bradford Pakistani mothers, and rely on the self-reported subgroup information. The y-axis indicates an estimate of the fraction of couples that would be ‘at risk’ of having a child with a recessive developmental disorder since both parents are carriers of P/LP variants in the same gene. Note that this is undoubtedly an under-estimate of absolute risk (see ‘Discussion’), but one can still compare the relative risk for intra- versus inter-biraderi unions. Error bars indicate standard deviations of the risk estimates and one-sided p-values are from permutation tests (see ‘Methods’). The ‘first cousins’ were sampled to match the IBD distribution observed for self-reported first cousin couples (see ‘Methods’), regardless of biraderi group.