Evidence of differential selective pressures in LRP5 within AMH
First, using ordinal SMACOF, we projected in two dimensions the relationships between the 1000G populations using ascertained statistics of positive selection computed at the LRP5 gene (Fig. 1). 1000G populations tend to cluster according to their continental origin, particularly for the AFR populations. The second dimension tends to distinguish CHB (Han Chinese) and STU (Sri Lankan Tamil in the UK). Overall, the presence of geographic population substructure for summary statistics accounting for positive selection suggests that this gene could have been under different selective pressures among human populations.
Supporting this interpretation, popHumanScan reported evidence of genomic positive selection in Sub-Saharan populations for the statistic (Supplementary Fig. S2A), accounting for the proportion of substitutions that are adaptive, and elevated values of iHS, summarizing (recent) departures in the allelic frequency given the observed haplotype length, in East Asian and, particularly, South Asian populations (Supplementary Fig. S2B).
Identification of Archaic introgression in AMH at the LRP5 locus
Next, we analysed whether the presence of differential selective pressures in AMH could be explained by archaic introgression. First, we checked reported maps of archaic introgression in AMH. For the first map of introgression, generated from 27,566 Icelandic genomes (21), LRP5 falls within a region of depletion of archaic introgression of 2.47 Mb, being one of the largest archaic-introgressed-free regions of the chromosome 11 (p-value = 0.0001). Analysis of a map of introgression of 1000G based on SPrime (20) supports the absence of signals of archaic introgression in populations out of Africa, with the exception of one CHS and two KHV haplotypes. In order to study this effect, we visualized the relationship between the introgressed haplotypes and the archaic populations. We constructed a genetic distance matrix between pairs of individuals using IBS. A weighted multidimensional scaling (wMDS) was run with this distance matrix by assigning the same weight to each of the four archaic samples, and dividing between all the 1000G samples the remaining weight (Fig. 2).
The first dimension (49.52% of explained variance) of the wMDS distinguishes the Denisovan sample against AMH and Neanderthal samples. The second dimension (36.42% of explained variance) distinguishes Altai Neanderthal against the rest. Interestingly, Vindija and Chagyrskaya cluster together and appear between AMH and Altai. Moreover, three haplotypes corresponding to CHS (Southern Han Chinese) and KHV (Kinh Vietnamese) populations appear as outliers from the AMH points, and closer to Vindija and Chagyrskaya.
Evidences of different selective pressures in LRP5 in AMH and archaic populations
Given the previous results, we wondered to which extent archaic populations and AMH populations showed evidence of different selective pressures. We used the map of nucleotide conservation among mammals (PhyloP) and the information of the ancestral allele of each SNP identified in AMH and archaic populations as defined in the 1000G to estimate the amount of conservation of SNPs that had appeared in the AMH genome compared to SNPs that appeared in the archaic populations. Our results (Fig. 3) show that SNPs that appeared (i.e., the derived allele is found) in the AMH lineage tend to occur in more conserved regions than SNPs from the archaic (ANC) lineage for all continental groups, except SAS (P(mean conserved D_AMH > mean conserved D_ANC) in AFR < 0.005, EUR = 0.003, EAS = 0.014, AMR = 0.014 and SAS = 0.102).
Given that highly conserved regions tend to be associated with deleterious effects (43), these results would support the presence of different pressures acting on archaic populations compared to AMH populations. In particular, stronger purifying selection would affect more archaic populations and/or relaxation on the selective pressures in the AMH lineage.
Overall, all these results support a complex recent evolution of LRP5, with different selective pressures acting on archaic and AMH populations.
Identification of putatively functional variants in LRP5 in Neanderthals and Denisovans
Considering the archaeological results together with the divergence between the archaic and AMH, BMD-increasing variants in the archaic genomes and absent in AMH might be expected. We searched available archaic LRP5 genomic sequences and identified four missense variants in Neanderthals and one in the Denisovan individual (p.R291Q), all having a suggestive evidence of functionality (Table 1, Fig. 4): all of them are located in the first β-propeller; two of them (p.A67T and p.A67V) result in a change of the same protein residue but towards a different amino acid; a third one ( p.R186Q) was found in two different Neanderthal individuals; and a fourth (p.M282R) affects a protein residue also mutated in human HBM cases.
Table 1
Archaic LRP5 variants analyzed in this work
Genomic position (GRCh37)
|
Variant
|
Protein effect
|
gnomAD frequency
|
SIFT
|
Polyphen
|
Individual
|
chr11:68115422
|
G > A
|
p.A67T
|
7.08·10− 6
|
0.0256
|
1.000
|
Vi33.26
|
chr11:68115423
|
C > T
|
p.A67V
|
|
0.001
|
1.000
|
Vi33.16
|
chr11:68125186
|
G > A
|
p.R186Q
|
4.60·10− 5
|
0.000
|
1.000
|
Vi33.16, Vi33.25
|
chr11:68131373
|
T > G
|
p.M282R
|
|
0.039
|
0.999
|
Mez1
|
chr11:68131400
|
G > A
|
p.R291Q
|
1.22·10− 5
|
0.012
|
0.995
|
Denisovan
|
Structure-based functional analyses of the impact of LRP5 variants
To determine the possible effect of the identified variants affecting residues p.A67, p.R186, p.M282 and p.R291, a protein homology model of the first β-propeller of LRP5 in interaction with DKK1 was generated (Fig. 5A). Interestingly, the p.A67 and p.M282 residues are located in the interaction region with DDK1 (Fig. 5A). For residue 67, we observe changes in the distances to p.D283, p.T80 and p.L113 in the mutated residues (Val or Thr), compared to the wild type (Ala),greater for Val than for Thr (Fig. 5B). In addition, the variants affect the structure of the β-sheets due to steric hindrance and cause changes in the stability of the protein [ΔΔG = 0.91 ± 0.02 kcal/mol (below the threshold of 1.6 kcal/mol, see methods) for the p.A67T and ΔΔG = 5.77 ± 0.10 kcal/mol for p.A67V]. In either case, no significant changes in interaction with DKK1 are observed (Supplementary Table S2).
The substitution of Met by Arg at position 282 causes three possible effects. On the one hand, a change in the surface electrostatic charge (Fig. 5C). Secondly, the atomic distances between the 282 residue of LRP5 and the I42 of DKK1 are longer with Arg than with Met, and the LRP5-DKK1 interaction ΔΔG is 2.7 kcal/mol (Fig. 5C, Supplementary Table S2). Finally, the protein stability ΔΔG is 6.16 ± 2.5 kcal/mol for p.M282R.
Variants p.R186Q and p.R291Q cause a change in the surface electrostatic charge (Fig. 5D, 5E) but do not affect protein stability (Supplementary Table S2).
Finally, we compared the protein stability changes (ΔΔG) of archaic variants with those of the HBM variants described in AMH and we did not observe any statistically significant difference between them (Supplementary Fig. S3).
In vitro functional analysis of the impact of archaic LRP5 variants
In order to evaluate the effect of the variants on the canonical Wnt pathway activity, we performed a luciferase reporter assay. Four of the variants (p.A67T, p.A67V, p.R186Q, and p.R291Q) displayed significantly greater Wnt pathway stimulation, compared to WT, with fold changes of 1.26 (p-value = 0.0013), 1.78 (p-value = 2.59·10− 10), 1.55 (p-value = 2.59·10− 10), and 1.18 (p-value = 0.0284), respectively, similar to the p.G171V variant (FC: 2.07; p-value = 2.59·10− 10), used as positive control (Fig. 6). Moreover, for two of these variants (p.A67T and p.A67V) DKK1 failed to significantly inhibit Wnt pathway activation, similarly to p.G171V. No significant differences were observed between p.M282R and WT, either in the Wnt pathway activation or in the DKK1 inhibition.