Evolutionary and functional analyses of LRP5 in Neanderthals, Denisovans and anatomically modern humans

Evolutionary and functional analyses of LRP5 in Neanderthals, Denisovans and anatomically modern humans | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Evolutionary and functional analyses of LRP5 in Neanderthals, Denisovans and anatomically modern humans Neus Roca-Ayats, Iago Maceda, Carlos David Bruque, Núria Martínez-Gil, and 7 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-3921272/v1 This work is licensed under a CC BY 4.0 License Status: Under Review Version 1 posted 7 You are reading this latest preprint version Abstract Background The human lineage has suffered a skeleton gracilization compared to other primates and archaic populations such as the Neanderthals. This gracilization has been traditionally explained by differences in the mechanical load that our ancestors exercised. However, there is growing evidence that gracilization could be also genetically determined. Results We have analyzed the LRP5 gene from an evolutionary and functional point of view, taking advantage of the published genomes of archaic Homo populations. Mutations in LRP5 are involved in high bone mineral density conditions. Our results suggest that this gene has a complex evolutionary history both between archaic and anatomically modern humans and within the anatomically modern human populations. In particular, we identified the presence of different selective pressures in archaics and anatomically modern humans, as well as evidence of positive selection in the African and South East Asian populations from the 1000G. Furthermore, we observed limited evidence of archaic introgression in this gene at haplotypes of East Asian ancestry, compatible with a general clearing of the archaic introgression due to functional differences in archaics compared to anatomically modern humans. In agreement with this hypothesis, we observed private mutations in the archaic genomes that we experimentally validated as putatively increasing high bone mineral density. In particular, four of five archaic missense mutations affecting the first β-propeller of LRP5 displayed enhanced Wnt pathway activation, of which two also displayed reduced negative regulation. Conclusions In summary, these data suggest a genetic component contributing to the understanding of skeletal differences between anatomically modern humans and archaic Homo populations. LRP5 bone mineral density Neanderthal Denisovan human evolution archaic introgression Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Introduction During the last three million years, the average bone strength has decayed on the Homo lineage, leading to a gracilization of the human skeleton compared to other primates ( 1 ) and to the Neanderthals. Bone strength is determined, in part, by bone mineral density (BMD) ( 2 ), which shows constitutive differences among human populations ( 3 ). In modern humans, BMD is a highly heritable trait, with up to 80% of the variance explained by genetic factors ( 4 ). So far, genome wide association studies (GWAS) have identified a large number of genomic regions associated with BMD variability, overall explaining 20% of the total estimated genetic variance ( 5 ). Rare monogenic forms of osteoporosis and high BMD provide another source for understanding the molecular pathways of BMD determination, highlighting genes and genetic variants with a significant impact on the BMD phenotype ( 6 ). Furthermore, several monogenic conditions characterized by high bone mass (HBM) are associated with mutations in the Wnt/β-catenin pathway ( 7 ), a major bone anabolic pathway. Low-density lipoprotein (LDL) receptor-related protein 5 gene ( LRP5 ), encoding its co-receptor, can be considered as one of the key genes regulating bone mass. It was one of the first to show association with BMD in the general population ( 8 ) and is always one of the top hits in GWAS. Additionally, and significantly, it bears rare variants producing extreme BMD phenotypes. The first LRP5 mutation causing HBM (p.G171V) was described in 2002 ( 9 ). Since then, several other heterozygous missense mutations have been described leading to the same phenotype. These are gain-of-function mutations, which result in a stimulation of osteoblastic bone formation ( 7 ). All HBM-associated LRP5 mutations identified are located in exons 2, 3 and 4, which collectively code for the first β-propeller domain of the protein, and reduce LRP5 binding affinity for the inhibitors sclerostin and DKK1 protein ( 10 ). In contrast, LRP5 mutations causing the osteoporosis pseudoglioma syndrome are scattered throughout the gene and are loss-of-function variants. Whereas the genetic architecture of BMD in Homo sapiens is becoming unraveled, little is known in archaic species such as Neanderthals, and the impact of introgressed BMD variants in modern populations. Analyses of the partition of Neanderthal heritability in BMD-related traits have shown an enrichment of introgressed variants in individuals of European ancestry. Nevertheless, introgressed variants are not directed toward increasing or decreasing BMD variants ( 11 ). It has also been reported that in modern human ancient samples, archaic ancestry decreased over time, particularly in areas near genes, and this observation has been diversely interpreted as evidence of hybrid sterility or a consequence of differences in effective population sizes between modern humans and Neanderthals (see ( 12 )). In this work we studied the genetic variation of LRP5 in anatomically modern humans (AMH) as well as in currently known archaic populations. We show that there is little evidence of archaic introgression in LRP5 and we identified the presence of different selective pressures in archaics and AMH. In addition, we studied five LRP5 variants found in archaic genomes, both structurally and functionally, showing that some of them might confer HBM since they cause changes in theoretical electrostatic charges and protein stability, they increase Wnt pathway activity and are less inhibited by Dkk1 protein, compared to wild-type modern human LRP5. Material and methods Databases 1000 genomes project Polymorphism information of AMH was retrieved from the publicly available VCF files from 1000 Genomes Project phase 3 ( http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ ) ( 13 ). This dataset was filtered to obtain polymorphic biallelic SNVs inside the region of the LRP5 gene according to Gencode v39 ± 500 kb. The ancestral allele was extracted from the 1000G VCF files (AA flag in the INFO field). The Altai ( 14 ), Denisovan ( 15 ), Vindija ( 16 ) and Chagyrskaya 8 ( 17 ) genomes were downloaded in VCF format from the Max Planck Institute for Evolutionary Anthropology ftp site ( http://cdna.eva.mpg.de/neandertal ). We selected polymorphisms that fall inside the region previously described and performed the same data filtering as before. The final dataset was created by merging the five previous datasets and we dropped out all SNVs that were not genotyped in all present samples. All filterings and the final merging were performed using the bcftools program ( 18 ). Base Conservation phyloP ( http://compgen.cshl.edu/phast/ ) ( 19 ) conservation scores compared to the expected under neutral drift for each nucleotide in LRP5 were retrieved from the Vertebrate Multiz Alignment & Conservation (46 Species) in UCSC ( http://hgdownload.soe.ucsc.edu/goldenPath/hg19/phyloP46way/ ). Assessment of Archaic introgressed regions in the LRP5 locus The presence of archaic introgressed regions in the LRP5 region were tested by two different approaches. First, we used the Sprime scores calculated for the 1000G individuals available in Browning, S. ( 20 ). We selected those SNVs present inside the LRP5 gene boundaries. For every individual we calculated the number of alleles predicted as introgressed by Sprime. We used all the individuals present in the CEU, CHS and CHB populations. We computed the number of differences between each haplotype from the 1,000 genomes project and the Altai Neanderthal and Denisovan genomes using the phased data from the VCF files from the 1000G ( 13 ). Second, we analyzed the presence of archaic introgressed SNPs identified in genomes from 27,566 Icelanders ( 21 ). Evidences of positive selection in LRP5 in AMH Statistics of positive selection were retrieved from PopHumanScan ( 22 ): π, the average number of nucleotide differences per site ( 23 ); Fay and Wu’s H ( 24 ); α, the proportion of substitutions that are adaptive ( 25 ); and iHS, based on the frequency of alleles in regions of high LD ( 26 ). Each of them was computed at PopHumanScan on windows of 100 kb for each population of the 1000G, with the exception of iHS, that was computed using windows of 10 kb. We used BioMart ( 27 ) to retrieve autosomal genes of the human genome. In order to estimate a single value of each statistic for each gene and population, we computed the amount of shared fragments between each gene and the PopHumanScan database, and estimated a weighted mean. Standardization of each statistic for each population at the LRP5 gene was conducted by constructing an empirical gene distribution estimated from a set of 1,098 autosomal genes with a similar length as LRP5 (136.69 ± 20 kb). Statistical analyses SMACOF Analysis of the patterns of positive selection in 1000G at LRP5 gene A Euclidean distance matrix between pairs of populations from the 1000G was computed using the standardized values of each statistic of positive selection. The relationship between the different populations from the 1000G using the patterns of selection from the considered statistics was projected in two dimensions using ordinal SMACOF ( 28 ). Weighted Multidimensional Scaling (wMDS) on AMH haplotypes and archaic genotypes at LRP5 In order to visualize the relationship between sequenced archaic individuals and haplotypes from individuals from the 1000G at the LRP5 gene, we computed an identical by state (IBS) distance between pair of haplotypes (in the case of comparing two haplotypes from the 1000G), haplotype and scaled genotype (in the case of comparing 1000G individuals and archaic) and scaled genotypes (in the case of comparing two archaic individuals). Given the unequal sample size of each continent and the archaic samples, in order to prevent biases in the estimated relationships by the MDS by sample size, we weighted the MDS using the function wcmdscale from the R package vegan ( 29 ), so the cluster of all AMH had the same weight as each archaic sample. Analysis of the degree of conservation of derived alleles present either in the AMH lineage or the archaic lineage We extracted the ancestral state of each SNP present at the LRP5 gene from the 1000G VCF. In order to test whether SNPs in the AMH lineage occurred more often at evolutionary conserved genomic regions than SNPs that occurred in the archaic lineage, for each continent we sampled at random without replacement 1000 sets of four AMH individuals, matching the number of archaic samples. For each set we identified the polymorphic SNPs where the derived allele was present only in AMH (D_AMH) and the ones present only in archaic individuals (D_ARC). For each SNP we retrieved the PhyloP score and computed the difference between the mean amount of conservation in D_AMH with regards to D_ARC. Identification and selection of Neanderthal and Denisovan LRP5 exonic variants Neanderthal and Denisovan publicly available sequencing data (UCSC Genome Browser) were used to retrieve missense variants in LRP5 , with a base quality score ≥ 23 and in a read with an alignment quality ≥ 150. Variants were filtered according to: 1) highly conserved positions; 2) damaging, according to SIFT ( 30 ) and Polyphen ( 31 ); 3) located in the HBM region of LRP5 (i.e. first β-propeller). Finally, putatively functional variants (i.e., present in more than one Neanderthal individual, affecting the same protein residue, affecting a protein residue mutated in reported human HBM cases) were selected for further analyses. The presence and frequency of the variants in AMH were assessed using the gnomAD database. Model building and assessment The sequence of the β-propeller region and EGF-LIKE 1 domain of LRP5 from UniProt (O75197-1 (NP_002326.2)) was used to perform a sequence identity search in the PDB database. Five LRP6 templates were evaluated to make the model (PDB IDs: 3S94, 3SOB, 3SOQ, 3SOV, 4DG6 ( 32 – 34 ); (Supplementary Table S1 ). The X-Ray crystallography templates with a resolution of less than 2 Å (3SOB, 3SOQ, 3SOV) were selected and the molecular homology model (MHM) was generated using MODELLER version 9.22 ( 35 ). The alignment of crystals and the LRP5 protein sequences was performed with Modeller alignment program and hand-curated in MEGA 5 software ( 36 ), (Supplementary Fig. S1 ). In addition, the model was generated with a region of 7 residues of the DKK1 protein from 3SOQ X-Ray crystallography. Model’s quality was assessed by DOPE ( 37 ), QMEANDisCo ( 38 ), and Ramachandran plots ( 39 ). Protein model is available in the Model Archive (DOI: 10.5452/ma-1smp3 ). UCSF Chimera program ( 40 ) was used for structural visualization and interpretation of the variants. In silico mutagenesis and stability calculations Protein variants were generated using FoldX 3.0 Beta 5.1 (foldx.crg.es) ( 41 ). Repair PDB command was used to optimize the total energy of the protein to FoldX's force field before residue changes were done. In silico mutagenesis was carried out using the BuildModel command, and each mutation was calculated five times. Protein interaction between LRP5 and DKK1 was calculated using the interaction command and protein stabilities, using Stability command. ∆∆G values were estimated as the difference between the energy of the wild type protein and the average of five replicas for each protein variant. A threshold of 1.6 kcal/mol was considered, as it corresponds to twice the standard deviation calculated with FoldX. Cell culture The Saos-2 cell line was used for luciferase reporter assays. It was obtained from the American Type Culture Collection (ATCC® htb-85™) and grown in Dulbecco’s Modified Eagle Medium (DMEM; Sigma-Aldrich), with 10% Fetal Bovine Serum (Gibco, Life Technologies) and 1% penicillin/streptomycin (Gibco, Life Technologies), at 37ºC and 5% of CO 2 . Plasmids and site-directed mutagenesis The pGL3-OT luciferase reporter construct, the Wnt1 -V5, mesdc2, LRP5 and DKK1 -FLAG expression vectors ( 42 ) were used. The LRP5 mutations p.A67T, p.A67V, p.G171V (positive control), p.R186Q, p.M282R, and p.R291Q were introduced with the QuikChange Lightning Site-Directed Mutagenesis Kit (Agilent), following the manufacturer instructions. All the plasmids were validated by Sanger sequencing. In vitro luciferase reporter assay Cells were seeded at a density of 1.5x10 5 cells per well in 12-well plates. After 24h, they were transfected with 1.072 µg of total DNA per well using the FuGENE HD reagent, according to manufacturer instructions (Promega): pGL3-OT (800 ng), pRL-TK (80 ng), containing the Renilla Luciferase gene, Wnt1 -V5 (32 ng), mesdc2 (64 ng), WT or mutated LRP5 (64 ng) and, depending on the experiment, DKK1 -FLAG (32 ng). When necessary, the empty pcDNA3 vector was used to adjust the total amount of DNA transfected. Forty-eight hours after transfection, cells were rinsed with PBS and lysed. The luciferase activity was measured using a Glomax Multi + luminometer (Promega), with the Dual-Luciferase® Reporter Assay System reagents (Promega). Each experiment was performed in triplicate and was repeated 3 times. Relative luciferase units (RLU, i.e., the ratio of the firefly luciferase activity over the Renilla luciferase activity) were calculated for each individual measurement and a one-way blocked ANOVA with Tukey HSD multiple comparisons tests were performed using R software version 3.4.1 and p-values < 0.05 were considered significant. All the data was ascertained for normality, homoscedasticity and atypical data points. Results 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. 5 A). Interestingly, the p.A67 and p.M282 residues are located in the interaction region with DDK1 (Fig. 5 A). 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. 5 B). 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. 5 C). 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. 5 C, 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. 5 D, 5 E) 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. Discussion It is well established that the Homo genus has undergone a skeletal gracilization, and particularly in AMH ( 1 ). Such gracilization has been mainly explained by changes in the mechanical load in AMH ( 44 ). However, GWAS studies using loci associated with BMD report differences between current human populations both in phenotype and genetics ( 5 ), raising new questions regarding which evolutionary processes, in terms of selective pressures and archaic introgression, have been at play. In the present study we focused on LRP5 , one of the key genes regulating bone mass, which has been found mutated both in high and low bone mass phenotypes in AMH ( 45 ). We analyzed it from an evolutionary point of view and studied the structure and activity of some archaic variants. The LRP5 gene appears as one of the top genes in the popHumScan showing evidence of positive selection in populations from the 1000G. In particular, Sub-Saharan African and South Asian populations show evidence of positive selective events acting on different types of genetic variation. In the case of Sub-Saharan African populations, the popHumanScan database identifies an excess of adaptive variants. In the case of South Asian populations, signals of positive selection are derived from analyzing the patterns of linkage disequilibrium. Additional evidence of positive selection has been identified in the LRP5 gene out of genes regulated by Vitamin D in East Asian populations from 1000G using frequency-spectrum-based tests ( 46 ). Overall, these results suggest that the LRP5 gene has a complex evolutionary history in human populations. It has been previously suggested that archaic introgression in allochthonous populations allows the introduction of genetic variants that have been positively selected in the archaic populations. Conversely, purifying selection could be more effectively acted against hybridization ( 12 ). This effect could be particularly important against phenotypes that are highly differentiated between AMH and archaic populations, such as BMD ( 1 ). From a genetic point of view, it has been suggested that populations out of Africa are enriched for derived low BMD-associated alleles at SNPs ascertained from GWAS compared to sub-Saharan populations, and that population phenotypic heterogeneity is the result of differential selective pressures in non-African versus Sub-Saharan populations ( 3 ). Therefore, if the LRP5 gene plays a main role in the BMD phenotype, we would expect to observe a depletion of archaic introgression in the LRP5 gene in populations out of Africa, as well as genetic variants present in archaic populations increasing the BMD. In our analyses, only three haplotypes in East Asian populations from the 1000G are suggestive of archaic introgression. This result agrees with the map of introgression based on SPrime on the same samples ( 20 ). Furthermore, Europeans from Iceland show an island of archaic introgression depletion at the LRP5 gene region, thus supporting that hybridization has not been tolerated in this genomic region. Moreover, when analyzing the sites where mutations specific to each lineage occur, we observe that AMH populations tend to accumulate mutations at positions that are highly conserved in the primate lineage as defined by phyloP statistic ( 19 ), compared to variants present in the archaic lineage. This suggests that AMH and archaic populations have been under different selective pressures. Under our hypothesis, given that high BMD is the ancestral phenotype, genetic variants modifying the function of LRP5 towards decreasing BMD would have been selectively ascertained in AMH. In contrast, genetic variants maintaining high BMD would have been selected in high BMD archaic populations. However, further directly testing this hypothesis in archaic individuals is not feasible. Nevertheless, in silico and in vitro analyses can be devised on mutations specific to the archaic lineage. Considering that several heterozygous missense variants in the first β-propeller domain of LRP5 are described to cause HBM ( 10 ), we specifically looked for mutations in this region in archaic genomes and identified 5 potential mutations that met the selection criteria. Human HBM mutations are gain-of-function changes that stimulate the Wnt pathway and reduce sclerostin and Dkk1 protein binding affinity and cell-based luciferase reporter systems have been extensively used to test them ( 42 , 47 – 49 ). Here, we took advantage of this system to find archaic mutations that similarly stimulate the Wnt pathway activity and gather in silico evidences supporting this, by creating a protein model. Two of the selected mutations (p.A67V and p.A67T) displayed the same in vitro effect as the well-known HBM p.G171V in agreement with in silico data, showing that they affect LRP5 stability. However, a loss of LRP5-DKK1 interaction was not observed in our structural model. Regarding p.R186Q and p.R291Q, our models displayed changes in surface electrostatic charges, which correlate with the luciferase results of higher pathway activation, but no effect on the DKK1 inhibition. However, we did not observe any significant change in Wnt pathway activity and DKK1 inhibition for the p.M282R mutation even though the protein model predicted a triple effect changing the electrostatic charge, destabilizing the protein and the LRP5-DKK1 interaction. Interestingly, the comparison of the change in protein stability caused by archaic variants and modern human HBM mutations does not show any statistically significant difference, which might suggest that they have similar functional consequences. Since we are modelling only small portions of LRP5 and DKK1, the discordance observed in some mutations between in vitro and in silico analyses may be explained by the fact that other LRP5 domains are involved in overall activity and DKK1 interaction, such as the third β-propeller ( 33 , 50 ). In addition, our static model might not fully represent the dynamic nature of LRP5 function. In this sense, the luciferase assay might better reflect the physiological context. On the other hand, the complexity of the assay used in this study, involving the cotransfection of several vectors, might have hindered differences in Wnt pathway activation below its sensitivity. Conclusions In conclusion, we provide data showing that LRP5 , a gene with an important role in BMD determination, follows a complex evolutionary history both within AMH and between AMHs and archaic Homo species. This evolutionary history agrees with the complexity of the evolution of the skeletal phenotype. Our structural and in vitro analyses of archaic LRP5 variants show that they resemble those causing HBM in AMH. Altogether, these data point to a genetic component that contributes to explain the skeletal differences between AMH and archaics. We might envision a set of archaic LRP5 variants which would explain their robust skeletons, that did not introgress into AMHs. Declarations Ethical Approval and consent: Ethical approval was not required for this study since it uses publicly available data and in vitro studies using commercially available cell lines. Data availability: The protein model underlying this article is available in Model Archive (DOI:10.5452/ma-1smp3). Acknowledgements and funding: NRA, NMG, MC, DG and SB acknowledge the financial support from Spanish Ministry of Science and Innovation (SAF 2016-75948R and PID2019-107188RB-C21) and Catalan Government (2017SGR:00738). NRA was a recipient of a FPU predoctoral fellowship from the Spanish Minsterio de Educación Cultura y Deporte. NMG was a recipient of a FI predoctoral fellowship from AGAUR (Catalan Agency for Management of University and Research Grants). OL and IM acknowledge the support from Spanish Ministry of Science and Innovation to the EMBL partnership, the Centro de Excelencia Severo Ochoa, CERCA Program/Generalitat de Catalunya, Spanish Ministry of Science and Innovation through the Instituto de Salud Carlos III, Generalitat de Catalunya through Departament de Salut and Departament d’Empresa i Coneixement, Co-financing with funds from the European Regional Development Fund by the Spanish Ministry of Science and Innovation corresponding to the Programa Operativo FEDER Plurirregional de España (POPE) 2014–2020 and by the Secretaria d’Universitats i Recerca, Departament d’Empresa i Coneixement of the Generalitat de Catalunya corresponding to the Programa Operatiu FEDER de Catalunya 2014–2020. OL gratefully acknowledges the financial support from Ministerio de Economía y Competitividad (Ministry of Economy and Competitiveness) PGC2018-098574-B-I00 and Generalitat de Catalunya (Government of Catalonia)—GRC 2017 SGR 937. IM gratefully acknowledges the financial support from the Government of Catalonia | Agència de Gestió d’Ajuts Universitaris i de Recerca (Agency for Management of University and Research Grants)—GRC 2014 SGR 615. Authors’ contributions: Conception: NRA, LM, OL, DG, SB Design of work: NRA, NGG, LM, WvH, OL, DG, SB Acquisition of data: NRA, IM, CDB, NMG, NGG, MC, OL Analysis of data: NRA, IM, CDB, OL Interpretation of data: NRA, CDB, LM, OL, DG, SB Writing—drafting: NRA, CDB, LM, OL, DB, SB Writing - Review & Editing: all authors Funding acquisition: OL, DG, SB All authors have read and agreed to the manuscript. Competing Interest: The authors declare no conflict of interest. References Chirchir H, Kivell TL, Ruff CB, Hublin J-J, Carlson KJ, Zipfel B, et al. Recent origin of low trabecular bone density in modern humans. Proc Natl Acad Sci USA. 2015;112(2):366–71. Kralick AE, Zemel BS. Evolutionary Perspectives on the Developing Skeleton and Implications for Lifelong Health. Front Endocrinol (Lausanne). 2020;11:513066. Medina-Gómez C, Chesi A, Heppe DHM, Zemel BS, Yin JL, Kalkwarf HJ, et al. BMD loci contribute to ethnic and developmental differences in skeletal fragility across populations: Assessment of evolutionary selection pressures. Mol Biol Evol. 2015;32(11):2961–72. Rivadeneira F, Mäkitie O. Osteoporosis and Bone Mass Disorders: From Gene Pathways to Treatments. Trends Endocrinol Metab. 2016;27(5):262–81. Morris JA, Kemp JP, Youlten SE, Laurent L, Logan JG, Chai R, et al. An atlas of genetic influences on osteoporosis in humans and mice. Nat Genet. 2019;51(2):258–66. Chen X, Yu H, Yu X. A review of the clinical, radiological and biochemical characteristics and genetic causes of high bone mass disorders. Curr Drug Targets. 2018;19. Gregson CL, Duncan EL. The Genetic Architecture of High Bone Mass. Front Endocrinol (Lausanne). 2020;11:595653. van Meurs JBJ, Trikalinos TA, Ralston SH, Balcells S, Brandi ML, Brixen K, et al. Large-Scale Analysis of Association Between LRP5 and LRP6 Variants and Osteoporosis. JAMA - J Am Med Assoc. 2008;299(11):1277–90. Little RD, Carulli JP, Del Mastro RG, Dupuis J, Osborne M, Folz C, et al. A mutation in the LDL receptor-related protein 5 gene results in the autosomal dominant high-bone-mass trait. Am J Hum Genet. 2002;70(1):11–9. Martínez-Gil N, Ugartondo N, Grinberg D, Balcells S. Wnt Pathway Extracellular Components and Their Essential Roles in Bone Homeostasis. Genes (Basel). 2022;13(1):138. McArthur E, Rinker DC, Capra JA. Quantifying the contribution of Neanderthal introgression to the heritability of complex traits. Nat Commun. 2021;12(1):4481. Ahlquist KD, Bañuelos MM, Funk A, Lai J, Rong S, Villanea FA, et al. Our Tangled Family Tree: New Genomic Methods Offer Insight into the Legacy of Archaic Admixture. Genome Biol Evol. 2021;13(7). Auton A, Abecasis GR, Altshuler DM, Durbin RM, Bentley DR, Chakravarti A, et al. A global reference for human genetic variation. Nature. 2015;526(7571):68–74. Prüfer K, Racimo F, Patterson N, Jay F, Sankararaman S, Sawyer S, et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Nature. 2014;505(7481):43–9. Meyer M, Kircher M, Gansauge MT, Li H, Racimo F, Mallick S, et al. A high-coverage genome sequence from an archaic Denisovan individual. Science. 2012;338(6104):222–6. Prüfer K, De Filippo C, Grote S, Mafessoni F, Korlević P, Hajdinjak M, et al. A high-coverage Neandertal genome from Vindija Cave in Croatia. Science. 2017;358(6363):655–8. Mafessoni F, Grote S, Filippo C De, Slon V, Kolobova KA, Viola B, et al. A high-coverage neandertal genome from chagyrskaya cave. Proc Natl Acad Sci USA. 2020;117(26):15132–6. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10(2):1–4. Pollard KS, Hubisz MJ, Rosenbloom KR, Siepel A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 2010;20(1):110–21. Browning SR, Browning BL, Zhou Y, Tucci S, Akey JM. Analysis of Human Sequence Data Reveals Two Pulses of Archaic Denisovan Admixture. Cell. 2018;173(1):53-61.e9. Skov L, Coll Macià M, Sveinbjörnsson G, Mafessoni F, Lucotte EA, Einarsdóttir MS, et al. The nature of Neanderthal introgression revealed by 27,566 Icelandic genomes. Nature. 2020;582(7810):78–83. Murga-Moreno J, Coronado-Zamora M, Bodelón A, Barbadilla A, Casillas S. PopHumanScan: the online catalog of human genome adaptation. Nucleic Acids Res. 2019;47(D1):D1080–9. Nei M, Gojobori T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986;3(5):418–26. Fay JC, Wu CI. Hitchhiking Under Positive Darwinian Selection. Genetics. 2000;155(3):1405–13. Smith NGC, Eyre-Walker A. Adaptive protein evolution in Drosophila. Nature. 2002;415(6875):1022–4. Voight BF, Kudaravalli S, Wen X, Pritchard JK. A Map of Recent Positive Selection in the Human Genome. PLOS Biol. 2006;4(3):e72. Smedley D, Haider S, Ballester B, Holland R, London D, Thorisson G, et al. BioMart - Biological queries made easy. BMC Genomics. 2009;10(1):1–12. de Leeuw J, Mair P. Multidimensional Scaling Using Majorization: SMACOF in R. J Stat Softw. 2009;31(3):1–30. Community Ecology Package [R package vegan version 2.6-4]. 2022 Kumar P, Henikoff S, Ng PC. Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat Protoc. 2009;4(7):1073–81. Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, et al. A method and server for predicting damaging missense mutations. Nat Methods. 2010;7(4):248–9. Bourhis E, Wang W, Tam C, Hwang J, Zhang Y, Spittler D, et al. Wnt antagonists bind through a short peptide to the first β-propeller domain of LRP5/6. Structure. 2011;19(10):1433–42. Cheng Z, Biechele T, Wei Z, Morrone S, Moon RT, Wang L, et al. Crystal structures of the extracellular domain of LRP6 and its complex with DKK1. Nat Struct Mol Biol. 2011;18(11):1204–10. Holdsworth G, Slocombe P, Doyle C, Sweeney B, Veverka V, Le Riche K, et al. Characterization of the interaction of sclerostin with the low density lipoprotein receptor-related protein (LRP) family of wnt co-receptors. J Biol Chem. 2012;287(32):26464–77. Šali A, Blundell TL. Comparative protein modelling by satisfaction of spatial restraints. J Mol Biol. 1993;234(3):779–815. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: Molecular Evolutionary Genetics Analysis Using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Mol Biol Evol. 2011;28(10):2731–9. Shen M, Sali A. Statistical potential for assessment and prediction of protein structures. Protein Sci. 2006;15(11):2507–24. Waterhouse A, Bertoni M, Bienert S, Studer G, Tauriello G, Gumienny R, et al. SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 2018;46(W1):W296–303. Ramachandran GN, Ramakrishnan C, Sasisekharan V. Stereochemistry of polypeptide chain configurations. J Mol Biol. 1963;7(1):95–9. Pettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera—A visualization system for exploratory research and analysis. J Comput Chem. 2004;25(13):1605–12. Schymkowitz J, Borg J, Stricher F, Nys R, Rousseau F, Serrano L. The FoldX web server: an online force field. Nucleic Acids Res. 2005;33(suppl_2):W382–8. Balemans W, Piters E, Cleiren E, Ai M, Van Wesenbeeck L, Warman ML, et al. The binding between sclerostin and LRP5 is altered by DKK1 and by high-bone mass LRP5 mutations. Calcif Tissue Int. 2008;82(6):445–53. Quintana-Murci L. Understanding rare and common diseases in the context of human evolution. Genome Biol. 2016;17(1):1–14. Ryan TM, Shaw CN. Gracility of the modern Homo sapiens skeleton is the result of decreased biomechanical loading. Proc Natl Acad Sci USA. 2015;112(2):372–7. Baron R, Kneissel M. WNT signaling in bone homeostasis and disease: from human mutations to treatments. Nat Med. 2013;19(2):179–92. Arciero E, Biagini SA, Chen Y, Xue Y, Luiselli D, Tyler-Smith C, et al. Genes Regulated by Vitamin D in Bone Cells Are Positively Selected in East Asians. PLoS One. 2015;10(12):e0146072. Patel MS, Karsenty G. Regulation of Bone Formation and Vision by LRP5. N Engl J Med. 2002;346(20):1572–4. Fenderico N, van Scherpenzeel RC, Goldflam M, Proverbio D, Jordens I, Kralj T, et al. Anti-LRP5/6 VHHs promote differentiation of Wnt-hypersensitive intestinal stem cells. Nat Commun 2019 101. 2019;10(1):1–13. Martínez-Gil N, Roca-Ayats N, Atalay N, Pineda-Moncusí M, Garcia-Giralt N, Van Hul W, et al. Functional Assessment of Coding and Regulatory Variants From the DKK1 Locus. JBMR Plus. 2020;4(12):e10423. Bourhis E, Tam C, Franke Y, Bazan JF, Ernst J, Hwang J, et al. Reconstitution of a Frizzled8·Wnt3a·LRP6 signaling complex reveals multiple Wnt and Dkk1 binding sites on LRP6. J Biol Chem. 2010;285(12):9172–9. Additional Declarations No competing interests reported. Supplementary Files SupplementaryInformationFINALHumanGenomics.docx Cite Share Download PDF Status: Under Review Version 1 posted Editorial decision: Revision requested 26 Feb, 2024 Reviews received at journal 13 Feb, 2024 Reviewers agreed at journal 09 Feb, 2024 Reviewers invited by journal 09 Feb, 2024 Editor assigned by journal 08 Feb, 2024 Submission checks completed at journal 08 Feb, 2024 First submitted to journal 02 Feb, 2024 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-3921272","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":272319266,"identity":"0b6a2020-cc11-4d09-b021-55c1cd95e1d9","order_by":0,"name":"Neus Roca-Ayats","email":"","orcid":"","institution":"IBUB, Universitat de Barcelona; Centro de Investigación Biomédica en Red de Enfermedades Raras (CIBERER) ISCIII, Barcelona, Spain; Institut de Recerca Sant Joan de Déu (IRSJD), Barcelona","correspondingAuthor":false,"prefix":"","firstName":"Neus","middleName":"","lastName":"Roca-Ayats","suffix":""},{"id":272319267,"identity":"872c9538-91ba-4caf-a59f-12f6b785c75e","order_by":1,"name":"Iago Maceda","email":"","orcid":"","institution":"CNAG, Centre Nacional d’Analisi Genòmic","correspondingAuthor":false,"prefix":"","firstName":"Iago","middleName":"","lastName":"Maceda","suffix":""},{"id":272319268,"identity":"b62dafd2-f374-4025-acd1-7d6170d962de","order_by":2,"name":"Carlos David Bruque","email":"","orcid":"","institution":"Hospital de Alta Complejidad El Calafate - S.A.M.I.C","correspondingAuthor":false,"prefix":"","firstName":"Carlos","middleName":"David","lastName":"Bruque","suffix":""},{"id":272319269,"identity":"7e8f99d6-9d92-4927-a456-5b91cdda13c1","order_by":3,"name":"Núria Martínez-Gil","email":"","orcid":"","institution":"IBUB, Universitat de Barcelona; Centro de Investigación Biomédica en Red de Enfermedades Raras (CIBERER) ISCIII, Barcelona, Spain; Institut de Recerca Sant Joan de Déu (IRSJD), Barcelona","correspondingAuthor":false,"prefix":"","firstName":"Núria","middleName":"","lastName":"Martínez-Gil","suffix":""},{"id":272319270,"identity":"1cc2cb99-0657-4f7c-b5ce-2857206b0079","order_by":4,"name":"Natàlia Garcia-Giralt","email":"","orcid":"","institution":"IMIM (Hospital del Mar Medical Research Institute), Centro de Investigación Biomédica en Red en Fragilidad y Envejecimiento Saludable (CIBERFES), ISCIII","correspondingAuthor":false,"prefix":"","firstName":"Natàlia","middleName":"","lastName":"Garcia-Giralt","suffix":""},{"id":272319271,"identity":"b2c20289-f2cc-4c34-a6aa-f36c0ce90846","order_by":5,"name":"Mónica Cozar","email":"","orcid":"","institution":"IBUB, Universitat de Barcelona; Centro de Investigación Biomédica en Red de Enfermedades Raras (CIBERER) ISCIII, Barcelona, Spain; Institut de Recerca Sant Joan de Déu (IRSJD), Barcelona","correspondingAuthor":false,"prefix":"","firstName":"Mónica","middleName":"","lastName":"Cozar","suffix":""},{"id":272319272,"identity":"57d751b8-7272-4041-bda2-dce6e260583c","order_by":6,"name":"Leonardo Mellibovsky","email":"","orcid":"","institution":"IMIM (Hospital del Mar Medical Research Institute), Centro de Investigación Biomédica en Red en Fragilidad y Envejecimiento Saludable (CIBERFES), ISCIII","correspondingAuthor":false,"prefix":"","firstName":"Leonardo","middleName":"","lastName":"Mellibovsky","suffix":""},{"id":272319274,"identity":"c9470b95-a631-4009-8f92-2f374078ee23","order_by":7,"name":"Wim Hul","email":"","orcid":"","institution":"University of Antwerp","correspondingAuthor":false,"prefix":"","firstName":"Wim","middleName":"","lastName":"Hul","suffix":""},{"id":272319276,"identity":"d314ab63-bbb3-4f46-92f8-39a8ef13be3e","order_by":8,"name":"Oscar Lao","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABDUlEQVRIiWNgGAWjYPACCwYG9sY2MJO9AUgkENYiwcDAcxCihecA0VokEtgQWvABfv7DDz9XMEjI80s+bnvwgeFwHg8D8+MPD2oY5OQbsGuRnJFmLHmGQcJw5uzEdsMZDIeLeRjYzCQSjjEYG+CwzuAGD4NkA4ME44bbiW3SPAyHE/fLv2FjSGxgSNyAw2EG588w/wRqsd9w82Cb9B+glh4GHuYPQC3183E4zOBADhvIlsQNNxjbpBkgWoDcBmCg4XAY0C9mlg0GEskze4B+6TFIh/lFwnADDi3AEHt8s6HCxraf/fizBz8qrMEh9vFHjY08rhCDOg/BSICyJPCpRwUJBFWMglEwCkbBiAMArx9SVsbkU1AAAAAASUVORK5CYII=","orcid":"","institution":"CSIC-Universitat Pompeu Fabra","correspondingAuthor":true,"prefix":"","firstName":"Oscar","middleName":"","lastName":"Lao","suffix":""},{"id":272319277,"identity":"04e1ff4d-783b-45bf-96c4-60fee24a798e","order_by":9,"name":"Daniel Grinberg","email":"","orcid":"","institution":"IBUB, Universitat de Barcelona; Centro de Investigación Biomédica en Red de Enfermedades Raras (CIBERER) ISCIII, Barcelona, Spain; Institut de Recerca Sant Joan de Déu (IRSJD), Barcelona","correspondingAuthor":false,"prefix":"","firstName":"Daniel","middleName":"","lastName":"Grinberg","suffix":""},{"id":272319278,"identity":"39e76c1b-d5ef-4687-a091-069b8cdd2882","order_by":10,"name":"Susanna Balcells","email":"","orcid":"","institution":"IBUB, Universitat de Barcelona; Centro de Investigación Biomédica en Red de Enfermedades Raras (CIBERER) ISCIII, Barcelona, Spain; Institut de Recerca Sant Joan de Déu (IRSJD), Barcelona","correspondingAuthor":false,"prefix":"","firstName":"Susanna","middleName":"","lastName":"Balcells","suffix":""}],"badges":[],"createdAt":"2024-02-02 14:47:20","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-3921272/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-3921272/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":51007074,"identity":"f91126dc-da09-4d4e-aa99-6ce046173e7d","added_by":"auto","created_at":"2024-02-12 15:29:05","extension":"jpg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":252536,"visible":true,"origin":"","legend":"\u003cp\u003eRelationship between the 1000G populations established by means of a non-metric ordinal SMACOF analysis using statistics of positive selection computed at the \u003cem\u003eLRP5\u003c/em\u003e gene.\u003c/p\u003e","description":"","filename":"Figure1.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/a1753d3525ed861546401b5f.jpg"},{"id":51007644,"identity":"b9214599-67f0-4f9b-98f6-8da4170b8a05","added_by":"auto","created_at":"2024-02-12 15:37:05","extension":"jpg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":449759,"visible":true,"origin":"","legend":"\u003cp\u003eWeighted multidimensional scaling of archaic and AMH samples using the genetic variation present in \u003cem\u003eLRP5\u003c/em\u003e. AMH samples have been weighted so all account for one fifth of the total weight. Each archaic sample accounts for one fifth of the total weight.\u003c/p\u003e","description":"","filename":"Figure2.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/3a0b419eab93a18b82d0cbbd.jpg"},{"id":51007079,"identity":"7c96c890-3a0a-49e7-8332-eff6ce4798f1","added_by":"auto","created_at":"2024-02-12 15:29:05","extension":"jpg","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":512842,"visible":true,"origin":"","legend":"\u003cp\u003eViolin plot of the distribution of the difference in the mean amount of PhyloP conservation at SNPs occurring in the AMH lineage of each continent compared to SNPs occurring in the ANC lineage. Each distribution for each continent was generated from 1000 datasets. Each dataset was obtained by sampling at random without replacement four individuals from the considered continent to match the number of archaic individuals, estimating the SNPs that occurred in the AMH lineage or in the ANC lineage, and computing the average PhyloP level of conservation. The red line indicates the expected value if SNPs occurred at each lineage on genomic positions with the same level of conservation. A value above 0 indicates that SNPs that occurred at the AMH lineage tend to happen in more conserved regions compared to ANC.\u003c/p\u003e","description":"","filename":"Figure3.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/5ef8806b6da366bb55670310.jpg"},{"id":51007643,"identity":"3933ce55-51d5-46f8-afd3-13b8682b0a44","added_by":"auto","created_at":"2024-02-12 15:37:05","extension":"jpg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":624277,"visible":true,"origin":"","legend":"\u003cp\u003e(A) Domain structure of the LRP5 co-receptor and localization of the missense variants associated with different human skeletal diseases (i.e., Osteoporosis pseudoglioma [OPPG], osteoporosis and HBM) according to the human gene mutation database (HGMD 2023.1), together with the archaic missense variants studied here. The size of the points indicates the number of variants described in each domain: small points represent 1 variant, medium points represent 2 variants and large points represent 3 or more variants. (B) Zoom of the first β-propeller domain of the LRP5 protein (adapted from Martínez-Gil et al. 2022).\u003c/p\u003e","description":"","filename":"Figure4.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/6f7b15031d2e11999b8ffd44.jpg"},{"id":51007077,"identity":"1ba875c0-ee5a-4b8d-8744-dd0afdf325bd","added_by":"auto","created_at":"2024-02-12 15:29:05","extension":"jpg","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1091567,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eMolecular structure of LRP5.\u003c/strong\u003e A. Molecular homology model of the first β-Propeller domain (at the top, top view, and at the bottom, side view). The residues mutated in Neanderthals and Denisovan are displayed in ball model and DKK1 is shown as a blue wire. B. The substitution of Ala 67 by either Val or Thr affects the structure of the β-sheets due to steric hindrance. Residue 67 of LRP5 interacts with Asp 283 which is a key residue interacting with DKK1. C. Detail of residue 282 of LRP5 interacting with isoleucine 42 of DKK1 (top: Met 282; bottom: Arg: 282; left: Solid electrostatic surface coloring; right: ribbon and wire display). D. and E. Evaluation of the substitutions Arg 186 by Gln and Arg 291 Gln by solid electrostatic surface coloring displays a decrease in the electrostatic charge in both cases (blue corresponds to positive charge and red to negative charge).\u003c/p\u003e","description":"","filename":"Figure5.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/b266a4824a3dac9505c9c9f1.jpg"},{"id":51007823,"identity":"829a28b7-64d6-477f-aa2d-e85cde02f170","added_by":"auto","created_at":"2024-02-12 15:45:05","extension":"jpg","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":307534,"visible":true,"origin":"","legend":"\u003cp\u003eRelative luciferase activity of Wnt pathway for the endogenous pathway (empty vectors), the WT or LRP5-mutated active pathway (Wnt1, LRP5, mesd2), and the WT or LRP5-mutated inhibited pathway (Wnt1, LRP5, mesd2, DKK1), in Saos-2 cells. The white bar corresponds to the endogenous pathway, black bars correspond to WT LRP5, the dark grey bar correspond to the HBM-causing LRP5 mutation (used as positive control), light grey bars correspond to the LRP5 mutations identified in Neanderthals, blue bars correspond to the LRP5 mutation identified in Denisovan. Results are expressed as mean±SD. *p\u0026lt;0.05; **p\u0026lt;0.01, ***p\u0026lt;0.001.\u003c/p\u003e","description":"","filename":"Figure6.jpg","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/fdb723ccbbf8ebdce1b69957.jpg"},{"id":51008385,"identity":"8910879a-493b-4a2a-b648-76f7504146f0","added_by":"auto","created_at":"2024-02-12 15:53:05","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1044991,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/fd607e71-c788-44f7-8d71-a6c00444171f.pdf"},{"id":51007080,"identity":"34850960-172c-40cb-8a1f-86f49afebf75","added_by":"auto","created_at":"2024-02-12 15:29:05","extension":"docx","order_by":8,"title":"","display":"","copyAsset":false,"role":"supplement","size":689305,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryInformationFINALHumanGenomics.docx","url":"https://assets-eu.researchsquare.com/files/rs-3921272/v1/8fc912eaa983044f85f76c0d.docx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Evolutionary and functional analyses of LRP5 in Neanderthals, Denisovans and anatomically modern humans","fulltext":[{"header":"Introduction","content":"\u003cp\u003eDuring the last three million years, the average bone strength has decayed on the Homo lineage, leading to a gracilization of the human skeleton compared to other primates (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e) and to the Neanderthals. Bone strength is determined, in part, by bone mineral density (BMD) (\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e), which shows constitutive differences among human populations (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eIn modern humans, BMD is a highly heritable trait, with up to 80% of the variance explained by genetic factors (\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e). So far, genome wide association studies (GWAS) have identified a large number of genomic regions associated with BMD variability, overall explaining 20% of the total estimated genetic variance (\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e). Rare monogenic forms of osteoporosis and high BMD provide another source for understanding the molecular pathways of BMD determination, highlighting genes and genetic variants with a significant impact on the BMD phenotype (\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e). Furthermore, several monogenic conditions characterized by high bone mass (HBM) are associated with mutations in the Wnt/β-catenin pathway (\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e), a major bone anabolic pathway. Low-density lipoprotein (LDL) receptor-related protein 5 gene (\u003cem\u003eLRP5\u003c/em\u003e), encoding its co-receptor, can be considered as one of the key genes regulating bone mass. It was one of the first to show association with BMD in the general population (\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e) and is always one of the top hits in GWAS. Additionally, and significantly, it bears rare variants producing extreme BMD phenotypes. The first LRP5 mutation causing HBM (p.G171V) was described in 2002 (\u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e). Since then, several other heterozygous missense mutations have been described leading to the same phenotype. These are gain-of-function mutations, which result in a stimulation of osteoblastic bone formation (\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e). All HBM-associated \u003cem\u003eLRP5\u003c/em\u003e mutations identified are located in exons 2, 3 and 4, which collectively code for the first β-propeller domain of the protein, and reduce LRP5 binding affinity for the inhibitors sclerostin and DKK1 protein (\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e). In contrast, \u003cem\u003eLRP5\u003c/em\u003e mutations causing the osteoporosis pseudoglioma syndrome are scattered throughout the gene and are loss-of-function variants.\u003c/p\u003e \u003cp\u003eWhereas the genetic architecture of BMD in \u003cem\u003eHomo sapiens\u003c/em\u003e is becoming unraveled, little is known in archaic species such as Neanderthals, and the impact of introgressed BMD variants in modern populations. Analyses of the partition of Neanderthal heritability in BMD-related traits have shown an enrichment of introgressed variants in individuals of European ancestry. Nevertheless, introgressed variants are not directed toward increasing or decreasing BMD variants (\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e). It has also been reported that in modern human ancient samples, archaic ancestry decreased over time, particularly in areas near genes, and this observation has been diversely interpreted as evidence of hybrid sterility or a consequence of differences in effective population sizes between modern humans and Neanderthals (see (\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e)).\u003c/p\u003e \u003cp\u003eIn this work we studied the genetic variation of \u003cem\u003eLRP5\u003c/em\u003e in anatomically modern humans (AMH) as well as in currently known archaic populations. We show that there is little evidence of archaic introgression in \u003cem\u003eLRP5\u003c/em\u003e and we identified the presence of different selective pressures in archaics and AMH. In addition, we studied five \u003cem\u003eLRP5\u003c/em\u003e variants found in archaic genomes, both structurally and functionally, showing that some of them might confer HBM since they cause changes in theoretical electrostatic charges and protein stability, they increase Wnt pathway activity and are less inhibited by Dkk1 protein, compared to wild-type modern human LRP5.\u003c/p\u003e"},{"header":"Material and methods","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eDatabases\u003c/h2\u003e \u003cdiv id=\"Sec4\" class=\"Section3\"\u003e \u003ch2\u003e1000 genomes project\u003c/h2\u003e \u003cp\u003ePolymorphism information of AMH was retrieved from the publicly available VCF files from 1000 Genomes Project phase 3 (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/\u003c/span\u003e\u003cspan address=\"http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) (\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e). This dataset was filtered to obtain polymorphic biallelic SNVs inside the region of the \u003cem\u003eLRP5\u003c/em\u003e gene according to Gencode v39\u0026thinsp;\u0026plusmn;\u0026thinsp;500 kb. The ancestral allele was extracted from the 1000G VCF files (AA flag in the INFO field). The Altai (\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e), Denisovan (\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e), Vindija (\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e) and Chagyrskaya 8 (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e) genomes were downloaded in VCF format from the Max Planck Institute for Evolutionary Anthropology ftp site (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://cdna.eva.mpg.de/neandertal\u003c/span\u003e\u003cspan address=\"http://cdna.eva.mpg.de/neandertal\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). We selected polymorphisms that fall inside the region previously described and performed the same data filtering as before. The final dataset was created by merging the five previous datasets and we dropped out all SNVs that were not genotyped in all present samples. All filterings and the final merging were performed using the bcftools program (\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003eBase Conservation\u003c/h2\u003e \u003cp\u003ephyloP (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://compgen.cshl.edu/phast/\u003c/span\u003e\u003cspan address=\"http://compgen.cshl.edu/phast/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) (\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e) conservation scores compared to the expected under neutral drift for each nucleotide in \u003cem\u003eLRP5\u003c/em\u003e were retrieved from the Vertebrate Multiz Alignment \u0026amp; Conservation (46 Species) in UCSC (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttp://hgdownload.soe.ucsc.edu/goldenPath/hg19/phyloP46way/\u003c/span\u003e\u003cspan address=\"http://hgdownload.soe.ucsc.edu/goldenPath/hg19/phyloP46way/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003cb\u003eAssessment of Archaic introgressed regions in the\u003c/b\u003e \u003cb\u003eLRP5\u003c/b\u003e \u003cb\u003elocus\u003c/b\u003e\u003c/p\u003e \u003cp\u003eThe presence of archaic introgressed regions in the \u003cem\u003eLRP5\u003c/em\u003e region were tested by two different approaches. First, we used the Sprime scores calculated for the 1000G individuals available in Browning, S. (\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e). We selected those SNVs present inside the \u003cem\u003eLRP5\u003c/em\u003e gene boundaries. For every individual we calculated the number of alleles predicted as introgressed by Sprime. We used all the individuals present in the CEU, CHS and CHB populations. We computed the number of differences between each haplotype from the 1,000 genomes project and the Altai Neanderthal and Denisovan genomes using the phased data from the VCF files from the 1000G (\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e). Second, we analyzed the presence of archaic introgressed SNPs identified in genomes from 27,566 Icelanders (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003cb\u003eEvidences of positive selection in\u003c/b\u003e \u003cb\u003eLRP5\u003c/b\u003e \u003cb\u003ein AMH\u003c/b\u003e\u003c/p\u003e \u003cp\u003eStatistics of positive selection were retrieved from PopHumanScan (\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e): π, the average number of nucleotide differences per site (\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e); Fay and Wu\u0026rsquo;s H (\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e); α, the proportion of substitutions that are adaptive (\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e); and iHS, based on the frequency of alleles in regions of high LD (\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e). Each of them was computed at PopHumanScan on windows of 100 kb for each population of the 1000G, with the exception of iHS, that was computed using windows of 10 kb. We used BioMart (\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e) to retrieve autosomal genes of the human genome. In order to estimate a single value of each statistic for each gene and population, we computed the amount of shared fragments between each gene and the PopHumanScan database, and estimated a weighted mean. Standardization of each statistic for each population at the \u003cem\u003eLRP5\u003c/em\u003e gene was conducted by constructing an empirical gene distribution estimated from a set of 1,098 autosomal genes with a similar length as \u003cem\u003eLRP5\u003c/em\u003e (136.69\u0026thinsp;\u0026plusmn;\u0026thinsp;20 kb).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003eStatistical analyses\u003c/h2\u003e \u003cp\u003e \u003cb\u003eSMACOF Analysis of the patterns of positive selection in 1000G at\u003c/b\u003e \u003cb\u003eLRP5\u003c/b\u003e \u003cb\u003egene\u003c/b\u003e\u003c/p\u003e \u003cp\u003eA Euclidean distance matrix between pairs of populations from the 1000G was computed using the standardized values of each statistic of positive selection. The relationship between the different populations from the 1000G using the patterns of selection from the considered statistics was projected in two dimensions using ordinal SMACOF (\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003cb\u003eWeighted Multidimensional Scaling (wMDS) on AMH haplotypes and archaic genotypes at\u003c/b\u003e \u003cb\u003eLRP5\u003c/b\u003e\u003c/p\u003e \u003cp\u003eIn order to visualize the relationship between sequenced archaic individuals and haplotypes from individuals from the 1000G at the \u003cem\u003eLRP5\u003c/em\u003e gene, we computed an identical by state (IBS) distance between pair of haplotypes (in the case of comparing two haplotypes from the 1000G), haplotype and scaled genotype (in the case of comparing 1000G individuals and archaic) and scaled genotypes (in the case of comparing two archaic individuals).\u003c/p\u003e \u003cp\u003eGiven the unequal sample size of each continent and the archaic samples, in order to prevent biases in the estimated relationships by the MDS by sample size, we weighted the MDS using the function \u003cem\u003ewcmdscale\u003c/em\u003e from the R package vegan (\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e), so the cluster of all AMH had the same weight as each archaic sample.\u003c/p\u003e \u003cp\u003e \u003cb\u003eAnalysis of the degree of conservation of derived alleles present either in the AMH lineage or the archaic lineage\u003c/b\u003e \u003c/p\u003e \u003cp\u003eWe extracted the ancestral state of each SNP present at the \u003cem\u003eLRP5\u003c/em\u003e gene from the 1000G VCF. In order to test whether SNPs in the AMH lineage occurred more often at evolutionary conserved genomic regions than SNPs that occurred in the archaic lineage, for each continent we sampled at random without replacement 1000 sets of four AMH individuals, matching the number of archaic samples. For each set we identified the polymorphic SNPs where the derived allele was present only in AMH (D_AMH) and the ones present only in archaic individuals (D_ARC). For each SNP we retrieved the PhyloP score and computed the difference between the mean amount of conservation in D_AMH with regards to D_ARC.\u003c/p\u003e \u003cp\u003e \u003cb\u003eIdentification and selection of Neanderthal and Denisovan\u003c/b\u003e \u003cb\u003eLRP5\u003c/b\u003e \u003cb\u003eexonic variants\u003c/b\u003e\u003c/p\u003e \u003cp\u003eNeanderthal and Denisovan publicly available sequencing data (UCSC Genome Browser) were used to retrieve missense variants in \u003cem\u003eLRP5\u003c/em\u003e, with a base quality score\u0026thinsp;\u0026ge;\u0026thinsp;23 and in a read with an alignment quality\u0026thinsp;\u0026ge;\u0026thinsp;150. Variants were filtered according to: 1) highly conserved positions; 2) damaging, according to SIFT (\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e) and Polyphen (\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e); 3) located in the HBM region of \u003cem\u003eLRP5\u003c/em\u003e (i.e. first β-propeller). Finally, putatively functional variants (i.e., present in more than one Neanderthal individual, affecting the same protein residue, affecting a protein residue mutated in reported human HBM cases) were selected for further analyses. The presence and frequency of the variants in AMH were assessed using the gnomAD database.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec7\" class=\"Section2\"\u003e \u003ch2\u003eModel building and assessment\u003c/h2\u003e \u003cp\u003eThe sequence of the β-propeller region and EGF-LIKE 1 domain of LRP5 from UniProt (O75197-1 (NP_002326.2)) was used to perform a sequence identity search in the PDB database. Five LRP6 templates were evaluated to make the model (PDB IDs: 3S94, 3SOB, 3SOQ, 3SOV, 4DG6 (\u003cspan additionalcitationids=\"CR33\" citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e); (Supplementary Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). The X-Ray crystallography templates with a resolution of less than 2 \u0026Aring; (3SOB, 3SOQ, 3SOV) were selected and the molecular homology model (MHM) was generated using MODELLER version 9.22 (\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e). The alignment of crystals and the LRP5 protein sequences was performed with Modeller alignment program and hand-curated in MEGA 5 software (\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e), (Supplementary Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). In addition, the model was generated with a region of 7 residues of the DKK1 protein from 3SOQ X-Ray crystallography. Model\u0026rsquo;s quality was assessed by DOPE (\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e), QMEANDisCo (\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e), and Ramachandran plots (\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e). Protein model is available in the Model Archive (DOI:\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.5452/ma-1smp3\u003c/span\u003e\u003cspan address=\"10.5452/ma-1smp3\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). UCSF Chimera program (\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e) was used for structural visualization and interpretation of the variants.\u003c/p\u003e \u003cp\u003e \u003cb\u003eIn silico\u003c/b\u003e \u003cb\u003emutagenesis and stability calculations\u003c/b\u003e\u003c/p\u003e \u003cp\u003eProtein variants were generated using FoldX 3.0 Beta 5.1 (foldx.crg.es) (\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e). Repair PDB command was used to optimize the total energy of the protein to FoldX's force field before residue changes were done. \u003cem\u003eIn silico\u003c/em\u003e mutagenesis was carried out using the BuildModel command, and each mutation was calculated five times. Protein interaction between LRP5 and DKK1 was calculated using the interaction command and protein stabilities, using Stability command. ∆∆G values were estimated as the difference between the energy of the wild type protein and the average of five replicas for each protein variant. A threshold of 1.6 kcal/mol was considered, as it corresponds to twice the standard deviation calculated with FoldX.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eCell culture\u003c/h2\u003e \u003cp\u003eThe Saos-2 cell line was used for luciferase reporter assays. It was obtained from the American Type Culture Collection (ATCC\u0026reg; htb-85\u0026trade;) and grown in Dulbecco\u0026rsquo;s Modified Eagle Medium (DMEM; Sigma-Aldrich), with 10% Fetal Bovine Serum (Gibco, Life Technologies) and 1% penicillin/streptomycin (Gibco, Life Technologies), at 37\u0026ordm;C and 5% of CO\u003csub\u003e2\u003c/sub\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003ePlasmids and site-directed mutagenesis\u003c/h2\u003e \u003cp\u003eThe pGL3-OT luciferase reporter construct, the \u003cem\u003eWnt1\u003c/em\u003e-V5, \u003cem\u003emesdc2, LRP5\u003c/em\u003e and \u003cem\u003eDKK1\u003c/em\u003e-FLAG expression vectors (\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e) were used. The \u003cem\u003eLRP5\u003c/em\u003e mutations p.A67T, p.A67V, p.G171V (positive control), p.R186Q, p.M282R, and p.R291Q were introduced with the QuikChange Lightning Site-Directed Mutagenesis Kit (Agilent), following the manufacturer instructions. All the plasmids were validated by Sanger sequencing.\u003c/p\u003e \u003cp\u003e \u003cb\u003eIn vitro\u003c/b\u003e \u003cb\u003eluciferase reporter assay\u003c/b\u003e\u003c/p\u003e \u003cp\u003eCells were seeded at a density of 1.5x10\u003csup\u003e5\u003c/sup\u003e cells per well in 12-well plates. After 24h, they were transfected with 1.072 \u0026micro;g of total DNA per well using the FuGENE HD reagent, according to manufacturer instructions (Promega): pGL3-OT (800 ng), pRL-TK (80 ng), containing the Renilla Luciferase gene, \u003cem\u003eWnt1\u003c/em\u003e-V5 (32 ng), \u003cem\u003emesdc2\u003c/em\u003e (64 ng), WT or mutated \u003cem\u003eLRP5\u003c/em\u003e (64 ng) and, depending on the experiment, \u003cem\u003eDKK1\u003c/em\u003e-FLAG (32 ng). When necessary, the empty pcDNA3 vector was used to adjust the total amount of DNA transfected. Forty-eight hours after transfection, cells were rinsed with PBS and lysed. The luciferase activity was measured using a Glomax Multi\u0026thinsp;+\u0026thinsp;luminometer (Promega), with the Dual-Luciferase\u0026reg; Reporter Assay System reagents (Promega). Each experiment was performed in triplicate and was repeated 3 times. Relative luciferase units (RLU, i.e., the ratio of the firefly luciferase activity over the Renilla luciferase activity) were calculated for each individual measurement and a one-way blocked ANOVA with Tukey HSD multiple comparisons tests were performed using R software version 3.4.1 and p-values\u0026thinsp;\u0026lt;\u0026thinsp;0.05 were considered significant. All the data was ascertained for normality, homoscedasticity and atypical data points.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cstrong\u003eEvidence of differential selective pressures in\u003c/strong\u003e \u003cstrong\u003eLRP5\u003c/strong\u003e \u003cstrong\u003ewithin AMH\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFirst, using ordinal SMACOF, we projected in two dimensions the relationships between the 1000G populations using ascertained statistics of positive selection computed at the \u003cem\u003eLRP5\u003c/em\u003e gene (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e). 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.\u003c/p\u003e\n\u003cp\u003eSupporting 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).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIdentification of Archaic introgression in AMH at the\u003c/strong\u003e \u003cstrong\u003eLRP5\u003c/strong\u003e \u003cstrong\u003elocus\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNext, 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 (\u003cspan class=\"CitationRef\"\u003e21\u003c/span\u003e), \u003cem\u003eLRP5\u003c/em\u003e 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\u0026thinsp;=\u0026thinsp;0.0001). Analysis of a map of introgression of 1000G based on SPrime (\u003cspan class=\"CitationRef\"\u003e20\u003c/span\u003e) 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.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e).\u003c/p\u003e\n\u003cp\u003eThe 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.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEvidences of different selective pressures in\u003c/strong\u003e \u003cstrong\u003eLRP5\u003c/strong\u003e \u003cstrong\u003ein AMH and archaic populations\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eGiven 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.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003e) 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\u0026thinsp;\u0026gt;\u0026thinsp;mean conserved D_ANC) in AFR\u0026thinsp;\u0026lt;\u0026thinsp;0.005, EUR\u0026thinsp;=\u0026thinsp;0.003, EAS\u0026thinsp;=\u0026thinsp;0.014, AMR\u0026thinsp;=\u0026thinsp;0.014 and SAS\u0026thinsp;=\u0026thinsp;0.102).\u003c/p\u003e\n\u003cp\u003eGiven that highly conserved regions tend to be associated with deleterious effects (\u003cspan class=\"CitationRef\"\u003e43\u003c/span\u003e), 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.\u003c/p\u003e\n\u003cp\u003eOverall, all these results support a complex recent evolution of \u003cem\u003eLRP5\u003c/em\u003e, with different selective pressures acting on archaic and AMH populations.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIdentification of putatively functional variants in\u003c/strong\u003e \u003cstrong\u003eLRP5\u003c/strong\u003e \u003cstrong\u003ein Neanderthals and Denisovans\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConsidering 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 \u003cem\u003eLRP5\u003c/em\u003e 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\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003e, Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e4\u003c/span\u003e): all of them are located in the first \u0026beta;-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.\u003c/p\u003e\n\u003cdiv class=\"gridtable\"\u003e\n\u003ctable id=\"Tab1\" border=\"1\"\u003e\u003ccaption\u003e\n\u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e\n\u003cdiv class=\"CaptionContent\"\u003e\n\u003cp\u003eArchaic \u003cem\u003eLRP5\u003c/em\u003e variants analyzed in this work\u003c/p\u003e\n\u003c/div\u003e\n\u003c/caption\u003e\n\u003cthead\u003e\n\u003ctr\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eGenomic position (GRCh37)\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eVariant\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eProtein effect\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003egnomAD frequency\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eSIFT\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003ePolyphen\u003c/p\u003e\n\u003c/th\u003e\n\u003cth align=\"left\"\u003e\n\u003cp\u003eIndividual\u003c/p\u003e\n\u003c/th\u003e\n\u003c/tr\u003e\n\u003c/thead\u003e\n\u003ctbody\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003echr11:68115422\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eG\u0026thinsp;\u0026gt;\u0026thinsp;A\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003e\u003cstrong\u003ep.A67T\u003c/strong\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\"\u0026minus;\"\u003e\n\u003cp\u003e7.08\u0026middot;10\u003csup\u003e\u0026minus;\u0026thinsp;6\u003c/sup\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.0256\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e1.000\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eVi33.26\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003echr11:68115423\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eC\u0026thinsp;\u0026gt;\u0026thinsp;T\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003e\u003cstrong\u003ep.A67V\u003c/strong\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.001\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e1.000\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eVi33.16\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003echr11:68125186\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eG\u0026thinsp;\u0026gt;\u0026thinsp;A\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003e\u003cstrong\u003ep.R186Q\u003c/strong\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\"\u0026minus;\"\u003e\n\u003cp\u003e4.60\u0026middot;10\u003csup\u003e\u0026minus;\u0026thinsp;5\u003c/sup\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.000\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e1.000\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eVi33.16, Vi33.25\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003echr11:68131373\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eT\u0026thinsp;\u0026gt;\u0026thinsp;G\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003e\u003cstrong\u003ep.M282R\u003c/strong\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\u0026nbsp;\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.039\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.999\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eMez1\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003ctr\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003echr11:68131400\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eG\u0026thinsp;\u0026gt;\u0026thinsp;A\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003e\u003cstrong\u003ep.R291Q\u003c/strong\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\"\u0026minus;\"\u003e\n\u003cp\u003e1.22\u0026middot;10\u003csup\u003e\u0026minus;\u0026thinsp;5\u003c/sup\u003e\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.012\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"char\" char=\".\"\u003e\n\u003cp\u003e0.995\u003c/p\u003e\n\u003c/td\u003e\n\u003ctd align=\"left\"\u003e\n\u003cp\u003eDenisovan\u003c/p\u003e\n\u003c/td\u003e\n\u003c/tr\u003e\n\u003c/tbody\u003e\n\u003c/table\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec11\" class=\"Section2\"\u003e\n\u003cp\u003e\u003cstrong\u003eStructure-based functional analyses of the impact of LRP5 variants\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo 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 \u0026beta;-propeller of LRP5 in interaction with DKK1 was generated (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eA). Interestingly, the p.A67 and p.M282 residues are located in the interaction region with DDK1 (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eA). 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.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eB). In addition, the variants affect the structure of the \u0026beta;-sheets due to steric hindrance and cause changes in the stability of the protein [\u0026Delta;\u0026Delta;G\u0026thinsp;=\u0026thinsp;0.91\u0026thinsp;\u0026plusmn;\u0026thinsp;0.02 kcal/mol (below the threshold of 1.6 kcal/mol, see methods) for the p.A67T and \u0026Delta;\u0026Delta;G\u0026thinsp;=\u0026thinsp;5.77\u0026thinsp;\u0026plusmn;\u0026thinsp;0.10 kcal/mol for p.A67V]. In either case, no significant changes in interaction with DKK1 are observed (Supplementary Table S2).\u003c/p\u003e\n\u003cp\u003eThe substitution of Met by Arg at position 282 causes three possible effects. On the one hand, a change in the surface electrostatic charge (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eC). 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 \u0026Delta;\u0026Delta;G is 2.7 kcal/mol (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eC, Supplementary Table S2). Finally, the protein stability \u0026Delta;\u0026Delta;G is 6.16\u0026thinsp;\u0026plusmn;\u0026thinsp;2.5 kcal/mol for p.M282R.\u003c/p\u003e\n\u003cp\u003eVariants p.R186Q and p.R291Q cause a change in the surface electrostatic charge (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eD, \u003cspan class=\"InternalRef\"\u003e5\u003c/span\u003eE) but do not affect protein stability (Supplementary Table S2).\u003c/p\u003e\n\u003cp\u003eFinally, we compared the protein stability changes (\u0026Delta;\u0026Delta;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).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIn vitro\u003c/strong\u003e \u003cstrong\u003efunctional analysis of the impact of archaic LRP5 variants\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eIn 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\u0026thinsp;=\u0026thinsp;0.0013), 1.78 (p-value\u0026thinsp;=\u0026thinsp;2.59\u0026middot;10\u003csup\u003e\u0026minus;\u0026thinsp;10\u003c/sup\u003e), 1.55 (p-value\u0026thinsp;=\u0026thinsp;2.59\u0026middot;10\u003csup\u003e\u0026minus;\u0026thinsp;10\u003c/sup\u003e), and 1.18 (p-value\u0026thinsp;=\u0026thinsp;0.0284), respectively, similar to the p.G171V variant (FC: 2.07; p-value\u0026thinsp;=\u0026thinsp;2.59\u0026middot;10\u003csup\u003e\u0026minus;\u0026thinsp;10\u003c/sup\u003e), used as positive control (Fig.\u0026nbsp;\u003cspan class=\"InternalRef\"\u003e6\u003c/span\u003e). 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.\u003c/p\u003e\n\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eIt is well established that the \u003cem\u003eHomo\u003c/em\u003e genus has undergone a skeletal gracilization, and particularly in AMH (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e). Such gracilization has been mainly explained by changes in the mechanical load in AMH (\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e). However, GWAS studies using loci associated with BMD report differences between current human populations both in phenotype and genetics (\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e), raising new questions regarding which evolutionary processes, in terms of selective pressures and archaic introgression, have been at play.\u003c/p\u003e \u003cp\u003eIn the present study we focused on \u003cem\u003eLRP5\u003c/em\u003e, one of the key genes regulating bone mass, which has been found mutated both in high and low bone mass phenotypes in AMH (\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e). We analyzed it from an evolutionary point of view and studied the structure and activity of some archaic variants.\u003c/p\u003e \u003cp\u003eThe \u003cem\u003eLRP5\u003c/em\u003e gene appears as one of the top genes in the popHumScan showing evidence of positive selection in populations from the 1000G. In particular, Sub-Saharan African and South Asian populations show evidence of positive selective events acting on different types of genetic variation. In the case of Sub-Saharan African populations, the popHumanScan database identifies an excess of adaptive variants. In the case of South Asian populations, signals of positive selection are derived from analyzing the patterns of linkage disequilibrium. Additional evidence of positive selection has been identified in the \u003cem\u003eLRP5\u003c/em\u003e gene out of genes regulated by Vitamin D in East Asian populations from 1000G using frequency-spectrum-based tests (\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e). Overall, these results suggest that the \u003cem\u003eLRP5\u003c/em\u003e gene has a complex evolutionary history in human populations. It has been previously suggested that archaic introgression in allochthonous populations allows the introduction of genetic variants that have been positively selected in the archaic populations. Conversely, purifying selection could be more effectively acted against hybridization (\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e). This effect could be particularly important against phenotypes that are highly differentiated between AMH and archaic populations, such as BMD (\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e). From a genetic point of view, it has been suggested that populations out of Africa are enriched for derived low BMD-associated alleles at SNPs ascertained from GWAS compared to sub-Saharan populations, and that population phenotypic heterogeneity is the result of differential selective pressures in non-African versus Sub-Saharan populations (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e). Therefore, if the \u003cem\u003eLRP5\u003c/em\u003e gene plays a main role in the BMD phenotype, we would expect to observe a depletion of archaic introgression in the \u003cem\u003eLRP5\u003c/em\u003e gene in populations out of Africa, as well as genetic variants present in archaic populations increasing the BMD.\u003c/p\u003e \u003cp\u003eIn our analyses, only three haplotypes in East Asian populations from the 1000G are suggestive of archaic introgression. This result agrees with the map of introgression based on SPrime on the same samples (\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e). Furthermore, Europeans from Iceland show an island of archaic introgression depletion at the \u003cem\u003eLRP5\u003c/em\u003e gene region, thus supporting that hybridization has not been tolerated in this genomic region. Moreover, when analyzing the sites where mutations specific to each lineage occur, we observe that AMH populations tend to accumulate mutations at positions that are highly conserved in the primate lineage as defined by phyloP statistic (\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e), compared to variants present in the archaic lineage. This suggests that AMH and archaic populations have been under different selective pressures. Under our hypothesis, given that high BMD is the ancestral phenotype, genetic variants modifying the function of \u003cem\u003eLRP5\u003c/em\u003e towards decreasing BMD would have been selectively ascertained in AMH. In contrast, genetic variants maintaining high BMD would have been selected in high BMD archaic populations. However, further directly testing this hypothesis in archaic individuals is not feasible. Nevertheless, \u003cem\u003ein silico\u003c/em\u003e and \u003cem\u003ein vitro\u003c/em\u003e analyses can be devised on mutations specific to the archaic lineage. Considering that several heterozygous missense variants in the first β-propeller domain of LRP5 are described to cause HBM (\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e), we specifically looked for mutations in this region in archaic genomes and identified 5 potential mutations that met the selection criteria.\u003c/p\u003e \u003cp\u003eHuman HBM mutations are gain-of-function changes that stimulate the Wnt pathway and reduce sclerostin and Dkk1 protein binding affinity and cell-based luciferase reporter systems have been extensively used to test them (\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e, \u003cspan additionalcitationids=\"CR48\" citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e). Here, we took advantage of this system to find archaic mutations that similarly stimulate the Wnt pathway activity and gather \u003cem\u003ein silico\u003c/em\u003e evidences supporting this, by creating a protein model. Two of the selected mutations (p.A67V and p.A67T) displayed the same \u003cem\u003ein vitro\u003c/em\u003e effect as the well-known HBM p.G171V in agreement with \u003cem\u003ein silico\u003c/em\u003e data, showing that they affect LRP5 stability. However, a loss of LRP5-DKK1 interaction was not observed in our structural model. Regarding p.R186Q and p.R291Q, our models displayed changes in surface electrostatic charges, which correlate with the luciferase results of higher pathway activation, but no effect on the DKK1 inhibition. However, we did not observe any significant change in Wnt pathway activity and DKK1 inhibition for the p.M282R mutation even though the protein model predicted a triple effect changing the electrostatic charge, destabilizing the protein and the LRP5-DKK1 interaction. Interestingly, the comparison of the change in protein stability caused by archaic variants and modern human HBM mutations does not show any statistically significant difference, which might suggest that they have similar functional consequences.\u003c/p\u003e \u003cp\u003eSince we are modelling only small portions of LRP5 and DKK1, the discordance observed in some mutations between \u003cem\u003ein vitro\u003c/em\u003e and \u003cem\u003ein silico\u003c/em\u003e analyses may be explained by the fact that other LRP5 domains are involved in overall activity and DKK1 interaction, such as the third β-propeller (\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e, \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e). In addition, our static model might not fully represent the dynamic nature of LRP5 function. In this sense, the luciferase assay might better reflect the physiological context. On the other hand, the complexity of the assay used in this study, involving the cotransfection of several vectors, might have hindered differences in Wnt pathway activation below its sensitivity.\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eIn conclusion, we provide data showing that \u003cem\u003eLRP5\u003c/em\u003e, a gene with an important role in BMD determination, follows a complex evolutionary history both within AMH and between AMHs and archaic \u003cem\u003eHomo\u003c/em\u003e species. This evolutionary history agrees with the complexity of the evolution of the skeletal phenotype. Our structural and \u003cem\u003ein vitro\u003c/em\u003e analyses of archaic \u003cem\u003eLRP5\u003c/em\u003e variants show that they resemble those causing HBM in AMH. Altogether, these data point to a genetic component that contributes to explain the skeletal differences between AMH and archaics. We might envision a set of archaic \u003cem\u003eLRP5\u003c/em\u003e variants which would explain their robust skeletons, that did not introgress into AMHs.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eEthical Approval and consent:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eEthical approval was not required for this study since it uses publicly available data and \u003cem\u003ein vitro\u003c/em\u003e studies using commercially available cell lines.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eData availability:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe protein model underlying this article is available in Model Archive (DOI:10.5452/ma-1smp3).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgements and funding:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eNRA, NMG, MC, DG and SB acknowledge the financial support from Spanish Ministry of Science and Innovation (SAF 2016-75948R and PID2019-107188RB-C21) and Catalan Government (2017SGR:00738). NRA was a recipient of a FPU predoctoral fellowship from the Spanish Minsterio de Educaci\u0026oacute;n Cultura y Deporte. NMG was a recipient of a FI predoctoral fellowship from AGAUR (Catalan Agency for Management of University and Research Grants). OL and IM acknowledge the support from Spanish Ministry of Science and Innovation to the EMBL partnership, the Centro de Excelencia Severo Ochoa, CERCA Program/Generalitat de Catalunya, Spanish Ministry of Science and Innovation through the Instituto de Salud Carlos III, Generalitat de Catalunya through Departament de Salut and Departament d\u0026rsquo;Empresa i Coneixement, Co-financing with funds from the European Regional Development Fund by the Spanish Ministry of Science and Innovation corresponding to the Programa Operativo FEDER Plurirregional de Espa\u0026ntilde;a (POPE) 2014\u0026ndash;2020 and by the Secretaria d\u0026rsquo;Universitats i Recerca, Departament d\u0026rsquo;Empresa i Coneixement of the Generalitat de Catalunya corresponding to the Programa Operatiu FEDER de Catalunya 2014\u0026ndash;2020. OL gratefully acknowledges the financial support from Ministerio de Econom\u0026iacute;a y Competitividad (Ministry of Economy and Competitiveness) PGC2018-098574-B-I00 and Generalitat de Catalunya (Government of Catalonia)\u0026mdash;GRC 2017 SGR 937. IM gratefully acknowledges the financial support from the Government of Catalonia | Ag\u0026egrave;ncia de Gesti\u0026oacute; d\u0026rsquo;Ajuts Universitaris i de Recerca (Agency for Management of University and Research Grants)\u0026mdash;GRC 2014 SGR 615.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026rsquo; contributions:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eConception: NRA, LM, OL, DG, SB\u003c/p\u003e\n\u003cp\u003eDesign of work: NRA, NGG, LM, WvH, OL, DG, SB\u003c/p\u003e\n\u003cp\u003eAcquisition of data: NRA, IM, CDB, NMG, NGG, MC, OL\u003c/p\u003e\n\u003cp\u003eAnalysis of data: NRA, IM, CDB, OL\u003c/p\u003e\n\u003cp\u003eInterpretation of data: NRA, CDB, LM, OL, DG, SB\u003c/p\u003e\n\u003cp\u003eWriting\u0026mdash;drafting: NRA, CDB, LM, OL, DB, SB\u003c/p\u003e\n\u003cp\u003eWriting - Review \u0026amp; Editing: all authors\u003c/p\u003e\n\u003cp\u003eFunding acquisition: OL, DG, SB\u003c/p\u003e\n\u003cp\u003eAll authors have read and agreed to the manuscript.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting Interest:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare no conflict of interest.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eChirchir H, Kivell TL, Ruff CB, Hublin J-J, Carlson KJ, Zipfel B, et al. Recent origin of low trabecular bone density in modern humans. Proc Natl Acad Sci USA. 2015;112(2):366\u0026ndash;71. \u003c/li\u003e\n\u003cli\u003eKralick AE, Zemel BS. Evolutionary Perspectives on the Developing Skeleton and Implications for Lifelong Health. Front Endocrinol (Lausanne). 2020;11:513066. \u003c/li\u003e\n\u003cli\u003eMedina-G\u0026oacute;mez C, Chesi A, Heppe DHM, Zemel BS, Yin JL, Kalkwarf HJ, et al. BMD loci contribute to ethnic and developmental differences in skeletal fragility across populations: Assessment of evolutionary selection pressures. Mol Biol Evol. 2015;32(11):2961\u0026ndash;72. \u003c/li\u003e\n\u003cli\u003eRivadeneira F, M\u0026auml;kitie O. Osteoporosis and Bone Mass Disorders: From Gene Pathways to Treatments. Trends Endocrinol Metab. 2016;27(5):262\u0026ndash;81. \u003c/li\u003e\n\u003cli\u003eMorris JA, Kemp JP, Youlten SE, Laurent L, Logan JG, Chai R, et al. An atlas of genetic influences on osteoporosis in humans and mice. Nat Genet. 2019;51(2):258\u0026ndash;66. \u003c/li\u003e\n\u003cli\u003eChen X, Yu H, Yu X. A review of the clinical, radiological and biochemical characteristics and genetic causes of high bone mass disorders. Curr Drug Targets. 2018;19. \u003c/li\u003e\n\u003cli\u003eGregson CL, Duncan EL. The Genetic Architecture of High Bone Mass. Front Endocrinol (Lausanne). 2020;11:595653. \u003c/li\u003e\n\u003cli\u003evan Meurs JBJ, Trikalinos TA, Ralston SH, Balcells S, Brandi ML, Brixen K, et al. Large-Scale Analysis of Association Between LRP5 and LRP6 Variants and Osteoporosis. JAMA - J Am Med Assoc. 2008;299(11):1277\u0026ndash;90. \u003c/li\u003e\n\u003cli\u003eLittle RD, Carulli JP, Del Mastro RG, Dupuis J, Osborne M, Folz C, et al. A mutation in the LDL receptor-related protein 5 gene results in the autosomal dominant high-bone-mass trait. Am J Hum Genet. 2002;70(1):11\u0026ndash;9. \u003c/li\u003e\n\u003cli\u003eMart\u0026iacute;nez-Gil N, Ugartondo N, Grinberg D, Balcells S. Wnt Pathway Extracellular Components and Their Essential Roles in Bone Homeostasis. Genes (Basel). 2022;13(1):138. \u003c/li\u003e\n\u003cli\u003eMcArthur E, Rinker DC, Capra JA. Quantifying the contribution of Neanderthal introgression to the heritability of complex traits. Nat Commun. 2021;12(1):4481. \u003c/li\u003e\n\u003cli\u003eAhlquist KD, Ba\u0026ntilde;uelos MM, Funk A, Lai J, Rong S, Villanea FA, et al. Our Tangled Family Tree: New Genomic Methods Offer Insight into the Legacy of Archaic Admixture. Genome Biol Evol. 2021;13(7). \u003c/li\u003e\n\u003cli\u003eAuton A, Abecasis GR, Altshuler DM, Durbin RM, Bentley DR, Chakravarti A, et al. A global reference for human genetic variation. Nature. 2015;526(7571):68\u0026ndash;74. \u003c/li\u003e\n\u003cli\u003ePr\u0026uuml;fer K, Racimo F, Patterson N, Jay F, Sankararaman S, Sawyer S, et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Nature. 2014;505(7481):43\u0026ndash;9. \u003c/li\u003e\n\u003cli\u003eMeyer M, Kircher M, Gansauge MT, Li H, Racimo F, Mallick S, et al. A high-coverage genome sequence from an archaic Denisovan individual. Science. 2012;338(6104):222\u0026ndash;6. \u003c/li\u003e\n\u003cli\u003ePr\u0026uuml;fer K, De Filippo C, Grote S, Mafessoni F, Korlević P, Hajdinjak M, et al. A high-coverage Neandertal genome from Vindija Cave in Croatia. Science. 2017;358(6363):655\u0026ndash;8. \u003c/li\u003e\n\u003cli\u003eMafessoni F, Grote S, Filippo C De, Slon V, Kolobova KA, Viola B, et al. A high-coverage neandertal genome from chagyrskaya cave. Proc Natl Acad Sci USA. 2020;117(26):15132\u0026ndash;6. \u003c/li\u003e\n\u003cli\u003eDanecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10(2):1\u0026ndash;4. \u003c/li\u003e\n\u003cli\u003ePollard KS, Hubisz MJ, Rosenbloom KR, Siepel A. Detection of nonneutral substitution rates on mammalian phylogenies. Genome Res. 2010;20(1):110\u0026ndash;21. \u003c/li\u003e\n\u003cli\u003eBrowning SR, Browning BL, Zhou Y, Tucci S, Akey JM. Analysis of Human Sequence Data Reveals Two Pulses of Archaic Denisovan Admixture. Cell. 2018;173(1):53-61.e9. \u003c/li\u003e\n\u003cli\u003eSkov L, Coll Maci\u0026agrave; M, Sveinbj\u0026ouml;rnsson G, Mafessoni F, Lucotte EA, Einarsd\u0026oacute;ttir MS, et al. The nature of Neanderthal introgression revealed by 27,566 Icelandic genomes. Nature. 2020;582(7810):78\u0026ndash;83. \u003c/li\u003e\n\u003cli\u003eMurga-Moreno J, Coronado-Zamora M, Bodel\u0026oacute;n A, Barbadilla A, Casillas S. PopHumanScan: the online catalog of human genome adaptation. Nucleic Acids Res. 2019;47(D1):D1080\u0026ndash;9. \u003c/li\u003e\n\u003cli\u003eNei M, Gojobori T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986;3(5):418\u0026ndash;26. \u003c/li\u003e\n\u003cli\u003eFay JC, Wu CI. Hitchhiking Under Positive Darwinian Selection. Genetics. 2000;155(3):1405\u0026ndash;13. \u003c/li\u003e\n\u003cli\u003eSmith NGC, Eyre-Walker A. Adaptive protein evolution in Drosophila. Nature. 2002;415(6875):1022\u0026ndash;4. \u003c/li\u003e\n\u003cli\u003eVoight BF, Kudaravalli S, Wen X, Pritchard JK. A Map of Recent Positive Selection in the Human Genome. PLOS Biol. 2006;4(3):e72. \u003c/li\u003e\n\u003cli\u003eSmedley D, Haider S, Ballester B, Holland R, London D, Thorisson G, et al. BioMart - Biological queries made easy. BMC Genomics. 2009;10(1):1\u0026ndash;12. \u003c/li\u003e\n\u003cli\u003ede Leeuw J, Mair P. Multidimensional Scaling Using Majorization: SMACOF in R. J Stat Softw. 2009;31(3):1\u0026ndash;30. \u003c/li\u003e\n\u003cli\u003eCommunity Ecology Package [R package vegan version 2.6-4]. 2022 \u003c/li\u003e\n\u003cli\u003eKumar P, Henikoff S, Ng PC. Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat Protoc. 2009;4(7):1073\u0026ndash;81. \u003c/li\u003e\n\u003cli\u003eAdzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, et al. A method and server for predicting damaging missense mutations. Nat Methods. 2010;7(4):248\u0026ndash;9. \u003c/li\u003e\n\u003cli\u003eBourhis E, Wang W, Tam C, Hwang J, Zhang Y, Spittler D, et al. Wnt antagonists bind through a short peptide to the first \u0026beta;-propeller domain of LRP5/6. Structure. 2011;19(10):1433\u0026ndash;42. \u003c/li\u003e\n\u003cli\u003eCheng Z, Biechele T, Wei Z, Morrone S, Moon RT, Wang L, et al. Crystal structures of the extracellular domain of LRP6 and its complex with DKK1. Nat Struct Mol Biol. 2011;18(11):1204\u0026ndash;10. \u003c/li\u003e\n\u003cli\u003eHoldsworth G, Slocombe P, Doyle C, Sweeney B, Veverka V, Le Riche K, et al. Characterization of the interaction of sclerostin with the low density lipoprotein receptor-related protein (LRP) family of wnt co-receptors. J Biol Chem. 2012;287(32):26464\u0026ndash;77. \u003c/li\u003e\n\u003cli\u003e\u0026Scaron;ali A, Blundell TL. Comparative protein modelling by satisfaction of spatial restraints. J Mol Biol. 1993;234(3):779\u0026ndash;815. \u003c/li\u003e\n\u003cli\u003eTamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: Molecular Evolutionary Genetics Analysis Using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Mol Biol Evol. 2011;28(10):2731\u0026ndash;9. \u003c/li\u003e\n\u003cli\u003eShen M, Sali A. Statistical potential for assessment and prediction of protein structures. Protein Sci. 2006;15(11):2507\u0026ndash;24. \u003c/li\u003e\n\u003cli\u003eWaterhouse A, Bertoni M, Bienert S, Studer G, Tauriello G, Gumienny R, et al. SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 2018;46(W1):W296\u0026ndash;303. \u003c/li\u003e\n\u003cli\u003eRamachandran GN, Ramakrishnan C, Sasisekharan V. Stereochemistry of polypeptide chain configurations. J Mol Biol. 1963;7(1):95\u0026ndash;9. \u003c/li\u003e\n\u003cli\u003ePettersen EF, Goddard TD, Huang CC, Couch GS, Greenblatt DM, Meng EC, et al. UCSF Chimera\u0026mdash;A visualization system for exploratory research and analysis. J Comput Chem. 2004;25(13):1605\u0026ndash;12. \u003c/li\u003e\n\u003cli\u003eSchymkowitz J, Borg J, Stricher F, Nys R, Rousseau F, Serrano L. The FoldX web server: an online force field. Nucleic Acids Res. 2005;33(suppl_2):W382\u0026ndash;8. \u003c/li\u003e\n\u003cli\u003eBalemans W, Piters E, Cleiren E, Ai M, Van Wesenbeeck L, Warman ML, et al. The binding between sclerostin and LRP5 is altered by DKK1 and by high-bone mass LRP5 mutations. Calcif Tissue Int. 2008;82(6):445\u0026ndash;53. \u003c/li\u003e\n\u003cli\u003eQuintana-Murci L. Understanding rare and common diseases in the context of human evolution. Genome Biol. 2016;17(1):1\u0026ndash;14. \u003c/li\u003e\n\u003cli\u003eRyan TM, Shaw CN. Gracility of the modern Homo sapiens skeleton is the result of decreased biomechanical loading. Proc Natl Acad Sci USA. 2015;112(2):372\u0026ndash;7. \u003c/li\u003e\n\u003cli\u003eBaron R, Kneissel M. WNT signaling in bone homeostasis and disease: from human mutations to treatments. Nat Med. 2013;19(2):179\u0026ndash;92. \u003c/li\u003e\n\u003cli\u003eArciero E, Biagini SA, Chen Y, Xue Y, Luiselli D, Tyler-Smith C, et al. Genes Regulated by Vitamin D in Bone Cells Are Positively Selected in East Asians. PLoS One. 2015;10(12):e0146072. \u003c/li\u003e\n\u003cli\u003ePatel MS, Karsenty G. Regulation of Bone Formation and Vision by LRP5. N Engl J Med. 2002;346(20):1572\u0026ndash;4. \u003c/li\u003e\n\u003cli\u003eFenderico N, van Scherpenzeel RC, Goldflam M, Proverbio D, Jordens I, Kralj T, et al. Anti-LRP5/6 VHHs promote differentiation of Wnt-hypersensitive intestinal stem cells. Nat Commun 2019 101. 2019;10(1):1\u0026ndash;13. \u003c/li\u003e\n\u003cli\u003eMart\u0026iacute;nez-Gil N, Roca-Ayats N, Atalay N, Pineda-Moncus\u0026iacute; M, Garcia-Giralt N, Van Hul W, et al. Functional Assessment of Coding and Regulatory Variants From the DKK1 Locus. JBMR Plus. 2020;4(12):e10423. \u003c/li\u003e\n\u003cli\u003eBourhis E, Tam C, Franke Y, Bazan JF, Ernst J, Hwang J, et al. Reconstitution of a Frizzled8\u0026middot;Wnt3a\u0026middot;LRP6 signaling complex reveals multiple Wnt and Dkk1 binding sites on LRP6. J Biol Chem. 2010;285(12):9172\u0026ndash;9. \u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"human-genomics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"hugm","sideBox":"Learn more about [Human Genomics](http://humgenomics.biomedcentral.com/)","snPcode":"40246","submissionUrl":"https://submission.nature.com/new-submission/40246/3","title":"Human Genomics","twitterHandle":"@OAgenetics","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true},"keywords":"LRP5, bone mineral density, Neanderthal, Denisovan, human evolution, archaic introgression","lastPublishedDoi":"10.21203/rs.3.rs-3921272/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-3921272/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003ch2\u003eBackground\u003c/h2\u003e \u003cp\u003eThe human lineage has suffered a skeleton gracilization compared to other primates and archaic populations such as the Neanderthals. This gracilization has been traditionally explained by differences in the mechanical load that our ancestors exercised. However, there is growing evidence that gracilization could be also genetically determined.\u003c/p\u003e\u003ch2\u003eResults\u003c/h2\u003e \u003cp\u003eWe have analyzed the \u003cem\u003eLRP5\u003c/em\u003e gene from an evolutionary and functional point of view, taking advantage of the published genomes of archaic \u003cem\u003eHomo\u003c/em\u003e populations. Mutations in \u003cem\u003eLRP5\u003c/em\u003e are involved in high bone mineral density conditions. Our results suggest that this gene has a complex evolutionary history both between archaic and anatomically modern humans and within the anatomically modern human populations. In particular, we identified the presence of different selective pressures in archaics and anatomically modern humans, as well as evidence of positive selection in the African and South East Asian populations from the 1000G. Furthermore, we observed limited evidence of archaic introgression in this gene at haplotypes of East Asian ancestry, compatible with a general clearing of the archaic introgression due to functional differences in archaics compared to anatomically modern humans. In agreement with this hypothesis, we observed private mutations in the archaic genomes that we experimentally validated as putatively increasing high bone mineral density. In particular, four of five archaic missense mutations affecting the first β-propeller of LRP5 displayed enhanced Wnt pathway activation, of which two also displayed reduced negative regulation.\u003c/p\u003e\u003ch2\u003eConclusions\u003c/h2\u003e \u003cp\u003eIn summary, these data suggest a genetic component contributing to the understanding of skeletal differences between anatomically modern humans and archaic \u003cem\u003eHomo\u003c/em\u003e populations.\u003c/p\u003e","manuscriptTitle":"Evolutionary and functional analyses of LRP5 in Neanderthals, Denisovans and anatomically modern humans","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-02-12 15:29:00","doi":"10.21203/rs.3.rs-3921272/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2024-02-26T05:55:37+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2024-02-14T00:25:52+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"87b97e1d-9a9a-4d42-8a93-e372c02ca220","date":"2024-02-09T07:27:11+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2024-02-09T05:04:41+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2024-02-09T03:40:54+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2024-02-09T03:40:53+00:00","index":"","fulltext":""},{"type":"submitted","content":"Human Genomics","date":"2024-02-02T14:45:36+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"human-genomics","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"hugm","sideBox":"Learn more about [Human Genomics](http://humgenomics.biomedcentral.com/)","snPcode":"40246","submissionUrl":"https://submission.nature.com/new-submission/40246/3","title":"Human Genomics","twitterHandle":"@OAgenetics","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"BMC/SO AJ","inReviewEnabled":true,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"dff72370-1d2f-414d-bbaf-c3d4e67d8fd6","owner":[],"postedDate":"February 12th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"under-review","subjectAreas":[],"tags":[],"updatedAt":"2024-05-07T07:11:45+00:00","versionOfRecord":[],"versionCreatedAt":"2024-02-12 15:29:00","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-3921272","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-3921272","identity":"rs-3921272","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}