Complex genetic architecture underlies human hand and foot evolution

Complex genetic architecture underlies human hand and foot evolution | 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 Article Complex genetic architecture underlies human hand and foot evolution Alexander Okamoto, Gayani Senevirathne, Pushpanathan Muthuirulan, and 4 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-7124496/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 11 May, 2026 Read the published version in Proceedings of the National Academy of Sciences → Version 1 posted You are reading this latest preprint version Abstract The transition to bipedal locomotion is a key event in human evolution, involving substantial changes to the skeleton, including the bones of the hands and feet (autopods). Hominins evolved a more muscular and opposable thumb while the other fingers are relatively shorter, enhancing manipulative capacity. The feet evolved robust first toes and short lateral toes to meet the challenges of bipedal walking and running. While adaptations in the hand and foot have often been considered separately, the fore- and hind limbs of primates are morphologically integrated, serially homologous structures, raising the possibility that natural selection on either autopod may have driven corresponding changes in the other. To explore the genetic architecture underlying human autopod evolution, we used functional genomics methods to identify regulatory elements and gene expression patterns in the developing phalanges and metacarpals of the human hand and foot. We find that gene expression and regulation differ along the proximal-distal axis and between timepoints but not between limb types or individual digits. We show that thousands of human-specific genomic features fall within autopod regulatory elements, some accessible in multiple tissues, others with tissue-specific accessibility. Our results highlight the complex genomic basis of human autopod evolution. Biological sciences/Evolution/Anthropology/Biological anthropology Biological sciences/Evolution/Evolutionary developmental biology Biological sciences/Genetics/Development Biological sciences/Genetics/Functional genomics Figures Figure 1 Figure 2 Figure 3 Figure 4 Main The functional divergence of the hand and foot for prehension and locomotion, respectively, has long been recognized as a critical event in human evolution 1 . In comparison to other apes, the human hand and foot have features in skeletal anatomy as well as the related musculature and innervation which are believed to have evolved adaptively in the context of their respective functional specializations 2 – 5 . Skeletally, the human hand evolved a large pollex (i.e., thumb) and a reduced length of the other digits, which improved manipulative capacity, a key factor in successful tool use for food acquisition and processing, and defense 6 – 10 . The human foot evolved a robust, adducted hallux (i.e., big toe) and short digits, creating a shorter, hinge-like forefoot with a stiff midfoot that improves performance in both walking and running 11 – 15 . These derived features contribute to the unique status of humans as the only extant primates that are obligate bipeds. Due to their differences in structure and function, the evolutionary histories of the human hand and foot have largely been considered separately 3 , 6 – 14 , 16 – 18 . However, data from both developmental genetics and hominoid anatomical studies provide evidence that the evolution of the hand and the foot should be considered together: hands and feet (hereafter collectively referred to as autopods) are serially homologous structures, share much of their genetic architecture 19 , and have been shown to covary strongly in size and shape 20 , 21 . Furthermore, when the anatomical changes between human and chimpanzee autopod skeletons are compared for both limb types, it becomes clear that some changes are remarkably similar – e.g., a larger, more robust first digit compared to the other digits, shorter phalanges compared to metapodials, loss of phalangeal curvature, and an overall size reduction (Fig. 1a-c). This raises the possibility that phenotypic changes in either the hand or foot might in part be explained by adaptive evolution in the other due to the potentially constraining effects of strong covariation between homologous limb elements of humans and other apes 20 – 22 . Simulations even suggest that selection for a modern human-like foot from a chimp-like ancestor is sufficient to produce a human-like hand purely as a byproduct due to anatomical integration of the autopods 23 . Despite the importance of the hand and foot in human evolution, very little is known about the genetic mechanisms patterning human-specific features of the hand and foot skeleton. Many of these features arise early in human development when the limb skeleton is first patterned in cartilage templates (Fig. 1d). While some studies have investigated the developmental genetics of human limbs during or slightly before this developmental window 24 , 25 , none have generated matched datasets for the hand and foot, nor for individual skeletal elements of each autopod. Therefore, to shed new light on the evolutionary history of the human hand and foot skeleton, we generated data on the gene expression and regulation of each metapodial and the pooled phalanges of each digit in both the hand and foot at two stages during this developmental window. Additionally, we performed corresponding experiments in stage-matched E15.5 mouse embryos on digits I, III, and V to explore evolutionary conservation of identified genes and regulatory elements. Using these datasets, we investigated the signature of natural selection on genomic regions patterning distinct elements of the human hand and foot skeleton. We tested (1) the prediction that there would be substantial overlap in gene expression and the genomic regulatory landscape between the hand and foot tissues (the genetic underpinnings of covariance and coevolution), and (2) that the foot would show stronger signatures of selection than the hand, in line with previous hypotheses and the fact the remodeling of the foot in the context of bipedalism was more extensive and predates the appearance of lithic technology 6 . Proportional changes in the autopod skeleton To quantitatively assess anatomical differences in proportion across the autopod skeleton, we first estimated the average percent difference in the length and width of each element in a large dataset of adult human, chimpanzee, and gorilla skeletons (See Methods) (Supplementary Table 1). As expected, the greatest amount of difference was observed in the more robust first digit in humans when compared with chimps (Fig. 1b-c). The other digits were generally shorter – especially the foot phalanges – with the exception of the metacarpals, which were similar in length between the two species (Figs. 1b-c). The same pattern is observed when humans are compared to gorillas (Extended Data Fig. 1). Patterns of gene expression and regulation across the autopod The identified proportional differences most likely reflect modifications to autopod development in each species, especially at the cartilage level. This is because the rudiments of each digital element are comprised of chondrocytes that organize to form growth plates governing longitudinal and transverse growth during ontogeny. We examined human foot and hand development and identified two gestational timepoints (E54 and E67) spanning a window where morphogenesis was occurring rapidly (Fig. 1d, Supplementary Note). For each digital ray of the human hand and foot, we generated matched transcriptomic (RNA-seq) and epigenomic (ATAC-seq) datasets for the phalangeal and metapodial regions separately (Supplementary Tables 2–3, Supplementary Results S1-3). For comparative purposes, we performed a similar approach on stage-matched mouse E15.5 forelimb and hind limb elements, focusing on digits I, III, and V (Extended Data Fig. 2, Tables S3-5, Supplementary Results S4-5). Analysis of our human RNA-seq dataset identified 1,263 differentially expressed genes (DEGs) between biologically relevant sets of autopod tissues (See Methods; Supplementary Table 2). These comparisons included adjacent tissues (along either the anterior-posterior or proximal-distal axes), homologous tissues between the hand and foot, and the same tissue across timepoints, as well as gene expression gradients across the autopod. Many DEGs separated the metapodials and phalanges but not individual digits (Fig. 2a; Extended Data Fig. 3). As predicted based on the patterns of covariance of coevolution, gene expression was overwhelmingly similar between the hand and foot, consistent with the expression patterns we observed in mouse (Extended Data Fig. 2) 26 , 27 . We also identified 3,031 DEGs between the two timepoints, with 1,453 genes upregulated early at E54 and 1,578 genes upregulated late at E67 (Fig. 2b). The later developmental time point was enriched for gene ontology (GO) terms related to “ossification,” which initiates around this stage. Using our human ATAC-seq data, we identified 60,150 regulatory elements accessible in the chondrocytes of at least one autopod skeletal element (Supplementary Table 3). Consistent with the RNA-seq data, principal component analysis (PCA) separated samples between the metapodials and phalanges and between timepoints but not by limb type (Fig. 2c). Indeed, no elements uniquely defined all hand or all foot tissues. Overall, regulatory elements are often accessible in many autopod tissues or tissue and timepoint specific (Fig. 2d). Using our previously collected stage-matched ATAC-seq datasets from human long bones, girdles, and the axial column 28 – 30 , we found that many of the identified regulatory elements show accessibility in at least one other developing skeletal tissue, albeit 12,488 regulatory elements have not been previously identified and appear unique to the autopod (Fig. 2d). To help elucidate broader regulatory patterns in functionally conserved sequences in other primates, we asked how many of the human autopod elements overlapped elements present in mouse, since elements conserved between human and mouse are likely present in other primates. While the number of individual rays examined in mouse was 3 versus 5 in humans, we found that 20,491 regulatory elements were accessible in both human and mouse (34.0% of all regulatory elements identified in human), only 2,654 of which were autopod-specific (21.3% of all human autopod-specific elements) (Extended Data Fig. 4). Genomic evolution in autopod regulatory elements Using our dataset, we sought to investigate the genomic basis of human-specific features of the autopod skeleton. To this end, we overlapped our autopod regulatory element sets with genomic regions marked by signatures of human-specific changes, including human accelerated regions (HARs) 31 – 38 , human conserved element deletions (hCONDELs) 39 , 40 , human ancestor quickly evolved regions (HAQERs), large human-specific inversions, and structurally divergent regions between humans and other apes (SDRs) 41 (Table 1, Extended Data Fig. 5–6). As structural changes can impact neighboring genomic loci by rearranging the local chromatin environment, overlaps were also determined for regions flanking each inversion and SDRs (Extended Data Fig. 7). We identified overlaps for these different types of genomic regions for elements with accessibility in multiple skeletal or autopod tissues, as well as those that are tissue and timepoint-specific (Fig. 3). Overall, 95% of the genes expressed in at least one autopod tissue fell within 500kb of one or more of these genomic regions. While all tissues exhibited accessibility for at least a few regulatory elements overlapping each type of genomic feature, this was not true for tissue-specific elements (Fig. 3f,i,j,o,r). Tissues with high levels of morphological change were not distinctly enriched for evolutionary signals, with relatively few overlaps detected in digit I but large numbers of overlaps seen in the foot phalanges of digits III and V (Fig. 3d-f,j-o). This analysis revealed that thousands of genomic loci with accessibility in the developing human autopod skeleton harbor evolutionary sequence alterations. To link these genomic changes to anatomical subdivisions of the autopod, we sought to partition our regulatory elements and overlaps according to their accessibility patterns. Based on the differences in gene expression and regulation between the phalanges and metapodials and between developmental stages, as well as the anatomical differences between the hand and foot, we explored the relative strength of the evolutionary signals along each of these axes (Fig. 4a-c). Given that evolutionary forces, such as natural selection, shape the patterns of genetic variation present within a species, we also compared the patterns of human, chimpanzee, and gorilla intraspecific variation within each regulatory set (See Methods). Enrichments within each comparison differed depending on the genomic feature and the specificity of the regulatory element set considered (i.e., all brain-filtered elements, autopod-specific elements, or autopod-specific, human-mouse conserved) (Fig. 4, Extended Data Figs. 8–9). Within the autopod-specific sets, the most significant differences in enrichments were identified for the hand versus foot comparison (Fig. 4c). For the autopod-specific elements, we sought to explore biological differences between these sets. To test whether different sets regulated the same or distinct target genes, we counted the number of regulatory elements from each set that fell within 500kb of genes expressed in at least one autopod tissue (Additional Data Fig. 10). As expected from the largely similar pattern of gene expression observed in all tissues, we found that the number of nearby regulatory elements from the shared sets were largely correlated with those of the more restricted ones. For the more specific sets, the numbers of hand and foot elements were not correlated (Pearson’s Correlation, R = 0.02, P = 0.365), while the numbers of phalangeal and metapodial elements were positively correlated (R = 0.11, P < 0.001) and the early and late elements negatively correlated, although given the strength of the correlation, this is unlikely to be biologically meaningful (R = -0.075, P < 0.001). Overall, many genes are near regulatory elements from multiple categories, consistent with complex regulatory control of gene expression. Accordingly, gene ontology enrichments for each set differed slightly but all included terms related to cartilage or bone development (Supplementary Table 3). Finally, we analyzed transcription factor (TF) binding profiles and found that the sets were differentially enriched for hundreds of TF motifs, suggesting that regulatory elements in each set preferentially interact with distinct sets of TFs (Supplementary Table 3). While hotspots of human-specific sequence changes are frequently prioritized, single nucleotide changes may also have substantial phenotypic consequences. Accordingly, we sought to estimate the total number of nucleotide changes that may have contributed to human autopod evolution. Briefly, we took our regulatory elements, removed polymorphic variant positions and misaligned regions, and counted the number of positions where the reference nucleotide in the human genome differed from that present in both the chimp and gorilla genomes (see Methods). While this approach ignores structurally divergent and misaligned regions, it gives a rough minimum estimate of the number of fixed substitutions that may contribute to human-specific differences. Using this approach, we identified 2,461,335 human-specific fixed substitutions that are accessible in the developing autopod skeleton, 291,867 of which are autopod-specific. Since humans and chimpanzees share the same last common ancestor to gorilla and therefore have evolved independently for the same length of time (i.e., the branch lengths are equal), we repeated this process using the chimpanzee branch. We identified 2,012,269 chimp-specific fixed substitutions that are accessible in the developing autopod skeleton, 236,518 of which are autopod-specific. Many of these changes likely resulted from neutral evolutionary processes and have minimal or no phenotypic effects, however, comparison of PhyloP scores between the human and chimp sets found a significant shift towards more positive scores in human (human mean 0.23, chimp mean 0.01, P < 2.2 × 10 − 16 , two-sided Wilcoxon rank-sum), suggesting that human fixed changes alter positions with deeper evolutionary constraint. This holds true for the autopod-specific set as well (human mean 0.24, chimp mean 0.05, P < 2.2 × 10 − 16 ). Overall, this analysis reveals that ~ 450,000 more fixed substitutions in autopod regulatory elements occurred on the human branch than the chimpanzee branch, with ~ 55,000 accessible only in the autopod. We partitioned these fixed substitutions along the same anatomical axes discussed above, finding that except for the elements unique to phalangeal and or late tissues, all sets showed more fixed substitutions along the human branch (Fig. 4d). Discussion The evolution of bipedalism is a critical event in hominin evolution and involved substantial remodeling of the postcranial skeleton. To further our understanding of the genetic mechanisms underlying these skeletal changes, we used functional genomics methods to investigate gene expression and regulation during development of the distal autopod skeleton in human. We identified thousands of regulatory elements and genes involved in the developing human hand and foot skeleton. These data revealed substantial differences in gene expression and regulation between the phalanges and metapodials and between the two timepoints but not between digits or limb types. While some of the identified regulatory elements are highly spatially and temporally restricted, most are accessible in multiple skeletal tissues, both within the autopod and across the postcranial skeleton. These results provide a genetic basis for the observed patterns of covariation across the autopod skeleton. We found that both multi-tissue and tissue-specific regulatory elements overlap thousands of genomic features that differ between humans and chimpanzees and could potentially underlie the evolution of human-specific autopod traits; these genomic features include human-specific nucleotide changes, losses or gains of sequence, as well as structural variation. Elements broadly accessible in autopod tissues provide a genetic basis for models of human limb evolution based on strong morphological integration (i.e. elevated phenotypic covariance) between homologous elements of the fore- and hind limbs 20 , 21 . While autopod covariation is relatively weaker in comparison to quadrupedal monkeys, homologous elements of the hands and feet of humans and other apes still covary substantially, suggesting the evolution in one limb could be driven at least in part by selection on the other 20 , 21 , 42 . Indeed, simulations suggest that selection for human-like foot proportions in a chimpanzee-like ancestor is sufficient to drive the evolution of the human hand, while the converse is not true 23 . The data presented in this study provide some support for this model. While we found many human-specific genomic features in regions accessible in multiple tissues, consistent with phenotypic covariation between elements, support for stronger selection on the foot depends on the evolutionary signal considered. In regulatory elements accessible only in autopod tissues, we found significantly more HARs, hCONDELs, and inversions but fewer SDRs within regulatory elements patterning the foot than the hand (Fig. 4b,c). We also found lower levels of intraspecific variation in the foot than the hand in human, chimpanzee, and gorilla, consistent with stronger stabilizing selection on the foot in each lineage. Finally, while there are overall more fixed substitutions found in foot elements than hand elements, both absolutely and per base pair, there is a greater excess of human compared to chimpanzee fixed substitutions in the hand (Fig. 4d). These patterns generally hold when considering the sets of all elements, autopod-specific elements, or human-mouse conserved autopod elements, but there are some differences, such as significantly more HAQERs overlapping hand than foot regulatory elements in the total set of brain-filtered elements. An important caveat of our approach is that as we rely on accessibility data from extant human samples, these patterns may not reflect the ancestral state. Altogether, we estimated that there are almost 2.5 million fixed nucleotide substitutions in the human genome compared to chimpanzee and gorilla that are accessible to transcriptional machinery in the developing autopod skeleton. This is ~ 450,000 more fixed nucleotides than are fixed in chimpanzee. Together with the fact that the human fixed nucleotides disrupt more evolutionarily conserved positions, this suggests a potentially greater extent of anatomical evolution along the human branch. ~55,000 of these fixed substitutions are uniquely accessible in autopod tissues, suggesting that autopod evolution likely has a complex underlying genetic architecture, even without considering integration with other elements of the skeleton. To the best of our knowledge, this is the first attempt to assign a minimum number of genome-wide fixed substitutions to a given set of human-specific phenotypic changes and highlights the extent of genomic changes that have potentially contributed to the evolution of the human autopod skeleton. Linking these genomic changes to observed evolutionary changes in anatomy remains a major challenge 43 , 44 . While regulatory element sets that are uniquely accessible in tissues with a particular trait might offer the most tempting candidate(s) for identifying the genomic changes underlying evolution in that trait, this approach to connecting genotype to phenotype excludes large numbers of regulatory elements that are biologically plausible by implicitly assuming that traits evolve independently 22 . In fact, our data demonstrate that many disparate skeletal elements share an underlying genetic architecture. These shared regulatory elements could be targeted by selection on a single skeletal trait so long as the off-target effects were minimal or could reflect more complex patterns of coevolution across skeletal elements. When we partition our data between the hands and the feet, between timepoints, or between phalanges and metapodials, in each case we find many shared regulatory elements that overlap human-specific changes just within the autopod-specific set. Many more regulatory elements are potentially evolutionarily relevant when elements shared with other postcranial bones are considered and even more would be identified if we removed the brain filtering step of our data processing pipeline. While substantially more challenging to interpret, elements accessible in multiple tissues are an important part of the genomic architecture patterning these skeletal elements. Of course, many of the human-specific sequence features highlighted in this study likely have no impact on gene expression. While estimating the number and magnitude of effect on gene expression for all these sequences features is beyond the capacities of current technologies, limited experimental research has shown that for regulatory HARs for which both the human and chimp sequence have been tested in transgenic animals, a third resulted in differential enhancer activity 45 . Massively parallel reporter assay testing of hCONDELs found that 8% exhibited differential enhancer activity between human and chimps in at least one cell line 40 . If HARs and hCONDELs show similar levels of activity in the developing skeleton, that would suggest that 158 HARs and 71 hCONDELs have altered regulatory activity in the developing human autopod skeleton (considering elements with accessibility in both autopod and non-autopod tissues, 25 HARs and 13 hCONDELs with accessible only in the autopod). Even if HARs and hCONDELs are substantially more predictive of human-specific gene regulatory evolution than the other sequence features considered, an assumption that has not been rigorously tested, this still suggests that the genomic basis of autopod evolution involved hundreds of genomic loci and should be considered highly polygenic. Under such a model, each individual genomic change is expected to have only a small phenotypic contribution, potentially complicating attempts to validate their effects in humanized mouse models 46 . To fully understand the genetic underpinning of human autopod evolution, future studies are necessary to generate comparable data on the bones of the wrist, midfoot, and hindfoot. These proximal autopod bones have also undergone substantial evolution in the human lineage 5 , 14 , 47 , 48 and are undoubtedly directly relevant to the evolution of the metapodials and structure of the autopod in general. Similarly, future studies should investigate the genomic changes underlying evolution in the soft tissues of the autopod. As this study was based on bulk chondrocyte tissue dissections, future research would also benefit from techniques with greater spatial resolution, such as spatial multi-omics, to unravel the genetics basis of these human-specific features with greater anatomical precision (e.g., Senevirathne et al., under review). Such methods could investigate human-specific changes in shape related to load bearing and the articulation of bones with one another. Furthermore, these methods could explore the impacts of extrinsic signaling and mechanical impacts from neighboring tissues, such as mesenchyme and muscle, on autopod skeletal development 2 . Given the intricate functional, developmental, and evolutionary relationships between all parts of the autopod, a holistic approach that integrates datasets from multiple tissue types will ultimately be necessary to understand the evolution of human-specific features. While only part of the story, this study provides novel insights into the evolutionary history of human bipedalism by highlighting the complex genetic architecture underlying human autopod evolution in the phalanges and metapodials. Materials and Methods Estimated human-chimpanzee anatomical change Measurements were taken of the length and head width of metapodials, proximal, and middle phalanges of the hand and foot in adult human ( n = 48), chimpanzee ( n = 44), and gorilla ( n = 57) skeletons (Appendix, Table 5). Chimpanzee and gorilla measurements come from wild-shot populations housed in the Powell-Cotton Museum and Hamman Todd Osteological Collection at the Cleveland Museum of Natural History, while human measurements are from the Hamman Todd Collection only. See ref 20 for details. The proximal and middle phalanges measurements were combined for each digit. The human-chimpanzee percent change for the length and width of each element was calculated as the mean human measurement/the chimpanzee measurement * 100. These calculations were repeated using the gorilla measurements in place of chimpanzee and to compare chimpanzee to gorilla. Human developmental sample collection. Human developmental samples were collected from first-trimester termination through the Birth Defects Research Laboratory (BDRL) at the University of Washington in full compliance with the ethical guidelines of the National Institutes of Health (NIH) and with the approval of the University of Washington Institutional Review Boards (IRB) for the collection and distribution of human tissues for research and Harvard University for the receipt and use of such materials. The BDRL obtained written consent from all tissue donors. Harvard University IRB determined that these samples constitute Non-Human Subjects Determination Status (Capellini: IRB16-1504). The fresh human samples were briefly washed in Hanks’ balanced salt solution and shipped at 4˚C. Upon arrival, the samples were immediately dissected under a light dissection microscope and directly subjected to RNA-seq or ATAC-seq protocols described below, following approved Harvard University IRB (IRB16-1504) and Committee on Microbiological Safety (COMS) (18–103) protocols. The BDRL performs polymerase chain reaction (PCR) using SRY and Amelogenin primers to determine the biological sex of each sample. Human ATAC-seq data collection and processing . Developmental samples were microdissected under a light microscope in 5% fetal bovine serum (FBS) in Dulbecco’s Modified Eagle Medium (DMEM) on ice. The pooled phalanges (proximal and distal for digit I; proximal, intermediate, and distal for digits II-V) and metapodials were collected separately for all digits in both the fore- and hind limbs. Samples were collected at two timepoints: an early stage (E53 to E59, n ≥3) and a later stage (E67 to E74, n ≥ 3) to capture the window of time when the autopod skeleton is prepatterned in cartilage and is just beginning to ossify and to account for heterogeneity in the timing of forelimb or hind limb as compared to overall limb development 49 , 50 . The tissue in each sample was then digested using 0.5% collagenase II in 5% FBS/DMEM in a 37˚C water bath for one hour. Every 30 minutes, the samples were spun down and gently pipetted to break up large clumps of cells. Next, the samples were incubated at 37˚C for an additional hour in a shaking incubator. Following these incubation steps, the sample were immediately placed on ice and then filtered through a 70 µm nylon cell strainer into a 50 mL conical tube by gently pressing any residual tissue through the filter followed by rinsing with 5% FBS/DMEM. Sample were centrifuged for 5 minutes at 500 g at 4˚C and most of the media aspirated. The cells were then resuspended and transferred to a 1.5 mL tube before another 5-minute centrifuge at 500 g at 4˚C. At this stage, live and dead cells were counted to ensure that all tissues had ~ 50,000 cells and cell death rates below 10%. 50,000 cells per technical replicated were resuspended in 1x PBS and then lysed using the ATAC-seq lysis buffer and centrifuged for 10 minutes at 4˚C 51 . The supernatant was discarded and 50 µl transposition mix containing 2.5 µL transposase was added before the samples were incubated at 37˚C for 30 minutes. Lastly, the Zymo DNA Clean and Concentrator kit was used, and the resulting purified DNA was eluted in 13 µl of water heated to 60–70˚C. Samples were stored at -20˚C prior to PCR amplification and barcoding. Samples were amplified and given an 8 bp barcode via 11 cycles of PCR amplification using the NEBNext High-Fidelity 2x PCR master mix. Following amplification, fragments were size-selected using Mag-Bind® RxnPure Plus beads. The sample pool was sequenced on a NovaSeq S4 machine at the Harvard Bauer Core to generate ≥40 million reads per sample. FastQC version 11.9 was used to determine the quality of fastq read files. Because some samples had been run on multiple lanes to achieve the desired number of reads, the reads for each sample were concatenated into a single file for R1 and a second file for R2. NGmerge version 0.3 was used to trim adapters from all reads 52 . Reads were then aligned to the Illumina prebuilt hg38 human reference genome using Bowtie2 (version 2.3.4.1) 53 , 54 . Reads were then indexed using samtools index 55 and duplicates removed using picard MarkDuplicates (version 2.9.0). The resulting files were indexed as before and mitochondrial reads were filtered out ( https://github.com/harvardinformatics/ATAC-seq/blob/master/atacseq/removeChrom.py ). The resulting .bam files were then used for peak calling via MACS software (version 2.1.1.2), using BAMPE and the following flags: --nolambda –bdg –verbose 56 . Reproducible peaks across replicates were identified at an IDR threshold of < 0.05, as defined by the IDR statistical test (version 2.0.3). Finally, IDR-called peak sets for autopod skeletal elements were filtered by peaks previously identified in E54 human brain to remove regulatory regions likely associated with general cellular housekeeping processes 30 . Peak sets were merged, subtracted, overlapped, etc. using the appropriate bedtools (v2.27.1) functions 57 . Human RNA-seq data collection and processing. Human samples as described above were dissected in 10% FBS/DMEM and a minimum of n = 6 per tissue per timepoint was collected. Each element was stripped of soft tissue and collected in a 2-ml tube containing 200 µl of TRIzol and one 5-mm stainless steel bead. Each sample was then homogenized at 50 Hz for 2 min, followed by 1 min on ice, and a second homogenization at 50 Hz for 2 min. Samples were stored at -80˚C until RNA extraction was performed. For each RNA extraction, samples were incubated at room temperature for 5 min and then centrifuged at 4˚C for 5 minutes at 12,000 x g. The supernatant was collected and transferred to a clean tube, to which 200 µl of chloroform was added per 1 ml of TRIzol. Each sample was then vortexed briefly and incubated at room temperature for 2 minutes before being transferred to a MaXtract tube and centrifuged at 4˚C for 5 minutes at 12,000 x g. Following centrifugation, the aqueous phase was removed and transferred to a new microcentrifuge Eppendorf tube and an equal volume of 100% ethanol was added. After completion of this phenol-chloroform step, the RNA extraction continued using the Direct-zol RNA MicroPrep kit following the manufacturer’s protocols and eluted in 15 µl of nuclease-free water. Samples were stored at -80˚C until library preparation. To maximize the number of tissues useable from each biological replicate, samples with RNA integrity number (RIN) scores higher than 6 were used for subsequent steps so long as the average RIN for all tissues from that sample was greater than 7. For cDNA library preparation and sequencing, samples were normalized to a single concentration and libraries were prepared using the Kapa mRNA HyperPrep kit with an input volume of 25 ng per sample following the manufacturer’s protocols. After an initial MiSeq Nano run to assess library quality, samples were repooled as necessary and then sequenced on an Illumina NovaSeq S4 six times to generate ~ 20 million paired-end reads per sample. See SI Appendix, Table S2 for detailed sequencing information for each sample, including index primer, read count, and other information. Six biological replicates per tissue per timepoint were sequenced. FastQC version 11.9 was used to assess the quality of each fastq read file 58 . Reads for each sample from multiples lanes were concatenated into a single file for R1 and a second file for R2. NGmerge version 0.3 was used to trim adapters from all reads 52 . Next, reads mapping to ribosomal RNA were removed using RiboDetector 59 version 0.2.7. STAR 60 version 2.7.1 was used to map reads to the human genome (hg38) with ≥80% of reads uniquely mapping for each sample. RSEM 61 version 1.3.3 was then used to generate read counts. For all samples, > 90% of reads were mapped uniquely to a gene. Mouse ATAC-seq data collection and processing. All mouse work was covered under the Capellini lab IACUC protocol (#13-04-161-3). E15.5 mouse embryos were collected from pregnant FVB/NJ females and transferred to cold 1x PBS. The fore- and hind limbs of three to five embryos were microdissected under a light microscope in 5% FBS/DMEM on ice. The pooled phalanges (proximal and distal for digit I, proximal, intermediate, and distal for digits III and V) and metapodials were collected separately for digits I, III, and V for both the forelimb and the hind limb. The right and left sides of each element were pooled together, resulting in a total pool of elements from six to ten individual limbs. All embryos in each biological replicate were from the same litter. Collagenase digestion, cell lysis, the transposase reaction, barcoding, and sequencing were all performed as described above for human samples. Data processing was also the same except that reads were aligned to the prebuild Illumina UCSC mm10 genome, E15.5 mouse brain accessibility data was used for the brain-filtering step 62 , and a different version of MACS (3.0.3) was used. See Supplementary Table 5 for detailed sequencing information for each sample, including index primer, read count, and other information. Mouse RNA-seq data collection and processing. E15.5 mouse embryos were collected from pregnant FVB/NJ females and transferred to cold 1x phosphate-buffer saline (PBS). Embryos were microdissected under a light microscope in 10% FBS/DMEM on ice and anatomically sexed 63 . The pooled phalanges (proximal and distal for digit I; proximal, intermediate, and distal for digits III and V) and metapodials were collected separately for digits I, III, and V for both the fore- and hind limbs. The right and left sides of each element were pooled together. All other elements of the RNA collection and extraction protocol were the same as described for human above. For cDNA library preparation and sequencing, samples were normalized to a single concentration and libraries were prepared using the Takara SMART-Seq v4 Ultra Low Input RNA Kit following the manufacturer’s protocols. Two samples with high concentrations were also prepared using the Kapa mRNA HyperPrep kit with an input volume of 25 ng per sample as per the human samples to detect any effects of library preparation method. After initial pooling, the relative concentrations of each sample were evaluated by running the pool on a MiSeq Nano and the pool concentrated adjusted accordingly. The final pool was then sequenced four times on an Illumina NovaSeq S4 to generate ≥40 million paired-end reads per sample. On average, samples were sequenced at 86 million reads per sample, within the recommended Encyclopedia of DNA Elements (ENCODE) guidelines ( www.encodeproject.org/about/experiment-guidelines/ ). See Supplementary Table 4 for detailed sequencing information for each sample, including index primer, read count, and other information. Six biological replicates of each tissue were sequenced. Differential expression and cluster analysis of RNA-seq. The large number of tissues investigated in this study (20 tissues at two timepoints in human, 12 in mouse) allows for many potentially pair-wise comparisons, however, only a subset of these comparisons is expected to be biologically meaningful. To minimize the number of factors that could cause gene expression differences, samples from a given tissue were only compared with directly adjacent tissues and with the homologous tissue in the other limb type. For humans, adjacent tissues came from sequential digits while in mouse, only digits I, III, V were collected so the nearest digit was used in place of the adjacent digit. Each human tissue was also compared across timepoints. This comparison scheme aims to target only a single variable for any given differential gene calculation, either spatial location within the autopod, limb type, or timepoint. To give an example, the gene expression profile of early samples of hand phalanges from digit III were compared separately with (1) the early hand phalanges of digit II, (2) the early hand phalanges of digit IV, (3) the early metacarpal of digit III, (4) the early foot phalanges of digit III, and (5) the late hand phalanges of digit III. All differentially expressed genes (DEGs) were identified using DESeq2 64 version 1.36.0 with a design including the sample type and the replicate ID except for comparisons between timepoints where only timepoint was included in the design. In addition, to identify genes expressed in a gradient across the digits, DESeq2 was rerun at both time points using a likelihood ratio test to compare a model of ~ replicate + digit with a baseline model of ~ replicate. Thresholds for calling DEGs were P adj 1.5. This approach was used for both human and mouse. Many of the DEGs identified above are expected to be DEGs in numerous comparisons, e.g., a gene highly expressed in a single tissue will be identified as a DEG in all comparisons involving that tissue. Similarly, many genes may be expressed in specific spatial domains encompassing multiple tissues and cannot be easily identified using a pair-wise approach. Therefore, to identify large patterns of gene expression present in our dataset, we performed a weighted gene co-expression network analysis (WGCNA) to identify clusters of genes with similar expression patterns 65 . The expression of DEGs was normalized using the variance transformation function vst in DESeq2. For comparison within a single timepoint, variation associated with biological replicates was removed using limma:removeBatchEffect) 66 . This was not used when comparing both human timepoints because this correction also removed timepoint differences that are of biological interest. The pickSoftThreshold function was used to estimate the lowest scale-free threshold for each dataset, which resulted in a minimum scale-free topology of 0.9. These soft power thresholds were used to construct signed correlation networks using the blockwiseModules function. Mouse data was processed following the same pipeline except that nearest digits were compared instead of neighboring digits (so digit III was compared to digits I and V instead of II and IV) due to reduced tissue sampling. DEG-calling was repeated for human with the reduced tissue set (digits I, III, V only) at each timepoint to facilitate direct comparison with mouse. Comparison of human and mouse ATAC/RNA-seq. Mouse orthologs of human genes (GRCh38.p14) were downloaded from the Ensembl BioMart database for comparison of gene sets between the two species. For comparison of orthologous regulatory elements, mouse regions (mm10) were lifted over to human genome coordinates (hg38), allowing for multiple matches and a minMatch = 0.1 threshold. Identification of accessibility sharing with other elements of the developing skeleton Regulatory elements accessible in the autopod were intersected with elements accessibility in at least one other developmental skeletal tissue for which data has been generated using the same ATAC-seq protocol 28 – 30 . This included the proximal and distal ends of the femur, tibia, humerus, and radius, three regions of the scapula (head and neck, blade, acromion), four regions of the pelvis (ilium, ischium, pubis, acetabulum), as well as the thoracic and the lumbar vertebrae. Data were available at both stages for all tissues except the vertebrae, for which only later stage was available. Identification of regulatory overlaps with genomic features To identify potential regulatory elements that were the targets of past positive selection in the human lineage, regulatory element sets were overlapped with human accelerated regions (HARs), which are non-coding regions of the genome marked by accelerated substitution rates in the human lineage compared to other vertebrates 31 – 38 , human conserved element deletions (hCONDELs) 39 , 40 , and human ancestor quickly evolved regions (HAQERs), the fasted evolving regions of the human genome 41 . Regulatory element sets were also overlapped with structural divergent regions of the human genome compared to other apes, including large chromosomal inversions and structural divergent regions (SDRs) 41 . As the HAQERs, inversions, and SDR sets were identified using the T2T-CHM13v2.0 genome, these sets were lifted over to hg38 using the liftover file available at https://github.com/marbl/CHM13 . Regions flanking inversions and SDRs were obtained for either 500kb or 1Mb windows on either side of these features. To investigate patterns of intraspecific variation in regulatory elements, the VCF file containing the positions of human SNPs in the dbSNP151 database were downloaded from the NIH website ( https://ftp.ncbi.nlm.nih.gov/snp/organisms/ ) and converted to bed format using the BEDOPS (v2.4.40) convert2bed function 67 . Chimp and gorilla SNP VCF files were downloaded from the Great Apes Genome Diversity Project 68 . We curated common (i.e., minor allele frequency > 0.05) single nucleotide polymorphisms (SNP) from each species. All intersections were performed using bedtools intersect 57 . To compare the number of regulatory elements from distinct sets falling near expressed genes, the coordinates of genes with a minimum normalized expression of 10 in at least one autopod tissue based on the UCSC hg38 refGene list were extended by 500kb in either direction and overlapped with the regulatory element sets using the bedtools intersect -c function. Statistical significance was assessed using a two-sided Pearson’s correlation test. Gene ontology enrichment analysis Gene ontology (GO) analysis of gene lists was performed using clusterProfiler with a background set of all genes expressed in at least one of the sampled tissues 69 . Enrichment analysis of genomic region sets was performed using rGREAT (v2.6.0) for five ontologies ("GO Molecular Function", "GO Biological Process", "GO Cellular Component", "Human Phenotype", "Mouse Phenotype") using the standard thresholds used in the online version of tool (binomial and hypergeometric P adj < 0.05, binomial fold enrichment ≥ 2) 70,71 . Transcription factor (TF) binding motif enrichment analysis TF binding motif enrichments were calculated in R using MonaLisa 72 with human PWM matrices from the JASPAR database 73 . Statistical significance was assessed using a negative log10 cut-off of 4. Regulatory element sizes were not standardized for these analyses. Calculation of human chimp sequence divergence Prior to performing this analysis, all elements in the genomic region sets were standardized to a length of 500 bps. To calculate the minimum number of single base pair differences between the human and chimpanzee genomes in each set of genomic regions, we first filtered out any positions that are variable within human, chimpanzees, or gorillas since interspecific differences in autopod morphology are much larger than intraspecific variation and therefore current allelic variants are insufficient to explain the evolutionary changes. VCF lists of variants for each species were obtained as described above and SNPs overlapping autopod regions were extracted using tabix (version 1.13) 74 and converted to bed format using bedtools and then lifted over to hg38 (for chimp and gorilla). Genomic regions sets were filtered by ENCODE blacklist regions 75 as well as the SNP sets from each species using bedtools to prioritize high quality genomic sequence without interspecific variants. The resulting regions were then lifted over to the chimpanzee (panTro6) and gorilla (gorGor6) genomes, and the DNA sequence extracted using the appropriate BSgenome package 76 . Regions had to be mappable between species along the entire length of the sequence and have at least 25% sequence identity between either target species (human or chimp) and its outgroups. This threshold was chosen to prioritize comparison of orthologous base pairs rather than faulty alignments or indel-related mismatches because 25% is the expected match for random sequences of the same length. A nucleotide was considered derived in a target lineage if it differed from the nucleotide observed (and likely conserved) in the other two species. This analysis was performed with both human and chimpanzee as the target species. The Zoonomia PhyloP scores for each genomic position were obtained for the human and chimpanzee sets 77 and compared using two-sided Wilcoxon rank-sum test. Segmentation of the adult autopod scans Raw TIFFs and mesh files for corresponding human and chimp autopods were downloaded from MorphoSource and imported into Amira v.2022.1. Slices showing ossification were segmented, and the selected voxels were assigned to a separate material. To visualize the autopod in 3D and generate a model, the “Generate Surface” and “Surface View” options were used in Amira. The final rendered models were exported as object files (.obj) and visualized using Blender 4.4. Declarations Data, Materials, and Software Availability All data needed to evaluate the conclusions in the paper are present in the paper and/or the SI Appendix. In addition, all ATAC-seq data raw sequencing fastq files and processed peak bed files and RNA-seq data raw sequencing fastq files and processed read count files have been deposited on National Center for Biotechnology Information Gene Expression Omnibus (GSE281502, GSE283854, GSE284187, GSE286924). The code used to perform these analyses and to generate these figures is available at: https://github.com/aokamoto-bio/Human_Autopod_Evolution. Acknowledgments We would like to thank the University of Washington Birth Defects Research Laboratory including Kimberly A. Aldinger, Dan Doherty, Ian G. Phelps, Jennifer C. Dempsey, Yasmeen Otaibi, Lucinda A. Cort; funded via NIH under NICHD Grant # R24HD000836, to I.A.G.; C. Reardon and the Harvard Bauer Core Facility for ATAC-seq and RNA-seq assistance and sequencing; C. Tabin, D. Lieberman, K. Cooper, D. Richard, and members of the Capellini laboratory (Harvard University) for critical insight and commentary on this project. T.D.C., A.S.O., G.S., and P.M., were supported by the Harvard University Dean’s Competitive Fund and Milton Fund for the functional genomics portion of this research. T.D.C. and A.S.O additionally received funds for mouse work from the National Science Foundation (BCS1518596 and BCS1847979). A.S.O and T.D.C. were supported by a National Science Foundation Doctoral Dissertation Improvement Grant (BSC-2337516). A.S.O. was supported by a National Science Foundation Graduate Research Fellowship (DGE-1745303). T.D.C. was supported by the American School of Paleontological Research (ASPR) via Harvard University. Author Contributions: A.S.O. performed human autopod sample dissection, ATAC-seq, RNA-seq, downstream data processing, all computational analyses, aided in study design, wrote the manuscript, and generated figures with input from all authors. G.S. and P.M. assisted with sample collection. C.R. provided comparative morphological data and helped write and revise the manuscript. I.G. and B.D.R.L. provided samples. T.D.C. conceived and supervised all aspects of the project, designed the study, and wrote and revised the manuscript with input from all authors. References Darwin, C. The Descent of Man, and Selection in Relation to Sex . (John Murray, London, 1871). Diogo, R., Siomava, N. & Gitton, Y. Development of human limb muscles based on whole-mount immunostaining and the links between ontogeny and evolution. Development 146 , (2019). Diogo, R., Richmond, B. G. & Wood, B. Evolution and homologies of primate and modern human hand and forearm muscles , with notes on thumb movements and tool use. J Hum Evol 63 , 64–78 (2012). Diogo, R., Molnar, J. L., Rolian, C. & Esteve-Altava, B. First anatomical network analysis of fore- and hindlimb musculoskeletal modularity in bonobos , common chimpanzees , and humans. Sci Rep 8 , 1–9 (2018). Kivell, T. L., Baraki, N., Lockwood, V., Williams‐Hatala, E. M. & Wood, B. A. Form, function and evolution of the human hand. American Journal of Biological Anthropology 181 , 6–57 (2023). Kivell, T. L. Evidence in hand: recent discoveries and the early evolution of human manual manipulation. Philosophical Transactions of the Royal Society B 370 , 20150105 (2015). Bardo, A., Vigouroux, L., Kivell, T. L. & Pouydebat, E. The impact of hand proportions on tool grip abilities in humans, great apes and fossil hominins: A biomechanical analysis using musculoskeletal simulation. J Hum Evol 125 , 106–121 (2018). Liu, M.-J., Xiong, C.-H. & Hu, D. Assessing the manipulative potentials of monkeys, apes and humans from hand proportions: implications for hand evolution. Proceedings of the Royal Society B: Biological Sciences 283 , 20161923 (2016). Young, R. W. Evolution of the human hand: the role of throwing and clubbing. J Anat 202 , 165–174 (2003). Marzke, M. W. & Marzke, R. F. Evolution of the human hand: Approaches to acquiring, analysing and interpreting the anatomical evidence. J Anat 197 , 121–140 (2000). Holowka, N. B. & Lieberman, D. E. Rethinking the evolution of the human foot: Insights from experimental research. Journal of Experimental Biology 221 , (2018). Lieberman, D. E. & Bramble, D. M. Endurance running and the evolution of Homo. Nature 432 , 345–352 (2004). Rolian, C., Lieberman, D. E., Hamill, J., Scott, J. W. & Werbel, W. Walking, running and the evolution of short toes in humans. Journal of Experimental Biology 212 , 713–721 (2009). DeSilva, J., McNutt, E., Benoit, J. & Zipfel, B. One small step: A review of Plio-Pleistocene hominin foot evolution. Am J Phys Anthropol 168 , 63–140 (2019). Hayafune, N., Hayafune, Y. & Jacob, H. A. C. Pressure and force distribution characteristics under the normal foot during the push-off phase in gait. The Foot 9 , 88–92 (1999). Almecijia, S., Smaers, J. B. & Jungers, W. L. The evolution of human and ape hand proportions. Nat Commun 6 , 1–11 (2015). Tocheri, M. W., Orr, C. M., Jacofsky, M. C. & Marzke, M. W. The evolutionary history of the hominin hand since the last common ancestor of Pan and Homo. J Anat 212 , 544–562 (2008). Williams-Hatala, E. M. et al. The manual pressures of stone tool behaviors and their implications for the evolution of the human hand. J Hum Evol 119 , 14–26 (2018). Shou, S., Scott, V., Reed, C., Hitzemann, R. & Stadler, H. S. Transcriptome analysis of the murine forelimb and hindlimb autopod. Dev Dyn 234 , 74–89 (2005). Rolian, C. Integration and evolvability in primate hands and feet. Evol Biol 36 , 100–117 (2009). Rolian, C., Lieberman, D. E. & Hallgrímsson, B. The coevolution of human hands and feet. Evolution (N Y) 64 , 1558–1568 (2010). Auerbach, B. M., Savell, K. R. R. & Agosto, E. R. Morphology, evolution, and the whole organism imperative: Why evolutionary questions need multi-trait evolutionary quantitative genetics. Yearbook Biological Anthropology 181 Suppl 76 , 180–211 (2023). Rolian, C. Re‐evaluating the correlated evolution of human hands and feet using viability selection modeling. American Journal of Biological Anthropology 183 , (2024). Cotney, J. et al. The Evolution of Lineage-Specific Regulatory Activities in the Human Embryonic Limb. Cell 154 , 185–196 (2013). Zhang, B. et al. A human embryonic limb cell atlas resolved in space and time. Nature (2023) doi:10.1038/s41586-023-06806-x. Taher, L. et al. Global gene expression analysis of murine limb development. PLoS One 6 , (2011). Cotney, J. et al. Chromatin state signatures associated with tissue-specific gene expression and enhancer activity in the embryonic limb. Genome Res 22 , 1069–1080 (2012). Richard, D. et al. Functional genomics of human skeletal development and the patterning of height heritability. Cell 188 , 15-32.e24 (2025). Young, M. et al. The developmental impacts of natural selection on human pelvic morphology. Sci Adv 8 , 1–24 (2022). Richard, D. et al. Evolutionary Selection and Constraint on Human Knee Chondrocyte Regulation Impacts Osteoarthritis Risk. Cell 181 , 362-381.e28 (2020). Gittelman, R. M. et al. Comprehensive identification and analysis of human accelerated regulatory DNA. Genome Res 25 , 1245–1255 (2015). Pollard, K. S. et al. Forces shaping the fastest evolving regions in the human genome. PLoS Genet 2 , 1599–1611 (2006). Bird, C. P. et al. Fast-evolving noncoding sequences in the human genome. Genome Biol 8 , 1–12 (2007). Bush, E. C. & Lahn, B. T. A genome-wide screen for noncoding elements important in primate evolution. BMC Evol Biol 8 , 1–10 (2008). Prabhakar, S., Noonan, J. P., Pääbo, S. & Rubin, E. M. Accelerated evolution of conserved noncoding sequences in humans. Science (1979) 314 , 786 (2006). Bi, X. et al. Lineage-specific accelerated sequences underlying primate evolution. Sci Adv 9 , (2023). Keough, K. C. et al. Three-dimensional genome rewiring in loci with human accelerated regions. Science (1979) 380 , (2023). Kostka, D., Holloway, A. K. & Pollard, K. S. Developmental Loci Harbor Clusters of Accelerated Regions That Evolved Independently in Ape Lineages. Mol Biol Evol 35 , 2034–2045 (2018). McLean, C. Y. et al. Human-specific loss of regulatory DNA and the evolution of human-specific traits. Nature 471 , 216–219 (2011). Xue, J. R. et al. The functional and evolutionary impacts of human-specific deletions in conserved elements. Science (1979) 380 , (2023). Yoo, D. et al. Complete sequencing of ape genomes. Nature 641 , 401–418 (2025). Young, N. M., Wagner, G. P. & Hallgrímsson, B. Development and the evolvability of human limbs. Proc Natl Acad Sci U S A 107 , 3400–3405 (2010). Cooper, K. L. The case against simplistic genetic explanations of evolution. Development 151 , dev203077 (2024). Varki, A. & Altheide, T. K. Comparing the human and chimpanzee genomes: Searching for needles in a haystack. Genome Res 15 , 1746–1758 (2005). Whalen, S. & Pollard, K. S. Enhancer Function and Evolutionary Roles of Human Accelerated Regions. Annu Rev Genet 56 , 423–439 (2022). Baumgartner, M., Ji, Y. & Noonan, J. P. Reconstructing human-specific regulatory functions in model systems. Curr Opin Genet Dev 89 , 102259 (2024). Aiello, L. & Dean, C. An Introduction to Human Evolutionary Anatomy . (Academic Press, London ; San Diego, 1990). Lewis, O. J. Functional Morphology of the Evolving Hand and Foot . (Clarendon Press, Oxford, 1989). Gardner, E., OʼRahilly, R. & Gray, D. J. The Prenatal Development of the Skeleton and Joints of the Human Foot. J Bone Joint Surg 41 , 847–876 (1959). Gray, D. J., Gardner, E. & O’Rahilly, R. The prenatal development of the skeleton and joints of the human hand. American Journal of Anatomy 101 , 169–223 (1957). Buenrostro, J. D., Wu, B., Chang, H. Y. & Greenleaf, W. J. ATAC‐seq: A Method for Assaying Chromatin Accessibility Genome‐Wide. Curr Protoc Mol Biol 109 , 139–148 (2015). Gaspar, J. M. NGmerge: merging paired-end reads via novel empirically-derived models of sequencing errors. BMC Bioinformatics 19 , 536 (2018). Langmead, B., Wilks, C., Antonescu, V. & Charles, R. Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics 35 , 421–432 (2019). Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat Methods 9 , 357–359 (2012). Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10 , (2021). Zhang, Y. et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol 9 , R137 (2008). Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26 , 841–842 (2010). Andrews, S., Krueger, F., Seconds-Pichon, A., Biggins, F. & Wingett, S. FastQC: A quality control tool for high throughput sequence data. Preprint at (2015). Deng, Z.-L., Unch, P. C. M. ¨, Mreches, R. & Mchardy, A. C. Rapid and accurate identification of ribosomal RNA sequences via deep learning. Nucleic Acids Res 50 , 60 (2022). Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29 , 15–21 (2013). Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12 , 323 (2011). Guo, M. et al. Epigenetic profiling of growth plate chondrocytes sheds insight into regulatory genetic variation influencing height. Elife 6 , 70–90 (2017). Staack, A., Donjacour, A. A., Brody, J., Cunha, G. R. & Carroll, P. Mouse urogenital development: a practical approach. Differentiation 71 , 402–13 (2003). Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15 , 550 (2014). Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9 , 559 (2008). Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43 , e47–e47 (2015). Neph, S. et al. BEDOPS: high-performance genomic feature operations. Bioinformatics 28 , 1919–1920 (2012). Prado-Martinez, J. et al. Great ape genetic diversity and population history. Nature 499 , (2013). Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation 2 , 100141 (2021). McLean, C. Y. et al. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol 28 , 495–501 (2010). Gu, Z. & Hübschmann, D. rGREAT : an R/bioconductor package for functional enrichment on genomic regions. Bioinformatics 39 , (2023). Machlab, D. et al. monaLisa: an R/Bioconductor package for identifying regulatory motifs. Bioinformatics 38 , 2624–2625 (2022). Fornes, O. et al. JASPAR 2020: update of the open-access database of transcription factor binding profiles. Nucleic Acids Res (2019) doi:10.1093/nar/gkz1001. Li, H. Tabix: fast retrieval of sequence features from generic TAB-delimited files. Bioinformatics 27 , 718–9 (2011). Amemiya, H. M., Kundaje, A. & Boyle, A. P. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. Sci Rep 9 , 9354 (2019). Pagès, H. BSgenome: Software infrastructure for efficient representation of full genomes and their SNPs. Preprint at https://bioconductor.org/packages/BSgenome (2024). Sullivan, P. F. et al. Leveraging base-pair mammalian constraint to understand genetic variation and human disease. Science (1979) 380 , (2023). Table 1 Genomic Feature All brain-filtered Autopod-specific Autopod-specific, human-mouse conserved Tissue and timepoint specific HARs 474 75 27 44 hCONDELs 886 157 59 93 HAQERs 215 35 9 17 Inversions 2511 488 116 365 SDRs 297 61 0 52 Table 1 Summary of regulatory element overlaps with genomic features HARs, human accelerated regions, hCONDELs, human conserved element deletions, HAQERs, human ancestor quickly evolved regions, and SDRs, structurally divergent regions between humans and other apes. Additional Declarations There is NO Competing Interest. Supplementary Files OkamotoAutopodS1HominoidMeasurements.xlsx Supplementary Table 1 OkamotoAutopodS2HumanRNA.xlsx Supplementary Table 2 OkamotoAutopodS3HumanATAC.xlsx Supplementary Table 3 OkamotoAutopodS4MouseRNA.xlsx Supplementary Table 4 OkamotoAutopodS5MouseATAC.xlsx Supplementary Table 5 AutopodEvolutionSupplement7132025.docx Supplemental Information Cite Share Download PDF Status: Published Journal Publication published 11 May, 2026 Read the published version in Proceedings of the National Academy of Sciences → Version 1 posted 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-7124496","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":486756022,"identity":"b8002df0-5c64-45fa-a85f-4510e6fcd38a","order_by":0,"name":"Alexander Okamoto","email":"","orcid":"https://orcid.org/0000-0002-4630-7548","institution":"Harvard University","correspondingAuthor":false,"prefix":"","firstName":"Alexander","middleName":"","lastName":"Okamoto","suffix":""},{"id":486756023,"identity":"d14a1dd6-2649-48d2-bd49-57812119ef3f","order_by":1,"name":"Gayani Senevirathne","email":"","orcid":"https://orcid.org/0000-0001-6704-9691","institution":"Harvard University","correspondingAuthor":false,"prefix":"","firstName":"Gayani","middleName":"","lastName":"Senevirathne","suffix":""},{"id":486756024,"identity":"a5b6df1f-8f41-4afe-a80b-22b1c6f4e664","order_by":2,"name":"Pushpanathan Muthuirulan","email":"","orcid":"","institution":"Harvard University","correspondingAuthor":false,"prefix":"","firstName":"Pushpanathan","middleName":"","lastName":"Muthuirulan","suffix":""},{"id":486756025,"identity":"bd25bd69-04ba-4665-9f70-94e86e38f1d8","order_by":3,"name":"Campbell Rolian","email":"","orcid":"","institution":"McGill University","correspondingAuthor":false,"prefix":"","firstName":"Campbell","middleName":"","lastName":"Rolian","suffix":""},{"id":486756026,"identity":"3b16022b-6bf7-40e2-8237-cd7f3a35fbcd","order_by":4,"name":"Ian Glass","email":"","orcid":"https://orcid.org/0000-0001-6762-8407","institution":"11.\tDivision of Genetic Medicine, Department of Pediatrics, University of Washington, Seattle, Washington","correspondingAuthor":false,"prefix":"","firstName":"Ian","middleName":"","lastName":"Glass","suffix":""},{"id":486756027,"identity":"d34a6272-a85f-42ef-b3fb-7742c4431971","order_by":5,"name":"Birth Defects Research Laboratory","email":"","orcid":"","institution":"University of Washington","correspondingAuthor":false,"prefix":"","firstName":"Birth","middleName":"Defects Research","lastName":"Laboratory","suffix":""},{"id":486756021,"identity":"0293936c-eab4-4ef5-b939-552963d22826","order_by":6,"name":"Terence Capellini","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABDklEQVRIiWNgGAWjYNCCAgseBmbmAx+AzASICBsQH8CumAdMGkgAtbAlziBJC4hpSJwWe/bDzx58MJCQMWfn+djw4w9DHv/sM4afK8oY5PhuJGC3hSfN3HAG0GGWzbwbG3vbGIolzuUYS545x2AsiUuLBIOZNA9Qi8Fh3u0PeBsYEhvOsCVINrYxJG7AqYX9m/QfsBaeh41//jAkzj/DlvwTqKUetxYeM2kGiBbGZh42oOFnmI+BbEkwwKXlTE65YQ9YC5ths2ybROJGoBbLhnMShjPPPMCqhb39+LYHPyps7A3OH37Y+OaPTeK8M4zNNxvKbOT5jmO3hQESBXAggcEgqGUUjIJRMApGASYAAFv/WomrN1TvAAAAAElFTkSuQmCC","orcid":"https://orcid.org/0000-0003-3842-8478","institution":"Harvard University","correspondingAuthor":true,"prefix":"","firstName":"Terence","middleName":"","lastName":"Capellini","suffix":""}],"badges":[],"createdAt":"2025-07-14 21:35:40","currentVersionCode":1,"declarations":{"humanSubjects":false,"vertebrateSubjects":false,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":false,"humanSubjectConsent":false,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":false},"doi":"10.21203/rs.3.rs-7124496/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-7124496/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1073/pnas.2603297123","type":"published","date":"2026-05-12T00:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":89279884,"identity":"25855db7-42e4-4444-bc02-79daea12bc27","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1575924,"visible":true,"origin":"","legend":"\u003cp\u003eAnatomy of the human and chimpanzee autopod skeleton\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea,\u003c/strong\u003e morphology of the adult hand (top) and foot (bottom) skeleton in human (left) and chimpanzee (right). \u003cstrong\u003eb, \u003c/strong\u003eaverage percent change in length and\u003cstrong\u003ec\u003c/strong\u003e, width between adult human and chimpanzee autopod skeletal elements. \u003cstrong\u003ed,\u003c/strong\u003e alizarin red/alcian blue bone/cartilage staining of developmental human autopod skeletons at gestational day E54 and E67. The intermediate and distal phalanges of digit III in the E50s sample were lost during sample preparation and have been drawn in. Metatarsal I in the E60-70s sample is broken and has been digitally corrected. Note that the phalanges are splayed due to the removal of soft tissue. Photos courtesy of Mariel Young.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/c7c25309863cec29ce13e657.png"},{"id":89280649,"identity":"62652ea1-8e47-435e-ab43-a18af576e8d5","added_by":"auto","created_at":"2025-08-18 10:27:33","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":560115,"visible":true,"origin":"","legend":"\u003cp\u003eGene expression and regulation in the human autopod\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea, \u003c/strong\u003ePrincipal component analysis (PCA) of normalized human gene count data, corrected for biological replicate. \u003cstrong\u003eb,\u003c/strong\u003e Volcano plot showing the differences in gene expression between early and late samples. Timepoint differences are not seen in \u003cstrong\u003ea, \u003c/strong\u003edue to correction for biological replicate in the Deseq2 model. \u003cstrong\u003ec, \u003c/strong\u003ePCA of normalized accessibility data corrected for variation between biological replicates. The top 10 most significantly differentially expressed genes early and late are labeled. \u003cstrong\u003ed,\u003c/strong\u003e Accessibility sharing of regulatory elements across timepoints. Dot sizes show the summation of all tissue combinations with the same total number of tissues with accessibility at each timepoint. Only tissue combinations with at least 5 regulatory elements were included in this calculation. See Methods for the full list of non-autopod skeletal tissues considered in this analysis.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/4e94b666849178297fc1d44a.png"},{"id":89279888,"identity":"1d581126-8bf3-49e3-bcce-47528aeb8046","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":828222,"visible":true,"origin":"","legend":"\u003cp\u003eEvolutionary signals overlapping autopod regulatory elements\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea\u003c/strong\u003e, total number of regulatory element accessible in each tissue after brain-filtering. \u003cstrong\u003eb\u003c/strong\u003e, number of regulatory elements in each tissue that are accessible only in autopod tissues, and \u003cstrong\u003ec\u003c/strong\u003e, number of regulatory elements that are accessible in only a single tissue at one or both timepoints. Number of overlaps for the brain-filtered (\u003cstrong\u003ed, g, j, m, p\u003c/strong\u003e), autopod specific (\u003cstrong\u003ee, h, k, n, q\u003c/strong\u003e), or tissue-specific (\u003cstrong\u003ef, i, l, o, r\u003c/strong\u003e) regulatory element sets with HARs (\u003cstrong\u003ed-f\u003c/strong\u003e), HAQERs (\u003cstrong\u003eg-i\u003c/strong\u003e), hCONDELs (\u003cstrong\u003ej-l\u003c/strong\u003e), major human inversions (\u003cstrong\u003em-o\u003c/strong\u003e), or human SDRs (\u003cstrong\u003ep-r\u003c/strong\u003e).\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/1f9e21ddab5657d408259f53.png"},{"id":89279892,"identity":"d6375092-ae59-4d8a-935f-9ba3abda9102","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":643165,"visible":true,"origin":"","legend":"\u003cp\u003ePartitioning of evolutionary signals across the autopods\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea\u003c/strong\u003e, division of regulatory elements accessible only in autopod tissues by the specificity of accessibility in limb type, proximal-distal region, or developmental stage. \u003cstrong\u003eb\u003c/strong\u003e, enrichment of each regulatory element set for HARs, HAQERs, hCONDELs, human-specific inversions, human SDRs, as well as common variants (MAF \u0026gt; 0.05) in human, chimpanzee, and gorilla per base pair of sequence. \u003cstrong\u003ec\u003c/strong\u003e, statistical significance of enrichment differences between set pairs. \u003cstrong\u003ed\u003c/strong\u003e, number of fixed substitutions along the human (colored bars) and chimpanzee (white bars) lineages per base pair for each set. Total number of fixed substitutions are listed above each bar. Color scheme in \u003cstrong\u003eb\u003c/strong\u003e,\u003cstrong\u003ed\u003c/strong\u003e, is the same as in \u003cstrong\u003ea\u003c/strong\u003e.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/d3adf8054174dd366416df7b.png"},{"id":109202421,"identity":"7ca1f130-2f4e-4fcf-b0f3-37daa4d90f96","added_by":"auto","created_at":"2026-05-13 14:02:27","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":4991585,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/176d5f5d-8889-47fb-b268-99e8fb591eac.pdf"},{"id":89279885,"identity":"767c2f82-26d2-4ff4-8b49-0ad27ecaa721","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"xlsx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":158830,"visible":true,"origin":"","legend":"Supplementary Table 1","description":"","filename":"OkamotoAutopodS1HominoidMeasurements.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/01b6429d8ff2b250236d972c.xlsx"},{"id":89279887,"identity":"89cc77d5-8b7d-4adb-acf8-dc1ddf4583d1","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":1335275,"visible":true,"origin":"","legend":"Supplementary Table 2","description":"","filename":"OkamotoAutopodS2HumanRNA.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/a2d6778ffd20746d68c54f9d.xlsx"},{"id":89280650,"identity":"0b197beb-8e69-41d2-b07a-a3efae920fe2","added_by":"auto","created_at":"2025-08-18 10:27:33","extension":"xlsx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":11441322,"visible":true,"origin":"","legend":"\u003cp\u003eSupplementary Table 3\u003c/p\u003e","description":"","filename":"OkamotoAutopodS3HumanATAC.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/8ca3b73741ee5a24ea24fc16.xlsx"},{"id":89279890,"identity":"ea6ab122-c10c-4afd-aac6-602f4b5974ff","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"xlsx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":8629617,"visible":true,"origin":"","legend":"Supplementary Table 4","description":"","filename":"OkamotoAutopodS4MouseRNA.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/5044450d65c4e63f479c249e.xlsx"},{"id":89281067,"identity":"e1dc8636-67fc-48bb-975a-85967e1e9fd3","added_by":"auto","created_at":"2025-08-18 10:35:33","extension":"xlsx","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":6865988,"visible":true,"origin":"","legend":"Supplementary Table 5","description":"","filename":"OkamotoAutopodS5MouseATAC.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/8c1df2781eddb21fffbaa3dc.xlsx"},{"id":89279889,"identity":"4e84ae45-8d9e-4e42-a403-20832ecda7cf","added_by":"auto","created_at":"2025-08-18 10:19:33","extension":"docx","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":3557088,"visible":true,"origin":"","legend":"Supplemental Information","description":"","filename":"AutopodEvolutionSupplement7132025.docx","url":"https://assets-eu.researchsquare.com/files/rs-7124496/v1/11a5eca0fc5dde03332e4bb9.docx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"Complex genetic architecture underlies human hand and foot evolution","fulltext":[{"header":"Main","content":"\u003cp\u003eThe functional divergence of the hand and foot for prehension and locomotion, respectively, has long been recognized as a critical event in human evolution\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u003c/sup\u003e. In comparison to other apes, the human hand and foot have features in skeletal anatomy as well as the related musculature and innervation which are believed to have evolved adaptively in the context of their respective functional specializations\u003csup\u003e\u003cspan additionalcitationids=\"CR3 CR4\" citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e–\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e. Skeletally, the human hand evolved a large pollex (i.e., thumb) and a reduced length of the other digits, which improved manipulative capacity, a key factor in successful tool use for food acquisition and processing, and defense\u003csup\u003e\u003cspan additionalcitationids=\"CR7 CR8 CR9\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e–\u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e10\u003c/span\u003e\u003c/sup\u003e. The human foot evolved a robust, adducted hallux (i.e., big toe) and short digits, creating a shorter, hinge-like forefoot with a stiff midfoot that improves performance in both walking and running\u003csup\u003e\u003cspan additionalcitationids=\"CR12 CR13 CR14\" citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e–\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e. These derived features contribute to the unique status of humans as the only extant primates that are obligate bipeds.\u003c/p\u003e\u003cp\u003eDue to their differences in structure and function, the evolutionary histories of the human hand and foot have largely been considered separately\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan additionalcitationids=\"CR7 CR8 CR9 CR10 CR11 CR12 CR13\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e–\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan additionalcitationids=\"CR17\" citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e–\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e. However, data from both developmental genetics and hominoid anatomical studies provide evidence that the evolution of the hand and the foot should be considered together: hands and feet (hereafter collectively referred to as autopods) are serially homologous structures, share much of their genetic architecture\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e, and have been shown to covary strongly in size and shape\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e,\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e. Furthermore, when the anatomical changes between human and chimpanzee autopod skeletons are compared for both limb types, it becomes clear that some changes are remarkably similar – e.g., a larger, more robust first digit compared to the other digits, shorter phalanges compared to metapodials, loss of phalangeal curvature, and an overall size reduction (Fig.\u0026nbsp;1a-c). This raises the possibility that phenotypic changes in either the hand or foot might in part be explained by adaptive evolution in the other due to the potentially constraining effects of strong covariation between homologous limb elements of humans and other apes\u003csup\u003e\u003cspan additionalcitationids=\"CR21\" citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e–\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e. Simulations even suggest that selection for a modern human-like foot from a chimp-like ancestor is sufficient to produce a human-like hand purely as a byproduct due to anatomical integration of the autopods\u003csup\u003e\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eDespite the importance of the hand and foot in human evolution, very little is known about the genetic mechanisms patterning human-specific features of the hand and foot skeleton. Many of these features arise early in human development when the limb skeleton is first patterned in cartilage templates (Fig.\u0026nbsp;1d). While some studies have investigated the developmental genetics of human limbs during or slightly before this developmental window\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e,\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e, none have generated matched datasets for the hand and foot, nor for individual skeletal elements of each autopod. Therefore, to shed new light on the evolutionary history of the human hand and foot skeleton, we generated data on the gene expression and regulation of each metapodial and the pooled phalanges of each digit in both the hand and foot at two stages during this developmental window. Additionally, we performed corresponding experiments in stage-matched E15.5 mouse embryos on digits I, III, and V to explore evolutionary conservation of identified genes and regulatory elements. Using these datasets, we investigated the signature of natural selection on genomic regions patterning distinct elements of the human hand and foot skeleton. We tested (1) the prediction that there would be substantial overlap in gene expression and the genomic regulatory landscape between the hand and foot tissues (the genetic underpinnings of covariance and coevolution), and (2) that the foot would show stronger signatures of selection than the hand, in line with previous hypotheses and the fact the remodeling of the foot in the context of bipedalism was more extensive and predates the appearance of lithic technology\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003e\u003cem\u003eProportional changes in the autopod skeleton\u003c/em\u003e\u003c/p\u003e\u003cp\u003eTo quantitatively assess anatomical differences in proportion across the autopod skeleton, we first estimated the average percent difference in the length and width of each element in a large dataset of adult human, chimpanzee, and gorilla skeletons (See Methods) (Supplementary Table\u0026nbsp;1). As expected, the greatest amount of difference was observed in the more robust first digit in humans when compared with chimps (Fig.\u0026nbsp;1b-c). The other digits were generally shorter – especially the foot phalanges – with the exception of the metacarpals, which were similar in length between the two species (Figs.\u0026nbsp;1b-c). The same pattern is observed when humans are compared to gorillas (Extended Data Fig.\u0026nbsp;1).\u003c/p\u003e\u003cp\u003e\u003cem\u003ePatterns of gene expression and regulation across the autopod\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe identified proportional differences most likely reflect modifications to autopod development in each species, especially at the cartilage level. This is because the rudiments of each digital element are comprised of chondrocytes that organize to form growth plates governing longitudinal and transverse growth during ontogeny. We examined human foot and hand development and identified two gestational timepoints (E54 and E67) spanning a window where morphogenesis was occurring rapidly (Fig.\u0026nbsp;1d, Supplementary Note). For each digital ray of the human hand and foot, we generated matched transcriptomic (RNA-seq) and epigenomic (ATAC-seq) datasets for the phalangeal and metapodial regions separately (Supplementary Tables\u0026nbsp;2–3, Supplementary Results S1-3). For comparative purposes, we performed a similar approach on stage-matched mouse E15.5 forelimb and hind limb elements, focusing on digits I, III, and V (Extended Data Fig.\u0026nbsp;2, Tables S3-5, Supplementary Results S4-5). Analysis of our human RNA-seq dataset identified 1,263 differentially expressed genes (DEGs) between biologically relevant sets of autopod tissues (See Methods; Supplementary Table\u0026nbsp;2). These comparisons included adjacent tissues (along either the anterior-posterior or proximal-distal axes), homologous tissues between the hand and foot, and the same tissue across timepoints, as well as gene expression gradients across the autopod. Many DEGs separated the metapodials and phalanges but not individual digits (Fig.\u0026nbsp;2a; Extended Data Fig.\u0026nbsp;3). As predicted based on the patterns of covariance of coevolution, gene expression was overwhelmingly similar between the hand and foot, consistent with the expression patterns we observed in mouse (Extended Data Fig.\u0026nbsp;2)\u003csup\u003e\u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e26\u003c/span\u003e,\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e. We also identified 3,031 DEGs between the two timepoints, with 1,453 genes upregulated early at E54 and 1,578 genes upregulated late at E67 (Fig.\u0026nbsp;2b). The later developmental time point was enriched for gene ontology (GO) terms related to “ossification,” which initiates around this stage.\u003c/p\u003e\u003cp\u003eUsing our human ATAC-seq data, we identified 60,150 regulatory elements accessible in the chondrocytes of at least one autopod skeletal element (Supplementary Table\u0026nbsp;3). Consistent with the RNA-seq data, principal component analysis (PCA) separated samples between the metapodials and phalanges and between timepoints but not by limb type (Fig.\u0026nbsp;2c). Indeed, no elements uniquely defined all hand or all foot tissues. Overall, regulatory elements are often accessible in many autopod tissues or tissue and timepoint specific (Fig.\u0026nbsp;2d). Using our previously collected stage-matched ATAC-seq datasets from human long bones, girdles, and the axial column\u003csup\u003e\u003cspan additionalcitationids=\"CR29\" citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e–\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e, we found that many of the identified regulatory elements show accessibility in at least one other developing skeletal tissue, albeit 12,488 regulatory elements have not been previously identified and appear unique to the autopod (Fig.\u0026nbsp;2d). To help elucidate broader regulatory patterns in functionally conserved sequences in other primates, we asked how many of the human autopod elements overlapped elements present in mouse, since elements conserved between human and mouse are likely present in other primates. While the number of individual rays examined in mouse was 3 versus 5 in humans, we found that 20,491 regulatory elements were accessible in both human and mouse (34.0% of all regulatory elements identified in human), only 2,654 of which were autopod-specific (21.3% of all human autopod-specific elements) (Extended Data Fig.\u0026nbsp;4).\u003c/p\u003e\u003cp\u003e\u003cem\u003eGenomic evolution in autopod regulatory elements\u003c/em\u003e\u003c/p\u003e\u003cp\u003eUsing our dataset, we sought to investigate the genomic basis of human-specific features of the autopod skeleton. To this end, we overlapped our autopod regulatory element sets with genomic regions marked by signatures of human-specific changes, including human accelerated regions (HARs)\u003csup\u003e\u003cspan additionalcitationids=\"CR32 CR33 CR34 CR35 CR36 CR37\" citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e–\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e, human conserved element deletions (hCONDELs)\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e,\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e, human ancestor quickly evolved regions (HAQERs), large human-specific inversions, and structurally divergent regions between humans and other apes (SDRs)\u003csup\u003e\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e (Table\u0026nbsp;1, Extended Data Fig.\u0026nbsp;5–6). As structural changes can impact neighboring genomic loci by rearranging the local chromatin environment, overlaps were also determined for regions flanking each inversion and SDRs (Extended Data Fig.\u0026nbsp;7). We identified overlaps for these different types of genomic regions for elements with accessibility in multiple skeletal or autopod tissues, as well as those that are tissue and timepoint-specific (Fig.\u0026nbsp;3). Overall, 95% of the genes expressed in at least one autopod tissue fell within 500kb of one or more of these genomic regions. While all tissues exhibited accessibility for at least a few regulatory elements overlapping each type of genomic feature, this was not true for tissue-specific elements (Fig.\u0026nbsp;3f,i,j,o,r). Tissues with high levels of morphological change were not distinctly enriched for evolutionary signals, with relatively few overlaps detected in digit I but large numbers of overlaps seen in the foot phalanges of digits III and V (Fig.\u0026nbsp;3d-f,j-o). This analysis revealed that thousands of genomic loci with accessibility in the developing human autopod skeleton harbor evolutionary sequence alterations.\u003c/p\u003e\u003cp\u003eTo link these genomic changes to anatomical subdivisions of the autopod, we sought to partition our regulatory elements and overlaps according to their accessibility patterns. Based on the differences in gene expression and regulation between the phalanges and metapodials and between developmental stages, as well as the anatomical differences between the hand and foot, we explored the relative strength of the evolutionary signals along each of these axes (Fig.\u0026nbsp;4a-c). Given that evolutionary forces, such as natural selection, shape the patterns of genetic variation present within a species, we also compared the patterns of human, chimpanzee, and gorilla intraspecific variation within each regulatory set (See Methods). Enrichments within each comparison differed depending on the genomic feature and the specificity of the regulatory element set considered (i.e., all brain-filtered elements, autopod-specific elements, or autopod-specific, human-mouse conserved) (Fig.\u0026nbsp;4, Extended Data Figs.\u0026nbsp;8–9). Within the autopod-specific sets, the most significant differences in enrichments were identified for the hand versus foot comparison (Fig.\u0026nbsp;4c).\u003c/p\u003e\u003cp\u003eFor the autopod-specific elements, we sought to explore biological differences between these sets. To test whether different sets regulated the same or distinct target genes, we counted the number of regulatory elements from each set that fell within 500kb of genes expressed in at least one autopod tissue (Additional Data Fig.\u0026nbsp;10). As expected from the largely similar pattern of gene expression observed in all tissues, we found that the number of nearby regulatory elements from the shared sets were largely correlated with those of the more restricted ones. For the more specific sets, the numbers of hand and foot elements were not correlated (Pearson’s Correlation, R = 0.02, P = 0.365), while the numbers of phalangeal and metapodial elements were positively correlated (R = 0.11, P \u0026lt; 0.001) and the early and late elements negatively correlated, although given the strength of the correlation, this is unlikely to be biologically meaningful (R = -0.075, P \u0026lt; 0.001). Overall, many genes are near regulatory elements from multiple categories, consistent with complex regulatory control of gene expression. Accordingly, gene ontology enrichments for each set differed slightly but all included terms related to cartilage or bone development (Supplementary Table\u0026nbsp;3). Finally, we analyzed transcription factor (TF) binding profiles and found that the sets were differentially enriched for hundreds of TF motifs, suggesting that regulatory elements in each set preferentially interact with distinct sets of TFs (Supplementary Table\u0026nbsp;3).\u003c/p\u003e\u003cp\u003eWhile hotspots of human-specific sequence changes are frequently prioritized, single nucleotide changes may also have substantial phenotypic consequences. Accordingly, we sought to estimate the total number of nucleotide changes that may have contributed to human autopod evolution. Briefly, we took our regulatory elements, removed polymorphic variant positions and misaligned regions, and counted the number of positions where the reference nucleotide in the human genome differed from that present in both the chimp and gorilla genomes (see Methods). While this approach ignores structurally divergent and misaligned regions, it gives a rough minimum estimate of the number of fixed substitutions that may contribute to human-specific differences. Using this approach, we identified 2,461,335 human-specific fixed substitutions that are accessible in the developing autopod skeleton, 291,867 of which are autopod-specific. Since humans and chimpanzees share the same last common ancestor to gorilla and therefore have evolved independently for the same length of time (i.e., the branch lengths are equal), we repeated this process using the chimpanzee branch. We identified 2,012,269 chimp-specific fixed substitutions that are accessible in the developing autopod skeleton, 236,518 of which are autopod-specific. Many of these changes likely resulted from neutral evolutionary processes and have minimal or no phenotypic effects, however, comparison of PhyloP scores between the human and chimp sets found a significant shift towards more positive scores in human (human mean 0.23, chimp mean 0.01, \u003cem\u003eP\u003c/em\u003e \u0026lt; 2.2 × 10\u003csup\u003e− 16\u003c/sup\u003e, two-sided Wilcoxon rank-sum), suggesting that human fixed changes alter positions with deeper evolutionary constraint. This holds true for the autopod-specific set as well (human mean 0.24, chimp mean 0.05, \u003cem\u003eP\u003c/em\u003e \u0026lt; 2.2 × 10\u003csup\u003e− 16\u003c/sup\u003e). Overall, this analysis reveals that ~ 450,000 more fixed substitutions in autopod regulatory elements occurred on the human branch than the chimpanzee branch, with ~ 55,000 accessible only in the autopod. We partitioned these fixed substitutions along the same anatomical axes discussed above, finding that except for the elements unique to phalangeal and or late tissues, all sets showed more fixed substitutions along the human branch (Fig.\u0026nbsp;4d).\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eThe evolution of bipedalism is a critical event in hominin evolution and involved substantial remodeling of the postcranial skeleton. To further our understanding of the genetic mechanisms underlying these skeletal changes, we used functional genomics methods to investigate gene expression and regulation during development of the distal autopod skeleton in human. We identified thousands of regulatory elements and genes involved in the developing human hand and foot skeleton. These data revealed substantial differences in gene expression and regulation between the phalanges and metapodials and between the two timepoints but not between digits or limb types. While some of the identified regulatory elements are highly spatially and temporally restricted, most are accessible in multiple skeletal tissues, both within the autopod and across the postcranial skeleton. These results provide a genetic basis for the observed patterns of covariation across the autopod skeleton.\u003c/p\u003e\u003cp\u003eWe found that both multi-tissue and tissue-specific regulatory elements overlap thousands of genomic features that differ between humans and chimpanzees and could potentially underlie the evolution of human-specific autopod traits; these genomic features include human-specific nucleotide changes, losses or gains of sequence, as well as structural variation. Elements broadly accessible in autopod tissues provide a genetic basis for models of human limb evolution based on strong morphological integration (i.e. elevated phenotypic covariance) between homologous elements of the fore- and hind limbs\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e,\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e. While autopod covariation is relatively weaker in comparison to quadrupedal monkeys, homologous elements of the hands and feet of humans and other apes still covary substantially, suggesting the evolution in one limb could be driven at least in part by selection on the other\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e,\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e,\u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e42\u003c/span\u003e\u003c/sup\u003e. Indeed, simulations suggest that selection for human-like foot proportions in a chimpanzee-like ancestor is sufficient to drive the evolution of the human hand, while the converse is not true\u003csup\u003e\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e. The data presented in this study provide some support for this model. While we found many human-specific genomic features in regions accessible in multiple tissues, consistent with phenotypic covariation between elements, support for stronger selection on the foot depends on the evolutionary signal considered. In regulatory elements accessible only in autopod tissues, we found significantly more HARs, hCONDELs, and inversions but fewer SDRs within regulatory elements patterning the foot than the hand (Fig.\u0026nbsp;4b,c). We also found lower levels of intraspecific variation in the foot than the hand in human, chimpanzee, and gorilla, consistent with stronger stabilizing selection on the foot in each lineage. Finally, while there are overall more fixed substitutions found in foot elements than hand elements, both absolutely and per base pair, there is a greater excess of human compared to chimpanzee fixed substitutions in the hand (Fig.\u0026nbsp;4d). These patterns generally hold when considering the sets of all elements, autopod-specific elements, or human-mouse conserved autopod elements, but there are some differences, such as significantly more HAQERs overlapping hand than foot regulatory elements in the total set of brain-filtered elements. An important caveat of our approach is that as we rely on accessibility data from extant human samples, these patterns may not reflect the ancestral state.\u003c/p\u003e\u003cp\u003eAltogether, we estimated that there are almost 2.5\u0026nbsp;million fixed nucleotide substitutions in the human genome compared to chimpanzee and gorilla that are accessible to transcriptional machinery in the developing autopod skeleton. This is ~ 450,000 more fixed nucleotides than are fixed in chimpanzee. Together with the fact that the human fixed nucleotides disrupt more evolutionarily conserved positions, this suggests a potentially greater extent of anatomical evolution along the human branch. ~55,000 of these fixed substitutions are uniquely accessible in autopod tissues, suggesting that autopod evolution likely has a complex underlying genetic architecture, even without considering integration with other elements of the skeleton. To the best of our knowledge, this is the first attempt to assign a minimum number of genome-wide fixed substitutions to a given set of human-specific phenotypic changes and highlights the extent of genomic changes that have potentially contributed to the evolution of the human autopod skeleton.\u003c/p\u003e\u003cp\u003eLinking these genomic changes to observed evolutionary changes in anatomy remains a major challenge\u003csup\u003e\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e,\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e44\u003c/span\u003e\u003c/sup\u003e. While regulatory element sets that are uniquely accessible in tissues with a particular trait might offer the most tempting candidate(s) for identifying the genomic changes underlying evolution in that trait, this approach to connecting genotype to phenotype excludes large numbers of regulatory elements that are biologically plausible by implicitly assuming that traits evolve independently\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e\u003c/sup\u003e. In fact, our data demonstrate that many disparate skeletal elements share an underlying genetic architecture. These shared regulatory elements could be targeted by selection on a single skeletal trait so long as the off-target effects were minimal or could reflect more complex patterns of coevolution across skeletal elements. When we partition our data between the hands and the feet, between timepoints, or between phalanges and metapodials, in each case we find many shared regulatory elements that overlap human-specific changes just within the autopod-specific set. Many more regulatory elements are potentially evolutionarily relevant when elements shared with other postcranial bones are considered and even more would be identified if we removed the brain filtering step of our data processing pipeline. While substantially more challenging to interpret, elements accessible in multiple tissues are an important part of the genomic architecture patterning these skeletal elements.\u003c/p\u003e\u003cp\u003eOf course, many of the human-specific sequence features highlighted in this study likely have no impact on gene expression. While estimating the number and magnitude of effect on gene expression for all these sequences features is beyond the capacities of current technologies, limited experimental research has shown that for regulatory HARs for which both the human and chimp sequence have been tested in transgenic animals, a third resulted in differential enhancer activity\u003csup\u003e\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e\u003c/sup\u003e. Massively parallel reporter assay testing of hCONDELs found that 8% exhibited differential enhancer activity between human and chimps in at least one cell line\u003csup\u003e\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e. If HARs and hCONDELs show similar levels of activity in the developing skeleton, that would suggest that 158 HARs and 71 hCONDELs have altered regulatory activity in the developing human autopod skeleton (considering elements with accessibility in both autopod and non-autopod tissues, 25 HARs and 13 hCONDELs with accessible only in the autopod). Even if HARs and hCONDELs are substantially more predictive of human-specific gene regulatory evolution than the other sequence features considered, an assumption that has not been rigorously tested, this still suggests that the genomic basis of autopod evolution involved hundreds of genomic loci and should be considered highly polygenic. Under such a model, each individual genomic change is expected to have only a small phenotypic contribution, potentially complicating attempts to validate their effects in humanized mouse models\u003csup\u003e\u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eTo fully understand the genetic underpinning of human autopod evolution, future studies are necessary to generate comparable data on the bones of the wrist, midfoot, and hindfoot. These proximal autopod bones have also undergone substantial evolution in the human lineage\u003csup\u003e\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e,\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e47\u003c/span\u003e,\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e and are undoubtedly directly relevant to the evolution of the metapodials and structure of the autopod in general. Similarly, future studies should investigate the genomic changes underlying evolution in the soft tissues of the autopod. As this study was based on bulk chondrocyte tissue dissections, future research would also benefit from techniques with greater spatial resolution, such as spatial multi-omics, to unravel the genetics basis of these human-specific features with greater anatomical precision (e.g., Senevirathne et al., under review). Such methods could investigate human-specific changes in shape related to load bearing and the articulation of bones with one another. Furthermore, these methods could explore the impacts of extrinsic signaling and mechanical impacts from neighboring tissues, such as mesenchyme and muscle, on autopod skeletal development\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e\u003c/sup\u003e. Given the intricate functional, developmental, and evolutionary relationships between all parts of the autopod, a holistic approach that integrates datasets from multiple tissue types will ultimately be necessary to understand the evolution of human-specific features. While only part of the story, this study provides novel insights into the evolutionary history of human bipedalism by highlighting the complex genetic architecture underlying human autopod evolution in the phalanges and metapodials.\u003c/p\u003e"},{"header":"Materials and Methods","content":"\u003cp\u003e\u003cem\u003eEstimated human-chimpanzee anatomical change\u003c/em\u003e\u003c/p\u003e\u003cp\u003eMeasurements were taken of the length and head width of metapodials, proximal, and middle phalanges of the hand and foot in adult human (\u003cem\u003en\u003c/em\u003e = 48), chimpanzee (\u003cem\u003en\u003c/em\u003e = 44), and gorilla (\u003cem\u003en\u003c/em\u003e = 57) skeletons (Appendix, Table\u0026nbsp;5). Chimpanzee and gorilla measurements come from wild-shot populations housed in the Powell-Cotton Museum and Hamman Todd Osteological Collection at the Cleveland Museum of Natural History, while human measurements are from the Hamman Todd Collection only. See ref\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e for details. The proximal and middle phalanges measurements were combined for each digit. The human-chimpanzee percent change for the length and width of each element was calculated as the mean human measurement/the chimpanzee measurement * 100. These calculations were repeated using the gorilla measurements in place of chimpanzee and to compare chimpanzee to gorilla.\u003c/p\u003e\u003cp\u003e\u003cem\u003eHuman developmental sample collection.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eHuman developmental samples were collected from first-trimester termination through the Birth Defects Research Laboratory (BDRL) at the University of Washington in full compliance with the ethical guidelines of the National Institutes of Health (NIH) and with the approval of the University of Washington Institutional Review Boards (IRB) for the collection and distribution of human tissues for research and Harvard University for the receipt and use of such materials. The BDRL obtained written consent from all tissue donors. Harvard University IRB determined that these samples constitute Non-Human Subjects Determination Status (Capellini: IRB16-1504). The fresh human samples were briefly washed in Hanks’ balanced salt solution and shipped at 4˚C. Upon arrival, the samples were immediately dissected under a light dissection microscope and directly subjected to RNA-seq or ATAC-seq protocols described below, following approved Harvard University IRB (IRB16-1504) and Committee on Microbiological Safety (COMS) (18–103) protocols. The BDRL performs polymerase chain reaction (PCR) using \u003cem\u003eSRY\u003c/em\u003e and \u003cem\u003eAmelogenin\u003c/em\u003e primers to determine the biological sex of each sample.\u003c/p\u003e\u003cp\u003e\u003cem\u003eHuman ATAC-seq data collection and processing\u003c/em\u003e.\u003c/p\u003e\u003cp\u003eDevelopmental samples were microdissected under a light microscope in 5% fetal bovine serum (FBS) in Dulbecco’s Modified Eagle Medium (DMEM) on ice. The pooled phalanges (proximal and distal for digit I; proximal, intermediate, and distal for digits II-V) and metapodials were collected separately for all digits in both the fore- and hind limbs. Samples were collected at two timepoints: an early stage (E53 to E59, \u003cem\u003en\u003c/em\u003e ≥3) and a later stage (E67 to E74, \u003cem\u003en\u003c/em\u003e ≥ 3) to capture the window of time when the autopod skeleton is prepatterned in cartilage and is just beginning to ossify and to account for heterogeneity in the timing of forelimb or hind limb as compared to overall limb development \u003csup\u003e\u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e49\u003c/span\u003e,\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eThe tissue in each sample was then digested using 0.5% collagenase II in 5% FBS/DMEM in a 37˚C water bath for one hour. Every 30 minutes, the samples were spun down and gently pipetted to break up large clumps of cells. Next, the samples were incubated at 37˚C for an additional hour in a shaking incubator. Following these incubation steps, the sample were immediately placed on ice and then filtered through a 70 µm nylon cell strainer into a 50 mL conical tube by gently pressing any residual tissue through the filter followed by rinsing with 5% FBS/DMEM. Sample were centrifuged for 5 minutes at 500\u003cem\u003eg\u003c/em\u003e at 4˚C and most of the media aspirated. The cells were then resuspended and transferred to a 1.5 mL tube before another 5-minute centrifuge at 500\u003cem\u003eg\u003c/em\u003e at 4˚C. At this stage, live and dead cells were counted to ensure that all tissues had ~ 50,000 cells and cell death rates below 10%. 50,000 cells per technical replicated were resuspended in 1x PBS and then lysed using the ATAC-seq lysis buffer and centrifuged for 10 minutes at 4˚C \u003csup\u003e\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e\u003c/sup\u003e. The supernatant was discarded and 50 µl transposition mix containing 2.5 µL transposase was added before the samples were incubated at 37˚C for 30 minutes. Lastly, the Zymo DNA Clean and Concentrator kit was used, and the resulting purified DNA was eluted in 13 µl of water heated to 60–70˚C. Samples were stored at -20˚C prior to PCR amplification and barcoding. Samples were amplified and given an 8 bp barcode via 11 cycles of PCR amplification using the NEBNext High-Fidelity 2x PCR master mix. Following amplification, fragments were size-selected using Mag-Bind® RxnPure Plus beads. The sample pool was sequenced on a NovaSeq S4 machine at the Harvard Bauer Core to generate ≥40\u0026nbsp;million reads per sample.\u003c/p\u003e\u003cp\u003eFastQC version 11.9 was used to determine the quality of fastq read files. Because some samples had been run on multiple lanes to achieve the desired number of reads, the reads for each sample were concatenated into a single file for R1 and a second file for R2. NGmerge version 0.3 was used to trim adapters from all reads \u003csup\u003e\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e\u003c/sup\u003e. Reads were then aligned to the Illumina prebuilt hg38 human reference genome using Bowtie2 (version 2.3.4.1) \u003csup\u003e\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e,\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e\u003c/sup\u003e. Reads were then indexed using \u003cem\u003esamtools index\u003c/em\u003e \u003csup\u003e\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e and duplicates removed using \u003cem\u003epicard MarkDuplicates\u003c/em\u003e (version 2.9.0). The resulting files were indexed as before and mitochondrial reads were filtered out (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/harvardinformatics/ATAC-seq/blob/master/atacseq/removeChrom.py\u003c/span\u003e\u003cspan address=\"https://github.com/harvardinformatics/ATAC-seq/blob/master/atacseq/removeChrom.py\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). The resulting .bam files were then used for peak calling via \u003cem\u003eMACS\u003c/em\u003e software (version 2.1.1.2), using BAMPE and the following flags: \u003cem\u003e--nolambda –bdg –verbose\u003c/em\u003e\u003csup\u003e\u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e56\u003c/span\u003e\u003c/sup\u003e. Reproducible peaks across replicates were identified at an IDR threshold of \u0026lt; 0.05, as defined by the IDR statistical test (version 2.0.3). Finally, IDR-called peak sets for autopod skeletal elements were filtered by peaks previously identified in E54 human brain to remove regulatory regions likely associated with general cellular housekeeping processes\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. Peak sets were merged, subtracted, overlapped, etc. using the appropriate \u003cem\u003ebedtools\u003c/em\u003e (v2.27.1) functions\u003csup\u003e\u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003e\u003cem\u003eHuman RNA-seq data collection and processing.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eHuman samples as described above were dissected in 10% FBS/DMEM and a minimum of \u003cem\u003en\u003c/em\u003e = 6 per tissue per timepoint was collected. Each element was stripped of soft tissue and collected in a 2-ml tube containing 200 µl of TRIzol and one 5-mm stainless steel bead. Each sample was then homogenized at 50 Hz for 2 min, followed by 1 min on ice, and a second homogenization at 50 Hz for 2 min. Samples were stored at -80˚C until RNA extraction was performed. For each RNA extraction, samples were incubated at room temperature for 5 min and then centrifuged at 4˚C for 5 minutes at 12,000 x g. The supernatant was collected and transferred to a clean tube, to which 200 µl of chloroform was added per 1 ml of TRIzol. Each sample was then vortexed briefly and incubated at room temperature for 2 minutes before being transferred to a MaXtract tube and centrifuged at 4˚C for 5 minutes at 12,000 x g. Following centrifugation, the aqueous phase was removed and transferred to a new microcentrifuge Eppendorf tube and an equal volume of 100% ethanol was added. After completion of this phenol-chloroform step, the RNA extraction continued using the Direct-zol RNA MicroPrep kit following the manufacturer’s protocols and eluted in 15 µl of nuclease-free water. Samples were stored at -80˚C until library preparation. To maximize the number of tissues useable from each biological replicate, samples with RNA integrity number (RIN) scores higher than 6 were used for subsequent steps so long as the average RIN for all tissues from that sample was greater than 7.\u003c/p\u003e\u003cp\u003eFor cDNA library preparation and sequencing, samples were normalized to a single concentration and libraries were prepared using the Kapa mRNA HyperPrep kit with an input volume of 25 ng per sample following the manufacturer’s protocols. After an initial MiSeq Nano run to assess library quality, samples were repooled as necessary and then sequenced on an Illumina NovaSeq S4 six times to generate ~ 20\u0026nbsp;million paired-end reads per sample. See SI Appendix, Table S2 for detailed sequencing information for each sample, including index primer, read count, and other information. Six biological replicates per tissue per timepoint were sequenced.\u003c/p\u003e\u003cp\u003eFastQC version 11.9 was used to assess the quality of each fastq read file\u003csup\u003e\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e58\u003c/span\u003e\u003c/sup\u003e. Reads for each sample from multiples lanes were concatenated into a single file for R1 and a second file for R2. NGmerge version 0.3 was used to trim adapters from all reads\u003csup\u003e\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e\u003c/sup\u003e. Next, reads mapping to ribosomal RNA were removed using RiboDetector\u003csup\u003e\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e\u003c/sup\u003e version 0.2.7. STAR\u003csup\u003e\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e\u003c/sup\u003e version 2.7.1 was used to map reads to the human genome (hg38) with ≥80% of reads uniquely mapping for each sample. RSEM\u003csup\u003e\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e61\u003c/span\u003e\u003c/sup\u003e version 1.3.3 was then used to generate read counts. For all samples, \u0026gt; 90% of reads were mapped uniquely to a gene.\u003c/p\u003e\u003cp\u003e\u003cem\u003eMouse ATAC-seq data collection and processing.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eAll mouse work was covered under the Capellini lab IACUC protocol (#13-04-161-3). E15.5 mouse embryos were collected from pregnant FVB/NJ females and transferred to cold 1x PBS. The fore- and hind limbs of three to five embryos were microdissected under a light microscope in 5% FBS/DMEM on ice. The pooled phalanges (proximal and distal for digit I, proximal, intermediate, and distal for digits III and V) and metapodials were collected separately for digits I, III, and V for both the forelimb and the hind limb. The right and left sides of each element were pooled together, resulting in a total pool of elements from six to ten individual limbs. All embryos in each biological replicate were from the same litter. Collagenase digestion, cell lysis, the transposase reaction, barcoding, and sequencing were all performed as described above for human samples. Data processing was also the same except that reads were aligned to the prebuild Illumina UCSC mm10 genome, E15.5 mouse brain accessibility data was used for the brain-filtering step\u003csup\u003e\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e62\u003c/span\u003e\u003c/sup\u003e, and a different version of MACS (3.0.3) was used. See Supplementary Table\u0026nbsp;5 for detailed sequencing information for each sample, including index primer, read count, and other information.\u003c/p\u003e\u003cp\u003e\u003cem\u003eMouse RNA-seq data collection and processing.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eE15.5 mouse embryos were collected from pregnant FVB/NJ females and transferred to cold 1x phosphate-buffer saline (PBS). Embryos were microdissected under a light microscope in 10% FBS/DMEM on ice and anatomically sexed\u003csup\u003e\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e63\u003c/span\u003e\u003c/sup\u003e. The pooled phalanges (proximal and distal for digit I; proximal, intermediate, and distal for digits III and V) and metapodials were collected separately for digits I, III, and V for both the fore- and hind limbs. The right and left sides of each element were pooled together. All other elements of the RNA collection and extraction protocol were the same as described for human above.\u003c/p\u003e\u003cp\u003eFor cDNA library preparation and sequencing, samples were normalized to a single concentration and libraries were prepared using the Takara SMART-Seq v4 Ultra Low Input RNA Kit following the manufacturer’s protocols. Two samples with high concentrations were also prepared using the Kapa mRNA HyperPrep kit with an input volume of 25 ng per sample as per the human samples to detect any effects of library preparation method. After initial pooling, the relative concentrations of each sample were evaluated by running the pool on a MiSeq Nano and the pool concentrated adjusted accordingly. The final pool was then sequenced four times on an Illumina NovaSeq S4 to generate ≥40\u0026nbsp;million paired-end reads per sample. On average, samples were sequenced at 86\u0026nbsp;million reads per sample, within the recommended Encyclopedia of DNA Elements (ENCODE) guidelines (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ewww.encodeproject.org/about/experiment-guidelines/\u003c/span\u003e\u003c/span\u003e). See Supplementary Table\u0026nbsp;4 for detailed sequencing information for each sample, including index primer, read count, and other information. Six biological replicates of each tissue were sequenced.\u003c/p\u003e\u003cp\u003e\u003cem\u003eDifferential expression and cluster analysis of RNA-seq.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eThe large number of tissues investigated in this study (20 tissues at two timepoints in human, 12 in mouse) allows for many potentially pair-wise comparisons, however, only a subset of these comparisons is expected to be biologically meaningful. To minimize the number of factors that could cause gene expression differences, samples from a given tissue were only compared with directly adjacent tissues and with the homologous tissue in the other limb type. For humans, adjacent tissues came from sequential digits while in mouse, only digits I, III, V were collected so the nearest digit was used in place of the adjacent digit. Each human tissue was also compared across timepoints. This comparison scheme aims to target only a single variable for any given differential gene calculation, either spatial location within the autopod, limb type, or timepoint. To give an example, the gene expression profile of early samples of hand phalanges from digit III were compared separately with (1) the early hand phalanges of digit II, (2) the early hand phalanges of digit IV, (3) the early metacarpal of digit III, (4) the early foot phalanges of digit III, and (5) the late hand phalanges of digit III.\u003c/p\u003e\u003cp\u003eAll differentially expressed genes (DEGs) were identified using DESeq2 \u003csup\u003e64\u003c/sup\u003e version 1.36.0 with a design including the sample type and the replicate ID except for comparisons between timepoints where only timepoint was included in the design. In addition, to identify genes expressed in a gradient across the digits, DESeq2 was rerun at both time points using a likelihood ratio test to compare a model of ~ replicate + digit with a baseline model of ~ replicate. Thresholds for calling DEGs were P\u003csub\u003eadj\u003c/sub\u003e \u0026lt; 0.05 and |fold-change| \u0026gt;1.5. This approach was used for both human and mouse.\u003c/p\u003e\u003cp\u003eMany of the DEGs identified above are expected to be DEGs in numerous comparisons, e.g., a gene highly expressed in a single tissue will be identified as a DEG in all comparisons involving that tissue. Similarly, many genes may be expressed in specific spatial domains encompassing multiple tissues and cannot be easily identified using a pair-wise approach. Therefore, to identify large patterns of gene expression present in our dataset, we performed a weighted gene co-expression network analysis (WGCNA) to identify clusters of genes with similar expression patterns\u003csup\u003e\u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e65\u003c/span\u003e\u003c/sup\u003e. The expression of DEGs was normalized using the variance transformation function \u003cem\u003evst\u003c/em\u003e in DESeq2. For comparison within a single timepoint, variation associated with biological replicates was removed using \u003cem\u003elimma:removeBatchEffect)\u003c/em\u003e\u003csup\u003e\u003cem\u003e\u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e66\u003c/span\u003e\u003c/em\u003e\u003c/sup\u003e. This was not used when comparing both human timepoints because this correction also removed timepoint differences that are of biological interest. The \u003cem\u003epickSoftThreshold\u003c/em\u003e function was used to estimate the lowest scale-free threshold for each dataset, which resulted in a minimum scale-free topology of 0.9. These soft power thresholds were used to construct signed correlation networks using the \u003cem\u003eblockwiseModules\u003c/em\u003e function.\u003c/p\u003e\u003cp\u003eMouse data was processed following the same pipeline except that nearest digits were compared instead of neighboring digits (so digit III was compared to digits I and V instead of II and IV) due to reduced tissue sampling. DEG-calling was repeated for human with the reduced tissue set (digits I, III, V only) at each timepoint to facilitate direct comparison with mouse.\u003c/p\u003e\u003cp\u003e\u003cem\u003eComparison of human and mouse ATAC/RNA-seq.\u003c/em\u003e\u003c/p\u003e\u003cp\u003eMouse orthologs of human genes (GRCh38.p14) were downloaded from the Ensembl BioMart database for comparison of gene sets between the two species. For comparison of orthologous regulatory elements, mouse regions (mm10) were lifted over to human genome coordinates (hg38), allowing for multiple matches and a minMatch = 0.1 threshold.\u003c/p\u003e\u003cp\u003e\u003cem\u003eIdentification of accessibility sharing with other elements of the developing skeleton\u003c/em\u003e\u003c/p\u003e\u003cp\u003eRegulatory elements accessible in the autopod were intersected with elements accessibility in at least one other developmental skeletal tissue for which data has been generated using the same ATAC-seq protocol\u003csup\u003e\u003cspan additionalcitationids=\"CR29\" citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e–\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e. This included the proximal and distal ends of the femur, tibia, humerus, and radius, three regions of the scapula (head and neck, blade, acromion), four regions of the pelvis (ilium, ischium, pubis, acetabulum), as well as the thoracic and the lumbar vertebrae. Data were available at both stages for all tissues except the vertebrae, for which only later stage was available.\u003c/p\u003e\u003cp\u003e\u003cem\u003eIdentification of regulatory overlaps with genomic features\u003c/em\u003e\u003c/p\u003e\u003cp\u003eTo identify potential regulatory elements that were the targets of past positive selection in the human lineage, regulatory element sets were overlapped with human accelerated regions (HARs), which are non-coding regions of the genome marked by accelerated substitution rates in the human lineage compared to other vertebrates\u003csup\u003e\u003cspan additionalcitationids=\"CR32 CR33 CR34 CR35 CR36 CR37\" citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e–\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e, human conserved element deletions (hCONDELs)\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e,\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e, and human ancestor quickly evolved regions (HAQERs), the fasted evolving regions of the human genome\u003csup\u003e\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e. Regulatory element sets were also overlapped with structural divergent regions of the human genome compared to other apes, including large chromosomal inversions and structural divergent regions (SDRs)\u003csup\u003e\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e. As the HAQERs, inversions, and SDR sets were identified using the T2T-CHM13v2.0 genome, these sets were lifted over to hg38 using the liftover file available at \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://github.com/marbl/CHM13\u003c/span\u003e\u003cspan address=\"https://github.com/marbl/CHM13\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e. Regions flanking inversions and SDRs were obtained for either 500kb or 1Mb windows on either side of these features.\u003c/p\u003e\u003cp\u003eTo investigate patterns of intraspecific variation in regulatory elements, the VCF file containing the positions of human SNPs in the dbSNP151 database were downloaded from the NIH website (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://ftp.ncbi.nlm.nih.gov/snp/organisms/\u003c/span\u003e\u003cspan address=\"https://ftp.ncbi.nlm.nih.gov/snp/organisms/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e) and converted to bed format using the BEDOPS (v2.4.40) convert2bed function\u003csup\u003e\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e67\u003c/span\u003e\u003c/sup\u003e. Chimp and gorilla SNP VCF files were downloaded from the Great Apes Genome Diversity Project\u003csup\u003e\u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e68\u003c/span\u003e\u003c/sup\u003e. We curated common (i.e., minor allele frequency \u0026gt; 0.05) single nucleotide polymorphisms (SNP) from each species. All intersections were performed using \u003cem\u003ebedtools intersect\u003c/em\u003e\u003csup\u003e\u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e57\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003eTo compare the number of regulatory elements from distinct sets falling near expressed genes, the coordinates of genes with a minimum normalized expression of 10 in at least one autopod tissue based on the UCSC hg38 refGene list were extended by 500kb in either direction and overlapped with the regulatory element sets using the \u003cem\u003ebedtools intersect -c\u003c/em\u003e function. Statistical significance was assessed using a two-sided Pearson’s correlation test.\u003c/p\u003e\u003cp\u003e\u003cem\u003eGene ontology enrichment analysis\u003c/em\u003e\u003c/p\u003e\u003cp\u003eGene ontology (GO) analysis of gene lists was performed using clusterProfiler with a background set of all genes expressed in at least one of the sampled tissues\u003csup\u003e\u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e69\u003c/span\u003e\u003c/sup\u003e. Enrichment analysis of genomic region sets was performed using rGREAT (v2.6.0) for five ontologies (\"GO Molecular Function\", \"GO Biological Process\", \"GO Cellular Component\", \"Human Phenotype\", \"Mouse Phenotype\") using the standard thresholds used in the online version of tool (binomial and hypergeometric P\u003csub\u003eadj\u003c/sub\u003e \u0026lt; 0.05, binomial fold enrichment ≥ 2)\u003csup\u003e70,71\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003e\u003cem\u003eTranscription factor (TF) binding motif enrichment analysis\u003c/em\u003e\u003c/p\u003e\u003cp\u003eTF binding motif enrichments were calculated in R using MonaLisa\u003csup\u003e\u003cspan citationid=\"CR72\" class=\"CitationRef\"\u003e72\u003c/span\u003e\u003c/sup\u003e with human PWM matrices from the JASPAR database\u003csup\u003e\u003cspan citationid=\"CR73\" class=\"CitationRef\"\u003e73\u003c/span\u003e\u003c/sup\u003e. Statistical significance was assessed using a negative log10 cut-off of 4. Regulatory element sizes were not standardized for these analyses.\u003c/p\u003e\u003cp\u003e\u003cem\u003eCalculation of human chimp sequence divergence\u003c/em\u003e\u003c/p\u003e\u003cp\u003ePrior to performing this analysis, all elements in the genomic region sets were standardized to a length of 500 bps. To calculate the minimum number of single base pair differences between the human and chimpanzee genomes in each set of genomic regions, we first filtered out any positions that are variable within human, chimpanzees, or gorillas since interspecific differences in autopod morphology are much larger than intraspecific variation and therefore current allelic variants are insufficient to explain the evolutionary changes. VCF lists of variants for each species were obtained as described above and SNPs overlapping autopod regions were extracted using tabix (version 1.13)\u003csup\u003e\u003cspan citationid=\"CR74\" class=\"CitationRef\"\u003e74\u003c/span\u003e\u003c/sup\u003e and converted to bed format using bedtools and then lifted over to hg38 (for chimp and gorilla). Genomic regions sets were filtered by ENCODE blacklist regions\u003csup\u003e\u003cspan citationid=\"CR75\" class=\"CitationRef\"\u003e75\u003c/span\u003e\u003c/sup\u003e as well as the SNP sets from each species using bedtools to prioritize high quality genomic sequence without interspecific variants. The resulting regions were then lifted over to the chimpanzee (panTro6) and gorilla (gorGor6) genomes, and the DNA sequence extracted using the appropriate BSgenome package\u003csup\u003e\u003cspan citationid=\"CR76\" class=\"CitationRef\"\u003e76\u003c/span\u003e\u003c/sup\u003e. Regions had to be mappable between species along the entire length of the sequence and have at least 25% sequence identity between either target species (human or chimp) and its outgroups. This threshold was chosen to prioritize comparison of orthologous base pairs rather than faulty alignments or indel-related mismatches because 25% is the expected match for random sequences of the same length. A nucleotide was considered derived in a target lineage if it differed from the nucleotide observed (and likely conserved) in the other two species. This analysis was performed with both human and chimpanzee as the target species. The Zoonomia PhyloP scores for each genomic position were obtained for the human and chimpanzee sets\u003csup\u003e\u003cspan citationid=\"CR77\" class=\"CitationRef\"\u003e77\u003c/span\u003e\u003c/sup\u003e and compared using two-sided Wilcoxon rank-sum test.\u003c/p\u003e\u003cp\u003e\u003cem\u003eSegmentation of the adult autopod scans\u003c/em\u003e\u003c/p\u003e\u003cp\u003eRaw TIFFs and mesh files for corresponding human and chimp autopods were downloaded from MorphoSource and imported into Amira v.2022.1. Slices showing ossification were segmented, and the selected voxels were assigned to a separate material. To visualize the autopod in 3D and generate a model, the “Generate Surface” and “Surface View” options were used in Amira. The final rendered models were exported as object files (.obj) and visualized using Blender 4.4.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eData, Materials, and Software Availability\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAll data needed to evaluate the conclusions in the paper are present in the paper and/or the SI Appendix. In addition, all ATAC-seq data raw sequencing fastq files and processed peak bed files and RNA-seq data raw sequencing fastq files and processed read count files have been deposited on National Center for Biotechnology Information Gene Expression Omnibus (GSE281502, GSE283854, GSE284187, GSE286924). The code used to perform these analyses and to generate these figures is available at: https://github.com/aokamoto-bio/Human_Autopod_Evolution.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003cstrong\u003e\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe would like to thank the University of Washington Birth Defects Research Laboratory including Kimberly A. Aldinger, Dan Doherty, Ian G. Phelps, Jennifer C. Dempsey, Yasmeen Otaibi, Lucinda A. Cort; funded via NIH under NICHD Grant # R24HD000836, to I.A.G.; C. Reardon and the Harvard Bauer Core Facility for ATAC-seq and RNA-seq assistance and sequencing; C. Tabin, D. Lieberman, K. Cooper, D. Richard, and members of the Capellini laboratory (Harvard University) for critical insight and commentary on this project. T.D.C., A.S.O., G.S., and P.M., were supported by the Harvard University Dean\u0026rsquo;s Competitive Fund and Milton Fund for the functional genomics portion of this research. T.D.C. and A.S.O additionally received funds for mouse work from the National Science Foundation (BCS1518596 and BCS1847979). A.S.O and T.D.C. were supported by a National Science Foundation Doctoral Dissertation Improvement Grant (BSC-2337516). A.S.O. was supported by a National Science Foundation Graduate Research Fellowship (DGE-1745303). T.D.C. was supported by the American School of Paleontological Research (ASPR) via Harvard University.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;\u003cstrong\u003eAuthor Contributions:\u0026nbsp;\u003c/strong\u003eA.S.O. performed human autopod sample dissection, ATAC-seq, RNA-seq, downstream data processing, all computational analyses, aided in study design, wrote the manuscript, and generated figures with input from all authors. G.S. and P.M. assisted with sample collection. C.R. provided comparative morphological data and helped write and revise the manuscript. I.G. and B.D.R.L. provided samples. T.D.C. conceived and supervised all aspects of the project, designed the study, and wrote and revised the manuscript with input from all authors.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eDarwin, C. \u003cem\u003eThe Descent of Man, and Selection in Relation to Sex\u003c/em\u003e. (John Murray, London, 1871).\u003c/li\u003e\n\u003cli\u003eDiogo, R., Siomava, N. \u0026amp; Gitton, Y. Development of human limb muscles based on whole-mount immunostaining and the links between ontogeny and evolution. \u003cem\u003eDevelopment\u003c/em\u003e \u003cstrong\u003e146\u003c/strong\u003e, (2019).\u003c/li\u003e\n\u003cli\u003eDiogo, R., Richmond, B. G. \u0026amp; Wood, B. Evolution and homologies of primate and modern human hand and forearm muscles , with notes on thumb movements and tool use. \u003cem\u003eJ Hum Evol\u003c/em\u003e \u003cstrong\u003e63\u003c/strong\u003e, 64\u0026ndash;78 (2012).\u003c/li\u003e\n\u003cli\u003eDiogo, R., Molnar, J. L., Rolian, C. \u0026amp; Esteve-Altava, B. First anatomical network analysis of fore- and hindlimb musculoskeletal modularity in bonobos , common chimpanzees , and humans. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 1\u0026ndash;9 (2018).\u003c/li\u003e\n\u003cli\u003eKivell, T. L., Baraki, N., Lockwood, V., Williams‐Hatala, E. M. \u0026amp; Wood, B. A. Form, function and evolution of the human hand. \u003cem\u003eAmerican Journal of Biological Anthropology\u003c/em\u003e \u003cstrong\u003e181\u003c/strong\u003e, 6\u0026ndash;57 (2023).\u003c/li\u003e\n\u003cli\u003eKivell, T. L. Evidence in hand: recent discoveries and the early evolution of human manual manipulation. \u003cem\u003ePhilosophical Transactions of the Royal Society B\u003c/em\u003e \u003cstrong\u003e370\u003c/strong\u003e, 20150105 (2015).\u003c/li\u003e\n\u003cli\u003eBardo, A., Vigouroux, L., Kivell, T. L. \u0026amp; Pouydebat, E. The impact of hand proportions on tool grip abilities in humans, great apes and fossil hominins: A biomechanical analysis using musculoskeletal simulation. \u003cem\u003eJ Hum Evol\u003c/em\u003e \u003cstrong\u003e125\u003c/strong\u003e, 106\u0026ndash;121 (2018).\u003c/li\u003e\n\u003cli\u003eLiu, M.-J., Xiong, C.-H. \u0026amp; Hu, D. Assessing the manipulative potentials of monkeys, apes and humans from hand proportions: implications for hand evolution. \u003cem\u003eProceedings of the Royal Society B: Biological Sciences\u003c/em\u003e \u003cstrong\u003e283\u003c/strong\u003e, 20161923 (2016).\u003c/li\u003e\n\u003cli\u003eYoung, R. W. Evolution of the human hand: the role of throwing and clubbing. \u003cem\u003eJ Anat\u003c/em\u003e \u003cstrong\u003e202\u003c/strong\u003e, 165\u0026ndash;174 (2003).\u003c/li\u003e\n\u003cli\u003eMarzke, M. W. \u0026amp; Marzke, R. F. Evolution of the human hand: Approaches to acquiring, analysing and interpreting the anatomical evidence. \u003cem\u003eJ Anat\u003c/em\u003e \u003cstrong\u003e197\u003c/strong\u003e, 121\u0026ndash;140 (2000).\u003c/li\u003e\n\u003cli\u003eHolowka, N. B. \u0026amp; Lieberman, D. E. Rethinking the evolution of the human foot: Insights from experimental research. \u003cem\u003eJournal of Experimental Biology\u003c/em\u003e \u003cstrong\u003e221\u003c/strong\u003e, (2018).\u003c/li\u003e\n\u003cli\u003eLieberman, D. E. \u0026amp; Bramble, D. M. Endurance running and the evolution of Homo. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e432\u003c/strong\u003e, 345\u0026ndash;352 (2004).\u003c/li\u003e\n\u003cli\u003eRolian, C., Lieberman, D. E., Hamill, J., Scott, J. W. \u0026amp; Werbel, W. Walking, running and the evolution of short toes in humans. \u003cem\u003eJournal of Experimental Biology\u003c/em\u003e \u003cstrong\u003e212\u003c/strong\u003e, 713\u0026ndash;721 (2009).\u003c/li\u003e\n\u003cli\u003eDeSilva, J., McNutt, E., Benoit, J. \u0026amp; Zipfel, B. One small step: A review of Plio-Pleistocene hominin foot evolution. \u003cem\u003eAm J Phys Anthropol\u003c/em\u003e \u003cstrong\u003e168\u003c/strong\u003e, 63\u0026ndash;140 (2019).\u003c/li\u003e\n\u003cli\u003eHayafune, N., Hayafune, Y. \u0026amp; Jacob, H. A. C. Pressure and force distribution characteristics under the normal foot during the push-off phase in gait. \u003cem\u003eThe Foot\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 88\u0026ndash;92 (1999).\u003c/li\u003e\n\u003cli\u003eAlmecijia, S., Smaers, J. B. \u0026amp; Jungers, W. L. The evolution of human and ape hand proportions. \u003cem\u003eNat Commun\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 1\u0026ndash;11 (2015).\u003c/li\u003e\n\u003cli\u003eTocheri, M. W., Orr, C. M., Jacofsky, M. C. \u0026amp; Marzke, M. W. The evolutionary history of the hominin hand since the last common ancestor of Pan and Homo. \u003cem\u003eJ Anat\u003c/em\u003e \u003cstrong\u003e212\u003c/strong\u003e, 544\u0026ndash;562 (2008).\u003c/li\u003e\n\u003cli\u003eWilliams-Hatala, E. M. \u003cem\u003eet al.\u003c/em\u003e The manual pressures of stone tool behaviors and their implications for the evolution of the human hand. \u003cem\u003eJ Hum Evol\u003c/em\u003e \u003cstrong\u003e119\u003c/strong\u003e, 14\u0026ndash;26 (2018).\u003c/li\u003e\n\u003cli\u003eShou, S., Scott, V., Reed, C., Hitzemann, R. \u0026amp; Stadler, H. S. Transcriptome analysis of the murine forelimb and hindlimb autopod. \u003cem\u003eDev Dyn\u003c/em\u003e \u003cstrong\u003e234\u003c/strong\u003e, 74\u0026ndash;89 (2005).\u003c/li\u003e\n\u003cli\u003eRolian, C. Integration and evolvability in primate hands and feet. \u003cem\u003eEvol Biol\u003c/em\u003e \u003cstrong\u003e36\u003c/strong\u003e, 100\u0026ndash;117 (2009).\u003c/li\u003e\n\u003cli\u003eRolian, C., Lieberman, D. E. \u0026amp; Hallgr\u0026iacute;msson, B. The coevolution of human hands and feet. \u003cem\u003eEvolution (N Y)\u003c/em\u003e \u003cstrong\u003e64\u003c/strong\u003e, 1558\u0026ndash;1568 (2010).\u003c/li\u003e\n\u003cli\u003eAuerbach, B. M., Savell, K. R. R. \u0026amp; Agosto, E. R. Morphology, evolution, and the whole organism imperative: Why evolutionary questions need multi-trait evolutionary quantitative genetics. \u003cem\u003eYearbook Biological Anthropology\u003c/em\u003e \u003cstrong\u003e181 Suppl 76\u003c/strong\u003e, 180\u0026ndash;211 (2023).\u003c/li\u003e\n\u003cli\u003eRolian, C. Re‐evaluating the correlated evolution of human hands and feet using viability selection modeling. \u003cem\u003eAmerican Journal of Biological Anthropology\u003c/em\u003e \u003cstrong\u003e183\u003c/strong\u003e, (2024).\u003c/li\u003e\n\u003cli\u003eCotney, J. \u003cem\u003eet al.\u003c/em\u003e The Evolution of Lineage-Specific Regulatory Activities in the Human Embryonic Limb. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e154\u003c/strong\u003e, 185\u0026ndash;196 (2013).\u003c/li\u003e\n\u003cli\u003eZhang, B. \u003cem\u003eet al.\u003c/em\u003e A human embryonic limb cell atlas resolved in space and time. \u003cem\u003eNature\u003c/em\u003e (2023) doi:10.1038/s41586-023-06806-x.\u003c/li\u003e\n\u003cli\u003eTaher, L. \u003cem\u003eet al.\u003c/em\u003e Global gene expression analysis of murine limb development. \u003cem\u003ePLoS One\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, (2011).\u003c/li\u003e\n\u003cli\u003eCotney, J. \u003cem\u003eet al.\u003c/em\u003e Chromatin state signatures associated with tissue-specific gene expression and enhancer activity in the embryonic limb. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e22\u003c/strong\u003e, 1069\u0026ndash;1080 (2012).\u003c/li\u003e\n\u003cli\u003eRichard, D. \u003cem\u003eet al.\u003c/em\u003e Functional genomics of human skeletal development and the patterning of height heritability. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e188\u003c/strong\u003e, 15-32.e24 (2025).\u003c/li\u003e\n\u003cli\u003eYoung, M. \u003cem\u003eet al.\u003c/em\u003e The developmental impacts of natural selection on human pelvic morphology. \u003cem\u003eSci Adv\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 1\u0026ndash;24 (2022).\u003c/li\u003e\n\u003cli\u003eRichard, D. \u003cem\u003eet al.\u003c/em\u003e Evolutionary Selection and Constraint on Human Knee Chondrocyte Regulation Impacts Osteoarthritis Risk. \u003cem\u003eCell\u003c/em\u003e \u003cstrong\u003e181\u003c/strong\u003e, 362-381.e28 (2020).\u003c/li\u003e\n\u003cli\u003eGittelman, R. M. \u003cem\u003eet al.\u003c/em\u003e Comprehensive identification and analysis of human accelerated regulatory DNA. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e25\u003c/strong\u003e, 1245\u0026ndash;1255 (2015).\u003c/li\u003e\n\u003cli\u003ePollard, K. S. \u003cem\u003eet al.\u003c/em\u003e Forces shaping the fastest evolving regions in the human genome. \u003cem\u003ePLoS Genet\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, 1599\u0026ndash;1611 (2006).\u003c/li\u003e\n\u003cli\u003eBird, C. P. \u003cem\u003eet al.\u003c/em\u003e Fast-evolving noncoding sequences in the human genome. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 1\u0026ndash;12 (2007).\u003c/li\u003e\n\u003cli\u003eBush, E. C. \u0026amp; Lahn, B. T. A genome-wide screen for noncoding elements important in primate evolution. \u003cem\u003eBMC Evol Biol\u003c/em\u003e \u003cstrong\u003e8\u003c/strong\u003e, 1\u0026ndash;10 (2008).\u003c/li\u003e\n\u003cli\u003ePrabhakar, S., Noonan, J. P., P\u0026auml;\u0026auml;bo, S. \u0026amp; Rubin, E. M. Accelerated evolution of conserved noncoding sequences in humans. \u003cem\u003eScience (1979)\u003c/em\u003e \u003cstrong\u003e314\u003c/strong\u003e, 786 (2006).\u003c/li\u003e\n\u003cli\u003eBi, X. \u003cem\u003eet al.\u003c/em\u003e Lineage-specific accelerated sequences underlying primate evolution. \u003cem\u003eSci Adv\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, (2023).\u003c/li\u003e\n\u003cli\u003eKeough, K. C. \u003cem\u003eet al.\u003c/em\u003e Three-dimensional genome rewiring in loci with human accelerated regions. \u003cem\u003eScience (1979)\u003c/em\u003e \u003cstrong\u003e380\u003c/strong\u003e, (2023).\u003c/li\u003e\n\u003cli\u003eKostka, D., Holloway, A. K. \u0026amp; Pollard, K. S. Developmental Loci Harbor Clusters of Accelerated Regions That Evolved Independently in Ape Lineages. \u003cem\u003eMol Biol Evol\u003c/em\u003e \u003cstrong\u003e35\u003c/strong\u003e, 2034\u0026ndash;2045 (2018).\u003c/li\u003e\n\u003cli\u003eMcLean, C. Y. \u003cem\u003eet al.\u003c/em\u003e Human-specific loss of regulatory DNA and the evolution of human-specific traits. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e471\u003c/strong\u003e, 216\u0026ndash;219 (2011).\u003c/li\u003e\n\u003cli\u003eXue, J. R. \u003cem\u003eet al.\u003c/em\u003e The functional and evolutionary impacts of human-specific deletions in conserved elements. \u003cem\u003eScience (1979)\u003c/em\u003e \u003cstrong\u003e380\u003c/strong\u003e, (2023).\u003c/li\u003e\n\u003cli\u003eYoo, D. \u003cem\u003eet al.\u003c/em\u003e Complete sequencing of ape genomes. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e641\u003c/strong\u003e, 401\u0026ndash;418 (2025).\u003c/li\u003e\n\u003cli\u003eYoung, N. M., Wagner, G. P. \u0026amp; Hallgr\u0026iacute;msson, B. Development and the evolvability of human limbs. \u003cem\u003eProc Natl Acad Sci U S A\u003c/em\u003e \u003cstrong\u003e107\u003c/strong\u003e, 3400\u0026ndash;3405 (2010).\u003c/li\u003e\n\u003cli\u003eCooper, K. L. The case against simplistic genetic explanations of evolution. \u003cem\u003eDevelopment\u003c/em\u003e \u003cstrong\u003e151\u003c/strong\u003e, dev203077 (2024).\u003c/li\u003e\n\u003cli\u003eVarki, A. \u0026amp; Altheide, T. K. Comparing the human and chimpanzee genomes: Searching for needles in a haystack. \u003cem\u003eGenome Res\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 1746\u0026ndash;1758 (2005).\u003c/li\u003e\n\u003cli\u003eWhalen, S. \u0026amp; Pollard, K. S. Enhancer Function and Evolutionary Roles of Human Accelerated Regions. \u003cem\u003eAnnu Rev Genet\u003c/em\u003e \u003cstrong\u003e56\u003c/strong\u003e, 423\u0026ndash;439 (2022).\u003c/li\u003e\n\u003cli\u003eBaumgartner, M., Ji, Y. \u0026amp; Noonan, J. P. Reconstructing human-specific regulatory functions in model systems. \u003cem\u003eCurr Opin Genet Dev\u003c/em\u003e \u003cstrong\u003e89\u003c/strong\u003e, 102259 (2024).\u003c/li\u003e\n\u003cli\u003eAiello, L. \u0026amp; Dean, C. \u003cem\u003eAn Introduction to Human Evolutionary Anatomy\u003c/em\u003e. (Academic Press, London ; San Diego, 1990).\u003c/li\u003e\n\u003cli\u003eLewis, O. J. \u003cem\u003eFunctional Morphology of the Evolving Hand and Foot\u003c/em\u003e. (Clarendon Press, Oxford, 1989).\u003c/li\u003e\n\u003cli\u003eGardner, E., OʼRahilly, R. \u0026amp; Gray, D. J. The Prenatal Development of the Skeleton and Joints of the Human Foot. \u003cem\u003eJ Bone Joint Surg\u003c/em\u003e \u003cstrong\u003e41\u003c/strong\u003e, 847\u0026ndash;876 (1959).\u003c/li\u003e\n\u003cli\u003eGray, D. J., Gardner, E. \u0026amp; O\u0026rsquo;Rahilly, R. The prenatal development of the skeleton and joints of the human hand. \u003cem\u003eAmerican Journal of Anatomy\u003c/em\u003e \u003cstrong\u003e101\u003c/strong\u003e, 169\u0026ndash;223 (1957).\u003c/li\u003e\n\u003cli\u003eBuenrostro, J. D., Wu, B., Chang, H. Y. \u0026amp; Greenleaf, W. J. ATAC‐seq: A Method for Assaying Chromatin Accessibility Genome‐Wide. \u003cem\u003eCurr Protoc Mol Biol\u003c/em\u003e \u003cstrong\u003e109\u003c/strong\u003e, 139\u0026ndash;148 (2015).\u003c/li\u003e\n\u003cli\u003eGaspar, J. M. NGmerge: merging paired-end reads via novel empirically-derived models of sequencing errors. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e19\u003c/strong\u003e, 536 (2018).\u003c/li\u003e\n\u003cli\u003eLangmead, B., Wilks, C., Antonescu, V. \u0026amp; Charles, R. Scaling read aligners to hundreds of threads on general-purpose processors. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e35\u003c/strong\u003e, 421\u0026ndash;432 (2019).\u003c/li\u003e\n\u003cli\u003eLangmead, B. \u0026amp; Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. \u003cem\u003eNat Methods\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 357\u0026ndash;359 (2012).\u003c/li\u003e\n\u003cli\u003eDanecek, P. \u003cem\u003eet al.\u003c/em\u003e Twelve years of SAMtools and BCFtools. \u003cem\u003eGigascience\u003c/em\u003e \u003cstrong\u003e10\u003c/strong\u003e, (2021).\u003c/li\u003e\n\u003cli\u003eZhang, Y. \u003cem\u003eet al.\u003c/em\u003e Model-based Analysis of ChIP-Seq (MACS). \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, R137 (2008).\u003c/li\u003e\n\u003cli\u003eQuinlan, A. R. \u0026amp; Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e26\u003c/strong\u003e, 841\u0026ndash;842 (2010).\u003c/li\u003e\n\u003cli\u003eAndrews, S., Krueger, F., Seconds-Pichon, A., Biggins, F. \u0026amp; Wingett, S. FastQC: A quality control tool for high throughput sequence data. Preprint at (2015).\u003c/li\u003e\n\u003cli\u003eDeng, Z.-L., Unch, P. C. M. \u0026uml;, Mreches, R. \u0026amp; Mchardy, A. C. Rapid and accurate identification of ribosomal RNA sequences via deep learning. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e50\u003c/strong\u003e, 60 (2022).\u003c/li\u003e\n\u003cli\u003eDobin, A. \u003cem\u003eet al.\u003c/em\u003e STAR: ultrafast universal RNA-seq aligner. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e29\u003c/strong\u003e, 15\u0026ndash;21 (2013).\u003c/li\u003e\n\u003cli\u003eLi, B. \u0026amp; Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e12\u003c/strong\u003e, 323 (2011).\u003c/li\u003e\n\u003cli\u003eGuo, M. \u003cem\u003eet al.\u003c/em\u003e Epigenetic profiling of growth plate chondrocytes sheds insight into regulatory genetic variation influencing height. \u003cem\u003eElife\u003c/em\u003e \u003cstrong\u003e6\u003c/strong\u003e, 70\u0026ndash;90 (2017).\u003c/li\u003e\n\u003cli\u003eStaack, A., Donjacour, A. A., Brody, J., Cunha, G. R. \u0026amp; Carroll, P. Mouse urogenital development: a practical approach. \u003cem\u003eDifferentiation\u003c/em\u003e \u003cstrong\u003e71\u003c/strong\u003e, 402\u0026ndash;13 (2003).\u003c/li\u003e\n\u003cli\u003eLove, M. I., Huber, W. \u0026amp; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. \u003cem\u003eGenome Biol\u003c/em\u003e \u003cstrong\u003e15\u003c/strong\u003e, 550 (2014).\u003c/li\u003e\n\u003cli\u003eLangfelder, P. \u0026amp; Horvath, S. WGCNA: an R package for weighted correlation network analysis. \u003cem\u003eBMC Bioinformatics\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 559 (2008).\u003c/li\u003e\n\u003cli\u003eRitchie, M. E. \u003cem\u003eet al.\u003c/em\u003e limma powers differential expression analyses for RNA-sequencing and microarray studies. \u003cem\u003eNucleic Acids Res\u003c/em\u003e \u003cstrong\u003e43\u003c/strong\u003e, e47\u0026ndash;e47 (2015).\u003c/li\u003e\n\u003cli\u003eNeph, S. \u003cem\u003eet al.\u003c/em\u003e BEDOPS: high-performance genomic feature operations. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 1919\u0026ndash;1920 (2012).\u003c/li\u003e\n\u003cli\u003ePrado-Martinez, J. \u003cem\u003eet al.\u003c/em\u003e Great ape genetic diversity and population history. \u003cem\u003eNature\u003c/em\u003e \u003cstrong\u003e499\u003c/strong\u003e, (2013).\u003c/li\u003e\n\u003cli\u003eWu, T. \u003cem\u003eet al.\u003c/em\u003e clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. \u003cem\u003eThe Innovation\u003c/em\u003e \u003cstrong\u003e2\u003c/strong\u003e, 100141 (2021).\u003c/li\u003e\n\u003cli\u003eMcLean, C. Y. \u003cem\u003eet al.\u003c/em\u003e GREAT improves functional interpretation of cis-regulatory regions. \u003cem\u003eNat Biotechnol\u003c/em\u003e \u003cstrong\u003e28\u003c/strong\u003e, 495\u0026ndash;501 (2010).\u003c/li\u003e\n\u003cli\u003eGu, Z. \u0026amp; H\u0026uuml;bschmann, D. rGREAT : an R/bioconductor package for functional enrichment on genomic regions. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e39\u003c/strong\u003e, (2023).\u003c/li\u003e\n\u003cli\u003eMachlab, D. \u003cem\u003eet al.\u003c/em\u003e monaLisa: an R/Bioconductor package for identifying regulatory motifs. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e38\u003c/strong\u003e, 2624\u0026ndash;2625 (2022).\u003c/li\u003e\n\u003cli\u003eFornes, O. \u003cem\u003eet al.\u003c/em\u003e JASPAR 2020: update of the open-access database of transcription factor binding profiles. \u003cem\u003eNucleic Acids Res\u003c/em\u003e (2019) doi:10.1093/nar/gkz1001.\u003c/li\u003e\n\u003cli\u003eLi, H. Tabix: fast retrieval of sequence features from generic TAB-delimited files. \u003cem\u003eBioinformatics\u003c/em\u003e \u003cstrong\u003e27\u003c/strong\u003e, 718\u0026ndash;9 (2011).\u003c/li\u003e\n\u003cli\u003eAmemiya, H. M., Kundaje, A. \u0026amp; Boyle, A. P. The ENCODE Blacklist: Identification of Problematic Regions of the Genome. \u003cem\u003eSci Rep\u003c/em\u003e \u003cstrong\u003e9\u003c/strong\u003e, 9354 (2019).\u003c/li\u003e\n\u003cli\u003ePag\u0026egrave;s, H. BSgenome: Software infrastructure for efficient representation of full genomes and their SNPs. Preprint at https://bioconductor.org/packages/BSgenome (2024).\u003c/li\u003e\n\u003cli\u003eSullivan, P. F. \u003cem\u003eet al.\u003c/em\u003e Leveraging base-pair mammalian constraint to understand genetic variation and human disease. \u003cem\u003eScience (1979)\u003c/em\u003e \u003cstrong\u003e380\u003c/strong\u003e, (2023).\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Table 1","content":"\u003ctable border=\"1\" cellspacing=\"0\" cellpadding=\"0\" class=\"fr-table-selection-hover\"\u003e\n \u003ctbody\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 21.3913%;\"\u003e\n \u003cp\u003eGenomic Feature\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 19.6522%;\"\u003e\n \u003cp\u003eAll brain-filtered\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 20.5217%;\"\u003e\n \u003cp\u003eAutopod-specific\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 21.0435%;\"\u003e\n \u003cp\u003eAutopod-specific, human-mouse conserved\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 17.3913%;\"\u003e\n \u003cp\u003eTissue and timepoint specific\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 21.3913%;\"\u003e\n \u003cp\u003eHARs\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 19.6522%;\"\u003e\n \u003cp\u003e474\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 20.5217%;\"\u003e\n \u003cp\u003e75\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 21.0435%;\"\u003e\n \u003cp\u003e27\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 17.3913%;\"\u003e\n \u003cp\u003e44\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 21.3913%;\"\u003e\n \u003cp\u003ehCONDELs\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 19.6522%;\"\u003e\n \u003cp\u003e886\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 20.5217%;\"\u003e\n \u003cp\u003e157\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 21.0435%;\"\u003e\n \u003cp\u003e59\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 17.3913%;\"\u003e\n \u003cp\u003e93\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 21.3913%;\"\u003e\n \u003cp\u003eHAQERs\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 19.6522%;\"\u003e\n \u003cp\u003e215\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 20.5217%;\"\u003e\n \u003cp\u003e35\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 21.0435%;\"\u003e\n \u003cp\u003e9\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 17.3913%;\"\u003e\n \u003cp\u003e17\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 21.3913%;\"\u003e\n \u003cp\u003eInversions\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 19.6522%;\"\u003e\n \u003cp\u003e2511\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 20.5217%;\"\u003e\n \u003cp\u003e488\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 21.0435%;\"\u003e\n \u003cp\u003e116\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 17.3913%;\"\u003e\n \u003cp\u003e365\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003ctr\u003e\n \u003ctd valign=\"top\" style=\"width: 21.3913%;\"\u003e\n \u003cp\u003eSDRs\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 19.6522%;\"\u003e\n \u003cp\u003e297\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 20.5217%;\"\u003e\n \u003cp\u003e61\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 21.0435%;\"\u003e\n \u003cp\u003e0\u003c/p\u003e\n \u003c/td\u003e\n \u003ctd valign=\"top\" style=\"width: 17.3913%;\"\u003e\n \u003cp\u003e52\u003c/p\u003e\n \u003c/td\u003e\n \u003c/tr\u003e\n \u003c/tbody\u003e\n\u003c/table\u003e\n\u003cp\u003e\u0026nbsp;Table 1 Summary of regulatory element overlaps with genomic features\u003c/p\u003e\n\u003cp\u003eHARs, human accelerated regions, hCONDELs, human conserved element deletions, HAQERs, human ancestor quickly evolved regions, and SDRs, structurally divergent regions between humans and other apes.\u003c/p\u003e\n"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":true,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-7124496/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-7124496/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"The transition to bipedal locomotion is a key event in human evolution, involving substantial changes to the skeleton, including the bones of the hands and feet (autopods). Hominins evolved a more muscular and opposable thumb while the other fingers are relatively shorter, enhancing manipulative capacity. The feet evolved robust first toes and short lateral toes to meet the challenges of bipedal walking and running. While adaptations in the hand and foot have often been considered separately, the fore- and hind limbs of primates are morphologically integrated, serially homologous structures, raising the possibility that natural selection on either autopod may have driven corresponding changes in the other. To explore the genetic architecture underlying human autopod evolution, we used functional genomics methods to identify regulatory elements and gene expression patterns in the developing phalanges and metacarpals of the human hand and foot. We find that gene expression and regulation differ along the proximal-distal axis and between timepoints but not between limb types or individual digits. We show that thousands of human-specific genomic features fall within autopod regulatory elements, some accessible in multiple tissues, others with tissue-specific accessibility. Our results highlight the complex genomic basis of human autopod evolution.","manuscriptTitle":"Complex genetic architecture underlies human hand and foot evolution","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-08-18 10:19:29","doi":"10.21203/rs.3.rs-7124496/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"e0e5ca48-38d7-4777-a90f-de0ac325d39e","owner":[],"postedDate":"August 18th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[{"id":51675486,"name":"Biological sciences/Evolution/Anthropology/Biological anthropology"},{"id":51675487,"name":"Biological sciences/Evolution/Evolutionary developmental biology"},{"id":51675488,"name":"Biological sciences/Genetics/Development"},{"id":51675489,"name":"Biological sciences/Genetics/Functional genomics"}],"tags":[],"updatedAt":"2026-05-13T14:02:20+00:00","versionOfRecord":{"articleIdentity":"rs-7124496","link":"https://doi.org/10.1073/pnas.2603297123","journal":{"identity":"proceedings-of-the-national-academy-of-sciences","isVorOnly":true,"title":"Proceedings of the National Academy of Sciences"},"publishedOn":"2026-05-12 00:00:00","publishedOnDateReadable":"May 12th, 2026"},"versionCreatedAt":"2025-08-18 10:19:29","video":"","vorDoi":"10.1073/pnas.2603297123","vorDoiUrl":"https://doi.org/10.1073/pnas.2603297123","workflowStages":[]},"version":"v1","identity":"rs-7124496","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-7124496","identity":"rs-7124496","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}