Multiple selective sweeps of ancient polymorphisms in and around LTα located in the MHC class III region on chromosome 6

Abstract

Background

Lymphotoxin-α (LTα), located in the Major Histocompatibility Complex (MHC) class III region on chromosome 6, encodes a cytotoxic protein that mediates a variety of antiviral responses among other biological functions. Furthermore, several genotypes at this gene have been implicated in the onset of a number of complex diseases, including myocardial infarction, autoimmunity, and various types of cancer. However, little is known about levels of nucleotide variation and linkage disequilibrium (LD) in and near LTα, which could also influence phenotypic variance. To address this gap in knowledge, we examined sequence variation across ~ 10 kilobases (kbs), encompassing LTα and the upstream region, in 2039 individuals from the 1000 Genomes Project originating from 21 global populations.

Results

Here, we observed striking patterns of diversity, including an excess of intermediate-frequency alleles, the maintenance of multiple common haplotypes and a deep coalescence time for variation (dating > 1.0 million years ago), in global populations. While these results are generally consistent with a model of balancing selection, we also uncovered a signature of positive selection in the form of long-range LD on chromosomes with derived alleles primarily in Eurasian populations. To reconcile these findings, which appear to support different models of selection, we argue that selective sweeps (particularly, soft sweeps) of multiple derived alleles in and/or near LTα occurred in non-Africans after their ancestors left Africa. Furthermore, these targets of selection were predicted to alter transcription factor binding site affinity and protein stability, suggesting they play a role in gene function. Additionally, our data also showed that a subset of these functional adaptive variants are present in archaic hominin genomes.

Conclusions

Overall, this study identified candidate functional alleles in a biologically-relevant genomic region, and offers new insights into the evolutionary origins of these loci in modern human populations.

Background

Lymphotoxin-α (LTα) in the MHC class III region encodes a potent cytotoxic polypeptide that plays a key role in regulating a number of biological processes, including lipid metabolism, coagulation, neurotransmission, and immunological response [1,2,3,4,5,6,7,8,9,10,11,12]. Studies have also shown that common variants in and/or near LTα contribute to the onset of complex diseases. For example, a recent analysis demonstrated that individuals with the A-allele at rs909253— located in an intronic region of LTα— have a higher risk for nasal NK/T-cell lymphoma relative to individuals with the G-allele in a Chinese population [13]. Data have also indicated that variability at rs2239704, rs909253, rs1041981, and rs2229094 was associated with elevated risk for other types of cancer, such as gastric and breast cancers in East Asians as well as non-Hodgkin lymphoma in individuals of European ancestry [14,15,16,17,18,19,20]. In vivo analyses using a mouse model further showed that upregulation of LTα in hepatocytic cells— infected by the hepatitis B or C virus — contributed to apoptosis and/or cell transformation leading to the development of hepatocellular carcinoma [21]. In addition, variability at rs909253, rs1800683 and rs1041981 has been correlated with increased susceptibility to cardiovascular disease in European, East Asian and/or Brazilian populations [10, 22,23,24]. Lastly, allelic variation at rs909253 and/or at rs2229094 has been implicated in the onset of inflammatory/autoimmune disorders, such as chronic periodontitis [25], ankylosing spondylitis [26], rheumatoid arthritis [27], systemic lupus erythematosus [28], vitiligo [29] and Sjogren’s syndrome [30] in individuals of non-African ancestry.

Despite the important role that LTα plays in complex traits, little is still known about levels of nucleotide variation and LD in and/or near this gene. Equally as important, the evolutionary processes that have shaped patterns of diversity in this region are similarly not known. These pieces of information are critical for identifying additional alleles in the LTα region that might contribute to phenotypes, including disease susceptibility. To address this gap in knowledge, we analyzed sequence variation across ~ 10 kbs on chromosome 6, encompassing the LTα gene and the 5′ region (and more broadly across the entire chromosome in some cases), in 2039 individuals from 21 distinct populations in the 1000 Genomes Project. Here, we observed striking patterns of variation in global populations, including an excess of intermediate-frequency alleles, the maintenance of multiple common haplotypes, and a deep coalescent time for variation (dating > 1 million years ago). We also identified a number of common alleles in or near LTα that are present in Neandertal and Denisovan genomes, further supporting the inferred ancient age of nucleotide variation. While these findings are generally consistent with a model of long-term balancing selection, we also observed extensive haplotype homozygosity on chromosomes carrying derived alleles primarily in non-African populations, suggestive of recent selection. Furthermore, a subset of these adaptive alleles were predicted to alter transcription factor binding site affinity and protein stability, suggesting they play a role in gene function. This latter finding could be highly informative for biomedical studies focused on the development of therapeutic interventions that mitigate, mimic or magnify the effects of these functional sites to combat diseases. Overall, our study identified candidate alleles that contribute to phenotypic variation, and offers additional insights into the evolutionary origins of these loci in modern human populations.

Results

Patterns of nucleotide variation

We identified a total of 183 bi-allelic single nucleotide polymorphisms (SNPs) across ~ 10 kbs of sequence, encompassing LTα (2226 base pairs [bps]) and the adjacent 5′ non-coding region (7411 bps) on chromosome 6, in 21 global populations from the 1000 Genomes Project. Of the 183 polymorphisms, 52 SNPs were located within LTα, while the remaining 131 SNPs were found in the 5′ region (Additional file 2: Table S1). The LTα gene is comprised of four exons (Fig. 1); however, the mature LTα protein is encoded by exons 2 (from bps + 461 to + 559), 3 (from bps + 646 to + 751) and 4 (from bps + 999 to + 1411). In the LTα coding region, we identified 12 polymorphisms in global populations, eight of which were nonsynonymous changes. Two of these nonsynonymous polymorphisms (rs2229094 and rs1041981) occurred at relatively high frequency in all populations. More specifically, the minor allele frequency (MAF) at rs2229094 (C-allele) and rs1041981 (A-allele) ranged from 27.3 to 45.8% and from 20.2 to 33.8%, respectively, in South Asian populations. In Europeans, the rs2229094 minor C-allele and rs1041981 minor A-allele varied from 22.9 to 38.3% and from 24.8 to 30.9%, respectively, while the MAF at these loci ranged from 14.4 to 22.8% and from 38.4 to 54.8%, respectively, in East Asian populations. Likewise, the minor C-allele at rs2229094 and the minor A-allele at rs1041981 varied from 22.4 to 30.8% and from 36.3 to 61.5%, respectively, in populations of African descent (which include indigenous Africans, African Americans and African Caribbeans). Further comparative analysis also showed that allelic variation at rs2229094 was present in the genomes of two closely related archaic species to modern humans—Neanderthals and Denisovans (Additional file 2: Table S1).

Fig. 1

Gene structure and distribution of variants among major geographic regions. a Gene structure of LTα consists of four exons (yellow shapes) and intervening introns (black horizontal lines). Exons are labeled 1 through 4, and the start-end nucleotide positions for each exon are given in bps. The translation start site (ATG) begins at nucleotide position + 461 in exon 2; the entire LTα coding region is 618 bps in length. Immediately adjacent to the LTα gene is the 5′ regulatory region shown in purple. The combined length of LTα and the 5′ region is 9637 bps; b Venn diagram showing the number of polymorphisms in LTα found in four geographic regions— namely, Africa, South Asia, Europe and East Asia; c Venn diagram showing the number of polymorphisms in the 5′ regulatory region present in the above-mentioned geographic regions. For b and c, the number of polymorphisms shared among geographic regions is given in the interior portions of intersecting ovals, while the number of polymorphisms private to a particular geographic region is given outside of the intersecting ovals. We did not include African American and African Caribbean (Barbadian) populations with indigenous Africans when we examined the continental differences in the number of polymorphisms

In addition, we detected population-specific nonsynonymous variation in our modern human samples. For example, the minor alleles at rs538402044, rs562333039, rs538877791, rs566451995 occurred at < 2.5% frequency in non-Africans, while these alleles were absent in populations of African descent. Furthermore, the minor C-allele at rs2229092 was common in Europeans varying from 5.1 to 7.0%, while it occurred at lower frequency in South Asians (from 1.7 to 5.5%) and in East Asians (from < 1 to 2.4%). In contrast, the rs2229092 C-allele was absent in indigenous Africans, but was observed at relatively low frequency (~ 1–2.5%) in African American and African Caribbean populations. Additionally, we found that synonymous variants were either absent or occurred at very low frequency (< 1.0%) in all populations (Additional file 2: Table S1). Overall, we observed a striking deficit of synonymous SNPs relative to nonsynonymous polymorphisms in global populations.

In the intronic regions of LTα, we identified a total of 31 SNPs (Additional file 2: Table S1). Alleles at four of these polymorphic sites (rs1800683 rs909253, rs2239704, and rs746868) occurred at relatively high frequency in global populations (Additional file 2: Table S1). In particular, the MAF at rs1800683 and rs909253 both varied from ~ 20.2 to 33.8% frequency in South Asians, from 38.5 to 54.7% in East Asians, from 26.6 to 30.9% in Europeans, and from 38.9 to 64.1% in populations of African descent (Additional file 2: Table S1). Furthermore, the MAF at rs2239704 and rs746868 (Additional file 2: Table S1) ranged from ~ 31 to ~ 44.0% in South Asians, from 33.3 to 48.1% in Europeans, and from 24.1 to 47.1% in East Asians. Comparatively, the frequency of minor alleles at these sites was lower in African and African-descended populations, varying from ~ 16.0 to ~ 38.0%. Our analysis also uncovered a moderate level of population-specific variation; specifically, 25.8, 12.9, 9.7 and 22.6% of the intronic variants were private to populations of African, South Asian, European and East Asian descent, respectively (Additional file 2: Table S1; Additional file 1: Figure S1). Additionally, we found that human derived alleles at rs909253 and rs2239704 were shared with Neandertal and Denisovan hominins.

In the adjacent 5′ region, we identified a total of 131 SNPs, and uncovered extensive allele sharing among distinct populations (Additional file 2: Table S1). Notably, a subset of these polymorphisms (specifically, rs2009658, rs2844484 and rs915654) have previously been implicated in the onset of complex diseases [19, 31]. Like in the LTα gene, our analyses also showed that modern and archaic humans shared a number of derived alleles located in the upstream region (Additional file 2: Table S1).

Lastly, we calculated standard measures of nucleotide diversity (θ_π and θ_W) in LTα and the upstream region, separately, for each global population (Table 1). We found that θ_π (the mean number of pairwise differences per nucleotide) in the LTα gene was similar across diverse populations. Specifically, θ_π ranged from 1.18 × 10^− 3 to 1.39 × 10^− 3 in populations of African descent; from 1.19 × 10^− 3 to 1.33 × 10^− 3 in South Asians, from 1.34 × 10^− 3 to 1.48 × 10^− 3 in Europeans and from 1.19 × 10^− 3 to 1.33 × 10^− 3 in East Asians (Table 1). Likewise, estimates of θ_W (nucleotide diversity calculated based on the number of segregating sites) did not vary greatly among populations (Table 1). In the 5′ region, we detected a similar pattern of nucleotide diversity among geographically distinct populations (Additional file 2: Table S2). Overall, we observed little difference in the level of nucleotide diversity between African and non-African populations across this ~ 10-kb region.

Table 1 Summary statistics based on sequence variation in the LTα gene

Tests of neutrality

To determine if patterns of variation are consistent with neutral evolution, we calculated Tajima’s D (D_T) and the Fay and Wu’s H (H) statistics for the LTα gene and the upstream region, separately, in each population (Table 1; Additional file 2: Table S2). We also generated expected D_T and H values under varying models of demographic growth using the ms software [32]. We incorporated growth as a parameter in these coalescent simulations given the genetic evidence for past human population expansion in prior studies [33, 34]. Based on these analyses, we found that D_T for LTα was more positive than expected (P < 0.05) in South Asians, Europeans and East Asians (Table 1; Additional file 2: Table S3). Moreover, we observed a significant departure of H values (P < 0.05) in South Asian and European populations (Table 1; Additional file 2: Table S4). Our results also showed a general pattern of positive D_T in populations of African descent, with the largest departures from expected values occurring in the Mende from Sierra Leone and the Luhya from Kenya (Table 1; Additional file 2: Table S3). However, we did not observe significant H values in African, African American and African Caribbean populations (Barbadians) (Table 1; Additional file 2: Table S4).

Additionally, we calculated D_T and H statistics for the adjacent 5′ region in each population (Additional file 2: Tables S5 and S6), and found significantly positive D_T values (P < 0.05) in South Asians, Europeans, and East Asians (Additional file 2: Tables S2 and S5). Furthermore, H statistics were more positive than expected (P < 0.05) in these populations (Additional file 2: Tables S2 and Table S6). In comparison, we observed a mix of positive and slightly negative D_T values (Additional file 2: Tables S2 and S5) among populations of African descent. However, we did detect significantly positive H values (P < 0.05) in both indigenous and recently admixed Africans (Additional file 2: Tables S2 and S6).

To further assess whether or not variation at LTα is evolving neutrally, we applied the McDonald-Kreitman (M-K) test, which compares the ratio of synonymous and non-synonymous sites within and between species, to our sequence data (Table 1). Our results showed a significant excess of nonsynonymous variants (P = 0.047) in the Kinh population from Vietnam and a borderline significant excess of amino acid variation (0.1 ≥ P ≥ 0.05) in other populations of African (African Americans and the Mende), South Asian (Bengali, Gujarati, Telugu, Punjabi and Sri Lankan Tamil), European (Finnish, Great British, and Iberian), and East Asian (Dai and southern Han Chinese) ancestry (Table 1). Thus, we consistently observed a higher proportion of nonsynonymous polymorphisms relative to synonymous changes in all populations. Overall, the excess of nonsynonymous variants and the sharply positive D_T and H statistics in global populations are consistent with the action of long-term balancing selection [35,36,37].

Haplotype variation and inferred relationships

We analyzed phased haplotype data, encompassing the entire 9637 base pair (bp) region, and identified 223 distinct haplotypes in global populations (Additional file 2: Table S7A). Of these variants, 23 haplotypes (H6, H8, H12, H30, H36, H37, H80, H82, H86, H105, H119, H128, H132, H136, H142, H158, H160, H187, H190, H204, H214, H217 and H219) were common, and six of these haplotypes were found in all populations (Additional file 2: Table S7A); 16–18 common haplotypes accounted for 80.4–88.4%, 77.6–83.1% and 77.1–88.9% of the total number of haplotypes in South Asian, European and East Asian populations, respectively (Additional file 2: Table S8). Comparatively, a smaller number of common haplotypes (H6, H30, H37, H82, H86, H119, H128, H142, H187, H190, H217 and H219) comprised 66.7–74.2% of all haplotype lineages in African and African-descended populations (Additional file 2: Table S9). Further analysis of all haplotype variation also showed that Africans and non-Africans exhibited high levels of haplotype diversity (h-diversity). Specifically, h-diversity varied from 90.5 to 95.0% in populations of African descent; from 91.5 to 93.6% in South Asians; from 91.9 to 93.3% in Europeans; and from 85.9 to 89.9% in East Asians (Table 1).

We also performed a median-joining network analysis [38] to explore the phylogenetic relationships among haplotype lineages, and found that haplotype variation grouped into three distinct clusters (arbitrarily labeled I, II, III; Fig. 2). Each cluster consisted of low-frequency haplotypes radiating from common haplotypes in a “star-like” pattern which is the genetic pattern predicted under a model of positive or purifying selection [39]. Moreover, the highest-frequency haplotypes (specifically, H37, H82, H128, H142 and H190) in clusters I, II and III were shared among globally diverse populations, and these haplotype lineages differed from one another at 28 polymorphic sites. Intriguingly, we also found that 22 of the 28 polymorphic sites on these high-frequency haplotypes were present in the Neandertal and Denisovan genomes (Additional file 2: Table S7B).

Fig. 2

Haplotype network of global variation across ~ 10 kbs encompassing LTα and the 5′ region. Haplotype relationships were constructed using the median-joining algorithm implemented in Network 5.0 [38]. The resulting phylogeny consists of nodes (representing haplotypes), and links which connect the nodes in the shortest possible tree. The size of the nodes is proportional to the number of chromosomes with a given haplotype, and the colors within haplotypes (or nodes) indicate the frequency of that haplotype in pooled populations (such as, South Asians or East Asians). Distinct clusters of haplotypes are surrounded by circles/ovals with dashed lines, and the alleles that define the highest-frequency haplotype(s) are also given in each cluster. Underlined nucleotides indicate alleles that are inferred to be under selection and are shared with archaic hominins. Each haplotype cluster was arbitrarily labelled as Cluster I, Cluster II and Cluster III

Pairwise linkage disequilibrium in the LTα region and surrounding genes

To quantify the allelic associations between SNPs across ~ 10 kbs, encompassing LTα and the 5′ region, we estimated pairwise LD using the D′ statistic [40]. Our results generally showed a faster decay of LD between loci (D′ < 100) in populations of African ancestry compared to non-Africans. For example, in the Esan from Nigeria, we observed lower levels of allelic association (D′ < 100) within and between LTα and the upstream region (Fig. 3; Additional file 1: Figure S1; Additional file 2: Table S10), implying a history of recombination in this population. By contrast, our LD plots revealed a higher level of allelic association in South Asians, Europeans, and East Asians. In particular, we observed a greater number of instances of complete LD (D′ = 100) between loci in non-Africans compared to populations of African descent (Fig. 3; Additional file 1: Figure S1; Additional file 2: Table S10). The inference of recombination across LTα and the 5′ region is further supported by our phylogenetic network which displayed a small number of reticulations or loops among haplotype lineages (Fig. 2). It is known that these reticulate (non-bifurcating) relationships among haplotypes can arise due to historical recombination events [41].

Fig. 3

Pairwise LD across ~ 10 kbs encompassing LTα and the 5′ region in selected populations. We measured pairwise LD, using the D′ statistic, calculated with Haploview [40], which outputs a graphical representation of allelic associations. In the resulting plots, pairwise LD coefficients D′ × 100 are given in each cell. The color scheme of each cell signifies the strength of association between SNP alleles. Specifically, bright red squares indicate statistically significant LD between SNP pairs (D′ = 100; logarithm of odds (LOD) > 2), while shades of pink/white cells signify little evidence of LD (D′ < 100; LOD < 2) between loci and purple cells signal high LD but with little statistical support (low LOD) [40]. The bold triangles in the plots also indicate strong blocks of LD between SNP markers. The corresponding genomic coordinate for each marker can be found in Additional file 2: Table S10. The genomic region analyzed is given at the top of each plot (exons are labeled 1 through 4)

To explore possible gene-gene interactions, we also quantified pairwise LD between genetic markers in the LTα region and neighboring genes (specifically, NFKBIL1, TNFα, and LTβ) within a 35,538-bp region on chromosome 6 (Additional file 1: Figure S2). This analysis showed that African and African-descended populations exhibited less LD (D′ < 100) between SNPs in the LTα region and nearby genes (Additional file 1: Figure S2). In comparison, however, we observed more instances of complete LD (D′ = 100) across this ~ 35.5 kilobase (kb) region in Eurasian populations (Additional file 1: Figure S2).

Long-range haplotype structure

We further characterized long-range LD across the entire ∼170 megabase (Mb) region of chromosome 6 for each population using the iHS statistic [42]. In particular, we calculated the absolute standardized |iHS| scores by normalizing raw values with the selscan and norm software [43]. Based on this analysis, we observed outlier |iHS| statistics for multiple SNPs— namely, rs4947326, rs4947327, rs9267497, rs36221456, rs73396237, rs148407582, rs9267499, rs2857709, rs2857708 and rs2229092— in the coding and 5′ regions in South Asians (Fig. 4; Additional file 1: Figure S3; Additional file 2: Table S12). We also identified multiple SNPs with elevated |iHS| values in the upstream, intronic and exonic regions in European populations (Fig. 4; Additional file 1: Figure S3; Additional file 2: Table S12). In addition, we detected several outlier |iHS| scores for loci in the 5′ region only (rs73396237, rs915654, rs2857709) in populations from East Asia (Fig. 4; Additional file 1: Figure S3; Additional file 2: Table S12). We did not, however, observe outlier values for SNPs in the LTα gene or the 5′ region in populations of African descent, except in African Americans where we detected a single extreme |iHS| score for rs73396237 located in the upstream region (Fig. 4; Additional file 1: Figure S3; Additional file 2: Table S12).

Fig. 4

Integrated haplotype score (iHS) and extended haplotype homozygosity (EHH) plots for selected populations. a–d Manhattan plots of standardized |iHS| statistics for selected populations. The dashed horizontal lines represent the cut-off for extreme |iHS| scores (|iHS | > 2), representing the most extreme 5% of empirical |iHS| values. The red dots in plots represent outlier SNPs in and/or near LTα; (E–H) EHH plots for selected populations show the decay of identity of haplotypes on chromosomes with a derived allele (the blue line) versus an ancestral allele (the red line) at core SNPs as a function of distance. The distance from the core SNP (at zero) is displayed on the x-axis; the negative numbers indicate distance upstream from the core SNP, while positive values indicate distance downstream from the core SNP. The EHH values are shown on the y-axis. It is important to note that these EHH plots are selected examples from a larger set of EHH curves that were generated in the present study, which can be found in Additional file 1: Figure S5)

To complement our |iHS| scans, we applied the nSL statistic which is a test for detecting both soft and hard sweeps in genomic data [44]. One advantage of this statistic is that it does not depend on the recombination rate, which can lead to false positive signals of haplotype homozygosity in chromosomal regions with low levels of recombination [44]. Using this method, we observed extreme standardized |nSL| scores (defined as |nSL| > 2) for variation in the 5′ region (rs36221456, rs73396237, rs2857709, rs4947326, rs2857601, rs2857708) and in the LTα coding region (rs2229092) in South Asians (Additional file 1: Figure S4; Additional file 2: Table S11), generally overlapping with the |iHS| results. We also detected outlier |nSL| values for a similar set of loci upstream (rs73396237, rs2857709, rs3093540, rs2071590, rs4947326, rs2857601, rs2239704 and rs3093546) and within the LTα gene (rs3093542 and rs2229092) in European populations (Additional file 1: Figure S4; Additional file 2: Table S12). Furthermore, we found a single extreme |nSL| score for rs73396237 in the upstream region in East Asians. By contrast, we did not observe unusual |nSL| scores for loci in any population of African descent (Additional file 1: Figure S4; Additional file 2: Table S12).

We also plotted the decay of haplotype homozygosity from core SNPs (represented by loci with extreme |iHS| or |nSL| scores) using the EHH statistic [45]. This analysis revealed long-distance haplotype homozygosity, spanning up to 1.4 megabases (Mbs), on chromosomes carrying the derived allele relative to chromosomes with the ancestral allele mainly in Eurasian populations (Fig. 4; Additional file 1: Figure S5). For example, we detected unusually long haplotype structure around derived alleles at rs73396237, rs2857709, rs2857708, rs4947326, rs4947327, rs148407582, rs9267497, rs9267499, and rs2229092 in the South Asian Telugu (ITU) population (Fig. 4; Additional file 1: Figure S5). We also observed a similar pattern of EHH on chromosomes carrying derived alleles in other non-African populations. In addition, we detected long-range EHH around the derived G-allele at rs73396237 (extending no more than 80 kbs) in African Americans in agreement with our |iHS| results (Fig. 4; Additional file 1: Figure S5). We did not, however, observe a difference in EHH curves for chromosomes harboring the derived versus the ancestral allele in other populations of African descent (Additional file 1: Figure S5).

Population differentiation and structure

To examine patterns of genetic differentiation, we calculated mean F_ST among geographic regions (namely Africa, South Asia, East Asia, and Europe) using polymorphisms across the ~ 10 kbs encompassing LTα and the upstream region. Mean global F_ST was estimated to be 0.010132, much lower than F_ST estimates for autosomal loci, which typically range from 0.10 to 0.16 [46, 47]. Our analyses also revealed that per-site F_ST values for common SNPs were below 0.10 (Additional file 2: Table S1). For example, F_ST values for rs1800683, rs2239704, rs909253, rs746868, rs2229094, rs3093542, rs2229092 and rs1041981 within the LTα gene were estimated to be 0.069033, 0.020427, 0.068599, 0.021288, 0.029649, 0.026510, 0.034892 and 0.067877, respectively (Additional file 2: Table S1). Similarly, F_ST values for common SNPs in the 5′ region varied from 0.0055 to 0.0745 (Additional file 2: Table S1). Although the above estimates of F_ST were not outliers in the empirical distribution of F_ST, calculated for ~ 7.5 million genome-wide SNPs, they did fall within the 15th percentile of empirical values. Thus, we consistently observed low estimates of F_ST in both the LTα gene and the upstream region.

Inferred ages of variants in LTα and the adjacent 5′ region

Using a coalescent-based method [48], we inferred the mean time to the most recent common ancestor (TMRCA) for variation in both LTα and the upstream region, and calculated the mean ages of individual polymorphisms (in years ago ± standard deviation in years) (Additional file 2: Table S13). Our analysis revealed a deep coalescence time for the origin of variation, estimated to be ~ 1.0 million years ago (ya) ± 57,207.68 years. We also found that the ages of common polymorphisms in the 5′ region ranged from 200,571.46 ya (± 77,877.56 years) to 734,694 ya (± 140,326.55 years) (Additional file 2: Table S13). Furthermore, variants with the oldest inferred ages—namely, rs2857710 (734,694 ya ± 140,326.55 years), rs3131641 (453,306.20 ya ± 69,869.40 years), rs2857602 (718,530.73 ya ± 141,061.25 years), rs2844486 (626,693.98 ya ± 177,795.95 years) and rs1121800 (472,408.24 ya ± 106,530.63 years)—were also found in the Neanderthal and/or Denisovan genomes (Additional file 2: Tables S1 and S13). Additionally, several of the SNPs with deep coalescence times (i.e., rs4947326, rs4947327, rs148407582 and rs9267497) were inferred to be targets of recent selection based on our long-range haplotype analyses (Additional file 2: Tables S11 and S13).

Inferred functional consequences

Given the higher proportion of nonsynonymous polymorphisms relative to synonymous changes in the LTα coding region, we inferred the possible effects of nonsynonymous SNPs on protein function using SIFT [49]. SIFT classifies nonsynonymous substitutions as either “tolerated” or “deleterious” based on sequence homology and the properties of amino acids [49]. This method predicted that seven of the eight nonsynonymous polymorphisms (specifically, rs538402044, rs2229094, rs368539892, rs2229092, rs1041981, rs538877791 and rs566451995) were “tolerated”, while variation at rs562333039 was predicted to be “deleterious”. Moreover, three “tolerated” variants— rs2229094, rs2229092, rs1041981— were common (MAF ≥ 5%) in at least one population, while the other four “tolerated” polymorphisms were either rare (MAF < 5%) or absent in global populations. The single “deleterious” SNP was present only in the Kinh population from Vietnam at very low frequency (MAF < 0.01). We also determined if common missense variation was conserved across mammalian species using GERP++ [50], phyloP [51] and phastCons [52]. Our results showed that the conservation score for each polymorphic site did not reach the “deleterious” threshold level, in agreement with our SIFT results (Table 2).

Table 2 Predictions of amino acid changes on LTα protein function

In addition, we analyzed changes in LTα protein stability (i.e. the folded conformation) caused by missense mutations at different common SNPs with two bioinformatic tools: FoldX [53] and I-Mutant3 [54]. Our analyses revealed that the total energy change (ΔΔG) for rs1041981, determined by FoldX, was − 0.802 kcal/mol, suggesting that mutation (the A-allele) at this site slightly destabilizes protein structure (Table 2). By contrast, the FoldX ΔΔG for rs2229092 was 1.180 kcal/mol, implying that variability at this site largely stabilizes the folding conformation of the LTα protein (Table 2). The I-Mutant3 method yielded similar predictions of protein stability changes due to nucleotide variation at the rs1041981 and rs2229092 loci. The folding energy difference caused by mutation at rs2229094 —calculated using I-Mutant3— was 0.07 kcal/mol (Table 2), indicating that variability at this site has little effect on the folded structure of the LTα protein. Regrettably, we were unable to determine the FoldX ΔΔG for rs2229094 because the available template structure, used as input for this bioinformatic tool, did not include residues encoded by this SNP.

Lastly, we examined the potential regulatory effects of variants in or near the LTα gene using the SNP2TFBS database [55]. Based on this analysis, we identified 13 non-coding polymorphisms that map to transcription factor binding sites (TFBSs), and are predicted to alter transcription factor binding site affinity (Table 3) [56]. Our data also showed that the transcription factors (TFs) affected by nucleotide changes in the TFBSs were Pdx1, Prrx2, KLF1, ZNF263, SP2, SP1, KLF5, ELF-1, Zfx, Nr1h3/Rxra, Arnt/Ahr, EGR2, Foxd3, FOXI1, IRF1, SPIB, ARID3A. Furthermore, a subset of these 13 SNPs (namely, rs2857709, rs4947326, rs4947327 and rs3093542, which are predicted to affect the binding of TFs ELF-1, Zfx and ZNF263) had extreme |iHS| and/or |nSL| scores, suggesting these candidate functional alleles are or have been targets of recent selection (Table 3; Additional file 2: Table S12). In addition, derived alleles at these sites had coalescence times that ranged from 56,718.38 ya (± 32,106.13 years) to 668,571.54 ya (± 105,795.94 years) (Additional file 2: Table S13). Our analysis also revealed that one of the 13 candidate functional polymorphisms (specifically, rs909253)— not inferred to be under selection— is a known genetic risk factor for several complex diseases [10, 13, 57, 58] , and is predicted to alter the binding of TF SPIB (Table 3).

Table 3 The effect of allelic variation on transcription factor binding affinity

Discussion

Evidence for soft sweeps in global populations

Our analyses uncovered striking patterns of diversity in LTα and the 5′ region in African and non-African populations. In particular, we found significantly positive D_T and H statistics, indicating an excess of intermediate-frequency-derived alleles, within and outside of Africa. Furthermore, we observed the maintenance of multiple high-frequency haplotype lineages in global populations, and inferred a deep coalescent time of variation (> 1.0 million years ago). We also detected low levels of genetic differentiation and extensive sharing of variation among globally diverse populations. A comparative analysis further revealed that a number of polymorphisms were shared between archaic and modern humans (i.e. trans-species polymorphisms). Taken together, these results are largely consistent with a model of long-term balancing selection. However, we also detected long-range EHH on chromosomes carrying derived alleles primarily in non-African populations, and observed a “star-like” phylogeny of haplotypes radiating around common haplotypes in our network. These genetic patterns are consistent with the predictions for positive selection. To reconcile these findings, which appear to support different models of selection, we argue that long-term balancing selection has acted to maintain alleles/haplotypes in African populations, while patterns of diversity and LD in non-Africans were shaped by recent selection.

Within Africa, we found evidence for balancing selection based on the site frequency spectrum and the persistence of multiple common haplotypes in populations. These results are similar to the findings in a prior study that detected signatures of balancing selection at the TAS2R38 locus—which is associated with normal variation in bitter taste perception—in ethnically diverse African populations [59]. In addition, we detected shorter blocks of EHH in most populations of African descent, likely reflecting a history of recombination in and around LTα [60,61,62]. One exception to this general pattern of LD, however, was observed in African Americans where we found unusually long haplotype structure around the derived allele at rs73396237. This SNP is also inferred to be a target of selection and is common in European populations. One possible explanation for this pattern is that the extended LD surrounding this selected SNP could have been introduced into African Americans through gene flow from the ancestors of western Europeans to enslaved West Africans during the Trans-Atlantic slave trade [63]. Indeed, prior studies have shown that contemporary African Americans have varying levels of European ancestry due to this historical event which occurred between the 15th and 19th centuries [63,64,65]. Additionally, the presence of other derived variation in African Americans (for example, the C-allele at rs2229092 in the LTα coding region), common in Europeans but notably absent in West Africans, also supports this model of historical admixture. Alternatively, and more provocatively, the common derived allele at rs73396237 (inferred to be under selection) could have already been present in the enslaved ancestors of African Americans and then became selectively advantageous in the new environments (for example, novel pathogens) to which these individuals were exposed, leading to extended haplotype structure over time. However, further analysis of a larger set of populations will be needed to distinguish: 1) whether the above extended LD in African Americans was due to strong selection on pre-existing common variation at rs73396237 that became beneficial in novel environments, or 2) a selective event had occurred first in the ancestors of Europeans who then admixed with the West African ancestors of present-day African Americans. Interestingly, recent studies have reported instances of populations acquiring selected alleles through past admixture [66,67,68]. Thus, it is highly conceivable that a selected allele could have been introduced into African Americans through admixture between African and non-African parental populations.

After modern humans migrated from Africa ~ 60,000–80,000 ya, we contend that positive selection on pre-existing genetic variation occurred in regions of the world outside of the African continent. Indeed, this model of adaptation has been reported for several genes in non-African populations [69]. In particular, a recent study uncovered evidence for long-term balancing selection at several loci (CLCNKB, PKDREJ, SDR39U1, and ZNF473) in West and East Africans [69]. However, these genes showed an absence of some of the signals of balancing selection in European and/or East Asian populations; data also indicated that one of the two alleles under balancing selection in Africa underwent a soft sweep in non-Africans [69], likely in response to changes in selective pressures outside of Africa [69,70,71,72,73]. In the present study, we propose that a similar phenomenon occurred in the LTα region during the history of non-Africans, leading to the presence of long-range haplotypes in these populations.

Additionally, we suggest that soft sweeps of multiple advantageous mutations occurred— such that pre-existing alleles of similar benefit were selected for and increased in frequency simultaneously, or in short succession of one another (in response to different selective pressures). Under this model of evolution (i.e. a multiple-locus soft sweep model), none of the favored mutations rise rapidly to fixation [74], but rather beneficial alleles at different loci undergo an incomplete selective sweep. In this scenario, multiple mutations can co-exist at intermediate frequency in populations, mimicking balancing selection (resulting in strongly positive D_T and H values) [37, 75,76,77,78]. Furthermore, multiple haplotypes carrying selectively advantageous alleles will become frequent within populations, leading to an increase in haplotype diversity [37, 79, 80]. Explicitly, both of these genetic patterns were consistently observed in Eurasian populations. This proposed model of selection is distinct from a classic selective sweep in which a novel beneficial mutation arises once and quickly rises in frequency until it becomes fixed [7].

Lastly, we detected a significant and marginally significant excess of nonsynonymous polymorphisms in African and non-African populations. This bias towards amino acid variation in the LTα coding region could arise due to several evolutionary forces, such as positive selection, a relaxation of functional constraint, or balancing selection in diverse populations [60, 81,82,83]. Although the precise environmental/external pressures driving the different types of selection in and outside of Africa are currently unknown, it is clear that multiple loci in LTα and the upstream region are or have been functionally important during human evolution.

Identification of functional candidate loci

The human genome contains thousands of experimentally corroborated transcription factor binding sites (TFBSs)[84, 85]; precise knowledge of the genetic variants that disrupt TFBSs is critical for understanding the molecular changes that contribute to phenotypic variation [86]. Using a computational approach, we identified a number of SNPs that were predicted to alter the binding of multiple TFs— most notably Zfx, ELF-1, SPIB and ZNF263. Prior studies have also documented the role that these TFs play in gene regulation and immunity. For example, studies have reported that Zfx is essential for the survival of recirculating mature B-cells and embryonic stem cells. Moreover, in vitro experiments demonstrated that Zfx-deficient peripheral T-cells failed to proliferate and expand after microbial antigen stimulation [87,88,89]. Other studies have also indicated that the ELF-1 binding site is important for the initiation of transcription and that the knockdown of this TF inhibited gene expression [90,91,92]. Data further showed that TF SPIB plays a role in B-cell differentiation (via gene regulation), enabling these lymphocytes to appropriately respond to environmental cues [93, 94]. In addition, analyses focused on ZNF263 found that this TF can have both positive and negative effects on transcriptional regulation of gene targets; more specifically, in cells where ZNF263 levels were low, a subset (∼15–20%) of genes exhibited decreased transcription, while a larger proportion of genes showed elevated expression [95]. As a result, we suggest that the candidate loci, identified in the present study, likely contribute to phenotypic variation through the regulation of gene expression (though it is not completely clear from our data in which direction these variants affect expression; that is to say, we do not know if they inhibit or enhance expression). Intriguingly, we also found that rs909253, which is predicted to alter the binding of SPIB, is associated with increased risk for myocardial infarction, non-Hodgkin lymphoma, and psoriatic arthritis [10, 13, 58], suggesting that the onset of these complex diseases may be due in part to LTα expression levels. However, further studies are needed to elucidate the role that LTα expression plays in complex disease susceptibility.

Our analyses also revealed strong LD among alleles in and between LTα and the 5′ region, raising the possibility that our candidate SNPs are in LD with previously identified variants in the literature. Therefore, it is conceivable that our candidate loci are either causal or act epistatically with previously described alleles in association studies. Furthermore, given the evidence for strong LD between alleles located in different genes (mainly in non-African populations), it is possible that loci in the LTα region and in neighboring genes could also interact epistatically, influencing gene expression. Additionally, a subset of the candidate functional sites (i.e., rs2857709, rs4947326, rs4947327, and rs3093542) were inferred to be under selection, suggesting that adaptive evolution in the LTα region involved changes in the level of gene expression.

Although computational approaches have been beneficial for predicting the impact of nucleotide variation on TF-DNA binding [96,97,98,99,100,101,102], the ability to infer functional consequences of nucleotide changes in TFBSs using such methods has its limits. One important limitation to note is that only a small fraction of all existing TFs in the genome have been characterized in terms of their binding properties to date [86]. As a result, a number of TFBSs in the genome do not have known cognate TFs. Equally as important, current computational methods examine each variant in DNA binding sites independently, neglecting the potential epistatic interactions among distinct TFBSs. Therefore, future in vitro and in vivo TF binding experiments will be important for further clarifying how the novel loci in our study singly or jointly, influence function, leading to differences in gene expression.

Our data also indicated that the high levels of amino acid change in the LTα protein is generally tolerated (i.e. they do not disrupt protein function), suggesting that sequence diversity may have a functional purpose. For example, greater amino acid diversity could expand the ability of the LTα protein to recognize a broad range of antigens, enhancing immune surveillance. Equally as important, the accumulation of amino acid changes could be critical for the other key functions that the LTα protein performs. While the reasons underlying this bias towards amino acid polymorphisms are not clear, this pattern of diversity has been previously reported for a number of other MHC genes [81]. In addition, our ΔΔG results indicated that the “mutant” proline residue at amino acid position 51— encoded by rs2229092— has a stabilizing effect on LTα protein structure relative to the “wild-type” histidine residue at the same site (Table 2). Indeed, it is widely recognized that high stability is important for preserving protein function in a range of conditions over time (for example, temperature and pH) [103]. Furthermore, extra stability is often correlated with protein evolution; specifically, it increases the tolerance of proteins to amino acid substitutions while still folding to its native structure [104]. Thus, highly stable proteins can function across broader physiochemical environments and accept a greater number of amino acid changes, increasing their capacity to perform diverse biochemical tasks [104]. Given that rs2229092 was also inferred to be a target of adaptive evolution in non-Africans, we contend that there has been recent selection for nucleotide variability that enhances protein monomer stability, suggesting that the LTα protein may play functionally important and diverse roles in these populations.

Though our analyses of protein structure indicated that nucleotide changes at rs2229092, rs2229094 and rs1041981 were generally tolerated, prior studies have reported that missense mutations at rs2229094 and rs1041981 are associated with elevated risk for several complex diseases [14, 17, 29]. Given our findings in the present study, we argue that protein instability (leading to dysfunction of the LTα protein) is not likely a major contributor to disease outcome. Further investigation of other possible molecular mechanisms may be more informative for understanding the underlying cause(s) of disease susceptibility (for example, rs2229094 and rs1041981 may be in LD with other loci that alter the expression of LTα proteins leading to increased risk for disease).

Timing of common variants in African and non-African populations

Our GENETREE analysis revealed an ancient TMRCA for variation in and near LTα (dating to > 1.0 million ya), and a deep coalescence time for common alleles that predate the origin of anatomically modern humans ~ 200,000–300,000 ya [63, 105]. The presence of human polymorphisms in Neanderthal and Denisovan genomes further support the inferred ancient age of common nucleotide variation in our dataset. This sharing of alleles among closely related species suggests that variation in and/or near LTα likely arose in modern H. sapiens because it was either: 1) inherited from the last ancestor that modern humans shared with Neandertals and Denisovans more than 500,000 ya [106], or 2) introduced into modern humans through admixture with Neanderthal and Denisovan populations who overlapped in space and time with H. sapiens. Because we observed these ancient polymorphisms in indigenous African populations, which are known to have low levels of Neandertal/Denisovan ancestry [107, 108], we argue that it is likely these variants originated from the last common ancestor of H. sapiens, Neandertals, and Denisovans in Africa rather than through introgression. Regardless, the inferred deep TMRCA of polymorphisms and their presence in closely related hominins suggest that these variants have persisted in the Homo lineage for a relatively long period of time.

Lastly, although classical and soft selective sweeps are expected to increase allele frequency differences among populations [109, 110], we did not observe extreme genetic differentiation at common SNPs in and/or near LTα. One possible explanation for this pattern is that both ancestral and derived alleles were common in populations in Africa. When modern humans migrated out of Africa within the last 100,000 years [63, 111], they carried ancestral and derived alleles at common sites with them. Because these alleles occur at appreciable frequencies in global populations, it is not surprising that we would not observe large among-population differences in allele frequency and thus not detect unusually large F_ST values at common sites [112].

The function and evolution of MHC genes

The MHC is a key component of the adaptive immune system in all jawed vertebrates [113,114,115]. In humans, MHC is a gene-dense region that spans ∼4 Mbs on the short arm of chromosome 6 and contains over 200 genes [113]. The classical MHC class I and class II genes (also called the human leukocyte antigen, HLA, system) encode cell-surface glycoproteins that play a key role in adaptive immunity [114, 115]. In particular, in cells infected by intracellular parasites, MHC class I molecules present parasite-derived peptides to cytotoxic T lymphocytes [114, 115]. Similarly, the MHC class II molecules present antigens (mainly from extracellular pathogens) on the surface of antigen-presenting cells [114, 115]. These exposed peptides are then recognized by helper T lymphocytes, resulting in a series of immune responses [116, 117]. In comparison, proteins encoded by MHC class III genes have somewhat different functions, playing roles in inflammation and complementary immune responses among other immune-related activities [118].

The classical MHC genes are among the most polymorphic in the human genome [114]. Furthermore, alleles at these loci are ancient in age (often predating species divergence events) and can be retained across multiple species. Studies have also indicated that balancing selection has likely acted at these genes to enhance the level of nonsynonymous substitutions in codons and allelic diversity over extremely long time periods [81, 114, 119, 120]. Comparatively, little is known about patterns of variation at genes within the MHC class III region or the microevolutionary forces that have shaped these patterns. However, some studies have identified intergenic sequences in this region that contain cis-acting elements for transcriptional or post-transcriptional regulation of gene expression [121, 122]. For example, non-coding sequences upstream of TNF-α were found to contain transcription factor and enhancer binding sites that influence TNF-α expression in blood serum [122]. In the present study, we similarly identified polymorphisms in the 5′ region of LTα that are predicted to alter LTα expression levels. However, more in-depth studies of the MHC class III region are needed in order to identify variants with a similar regulatory effect on other structural genes in this chromosomal region, and to better understand the evolutionary history of these genes in the human lineage.

Concluding remarks

It has been estimated that more than 90% of the SNPs identified in prior association studies lie outside of protein coding genes and that a substantial fraction of these polymorphisms affect gene regulation [123]. Therefore, understanding whether or not variation within TFBSs disrupts or creates functional sites will be informative for elucidating the molecular basis of gene regulation and complex traits, including disease susceptibility. Our study identified SNPs in TFBSs that potentially influence gene expression and were also shared with archaic hominins. The biological consequences of polymorphisms (in modern H. sapiens) that originate from closely related species have been a focus of intense scientific inquiry in the field of human genomics [124, 125]. For example, a number of studies have indicated that introgressed Neandertal alleles, occurring at relatively high frequency in contemporary Eurasian populations, contribute to immunity, metabolism, height, hair color, and skin tone [124,125,126,127,128,129,130,131,132,133,134,135,136]. Analyses have also shown that introgressed archaic alleles have been targets of positive selection in modern populations. Thus, polymorphisms of non-modern origin are argued to have played important roles in human biology. In the current study, our data suggest that ancient shared alleles in the LTα region— likely inherited from the common ancestor of H. sapiens, Neandertals, and Denisovans— were functional and adaptive in modern humans. These findings raise new and intriguing questions about whether or not archaic and modern humans exhibited similar phenotypic variation due to their shared allelic variation. Indeed, interrogating the genomes of related species— such as Neandertals and Denisovans— can provide insights into the evolutionary origins of alleles associated with traits that may or may not be unique to modern humans.

Methods

Population samples

We analyzed the 9637 bps, encompassing the LTα gene (2226 bps) and the adjacent 5′ region (7411 bps) on chromosome 6 in 2039 unrelated individuals from the 1000 Genomes Project (Phase 3). These individuals originated from four geographic regions: [1] Africa includes 99 Esan (ESN, Nigeria), 113 Mandinka (GWD, The Gambia), 54 Mende (MSL, Sierra Leone), 107 Yoruba (YRI, Nigeria), 99 Luhya (LWK, Kenya), 60 African Americans (ASW, United States), 96 Barbadians (ACB, Barbados); [2] South Asia encompasses 84 Bengali (BEB, Bangladesh), 103 Gurjarati (GIH, India), 99 Indian Telugu (ITU, India), 93 Punjabi (PJL, Pakistan), 99 Tamil (STU, Sri Lanka); [3] East Asia consists of 93 Dai (CDX, China), 103 Han (CHB, Beijing), 105 Han (CHS, southern China), 104 Japanese (JPT, Japan), 99 Kinh (KHV, Vietnam); [4] Europe comprise 90 Great British (GBR, Great Britain), 99 Finnish (FIN, Finland), 107 Iberian (IBS, Spain), 104 Toscani Italians (TSI, Italy) [137]. Although we grouped African Americans and African Caribbeans (Barbadians) with indigenous Africans above, we did not combine these recently admixed population with Africans when we examined the geographic distribution of polymorphisms in Fig. 1b and c. As a result, the total number of variants in Fig. 1b and c does not include the population-specific variation present in African Americans and African Caribbeans.

Nucleotide variation

We extracted variant calls (in vcf files; build GRCh 37 human assembly) from the LTα gene and the adjacent 5′ region in the 1000 Genomes sequencing dataset using vcftools [138]. The start and end positions for this region (build GRCh 37) were obtained from NCBI [139]. The minor allele frequency (MAF) was defined as the second most frequent allele at a given site in pooled populations (i.e., 4078 chromosomes). SNPs were also broadly classified as common or rare based on the MAF at a given site. More explicitly, SNPs with an MAF ≥ 5% were classified as common, while SNPs with an MAF < 5% were categorized as rare. For comparative analysis, we examined sequence contigs for ~ 10 kbs, encompassing the LTα and 5′ region, in Neandertal and Denisovan samples from the UCSC genome browser [140].

Tests of neutrality

We calculated the Tajima’s D (D_T) [141] and Fay and Wu’s H (H) [142] statistics for LTα and the adjacent 5′ region, separately, in each population. Significance of the test statistics was assessed by comparing the observed values to expected values generated from 10,000 neutral coalescent simulations incorporating different models of growth using the ms software [32]. If a sample showed a statistically significant deviation from the expected theoretical distribution, the null hypothesis of no selection was rejected. For Africans, we calculated expected D_T and H statistics under a 2-, 4-, 6-, 8-, and 10-fold increase in population size (starting from 10,000 individuals) beginning 70,000 ya until the present [34, 143,144,145]. For non-Africans, we calculated the expected D_T and H values under a range of demographic scenarios that included a population bottleneck at 60,000 ya (from an initial population size of 10,000 individuals and decreasing to 2000 individuals), followed by 10-, 20-, and 40- fold growth beginning at 50,000 ya until the present [144, 145]. For all simulations, we used a generation time of 20 years. Because LTα and the 5′ region are different sizes, we also incorporated sequence length as a parameter in these simulations.

In addition, we applied the McDonald–Kreitman (M-K) test to the LTα coding region using DnaSP [146]. The M-K test compares the ratio of polymorphism to divergence at replacement and silent sites. Under neutrality, the ratio of replacement to silent substitutions between species (D_N/D_S) is expected to equal the ratio of replacement to silent polymorphisms within species (P_N/P_S).

Haplotype variation and inferred relationships

We extracted fully phased haplotype data from the 1000 Genomes Project for 2039 individuals using vcftools [138], and then applied a custom script to identify unique haplotypes along with the number of times that each haplotype appeared in the pooled global dataset. The genealogical relationships among haplotypes were inferred using the median-joining algorithm implemented in the Network 5.0 program [38]. The resulting phylogeny was a tree with the minimum number of changes among all possible trees [38]. Haplotype diversity (h-diversity) was also calculated for each population using the DnaSP software [146].

Linkage disequilibrium

We examined pairwise LD: 1) in the LTα region spanning ~ 10 kbs, and separately 2) in an ~ 35.5-kb genomic region, encompassing neighboring genes (NFKBIL1, LTα, TNFα, and LTβ), using the Haploview software [40]. LD was quantified using the D′ statistic [147], which indicates the magnitude of LD between SNP loci based on allele frequency. To maximize our power to detect a relationship between SNP loci, we filtered out SNPs with a MAF < 1% in our analyses. In the resulting LD plots, the color intensity of each square represents the strength of the relationship between SNP alleles. Specifically, bright red squares indicate complete LD between SNP pairs (D′ = 100; logarithm of odds (LOD) > 2); shades of pink/white squares signify little evidence of LD (D′ < 100; LOD < 2); purple squares denote high LD but with little statistical support (low LOD) [40]. The bold triangles in the plots also indicate strong blocks of LD between SNP markers.

Extended haplotype homozygosity

We characterized long-range LD on chromosome 6 for each population using the iHS statistic [42], which indicates the amount of haplotype homozygosity on chromosomes carrying the derived allele compared to chromosomes with the ancestral allele. To identify outlier values, the unstandardized scores for > 340,000 SNPs across the ~ 170-Mb region of chromosome 6 were normalized with the norm program implemented in the selscan package [43]. SNPs with a standardized |iHS| > 2 represent the most extreme 5% of scores. We also estimated the length of haplotype homozygosity with another haplotype-based statistic, nSL, using selscan [43]. This statistic measures haplotype lengths based on the number of segregating sites in a sample and does not depend on the recombination rate, unlike the iHS statistic, making it robust to recombination rate variation [44]. The output results were normalized following the same procedure used for iHS. To complement these scans for selection, we quantified the decay of identity of haplotypes with distance by calculating the EHH statistic using loci with extreme iHS and/or nSL scores as core SNPs [45].

Population differentiation and structure

To measure the degree of genetic divergence among global populations, we calculated average F_ST derived from genomic data using the Weir and Cockerham method implemented in vcftools [138, 148]. We also computed among-population F_ST at individual polymorphic sites across the ~ 10-kb region of interest and at ~ 7.5 million randomly selected SNPs from across the genome. The observed F_ST estimates for SNPs in and/or near LTα were then compared with the empirical distribution of F_ST values derived from genome-wide SNPs in order to identify outlier values (≤ 2.5th percentile or ≥ 95th percentile of the distribution). F_ST values at individual sites in or near LTα, and F_ST values at the ~ 7.5 million genome-wide polymorphisms were calculated using the same 21 populations.

Age estimates of mutations

We used a coalescent-based approach to estimate: 1) the expected TMRCA of the gene tree and 2) the expected ages of individual polymorphisms. This method requires an outgroup that provides the ancestral state at each polymorphic site. We applied the GENETREE software [48] to our sequence data to obtain the maximum likelihood estimate (MLE) of θ over 1,000,000 runs [48]. Using the MLE of θ and our estimate of μ, we calculated the effective population size (N_e) parameter, which was determined to be 18,367 based on the formula N_e = θ/4μ. Then, the TMRCA of the gene tree and the ages of individual polymorphisms were estimated from the weighted average of simulated ages over 100,000 independent runs [48]. GENETREE assumes an infinite alleles model and no recombination [48]. We removed haplotypes and/or sites that were not compatible with these underlying assumptions before applying the GENETREE algorithm.

Functional analysis

The effects of nonsynonymous SNPs at LTα on protein function were inferred using the SIFT algorithm implemented in the Ensembl Variant Effect Predictor toolset. SIFT predicts the functional impact of amino acid substitutions (i.e., “tolerated” or “deleterious”) based on sequence homology, the physical properties of amino acids and multiple alignment information [49]. In addition to SIFT, we executed GERP++ [49], phyloP [50] and phastCons [51] to evaluate whether or not common missense variation was conserved across 20 different mammalian species [149]. We also predicted changes in LTα protein stability caused by common missense mutations using two bioinformatic tools: FoldX [53] and I-Mutant3 [54] tools. ΔΔG is the difference in the Gibbs free energy for folding between the final state (the mutant) and the reference state (the wild-type). For each polymorphic site, ΔΔG was computed from the free energy of the wild type protein (encoded by the major allele at a given site) minus the Gibbs free energy of the mutated protein (encoded by the minor allele at a given site). In other words, ΔΔG = ΔG WT – ΔG MUT. Given this definition, mutations resulting in negative ΔΔG values were classified as “destabilizing”, while mutations resulting in positive ΔΔG values were categorized as “stabilizing”. We constructed the protein structure of LTα, which served as the input for FoldX, with the homology modelling server SWISS-MODEL [150]. The resulting LTα protein structure only consisted of amino acids from positions 48 to 80, which encompassed rs2229092 (H51P) and rs1041981 (T60N), but excluded rs2229094 (C13R). Therefore, the ΔΔG could not be calculated for rs2229094 using the FoldX tool.

We also analyzed polymorphisms in and near the LTα gene using the SNP2TFBS tool [55], which predicts if SNPs in TFBSs affect transcription factor binding in the human genome. The effect of a given SNP on transcription factor (TF) binding is estimated based on a position weight matrix (PWM) model for the binding specificity of the corresponding TF factor [55, 56]. In addition, using the NCBI ClinVar database [57], we searched for any reports of an association between variants (inferred to be functional in Table 3) and human health outcomes.