The developmental impacts of natural selection on human pelvic morphology.

Evolutionary responses to selection for bipedalism and childbirth have shaped the human pelvis, a structure that differs substantially from that in apes. Morphology related to these factors is present by birth, yet the developmental-genetic mechanisms governing pelvic shape remain largely unknown. Here, we pinpoint and characterize a key gestational window when human-specific pelvic morphology becomes recognizable, as the ilium and the entire pelvis acquire traits essential for human walking and birth. We next use functional genomics to molecularly characterize chondrocytes from different pelvic subelements during this window to reveal their developmental-genetic architectures. We then find notable evidence of ancient selection and genetic constraint on regulatory sequences involved in ilium expansion and growth, findings complemented by our phenotypic analyses showing that variation in iliac traits is reduced in humans compared to African apes. Our datasets provide important resources for musculoskeletal biology and begin to elucidate developmental mechanisms that shape human-specific morphology.

INTRODUCTION

Bipedalism has been viewed as a key adaptation that allowed for crucial behavioral shifts in human evolution, including freeing the arms to use tools, transporting food, carrying offspring, and hunting and gathering (, ). Fossil evidence suggests that the hominin skeleton evolved to accommodate bipedalism relatively quickly after the hominin lineage diverged from its last common ancestor with chimpanzees (, ). Skeletal adaptations for bipedalism were evident in early hominins (, ), and in Homo erectus, much of the hind limb approximated its modern configuration (). Bipedalism likely facilitated other derived hominin skeletal adaptations (e.g., changes in forelimb anatomy facilitated tool making and use) and subsequent increases in cranial capacities and behavioral changes in the genus Homo (–).

The human pelvic girdle is a complex anatomical structure that facilitates both bipedal locomotion and the passage of a relatively large fetal head during parturition (Fig. 1A). To accommodate the biomechanical demands unique to bipedalism, and possibly the effects of increased brain size on parturition () and other factors (, ), natural selection markedly sculpted human pelvic morphology to function differently from chimpanzees and other apes (the pelves of which have also evolved in response to selection on positional behavior) (, ). Evolutionary changes in the human pelvis are present in each of its subelements, i.e., the ischium, pubis, and ilium. For example, compared to chimpanzees and gorillas (African apes), the human ischium has been shifted posteriorly and substantially shortened (), allowing for the repositioning of the hamstrings, which increases efficiency during locomotion (Fig. 1A) (). Likewise, the pubic rami are ventrally located, which supports abdominal organs during upright posture and allows for birth canal (i.e., pelvic inlet) expansion (Fig. 1A) (, ). However, the human pelvic structure most distinct from other apes is the ilium (Fig. 1A), which has been reduced in height and expanded in breadth, curving around the human body in the parasagittal orientation rather than forming the blade-like coronally oriented structure characteristic of other apes (). This expansion created the hallmark basin-like human pelvis and shifted the function of the lesser gluteal muscles to become major hip stabilizers crucial for resisting pelvic tilt during the single leg support phases of bipedal walking and running (, ).

Fig. 1.

Pelvic morphology in humans and African apes.

(A) The morphology of the human pelvis contrasted to the pelvic morphologies of our closest living relatives, chimpanzees and gorillas, as viewed ventrally (top) and laterally (bottom). (B) Comparison of difference in sex standardized variation in pelvic traits among gorillas (Gorilla gorilla), chimpanzees (Pan troglodytes), and modern humans (Homo sapiens). Measurements correspond to those in (A) and colors to regions in (B) to (D) and Fig. 3A. (C) Conditional evolvabilities for individual pelvic traits among gorillas (G. gorilla), chimpanzees (P. troglodytes), and modern humans (H. sapiens). Measurements correspond to those in (A), colors to regions in Fig. 3A, and data correspond to table S1. (D) Comparisons of conditional evolvability of iliac breadth and acetabular height across modern human groups. Acetab., acetabulum; A.H., acetabular height; I.B., ilium breadth; I.H., ilium height; I.L., ischium length; P.L., pubis length.

Fig. 3.

Human pelvic subelement chondrocyte transcriptomics.

(A) Developing CS21 to CS23 (E53 to E59) human pelvic subelements (ilium, pubis, ischium, and acetabulum) chosen for transcriptomic (RNA-seq; N = 7 biological replicates per subelement) assays (note that CS21 is shown). (B) PCA of normalized count data from RNA-seq for all subelements and replicates. (C) Heatmap of 1536 significant DEGs determined by comparing the four subelements against each other. (D) Heatmap of 883 DEGs significantly and uniquely up-regulated or down-regulated in the ilium during development, as determined by comparing the ilium to all other subelements. See fig. S3 (B to D) for related content and table S2 for gene lists and expression values. (E) The top 20 significantly up- or down-regulated DEGs in the ilium, including known pelvic transcription factors EMX2 and HOXD9. (F) Significantly enriched GO terms for DEGs within the chondrogenic pelvis.

Underlying the suite of anatomical changes to the human pelvis are likely changes in the developmental regulation of genes, as found for other adaptive differences between closely related species (, ). This occurs because most gene regions exhibit functional pleiotropy, but regulatory sequences typically control gene expression and affect the phenotype in highly specific anatomical domains and at developmental time points (). For example, regulatory regions unique to human knee and brain development exhibit evidence of substantial sequence divergence and acceleration during human but not African ape evolution (, , ). Compared to the knee, the human pelvis is even more markedly derived, with its morphology evident before birth (, ) and likely under a complex regulatory architecture to accommodate different demands. Our understanding of the developmental genetics of the superior (ilium) and inferior (ischium and pubis) pelvis is, at best, rudimentary in mammals and even less in humans (). Currently, there is little understanding of the transcriptional and epigenetic mechanisms governing human skeletal, let alone pelvic development, and there is correspondingly even less known on the way natural selection has shaped human pelvic architecture, especially for the highly derived ilium.

Here, we first investigate the comparative morphology of the human and ape pelvis, revealing derived features of the human pelvis, including reduced morphological variation in the human ilium compared to other human pelvic structures and those of African apes. We next investigate the prenatal development of the human pelvis, pinpointing a key developmental window when the ilia expand and reorient in the parasagittal plane. We then reveal the transcriptomic and epigenomic architecture underlying this early pelvic developmental stage. In this context, we identify genomic regions, specifically those involved in iliac-specific growth and overall pelvic formation, which were likely targeted by natural selection to produce morphological differences observed between human and nonhuman primate pelves. Overall, we reveal that patterns of human pelvic morphological variation are recapitulated in its underlying regulatory sequence architecture and variation, likely reflecting complex polygenic selection on the pelvis and its specific subelements, especially the ilium and its growth control.

RESULTS

Morphometric insights into human pelvic morphology

One of the most consequential evolutionary changes in human pelvis shape is the reduction in height and breadth of the ilium in humans compared to chimpanzees (Fig. 1A) (). Likewise, the reorientation of the ilium from a coronal (in African apes) to a parasagittal orientation (in humans) has been critical to shifting and augmenting the functions of key lateral and ventral muscles required for human walking. These iliac changes are coupled with a decrease in ischial length compared to other African apes. We quantified pelvic variation using a large comparative morphological sample of humans, chimpanzees, and gorillas (table S1). We confirmed that iliac height and breadth changed more comprehensively in humans relative to other traits of the pubis and ischium, although the latter has become substantially shortened, likely because of selection for efficient hip extension during bipedalism (fig. S1 and table S1) (). Results from a large number of studies (–) [and see review by Love et al. ()] have suggested that natural selection can affect patterns and magnitudes of variation and covariation among traits, and following previous work (, ), we hypothesized that natural selection for parturition and bipedalism could affect patterns of morphological variation in the pelvis when humans are compared to African apes. While humans showed similar levels of variation in ischium length and acetabular diameter to other African apes, they showed less variation in iliac height and breadth, as well as pubis length, while gorillas and chimpanzees did not differ in those traits (Fig. 1B).

On the basis of the common assumption that the phenotypic variance/covariance matrix is proportional to the genetic variance/covariance matrix (see the Supplementary Materials and Methods) (, ), as was done in numerous previous studies (, , –), we made this substitution and quantified the ability of a population to evolve in the direction of selection given stabilizing selection on other traits [i.e., conditional evolvability (); see also the Supplementary Materials]. Our results show that humans have significantly lower levels of conditional evolvability than African apes in iliac breadth, pubis length, and acetabulum height, while ischium length displayed the opposite pattern (Fig. 1C and table S1) (see Discussion). Overall, these results suggest that, compared with other African apes, the human pelvis has less variation available in individual traits upon which selection can act and may be more evolutionarily constrained, possibly reflecting past selection pressures.

We additionally investigated whether natural selection affected among-trait differences in variation by comparing conditional evolvability across human groups. While iliac height and breadth have greater variation overall when compared to other pelvic traits after controlling for sex and group (fig. S1), they are less evolvable than other pelvic traits (e.g., acetabular height) across multiple human groups, and there is substantially less variation among those groups in these iliac traits (Fig. 1D and fig. S1). This could indicate greater constraint in the human ilium in response to selection compared with the other subelements of the pelvis. Together with the results of the interspecific analyses, these findings suggest that natural selection shaped patterns of (i.e., reduced) morphological variation in humans in a way that could reflect the derived iliac morphology associated with bipedalism.

Developmental insights into human pelvic morphology

Given that previous studies have demonstrated that the bony anatomy of the human pelvis arises prenatally (, ) and that we observed distinct human patterns of morphological variation in the adult ilium (e.g., breadth), we next sought to understand when in utero the human pelvis forms and importantly takes on the general appearance of adult structures. Following our previously described protocols (), we first reconstructed the three-dimensional (3D) cartilaginous skeletal morphologies of pelvic subelements. Here, we focused on specimens from Carnegie stage 18 (CS18) to CS23 (N = 2 per stage), when the chondrogenic anlage of most postcranial skeletal elements arise (Fig. 2). While early stages (CS18 to CS19) show rudimentary rod-like ilia and pelvic structures (Fig. 2A) (, , , ), we pinpointed CS20 to CS22 as the critical growth window when the ilium anlage markedly expands in breadth mediolaterally to form a blade shape (Fig. 2A). During this growth phase, the expanding iliac blade (superior portion of the ilium) also reorients from a coronal (African ape-like) plane to a more parasagittal plane (Fig. 2A and fig. S2). This reorientation, which continues through CS23, causes the attachment sites for stabilizing gluteal muscles to face laterally and the presumptive chondrocyte populations of the anterior superior iliac spine (ASIS) and anterior inferior iliac spine (AIIS), key anatomical sites of future muscle and ligament attachments, to now face more ventrally (Fig. 2B and fig. S2). This ilium orientation reflects the adult human morphology and differs from that observed in adult African apes.

Fig. 2.

Human embryonic pelvic development.

(A) 3D image reconstructions of the chondrogenic anlage of the pelvis from CS18 to CS22, showing ventral (left) and lateral (right) views [specimens from () and in the Supplementary Materials: 6524 (CS18), 2114 (CS19), 0426 (CS20), 7254 (CS21), 0895 (CS22), and 9226 (CS23)]. (B) 3D image reconstructions of the chondrogenic anlage of the pelvis at CS23, showing the ventral view of the entire pelvic girdle. (C) Alcian blue (cartilage) and Alizarin red (bone) staining of developing human pelvis, vertebral column, and upper hind limb at E54. (D) Alcian blue (cartilage/chondrocyte) and Alizarin red (bone/primary ossification) staining of developing human pelvis at E54 (left) and E67 (middle) and mouse E15.5 pelvis (right). ac, acetabulum; fe, femur; il, ilium; inf, inferior pelvis; is, ischium; Prim. Ossific., primary ossification center; pu, pubis; svc, sacral vertebral column.

For the pubis and ischium, at CS18, the inferior pelvis is poorly differentiated with an inferior projection possibly representing the early ischium. This suggestion is based on previous work in mouse and chick systems, suggesting that the ischium may be the earlier of the two inferior pelvic structures to form (). From CS19 to CS22, each inferior subelement then elongates inferiorly (caudally) but superiorly merges with the inferior portion of the iliac condensation to form the acetabulum or hip socket proper, while for the pubis, additional growth occurs inferoventrally through CS22 [i.e., embryonic day 55 (E55) to E56] (Fig. 2A). By CS23, the ischium faces posteriorly, while the ventral expansion of the pubis has occurred in such a manner that its superior ramus or branch resides adjacent to the midsagittal plane to meet the superior ramus of the contralateral pelvis (Fig. 2, B and C).

While the early chondrocyte model of the human pelvis prefigures adult morphology by CS23, patterns of primary ossification (i.e., bone formation) of the pelvis, notably ilium, provide additional insights into this important growth phase. At CS21/22, while the breadth of the chondrogenic anlage of the pelvis greatly expands, primary ossification or bone formation has not yet occurred (Fig. 2, C and D). However, at CS23, and through early fetal time points (e.g., E67), a primary ossification center forms within the iliac chondrocyte model, wherein osteoblasts begin to secrete hardened calcium-rich bone. We observed that, at this stage, ossification does not occur uniformly across the human iliac chondrocyte anlage (Fig. 2D). Notably, the iliac region immediately adjacent to the superior acetabular margin remains protected from ossification at CS23 until much later in fetal development when the ossification center extends across the ilium (, , ). This pattern is distinct from the mouse where the ossification zone forms across the entire midportion of the iliac blade by E14.5 and E15.5 (Fig. 2D) (). This human ventral chondrogenic region likely reflects continued iliac cartilage growth and serves as anlage for the presumptive ASIS and AIIS to form. Much later in fetal and postnatal development, secondary ossification centers appear at the AIIS and ASIS, helping to anchor muscle and connective tissue attachments important for bipedalism (, , ). Overall, the described patterns of chondrocyte growth and ossification suggest that natural selection acted on the early human endochondral program for a pelvic configuration that accommodated bipedal locomotion (e.g., the breadth of the ilium).

Transcriptomic insights into the developing human pelvis

To understand the emergence of the adult pelvic pattern at the transcriptomic level, we developed a novel protocol that permits the capture of high-quality RNA from individual microdissected pelvic chondrocyte subpopulations from single human developmental samples (see the Supplementary Materials) (Fig. 3A). We performed RNA sequencing (RNA-seq) on individual subelement chondrocyte populations spanning the critical pelvis, especially iliac, growth window, CS21 to CS23 (i.e., E53 to E59) (N = 7 per subelement) (table S2), which corresponds to mouse E15.5, an early stage of chondrogenesis in which chondrocyte populations are rather homogeneous (). Principal components analysis (PCA) revealed that sequenced libraries generally cluster by subelement type (PC1), with pubic and ischial samples showing some overlapping clustering likely due to their shared tissue origins as inferior pelvic structures (Fig. 3B) (). While we observed an effect of biological sex on gene expression (PC2), sex gene expression patterns have been observed in early human embryos (, ) and do not directly explain sex differences in human pelvic anatomy, as marked dimorphism arises much later during puberty and not early during the embryonic period (). Thus, downstream comparisons of differential gene expression between pelvic subelements were done while controlling for sex (see the Supplementary Materials). We identified 1536 differentially expressed genes (DEGs) (table S2) that exhibited unique expression in specific subelement(s) when compared to the others (Fig. 3C). By clustering significant DEGs across the pelvis, we then identified four clusters of genes that are uniquely up- or down-regulated in specific pelvis subelements, most of which have previously unknown roles in pelvis development (fig. S3A). We identified five additional clusters that reflect shared patterns of gene expression across two or more subelements likely reflecting patterns of complex cartilage growth and integration across the pelvis (see fig. S3A and table S2). When performing gene ontology (GO) analyses of genes differentially expressed across pelvic subelements, we found significant enrichments for skeletal-associated GO terms such as “skeletal system development” [adjusted P value (Padj) = 2.04 × 10−34], “cartilage development” (Padj = 2.07 × 10−19), and “chondrocyte differentiation” (Padj = 1.88 × 10−14) (Fig. 3F and table S2).

Given the marked expansion of the ilium during this developmental window, we next analyzed genes that were up- or down-regulated specifically in the ilium relative to the rest of the pelvis (Fig. 3D). We found that 371 genes were up-regulated and 512 genes were down-regulated in the ilium, with the top 20 shown in Fig. 3E [note that similar analyses were performed for the pubis, ischium, and acetabulum, with results described in fig. S3 (B to D) and table S2]. Up-regulated ilium genes had many associations with skeletal development (table S2). For example, EMX2, DLX5, MMP13, and IFITM5 are all likely to play a role in the growth and expansion of the iliac blade in humans. EMX2 is a hierarchical regulator of ilium development whose loss in mouse results in the absence of the ilium altogether (), DLX5 is associated with cell proliferation (), MMP13 is associated with cartilage development (), and IFITM5 has been associated with abnormal ilium morphology and increased width of the hypertrophic chondrocyte zone ().

We next performed a time course analysis across this key developmental window by comparing samples from CS21 (early E50s) to CS22/23 (late E50s) stages to identify DEGs specifically involved in the expansion of the iliac blade. Despite detecting up-regulated versus down-regulated DEG sets for the pubis, ischium, and acetabulum, we generally did not see enrichments for skeletal development terms for either early or late DEG sets or the combined DEG sets for these subelements (table S2). Conversely, we found that DEGs for the ilium in this analysis (261 genes) were highly enriched for skeletal GO terms observed previously across the pelvis (table S2). By next dividing the ilium set into those up-regulated early (CS21) versus late (CS22/23), 105 genes were highly expressed early in the ilium, but not as highly expressed 3 days later (table S2). These genes were enriched in unexpected pathways involved in neural development, albeit significant enrichments were seen for “cell-cell adhesion,” “ossification,” and “bone mineralization” (table S2). This latter finding suggests that bone formation/ossification processes may be delayed or down-regulated during this window. Conversely, those 133 ilium genes up-regulated later at CS22/23 in the chondrogenic zone of protection showed strong enrichments for GO terms including “cartilage development,” “connective tissue development,” “extracellular matrix organization,” and “chondrocyte differentiation” (table S2). Of these, we identified 36 genes as having been directly associated with skeletal and cartilage-related phenotypes in mice or humans and that additionally show expression in embryonic mouse pelves (table S2). These include genes such as NID2 (), FRZB (), COL2A1 (), and COL27A1 (), among many others. We then identified an additional 50 of 133 genes with general roles in skeletal development and/or skeletal disease (table S2). Therefore, approximately 65% (86 of 133) of genes up-regulated in the ilium during its continued growth and expansion phase have known roles in skeletal development. Overall, our time course analysis of DEGs in the ilium is in line with our skeletal staining analyses, revealing alterations to chondrogenic (i.e., up-regulations) and osteogenic pathways (i.e., down-regulations) that facilitate increased ilium growth and expansion.

Epigenetic insights into the developing human pelvis

The patterns of differential gene expression between the chondrogenic subelements of the human pelvis suggest underlying cis-regulatory differences in the control of gene expression, notably for growth and expansion of the ilium. Moreover, the morphological differences in the pelvis and its subelements between humans and African apes point toward cis-regulatory divergence between species. Using our chondrocyte-optimized ATAC-seq (assay for transposase-accessible chromatin using high-throughput sequencing) protocol (, ), we epigenetically profiled chondrocyte populations from each subelement anlage at human gestational day E54 (N = 3 per subelement; N = 12) (see the Supplementary Materials). Given ethical and practical concerns on collecting African ape embryonic samples, ATAC-seq was also used on chondrocytes from each pelvic subelement from stage-matched mouse E15.5 embryos (N = 3 per subelement; N = 9) (table S3). As mice have a quadrupedal pelvic structure, which is similar but not identical to quadrupedal primates, this approach allowed us to identify functionally conserved versus divergent pelvic regulatory regions in humans and mouse with potential ties to gross morphological differences in pelvic structure. We also performed ATAC-seq on stage-matched human E54 and mouse E15.5 brains (, ), and significantly called peaks (herein, regulatory regions) were then filtered to remove chondrocyte sequences that overlapped with brain and are likely involved in regulating more general cellular “housekeeping” processes (table S3). While reducing the overall regulatory region set sizes, removal of accessible sequences shared with the brain substantially strengthened the specificity of chondrocyte signals for each pelvic subelement set, as revealed using GREAT (genomic regions enrichment of annotations tool) (table S3 and the Supplementary Materials) (). While there is likely a subset of accessible regions shared between the brain and pelvis that is important to pelvic biology, here, we take a conservative approach and only further report on regions unique to the pelvis.

Accessible pelvic regions generated in each species were strongly localized in noncoding distal intergenic, intronic, and promoter regions, emphasizing the importance of differential cis-regulatory usage during skeletal development [e.g., fig. S4 (A and B)]. While most accessible regulatory regions were shared across all pelvic subelements, reflecting their general role in pelvic and chondrocyte biology (Fig. 4A and fig. S4C), a sizable portion of accessible regions was unique to individual pelvic subelements (e.g., in humans: 20.1% in ilium, 18.3% in pubis, 4.6% in ischium, and 5.3% in acetabulum), representing subelement-specific regulatory regions (Fig. 4A, fig. S4C, and table S3). Moreover, intersecting human chondrocyte pelvic subelement-specific sets with aggregated accessibility data from a number of distinct human fetal tissues (see the Supplementary Materials) further revealed their pelvic and pelvic subelement specificity (fig. S4D). Notably, ATAC-seq regions are located in proximity to genes that were uniquely and highly up-regulated in specific pelvic subelements [see fig. S5 (A to C) for examples at EMX2 (ilium) (), EYA1 (pubis) (), and PAPPA2 (ischium) ()].

Fig. 4.

Epigenomic, transcriptomic, and evolutionary sequence features of human pelvic chondrocyte development.

(A) ATAC-seq peak sets for pelvic subelements. Lines connecting dots depict shared peaks between subelements. (B) Partitioned heritability enrichments for region sets using GWAS birth weight, BMI-adjusted waist-hip ratio, and hip circumference data. Asterisks depict significant enrichments (Padj < 0.05). (C) Left; LiftOver analysis for ilium peaks (brain-filtered, top; subelement-filtered, bottom) between human and mouse. Right: Enrichment of overlapping human/mouse brain-filtered and subelement-specific peaks identifying overlaps. (D) Enrichments of subelement-specific and shared peaks with human accelerated regions (HARs). (E) Integrative approach: Gene A, with increased ilium expression versus other subelements, has enhancers within 100 kb. Each enhancer has scores: (1) differential accessibility in a subelement relative to other subelements, (2) pairwise sequence divergence between humans and chimpanzees (blue lines depict fixed mutations relative to ancestral state), and (3) primate sequence conservation (phyloP20ways; blue line) averaged over enhancer length. Each enhancer has a distance-based scaling factor (bottom left). Distance-adjusted scores are summed across enhancers (top right), yielding three values. All expressed genes are ranked by scores, with differential expression values (bottom right); sets of ranks are RRA-aggregated. (F) Gene set enrichment analysis for rank-aggregated genes with increased ilium expression. Enriched terms sorted by gene ratio with dot size and color corresponding to associated gene number and term’s adjusted P value. (G) Genes associated with identified term in ilium analysis had their average cis-regulatory human-chimpanzee similarity values calculated (red line), compared to 1000 sets of randomly sampled genes (black lines). Genes with this term have greater cis-regulatory similarity than expected (Padj = 0.0001). (H) Genes in the indicated term have significantly lower cis-regulatory similarity than expected (Padj = 0.02). See tables S3 and S4.

Given the expectation that regulatory regions active in pelvic development should affect measurable human pelvic phenotypes, we sought to identify any signals of heritability in pelvic morphology, which may be captured by our ATAC-seq datasets. While there are no genome-wide association studies (GWAS) explicitly studying fine-grained pelvic dimensions (e.g., iliac breadth), we instead used other related pelvic size/shape parameters including hip circumference () and body mass index (BMI)–adjusted waist-to-hip ratio () as broad phenotypes. In addition, given the role of pelvic structures in forming the birth canal, neonate features such as size (, ) and birth weight () may reflect aspects of pelvic morphology. Thus, we also considered neonate GWAS datasets from the Early Growth Genetics consortium (see the Supplementary Materials). Using partitioned heritability analysis [see the Supplementary Materials and ()] on our brain-filtered pelvic subelement sets [specific sets being too sparse to provide reliable heritability estimates ()], we found that all sets captured a significant proportion of heritability in hip circumference, while significant enrichments of waist-to-hip ratio signal were limited to the ischium and acetabulum sets (Fig. 4B and table S3). Significant enrichments were not observed for birth length or head circumference [although this may be due to smaller sample sizes in these sets ()], while both the maternal and fetal effects on birth weight were enriched for all sets, with the strongest enrichment observed for acetabular regulatory regions (table S3 and Discussion).

We next examined evidence of between-species functional (i.e., overlapping; shared) orthology and divergence at the chromatin accessibility level (table S3). In general, regulatory regions from the entire pelvis show the most functional orthology between humans and mice, but the number of regulatory region overlaps markedly decrease when comparing more-refined subelement and subelement-specific datasets between species. For example, for the ilium, we found that only 40.04% of human regulatory regions were functionally shared with mouse. The extent of overlap decreased to 11.5% of human regulatory regions when comparing brain-filtered sequences (Fig. 4C, upper left) and further decreased to just 1.1% of human ilium-specific regulatory regions when comparing subelement-specific sequences (Fig. 4C, lower left). Despite these decreases in overlap in more refined filtered sets, the overlapping/shared regions per set intersected more than expected by chance, with more general sets showing stronger enrichments as expected (Fig. 4C, right). Moreover, despite smaller set sizes, these overlapping sets for each tissue still generally retained strong signals of chondrocyte biology (table S3).

We also observed that the more refined regulatory region sets (i.e., those filtered first by brain and then subelement) were predominantly located in distal intergenic regions following with an expectation for greater regulatory divergence for more distal cis-regulatory sequences (fig. S6, A and B) (). This pattern was notable for unique or divergent regulatory region sets between humans and mice, where they were often located far away from the transcriptional start site (TSS) while still remaining relevant to chondrocyte and skeletal biology. For example, both the human-specific and mouse-specific (brain-filtered) ilium regulatory regions were located typically at least 5 kb from the promoter and showed enrichments for collagen biology, chondrocyte biology, Wnt signaling, and skeletal phenotypes in humans and mouse (table S3 and fig. S6B). For more refined species-specific (i.e., gained) ilium-specific regulatory regions, we observed that, while each displayed few regions located within 5 kb of a gene’s TSS and most regions located between 50 and 500 kb away (fig. S6B), only the human-specific ilium regions retained strong signals of association, with chondrocyte biology potentially reflecting unique continued chondrogenic growth and expansion of the iliac blade (table S3). Collectively, these findings indicate that, regardless of the potential regulatory sequence turnover during evolution, marked functional cis-regulatory conservation and functional divergence exist at the level of subelement-specific regulatory biology.

Evolutionary genetic insights into human pelvic regulation

Using comparative genomic methods on human unique regulatory sets and shared human/mouse regulatory sets, we next sought evidence of putatively functional nucleotide changes occurring along the hominin lineage. Both types of regulatory region sets are relevant, as it has been demonstrated that nucleotide substitutions in functionally conserved enhancers can underlie phenotypic differences in limb biology (), as can unique regulatory sequence gains/losses in species (). We first examined human accelerated regions (HARs), genomic sequences that have accumulated nucleotide substitutions along the human lineage at a much faster rate than expected by neutral evolution (see the Supplementary Materials). As HARs are predominantly present in both Neanderthal and Denisovan genomes (i.e., they predate their split from modern humans), their co-occurrence with pelvic regulatory sets likely reflects very ancient natural selection on these sequences (–). Table S4 shows HAR/ATAC-seq regulatory region intersections, identifying that both pelvic subelement shared regions and subelement-specific regions show overlap with specific HAR sequences. We identified 222 HARs intersecting human pelvic subelement regions, 152 of which remain after removing brain-overlapping regulatory regions. Of these 152, only 30 HARs fall within sequences that are functionally orthologous in mouse, indicating that human acceleration has often occurred on regulatory sequences that have functionally diverged because humans and mice shared a last common ancestor. These 152 brain-filtered HARs/pelvis ATAC regions are found throughout the pelvis, with 107 identified in ilium regulatory regions, 111 in pubis regions, 78 in ischium regions, and 81 in acetabulum regions. These can be further partitioned into human subelement-specific regulatory regions that overlap with HARs, of which we identified 10 ilium-specific overlaps, 10 pubis-specific overlaps, 2 acetabulum-specific overlaps, and 1 ischium-specific overlap. The rest overlap either in regions with chromatin accessibility across all four tissues or in some combination of two to three tissues. Figure S7 (A to C) provides examples of loci exhibiting subelement-specific/HAR intersections. We also identified HAR sequences intersecting ATAC-seq regions accessible in multiple subelements that may indicate complex roles in growth, development, and pelvic integration (see table S4).

To understand whether pelvic ATAC-seq sets contain more HARs than expected by chance, we also examined all subelement-specific sets, as well as two unrelated positive and negative control tissue sets [i.e., human brain and B lymphocytes, respectively ()] for whether each significantly overlapped with HARs (table S4). First, previously reported E54 brain ATAC-seq regions were enriched in HARs, while those regions from B lymphocytes were not (Fig. 4D). Second, each subelement-specific regulatory region set (e.g., the ilium-specific set) was specifically enriched for overlap with HAR intersections (Fig. 4D). Third, shared regulatory regions across all four pelvic subelements (i.e., the pelvic common set) also showed enrichment for HAR sequences (Fig. 4D), as did shared ischium/pubis and ischium-specific and pubis-specific regulatory region sets, indicating that the regulatory architecture of inferior pelvic (inlet) chondrocyte biology has been under ancient natural selection (fig. S8). These latter two findings also suggest that there has been highly complex polygenic selection on the developmental regulation of the pelvis at the chondrocyte level.

Integrating transcriptomic, epigenomic, and genetic data to study pelvic evolution

To further refine our understanding of the genes and pathways involved in human pelvic evolution, we integrated data on chromatin accessibility and putative target gene expression, with a more comprehensive assessment of sequence change between humans and African apes (details presented in Fig. 4E, fig. S9, table S4, and the Supplementary Materials). Briefly, for each gene, we consider all nearby ATAC-seq regulatory regions and calculate metrics on the basis of phylogenetic conservation, sequence divergence (between humans and chimpanzees), and differential accessibility of these regions, the latter of which significantly correlates with subelement-specific gene expression (fig. S10). These metrics are used to prioritize genes that are differentially expressed in each pelvic subelement and for which local cis-regulation may be affected by sequence modifications. Focusing here on the ilium subelement (but see table S4 for results on other subelements), GO enrichments on prioritized ilium gene sets revealed key terms such as “skeletal system development,” “ossification,” and “biological adhesion,” among others (Fig. 4F). While these gene-set enrichments were defined using an aggregation of multiple metrics (see the Supplementary Materials), different GO terms had slight biases toward some metrics (e.g., were driven more by changes to enhancer conservation than gene expression). Thus, we pursued terms biased toward regulatory sequence divergence (see the Supplementary Materials) and observed that ilium-enriched genes associated with “embryonic skeletal system morphogenesis,” along with several other developmentally related terms, were actually biased against human-chimp regulatory sequence divergence (Fig. 4G and table S4). This finding may generally speak to evolutionary constraints on regulatory activity around developmental genes or more likely the importance that single base pair modifications have in modifying normally highly conserved transcription factor binding site positions (, ). However, genes associated with “negative regulation of development” did not follow this trend and instead have nearby regulatory sequences with increased and marked human-chimpanzee divergence (Fig. 4H and table S4). This set includes BMP7, which has roles in mouse skeletogenesis () and exhibits strong expression in mouse and human pelvis (, ), as well as FOXA2, which plays an important role in chondrocyte hypertrophy and skeletal development () and is expressed moderately in mouse pelvis () but strongly in human pelvis (fig. S11, A and B). Both loci contain several human-derived sequence changes that may act to alter the regulatory activity of each element.

We next adapted this ranking method to prioritize DEGs detected during ilium growth and expansion (during CS21 to CS23) based on possible cis-regulatory modification (see the Supplementary Materials). We found nine genes for which expression was significantly increased later in ilium development (i.e., in the ventral chondrogenic region) and which had nearby regulatory regions having a greater-than-average rate of human-chimpanzee differences (table S4). These included several genes involved in chondrogenesis and skeletal (including pelvic) development, including the major marker of chondrocytes, COL2A1 (–), as well as COL9A1 (, ), SFRP5 (), and WWP2 (see fig. S11C) ().

Linking patterns of variation in genetics and morphology

Last, given the observed signals of adaptive sequence evolution for the pelvis, and as natural selection should affect within-species sequence diversity (, ), we sought evidence that the regulatory architecture of the pelvis and its subelements have been uniquely affected in modern humans compared to African apes. Analyzing variant data for humans (), we first found that sequence diversity in ilium-specific regulatory regions was markedly constrained compared to human genomic backgrounds, falling significantly to the left of promoter-TSS and intronic regions, while acetabulum-specific regions were instead enriched for variation (Fig. 5A, table S4, and the Supplementary Materials). This latter finding may have implications to hip joint shape and disease. Comparing pelvic subelement sets directly, we next found that only ilium-specific regulatory regions have significantly reduced diversity compared to all other pelvic subelement sets including those shared across the pelvis (Fig. 5C and table S4). These human ilium diversity patterns are unique; our investigation of ilium sequence diversity using chimpanzee () and gorilla () variant data revealed that ilium-specific regulatory regions had no significant reductions in sequence diversity in these species (Fig. 5C). This unique human iliac-specific pattern was also seen when comparing diversity patterns in each subelement set between species (fig. S12B). As further confirmation, we compared chimpanzee sequence diversity within pelvic subelement sets relative to chimpanzee genomic backgrounds and functional annotations (fig. S12A and the Supplementary Materials). We found a markedly different pattern from humans, noting higher sequence diversity levels in chimpanzees compared to humans (); specifically, acetabulum-specific regions were significantly depleted for variation, while ilium-specific regions (as a set) fell near the middle of the distribution (table S4). Overall, these data reveal a unique human pattern of constraint among iliac-specific regions.

Fig. 5.

Genetic diversity levels in human and African ape pelvic regulatory regions.

(A) Counts of common human genetic variants per base pair of sequence for ATAC-seq regulatory region sets compared to random region sets along with other genomic features; labels correspond to table S4. (B) Properties of common human variants within pelvic ATAC-seq regions: (Top) ARGweaver () estimated allele age of variants falling within the indicated subelement region set. (Bottom) LINSIGHT () predicted variant effect scores of variants within sets. See table S4 for statistics. (C) Common variants in human, chimpanzee, and gorilla intersecting a given ATAC-seq region were counted for all subelement-specific sets, expressed as “SNPs per sequence.” Mean and median values are shown in dashed and bold lines, respectively. In chimpanzees and gorillas, all pelvic subelements have comparable numbers of SNPs. TTS, transcription terminal site. Ilium-ischium, *P = 0.0039; ilium-pubis, **P = 0.00048; ilium-acetabulum, **P = 0.000099. See table S4, sheets titled Intraspecies, Human; Intra-species, Chimp; and Intraspecies, Gorilla for additional statistics on all non-significant (ns) findings.

These sequence findings complement our morphological results that show decreased morphological variation in humans relative to chimpanzees and gorillas and lower conditional evolvabilities for the human ilium across 17 human groups (Fig. 1D and fig. S1). These are patterns that were not observed in other pelvic subelements. As an additional means of assessing sequence constraint in pelvic subelement sets within humans, we considered the estimated allele age of variants falling within regulatory regions (see the Supplementary Materials) (). Comparing across pelvic sets, we observed that variants within ilium-specific regions tended to be younger than those in other subelement sets and variants not falling within ATAC-seq regions (Fig. 5B, top, and table S4). An increased prevalence of younger mutations is suggestive of stronger purifying selection in a locus, as older mutations are more efficiently eliminated (, ). This ilium-specific signal of constraint is likely not driven by a general increase in functional effects of variants falling within these regions, as when we compared predicted per–single nucleotide polymorphism (SNP) effects () across sets, we observed no significant differences (Fig. 5B, bottom, and table S4). These two findings, along with our observations at the inter- and intraspecies level, suggest that recent purifying selection on the backdrop of ancient (positive) selection has had a profound and continued impact in shaping the human ilium, possibly reflecting selection on an important iliac growth phase during prenatal life.

DISCUSSION

In this study, we found unique human patterns of morphological variation in the ilium, and pelvis in general, which are also reflected at the sequence level for the genome-wide set of regulatory regions involved in pelvic formation. These morphological and genetic patterns are not present in chimpanzees and gorillas, pointing to evolutionary changes in the genetic architecture of hominin and human pelvis development that possibly arose as human ancestors acquired and sustained habitual bipedal locomotion. In this context, we observed ilium growth, expansion, and reorientation and marked evidence of the effects of natural selection (i.e., HAR enrichments and depletions in genetic diversity in site-specific regulatory regions) on ilium and other pelvic regulatory sequences involved in a critical window of human pelvic chondrogenesis.

On the basis of fully ossified late fetal and neonatal samples, human ilium growth and expansion has been posited to occur via an anteriorly located growth zone [which remains undetected or underdeveloped in other primates, most notably chimpanzees and gorilla ()]. This zone has been hypothesized to constitute a unique hominin synapomorphy (). We pinpointed this chondrocyte zone in developing human pelves in utero and found that it remains protected from ossification, permitting enhanced growth of the ventral portion of the iliac blade and contributing to the early anlage of the AIIS and ASIS. Growth and expansion of this portion of the blade permits the attachment of key muscles (both agonistic and stabilizing) and soft tissues of human bipedalism, including the inguinal ligament, sartorius muscle, and tensor fasciae latae muscle to the ASIS, and the rectus femoris muscle and iliofemoral ligament to the AIIS (). We also found that, during this growth window, the superior portion of the iliac blade reorients along the parasagittal plane (see below), which causes the ASIS and AIIS to reside ventrally. While additional and difficult-to-execute functional studies are required to understand the root cause of this reorientation, this shift is a quintessential feature that underscores human bipedalism because it repositions key gluteal muscles (gluteus minimus and medius) to the lateral side of the pelvis to stabilize the pelvis during upright walking. This suite of pelvic morphological features is unique to hominins (), and some of these features, such as a prominent AIIS, were present 4.4 million years ago in the earliest bipedal hominin, Ardipithecus ramidus (, ).

This critical growth window is the subject of a complex polygenic architecture; thousands of genes and regulatory regions, many specific to ilium, likely contribute to ilium expansion and reorientation. For example, we identified a large subset of genes uniquely up-regulated in the chondrocyte anlage of the ilium relative to the rest of the pelvis, and iliac genes whose expression change over the iliac growth, expansion, and reorientation window (i.e., from CS21 to CS23). Those genes that are down-regulated during this developmental window correspond to functions involved in bone ossification, whereas those that are up-regulated instead are involved in chondrogenesis and chondrocyte biology and have a myriad of known roles in chondrogenesis and pelvic formation (e.g., at least 65% of the 133 up-regulated iliac growth genes have known skeletal roles). This is consistent with our skeletal findings that the anterior/ventral ilium is protected from ossification for prolonged chondrogenesis. Likewise, the myriad regulatory regions specific to the ilium also revealed the importance of pattern specification, chondrocyte differentiation, and metabolic pathways in underlying chondrogenic growth and expansion of the iliac blade. Overall, these findings specifically point to a series of chondrocyte regulatory mechanisms underlying human-specific ilium development, in line with prior studies noting the role of chondrogenesis in the development of anatomical specializations (, , , , ). Moreover, our novel approach has also now substantially added to what was previously a dearth of knowledge on the molecular control of mammalian pelvicogenesis () and provides important suites of transcription factors, signaling molecules, and biological pathways that influence general pelvic as well as pelvic subelement-specific anatomical development.

We previously reported for the human knee (another skeletal structure critical for bipedalism) that regulatory sequences shared between chondrocyte populations of the developing joint (i.e., shared between the distal femur and proximal tibia) did not significantly overlap HARs, whereas only sequences specific to each subelement (e.g., regulatory regions specific to the distal femur) showed such enrichment. In the case of the knee, ancient selection likely targeted anatomically specific regulatory sequences shaping growth and morphogenesis above versus below the knee (). In the context of polygenicity underlying human pelvic morphology, we found that not only are chondrocyte subelement-specific regulatory sequences (e.g., ilium-specific regions) enriched for HARs, but so too are regulatory regions accessible across all pelvic elements. Moreover, we observed strong HAR enrichments in sequences shared between inferior pelvic subelements (pubis and ischium) likely reflecting the role of development and evolution on shaping our unique birth canal and roles in bipedalism. These findings indicate that chondrocyte and cartilage biology has been a target of selection in the ancient human past, and this is supported by earlier findings on human cartilage biology ().

By integrating ATAC-seq (i.e., differential accessibility), RNA-seq (i.e., differential gene expression), and human and primate genetic sequence data, we also revealed finer-grain patterns of regulatory sequence divergence modifying pelvic subelement-specific expression programs. Specifically, by stratifying genes by ilium-specific expression and both chromatin accessibility and between-species sequence divergence of their nearby regulatory regions, our results further support the substantial polygenicity of ilium development and the functional modularity of regions regulating these genes, both of which may have been taken advantage of by ancient natural selection acting to shape the human ilium. That we observed an enrichment for the term “negative regulation of development” for the set of iliac up-regulated genes whose nearby regulatory regions exhibit marked human/chimpanzee sequence divergence supports the findings of delayed ossification/increased chondrogenesis during iliac growth. This is in line with the heterochronic alterations to human life history strategy in which humans have delayed growth and development in utero compared to apes (, ). The molecular targets of this evolutionary process are possibly modifications (via HARs and single base pair modifications) to conserved chondrocyte regulatory sequences, such as those near BMP7 (, ), COL2A1 (–), COL9A1 (, ), HDAC2 (), RUNX2 (, ), SFRP5 (), WWP2 (), and others, each of which may mediate shifts in gene expression during early human prenatal ilium patterning and growth. The cumulative effect of these shifts acting on developmental chondrogenesis pathways may underlie changes in ilium morphology.

Last, we found that ilium morphology exhibits reduced variation in modern humans when compared to other African apes and, within humans, reduced conditional evolvability compared to the acetabulum, pubis, or ischium. Note that between taxa, we observe similar levels of conditional evolvability in iliac height, potentially reflecting constraints related to body size and locomotor loading regimes in African apes (). However, iliac breadth in humans has lower conditional evolvability, implying that the evolutionary potential in the ilium has not been apportioned equally among developmental axes. Reasons for this will require further study, although we note that the delayed ossification and, thus, extended chondral growth on the anterior ilium in our human gestational samples notably coincides with the anteroposterior (breadth) axis. Across human populations, the lower levels of evolutionary constraint in the inferior pelvic subelements compared with the ilium (see Fig. 1D and fig. S1C) also await further exploration, although we provisionally hypothesize that the human ischium may have greater potential autonomous response to directional selection given its importance in both parturition () and bipedalism (, ).

Notably, we also observed the effects of such constraint on the genome-wide set of ilium-specific regulatory regions, especially compared to pubis-, ischium-, and acetabulum-specific sets; the set of regions shared across the pelvis; and several functional genomic annotations (e.g., promoter-TSS sites). The variants falling within ilium-specific regions also tended to be more recently manifest than those within other subelement regions. This trend does not appear to be explained by increased predicted functional effects of these variants but, instead, may suggest that these regions have been specifically placed under stronger purifying selection following ancient natural selection. Overall, our findings for the ilium on the phenotypic level are matched on the genome-wide (genotypic) level and suggest how selection for uniquely human pelvic characteristics has shaped pelvic regulatory variation at a key stage of ilium expansion. That we saw no overt comparative human signal of reduced genetic diversity within ischium and pubis regulatory regions in part mirrors a lack of a marked effect on morphological measures for these subelements and may suggest that modern drift processes or relaxed selection has been occurring.

We note that our study has a number of inherent weaknesses, including our inability to functionally test sequences using in vivo assays. For example, many of our findings directly on human iliac chondrocytes would not be easily functionally testable in the mouse model, considering that the mouse lacks the same anterior ilium structures found in the human ilium. We also note that, given the difficulty of acquiring rare human samples along with the large number of cells typically required for chromatin capture techniques such as high-throughput chromosome conformation capture (HiC), we unfortunately could not directly ascertain regulatory region (enhancer)–promoter interactions. Last, given issues with accessing earlier staged human samples, we caution that we cannot assess the potential yet limited roles that earlier patterning events may have on pelvic morphology.

In summary, our integrated morphological, transcriptional, and epigenetic analysis on the pelvic girdle has yielded substantial insight into a feature that has defined much of the evolution of human walking and running as well as childbirth. We present these rare human developmental and morphological datasets so that future studies can investigate human-specific pelvic development, the development of pelvic sexual dimorphism, and the genetic etiology of congenital pelvic deformities, such as hip dysplasia, and adult hip joint diseases, such as hip osteoarthritis.

MATERIALS AND METHODS

Morphological samples

For the interspecific analyses, 3D landmark data from the hip bone of adult Pan troglodytes, Gorilla, and Homo sapiens were collected using a MicroScribe digitizer (Immersion, San Jose) based on homologous bony landmarks from collections at the National Museum of Natural History (Washington, DC), the American Museum of Natural History (New York, NY), the Museum of Comparative Zoology (Cambridge, MA), the Cleveland Museum of Natural History (Cleveland, OH), the Anthropological Institute and Museum at the University of Zurich (Zurich, Switzerland), the Museum fur Naturkunde (Berlin, Germany), the Royal Museum for Central Africa (Tervuren, Belgium), and the Royal Belgian Institute for Natural Sciences (Brussels, Belgium). All individuals were adults, and sexes of most were known; in the few instances in which the sex was unknown, it was estimated using a discriminant function analysis. The right hip bone was used where possible. From these 3D coordinates, interlandmark distances were calculated [see also Grabowski et al. () and Grabowski and Roseman ()] and used here for the traits acetabulum height and pubis length. These interlandmark distances were complemented by distances taken using digital calipers (Mitutoyo Series 500-, 19X-, 20, Plymouth, MI), here the traits iliac height, iliac breadth, and ischium length (see table S1 for sample details).

For the analyses of recent human variation and evolvability statistics, 3D landmark data were collected from human os coxae using a G2X Immersion MicroScribe digitizer (Immersion, San Jose). Prior studies have demonstrated that these landmark data are repeatable, precise, and accurate (, ). Measurements in all cases were taken from the best-preserved side. All data used herein were obtained by L. Betti from collections at the National Museum of Natural History (Washington, DC); the American Museum of Natural History (New York, NY); the Cleveland Museum of Natural History (Cleveland, OH); Chiang Mai University (Chang Mai, Thailand); Coimbra University (Coimbra, Portugal); the Duckworth Collection, University of Cambridge (Cambridge, UK); Kyoto University (Kyoto, Japan); Musée Canadien des Civilisations (Gatineau, Canada); Musée de l’Homme (Paris, France); McGregor Museum (Kimberley, South Africa); Museo Nazionale di Antropologia e Etnologia (Florence, Italy); the San Diego Museum of Us (San Diego, CA); the Natural History Museum (London, UK); Naturhistorisches Museum (Vienna, Austria); the National Museum of Kenya (Nairobi, Kenya); Tokyo University (Tokyo, Japan); Phoebe Hearst Museum of Anthropology (Berkeley, CA); the University of Pennsylvania (Philadelphia, PA); the University of Rome “La Sapienza” (Rome, Italy); the University of Tennessee (Knoxville, TN); and the University of Witwatersrand (Johannesburg, South Africa). In addition, Arikara data were obtained by B. M. Auerbach and A. M. Mallard from the collections housed at the University of Tennessee, Knoxville, using the same methods as L. Betti; the morphological data obtained from the Arikara by L. Betti and by B. M. Auerbach and A. M. Mallard were compared and found to yield identical results. Using the same interlandmark calculation method as noted above, we obtained linear dimensions for five traits: iliac height, iliac breadth, acetabulum height, pubis length, and ischial length. All individuals selected were determined to be adult (based on complete fusion of all pelvic secondary ossification centers), and sex was estimated using visual sexing of the pelvis (, ). To improve the reliability of the calculation of variance and evolvability estimates, human groups used in this study were restricted to those with sample sizes greater than 20 of each sex.

Morphological analysis

Here, we substitute the phenotypic variance/covariance matrix (e.g., P) for the genetic variance/covariance matrix (e.g., G) based on the assumption that the two are proportional, a conjecture originally made by Cheverud (1988) () and later supported in a wide range of studies, particularly for morphological traits (), and which has been done in numerous previous studies (, , , , , ).

We calculated sex-corrected mean standardized covariance matrices (the P matrix) for each species using the residuals from a multivariate analysis of variance (MANOVA) with the five traits as the dependent variable and sex as the independent variable following previous works (, , , ), and the effects of group for modern humans using the same approach to ensure that our estimate was free of variation owing to geographic patterning (). We repeated our analyses, leaving in these effects, and found no change in our results. All covariance matrices were mean standardized () to allow for interspecific comparisons. Trait-level conditional evolvability is derived using the formula , where x and y are the individual traits and P is the phenotypic variance/covariance matrix (). This equation calculates the variance in each trait that remains after conditioning on all other traits in the covariance matrix.

To compare differences in variation across modern human groups for pelvic traits, mean sex-corrected mean standardized covariance matrices were calculated using the residuals from a MANOVA with the five traits as dependent variables and sex as the independent variable. Trait-level conditional evolvability was obtained from the mean-scaled P matrices using the same methods as described above.

All SEs and tests for significant differences were calculated using a parametric bootstrap technique (). All morphological analyses were conducted in R ().

Human developmental sample collection

The human products of conception at gestational days E53 to E59 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. The BDRL was supported by NIH award number 5R24HD000836 from the Eunice Kennedy Shriver National Institute of Child Health and Human Development. Harvard University IRB determined that these samples constitute Non-Human Subjects Determination Status (Capellini: IRB16-1504). The corresponding author (T.D.C.) received no federal funds (e.g., NIH or National Science Foundation) to acquire, receive, process, or use samples. The fresh human samples were briefly washed in Hanks’ balanced salt solution and transported at 4°C during shipment. 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. Samples used for skeletal preparation were treated differently as they were flash-frozen in liquid nitrogen at the BDRL, shipped on dry ice, stored at −80°C, and then microdissected and processed as described below.

Skeletal preparation

Skeletal preparation was performed on gestational samples that had been flash-frozen and stored at −80°C. Samples were initially thawed in a 37°C water bath before dissection in 1× phosphate-buffered saline (PBS). For skeletal preparation dissection, major muscle bellies and large pieces of soft tissue were removed, but joints were not disarticulated, and the pelvis was left intact rather than separated into individual subelements, as was done for RNA-seq and ATAC-seq dissections. The samples were then placed in 95% EtOH for 5 days and shielded from light. Samples were then placed in acetone for 2 to 3 days to remove fat and soft tissue. Samples were then stained for bone (Alizarin red) and cartilage (Alcian blue) in a staining solution containing 5 ml of 0.3% Alcian blue in 70% EtOH, 5 ml of 0.1% Alizarin red in 95% EtOH, 5 ml of glacial acetic acid, and 85 ml of 70% EtOH. Samples remained in staining solution and were shielded from light for 3 days while staining progress was monitored daily. Following the completion of staining, samples were washed in water and then cleared in 1% KOH for 2 days. They then went through a glycerol series, with samples placed in 25% glycerol for 1 to 2 days, then 50% glycerol for 1 to 2 days, and ultimately stored in 100% glycerol. Samples were imaged on a light microscope in 100% glycerol.

Human cartilage anlagen reconstruction

Embryonic samples from the 3D Atlas of Human Embryology were obtained and processed following protocols in de Bakker et al. (). Digital reconstructions of the pelvis and its subregions of each specimen between CS18 (E44 to E48) and CS23 (E56 to E60) were obtained as 3D surfaces in pdf format and were used for analysis and image presentation. At each stage shown, a minimum of N = 2 biological replicates were examined per stage.

Human RNA-seq data collection and analysis

Fresh pelvic samples were microdissected under a light microscope in 1× PBS on ice, and the component parts or subelements of the pelvis (ilium, ischium, pubis, and acetabulum) were separated, stripped clear of soft tissue, and collected in 2-ml tubes containing 200 μl of TRIzol and one 5-mm stainless steel bead. The right and left sides of each pelvic subelement were pooled from each donor sample. Each sample was then homogenized at 50 Hz for 2 min, followed by 1 min on ice, and a final 50-Hz homogenization for 2 min. Samples were then placed at −80°C until RNA extraction was performed.

For each RNA extraction, samples were first subjected to a phenol-chloroform reaction in which samples were moved to a new microcentrifuge Eppendorf tube with 200 μl of chloroform per 1 ml of TRIzol in the sample. They were then vigorously vortexed and incubated at room temperature for 2 min before centrifugation for 5 min at 12,000g at 4°C. Following centrifugation, the aqueous phase was removed and transferred to a new microcentrifuge Eppendorf tube. After completion of this phenol-chloroform step, the RNA extraction continued using the Zymo Direct-zol RNA MicroPrep kit following the manufacturer’s protocols and eluted in 15 μl of nuclease-free water. Samples were then quantified using a Qubit and following the manufacturer’s protocols. Samples were also nanodropped to determine 260/230 and 260/280 values and were then stored at −80°C. Aliquots of samples were also run on a TapeStation to determine their RNA integrity number (RIN) value. Only samples with RIN scores of 7 or higher were used for subsequent steps.

For cDNA library generation and sequencing, samples were normalized to a single concentration, and libraries were prepared using KAPA mRNA Directional Library Preparation methods following the manufacturer’s protocols. Twenty libraries were then run as quality control on a TapeStation to verify their quality. Quantitative polymerase chain reaction (PCR) was then performed on the library pool before sequencing of the library on NextSeq High 2 × 38. The library was sequenced repeatedly using paired-end sequencing on three lanes to obtain a minimum of 10 million reads for each sample, with some samples requiring only one or two lanes. On average, samples were sequenced at 23 million reads per sample (when including one extreme outlier at 80 million reads) or 20.7 million reads per sample (when removing this outlier) within the recommended Encyclopedia of DNA Elements (ENCODE) guidelines (www.encodeproject.org/about/experiment-guidelines/). See table S2 for sequencing sample, index primer, read count, and other information. Seven biological replicates of each pelvic subelement were sequenced, but one acetabulum sample failed in sequencing because of low concentration, and one ilium sample was removed in downstream processing because of highly outlying values. Therefore, N = 6 for acetabulum and ilium.

Computational analysis of RNA-seq data began with running Fast QC () version 11.5 on each Fastq file to determine the per-base sequence quality, GC content, etc. to ensure that all files met our standards. Because the samples had each been run on three 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. SortMeRNA version 2.1b () was then used to blacklist ribosomal RNA databases and filter samples for ribosomal RNA, which constituted 10 to 20% of most samples. STAR version 2.6.0 () was then used to map reads to the human genome, hg19, with about 90% of reads uniquely mapping for each sample. RSEM version 1.2.29 () was then used to generate read counts. For all samples, >90% of reads were mapped uniquely to a gene. DESeq2 version 1.26.0 () was then used to quantify differential expression and downstream comparisons between samples.

PCA of sample-level regularized and log-transformed normalized count data revealed that sequenced libraries were clustered by sex of the sample along PC2; we then incorporated this source of variation into our DESeq2 model and controlled for sex when assessing DEGs by pelvic subelement. We also performed two types of analyses: whole pelvis analyses in which the ilium, ischium, pubis, and acetabulum samples were compared against each other to identify genes differentially expressed across the entire structure and subelement-specific analyses in which we sought to identify genes that are up-regulated or down-regulated in one specific pelvic subelement relative to others (i.e., ilium versus everything non-ilium), or between time points (e.g., early ilium at E53 to late ilium at E57) (table S2).

GO analyses were performed using clusterProfiler (). In addition, genes demonstrating significant expression changes relative to other subelements (e.g., ilium-specific DEGs) were aggregated (table S2), with the final unique gene set used with the MGI Batch Query tool () to find associated Mammalian Phenotype (MP) and Human Disease Ontology terms for these genes (table S2).

Human ATAC-seq data collection and analysis

Fresh pelvic samples were microdissected under a microscope in 1× PBS on ice, and the subelements of the pelvis were separated, stripped clear of soft tissue, and collected in microcentrifuge Eppendorf tubes containing 200 μl of 5% fetal bovine serum/Dulbecco’s modified Eagle’s medium (FBS/DMEM). The right and left sides of each pelvic subelement were pooled from each donor sample. Each sample was then subjected to 1% collagenase II (VWR, 80056–222, Radnor, Pennsylvania) digestion for 2 hours at 37°C rocking, mixing every 30 min. After placing on ice, samples were next filtered using a microcentrifuge filter setup by gently mashing the residual tissues through the filter followed by rinsing with 5% FBS/DMEM. Samples were then spun down at 500g at 4°C for 5 min. All cell counting methods were performed using trypan blue and a hemocytometer under an inverted microscope, and subsequent ATAC-seq steps were performed on those samples that had cell death rates well below 10%. Replicates consistently yielded several hundred thousand cells. Next, cells were resuspended in concentrations of 50,000 cells in 1× PBS. Cell samples were then subjected to the ATAC-seq protocol as described previously (, ), modifying the protocol by using 2 μl of transposase per reaction [see Richard et al. () and Guo et al. ()]. The transposase reaction product was then purified using the Omega MicroElute DNA Clean Up Kit following the manufacturer’s protocols, eluted in 10 μl of warmed double-distilled water (ddH2O), and stored at 20°C.

All samples were next subjected to PCR amplification and barcoding following Buenrostro et al. (2013, 2015) (, ). Ten microliters of transposed DNA was placed in a reaction containing the NEBNext High-Fidelity PCR Master Mix, ddH2O, and barcoding primers. Following amplification, samples were transferred to new tubes and treated for fragment size selection using the OMEGA Bead Purification Protocol following the manufacturer’s instructions. The samples were eluted in 15 μl of tris EDTA (TE), nanodropped, diluted to 5 ng/μl, and run on a Bioanalyzer for quality control and visualization of fragment size distribution. Before sequencing, sample concentrations were determined using the KAPA Library Quantification Complete Kit (KK4824). Samples were then sent out to the Harvard University Bauer Core Facility for paired-end sequencing on one lane of the Illumina NextSeq 500. Each biological replicate was sequenced on average to 83 million reads per sample. Please see table S3 for sequencing sample, index primer, read count, and other information.

Computational analysis of ATAC-seq data began with checking read quality via FastQC. Adapters were then removed using NGmerge () version 0.2. Reads were subsequently aligned to human reference hg19 genome assembly with Bowtie2 v2.3.2 () using default parameters for paired-end alignment. Duplicate reads were removed using Picard’s MarkDuplicates function (version 2.9.0) (http://broadinstitute.github.io/picard), and mitochondrial reads were also filtered out. The resulting .bam files were then used for peak calling via MACS2 software (version 2.1.1.2) (), using BAMPE and the following flags: --nolambda –bdg --verbose. Peaks reproducible across biological replicates were screened using a conservative irreproducible discovery rate (IDR) threshold of <0.01, as defined by the IDR statistical test (version 2.0.3) (). The IDR method identifies overlaps in peak calls across pairs of sample replicates by comparing ranked peak lists (using MACS2 q value) to define a reproducibility score curve. These paired peaks are then assigned a pointwise score based on this curve. They are then sorted, and all peaks that fall below an IDR threshold (here defined as 0.01) are taken as our final reproducible peak set across replicates. This yielded peak calls, including those called using normal stringency IDR thresholds (0.05) within ENCODE guidelines of greater than 50,000 peaks (). Following the identification of the IDR set of peaks for each pelvic subelement, the peaks were then reduced by filtering for brain. ATAC-seq was also performed on stage-matched E54 brain samples, and called peaks were then filtered to remove chondrocyte sequences that overlapped with brain sequences and are likely involved in regulating more general cellular processes. Intersections between brain peaks and pelvic peaks were done using bedtools version 2.27.1 (), which was also used to intersect pelvic sets to identify which peaks were shared across pelvic subelements and which were specific to only one subelement. GREAT version 4.0.4 () was used to identify functional enrichments for peak sets.

As an additional means to assess tissue specificity of pelvic subelement ATAC-seq peak sets, we collected ENCODE () deoxyribonuclease I hypersensitivity datasets for eight different fetal tissues (adrenal gland, brain, heart, lung, muscle, skin, spleen, and stomach) retrieving sorted, duplicate-filtered mapped read files (.bam) via the ENCODE web portal (). To define reproducible hypersensitivity sites within each tissue, we applied the IDR statistical test (version 2.0.3) as described above. For each sample, peaks were called with MACS2 (version 2.1.1.2) using the following parameters: “-f BAMPE –nolambda” and “-f BAM –no-model –shift -100 –extsize 200” for paired-end and single-end experiments, respectively. An IDR threshold of 0.05 was applied, with resulting filtered peak sets combined using “bedtools merge” in those instances where both single-end and paired-end experiments for a given tissue were obtained and processed separately with MACS2/IDR. The resulting group of eight peak sets were then intersected with a given ATAC-seq peak set (e.g., the ilium-specific peaks), with the number of overlaps of different fetal tissues with each individual ATAC-seq peak counted per peak. Counts of the number of ATAC-seq peaks in a given set that intersected with one or more fetal tissue were visualized using barplots in base R. The following were the ENCODE file accessions used: ENCFF542NTO, ENCFF846NFV, ENCFF848SVJ, ENCFF542AQF, ENCFF602FMG, ENCFF537DKA, ENCFF948DOQ, ENCFF108NTA, ENCFF431RLT, ENCFF129ZBZ, ENCFF952QLZ, ENCFF436NOR, ENCFF808BII, ENCFF807SCV, ENCFF196CDZ, ENCFF019NIR, ENCFF376OPD, ENCFF562IIU, ENCFF544SGF, ENCFF152KTS, ENCFF900LLD, ENCFF124QRT, ENCFF352BKQ, ENCFF813XJE, ENCFF774DEI, ENCFF431JSL, ENCFF594WQM, ENCFF706ECO, ENCFF280OVS, ENCFF523OSH, ENCFF934CDN, ENCFF987PLB, ENCFF496KEC, ENCFF870CEG, ENCFF566BGJ, ENCFF519UYK, ENCFF968VWA, ENCFF349QPC, ENCFF667MRJ, ENCFF343CEI, ENCFF852QUY, ENCFF746WUF, ENCFF008VCJ, ENCFF890JQE, ENCFF349QPC, ENCFF667MRJ, ENCFF343CEI, ENCFF852QUY, ENCFF746WUF, ENCFF008VCJ, ENCFF890JQE, ENCFF767PZY, ENCFF113AJV, ENCFF333ZXH, ENCFF221DVC, ENCFF324KEB, ENCFF159KZC, ENCFF746QTI, ENCFF969UDM, ENCFF615FRI, ENCFF584WDA, ENCFF890KBJ, ENCFF949LCP, ENCFF717FJE, ENCFF661IGM, ENCFF085XYO, ENCFF374BUY, ENCFF682DZN, ENCFF215BUS, ENCFF121ZBR, ENCFF285BSG, ENCFF961FZQ, ENCFF698LQU, ENCFF654DRD, ENCFF859ERU, ENCFF663KHG, ENCFF630SFP, ENCFF565CIZ, ENCFF232UOO, ENCFF935OOV, ENCFF119VPN, ENCFF337MIR, ENCFF611AGL, ENCFF054IPB, ENCFF308HJB, ENCFF817QEH, ENCFF635UOX, ENCFF186BCI, ENCFF020GHN, ENCFF658ZEV, ENCFF672LIO, ENCFF951KSP, ENCFF833XQA, ENCFF151PSG, ENCFF032KCG, ENCFF959RWP, ENCFF678XNA, ENCFF429OVI, ENCFF663QIV, ENCFF245HUN, ENCFF799LOQ, ENCFF922HOC, ENCFF141FQR, ENCFF730UQK, ENCFF156XGI, ENCFF269RUB, ENCFF248VBW, ENCFF775QZL, ENCFF412NAG, ENCFF065ECG, ENCFF468GQN, ENCFF863QYT, ENCFF339UWC, ENCFF426BXR, ENCFF108JUG, ENCFF989EMJ, ENCFF032YRW, ENCFF783EIR, ENCFF845MVQ, ENCFF794BQM, ENCFF904SUV, ENCFF014ZPA, ENCFF127ZSV, ENCFF559JLL, ENCFF574VDU, ENCFF659MZH, ENCFF936FGG, ENCFF772QLS, ENCFF568IND, ENCFF110IZT, ENCFF229SBH, ENCFF560XJK, ENCFF216WEZ, ENCFF693BFS, ENCFF151JHW, ENCFF663SJQ, ENCFF467OYK, ENCFF376AJT, ENCFF141NAT, ENCFF077JZB, ENCFF213QEX, ENCFF039NMI, ENCFF999WJD, ENCFF561KHG, ENCFF679ZBQ, ENCFF591NPU, ENCFF146ADR, ENCFF085BXU, ENCFF371NCW, ENCFF636FHS, ENCFF792JOV, ENCFF264LLE, ENCFF148JEQ, ENCFF065WZW, ENCFF067LVL, ENCFF971APS, ENCFF056UYZ, and ENCFF127RJK.

Mouse ATAC-seq data collection and analysis

Fresh mouse pelvic samples at E15.5 were microdissected under a microscope in 1× PBS on ice, and the subelements of the pelvis were separated, stripped clear of soft tissue, and collected in microcentrifuge Eppendorf tubes containing 200 μl of 5% FBS/DMEM. Because of their diminutive sizes, the pubis and ischium, constituting the inferior pelvis, were collected together. The right and left sides of each pelvic subelement were pooled from each sample and pooled across a litter of eight specimens. Each pooled sample was then subjected to 1% collagenase II (VWR, 80056–222, Radnor, Pennsylvania) digestion for 2 hours at 37°C rocking, mixing every 30 min. After placing on ice, samples were next filtered using a microcentrifuge filter setup by gently mashing the residual tissues through the filter followed by rinsing with 5% FBS/DMEM. Samples were then spun down at 500g at 4°C for 5 min. Samples were then counted and yielded technical replicates of 50,000 chondrocyte cell populations per girdle subelement, and these concentrations of 50,000 cells were then subjected to the standard ATAC-seq protocol as described previously (, ), modifying the protocol by using 2 μl of transposase per reaction [see Richard et al. () and Guo et al. ()]. The transposase reaction product was then purified using the Omega MicroElute DNA Clean Up Kit following the manufacturer’s protocols, eluted in 10 μl of warmed ddH2O, and stored at 20°C.

All samples were next subjected to PCR amplification and barcoding following Buenrostro et al. (2013, 2015) (, ). Ten microliters of transposed DNA were placed in a reaction containing NEBNext High-Fidelity PCR Master Mix, ddH2O, and barcoding primers. Following amplification, samples were transferred to new tubes and treated for fragment size selection using the OMEGA Bead Purification Protocol following the manufacturer’s instructions. The samples were eluted in 15 μl of TE, nanodropped, diluted to 5 ng/μl, and run on a Bioanalyzer for quality control and visualization of fragment size distribution. Before sequencing, sample concentrations were determined using the KAPA Library Quantification Complete Kit (KK4824). Samples were then sent out to the Harvard University Bauer Core Facility for paired-end sequencing on one lane of the Illumina NextSeq 500. See table S3 for sequencing sample, index primer, read count, and other information.

Computational analysis of ATAC-seq data began with checking read quality via FastQC. Adapters were then removed using NGmerge () version 0.2. Reads were subsequently aligned to mouse reference mm10 genome assembly with Bowtie2 v2.3.2 () using default parameters for paired-end alignment. Duplicate reads were removed using Picard’s MarkDuplicates function (version 2.9.0), and mitochondrial reads were also filtered out. The resulting .bam files were then used for peak calling via MACS2 software (version 2.1.1.2) (), using BAMPE and the following flags: --nolambda –bdg --verbose. Peaks reproducible across biological replicates were screened using a conservative IDR threshold of <0.01, as defined by the IDR statistical test (version 2.0.3). The IDR method identifies overlaps in peak calls across pairs of sample replicates by comparing ranked peak lists (using MACS2 q value) to define a reproducibility score curve. These paired peaks are then assigned a pointwise score based on this curve. They are then sorted, and all peaks that fall below an IDR threshold (here defined as 0.01) are taken as our final reproducible peak set across replicates. This yielded peak calls, including those called using normal stringency IDR thresholds (0.05) within ENCODE guidelines of greater than 50,000 peaks (). Following the identification of the IDR set of peaks for each pelvic subelement, the peaks were then reduced by filtering for brain. ATAC-seq was also performed on stage-matched E15.5 brain samples, and called peaks were then filtered to remove chondrocyte sequences that overlapped with brain sequences and are likely involved in regulating more general cellular processes. Intersections between brain peaks and pelvic peaks were done using bedtools version 2.27.1, which was also used to intersect pelvic sets to identify which peaks were shared across pelvic subelements and which were specific to only one subelement. GREAT version 4.0.4 was used to identify functional enrichments for peak sets.

Cross-species ATAC-seq region orthology analysis

Mouse ATAC-seq IDR datasets for pelvic subelements were lifted over from mm10 to hg19 using the UCSC Genome Browser LiftOver tool. They were then directly compared to corresponding human pelvic ATAC-seq IDR datasets using bedtools version 2.27.1 to identify regions of orthology between the two species and regions specific to either species. The reciprocal LiftOver was also performed. For this comparison, the human pubis and ischium datasets were combined into a merged pubis-ischium dataset due to the fact that these two subelements were taken together in the mouse because of their small size. Enriched overlap between human and mouse pelvic subelement sets was then tested using the regioneR package (version 1.8.1) () using the “permTest” function, generating 1000 randomized region sets as a background using the “circularRandomizeRegions” option and the “count.once” flag, with all other options set to defaults. Significance was assessed at P < 0.05 (table S3).

GWAS phenotypic analysis

GWAS summary statistics for hip circumference (), BMI-adjusted waist-to-hip ratio (), birth length (), birth weight (both fetal and maternal effects) (), and head circumference () were obtained and processed for use with the LDSC software version 1.0.0 (https://github.com/bulik/ldsc) (). Summary data were converted to an appropriate format for the “munge_numstats.py” script using custom R code. The latest version of the “baseline-LD” model available from the Alkes group (version 2.2) (http://data.broadinstitute.org/alkesgroup/LDSCORE) was downloaded. For partitioned heritability analysis, custom scripts were used to generate a modified version of the baseline-LD model with a reduced set of features including our regulatory regions sets (see table S4 for full list of features), which was subsequently used to recalculate LD scores with the “ldsc.py” script with the following options: plink files for the EUR subset of 1000G (obtained from the Alkes group; see above link), --ld-wind-cm 1, constraining to the same set of SNPs used in the original baseline LD model. Partitioned heritability of features was done using the “ldsc.py --h2” mode, using recalculated LD scores, extracted weights, reference 1000 Genomes Project phase 3 (1KG3) frequencies, the “--overlap-annot” flag, and all other settings left to defaults. P values for significant proportion of heritability captured by particular features was corrected for the number of features tested (n = 26) using Benjamini-Hochberg (BH) correction. Results are shown in table S4, with statistical significance defined as Padj < 0.05. Analyses were performed on both brain-filtered pelvic subelement sets and the smaller specific pelvic subelement sets, although, given issues of sparsity, the results of the latter were not considered reliable ().

HAR analysis

A set of human sequences displaying evidence of nucleotide acceleration in a variety of different contexts was aggregated from multiple studies (see table S4) (–) and intersected with the subelement-specific sets (e.g., ilium-specific) as well as pelvic common regions and paired subelement sets (e.g., pubis-ilium shared sites) along with ATAC-seq data obtained from stage-matched E54 brain tissue and a previously published B lymphocyte dataset () using bedtools (version 2.29.1). Similar analyses were performed on mouse ATAC-seq datasets. To account for the larger average peak size in pelvic regions relative to our two control sets (brain and GM12878 ATAC-seq), all regulatory regions in all pelvic sets were padded to a fixed size of 500 base pairs (bp). Further accounting for differences in set size, intersections were calculated as intersections per base pair of sequence in a given set. A background distribution was made by generating 10,000 random region sets consisting of 8332 (total number of subelement-specific regions) regions from across the genome with a constant length of 500 bp via the bedtools “random” function, followed by set intersections with HARs. Intersections per base pair values for both target and background sets were plotted as a histogram (Fig. 4D and fig. S8); the distribution of these values was assessed using the “qqnorm” (R base; version 4.0.3) and “qqPlot” (car version 3.0.8) functions—these suggested possible deviations from a normal distribution. Accordingly, the “fitdistrplus” package (version 1.1.1) was used to determine an appropriate distribution to fit the data. Given limitations in fitting a number of distributions (e.g., gamma, beta, and log-normal) for positive values, background sets with no accelerated region overlaps were removed for curve fitting (n = 9986 after filtering). The “descdist” function was initially used to assess curve behavior; goodness-of-fit statistics (from the “gofstat” function) for gamma, beta, exponential, and log-normal distributions were subsequently compared, with the beta distribution subsequently selected [this choice also being appropriate given the fractional nature of the data points ()]. Beta distribution parameters (“shape1” and “shape2” in the R implementation of “pbeta”) were fit using a bootstrap method (“bootdist” from “fitdistrplus”), with the median parameter estimates from 1000 samples used to define the distribution for significance testing of target intersections per base pair values with “pbeta” (upper-tail P values). Results were subsequently adjusted using BH correction, along with those obtained using the normal cumulative distribution function (CDF) distribution (“pnorm” in base R) (table S4). Regulatory regions intersecting regions of acceleration were associated with the closest annotated TSS with the HOMER (version 4.11) () “annotatePeaks.pl” script (table S4). Genes associated with regulatory regions from multiple pelvic sets were collapsed, with the final unique gene set used with the MGI Batch Query tool to find associated MP for these genes (table S4).

Cross-species conservation

To examine interspecies conservation, phyloP () scores for 20 primate species (phyloP20ways) were obtained from University of California, Santa Cruz (UCSC) () and used to extract per–base pair conservation scores for subelement-specific sets as well as pelvic common regions and paired subelement sets. To aggregate per–base pair conservation scores over the length of different regulatory regions (see below), for these analyses, regions were fixed to a constant size of 1500 bp (centered on region centers), representing the average peak size for all the pelvic ATAC-seq region sets generated. For analysis of trends in phyloP score over the length of a peak (fig. S9A), conservation scores were averaged per base pair across all ATAC-seq regions in a set and plotted as a function of distance along the region. A set of 1500-bp regions (matched for set size) randomly sampled genome-wide was also generated, and average per–base pair conservation scores were calculated as a background. For analyses comparing primate conservation to human-chimpanzee divergence (here and below), aligned sequences for human and chimpanzee were extracted from a primate alignment (multiz20way) obtained from UCSC, with sequence identity (%ID) for a given ATAC-seq region calculated as the number of nucleotide matches divided by total (ungapped) sequence length. These values were compared to phyloP20ways scores averaged over the length of the same region and subsequently plotted in order of increasing %ID (as seen in fig. S9B). Regions sorted by %ID were sliced into the top and bottom 25% for each sequence set and were analyzed using GREAT (version 4.0.4) (table S4).

Note about GREAT: GREAT takes an input set of genomic regions along with a defined ontology of gene annotations; first, it defines regulatory domains for all genes genome-wide and then measures the fraction of the genome covered by the regulatory domains of genes associated with a particular annotation (e.g., “cartilage development”). These fractions are used as the expectation in a binomial test, counting the number of input genomic regions falling within a given set of regulatory domains, which results in the reported significance of association between an input region set and a particular GO term. GREAT also performs a more traditional gene-based hypergeometric test to test for significance of region set–ontology association. The program returns a set of enriched ontologies sorted by the joint rankings of false discovery rate–corrected binomial and hypergeometric tests, as reported here in tables S1 to S4. For this study, we chose terms in the GREAT output with relevance to chondrocyte biology—pelvic-specific anatomy is not well annotated, particularly at the genetic level, and it is well known that there is a bias toward GO annotations for other phenotypes (such as immune function), particularly for phenotypes involving better-studied genes that tend to have increased GO term representation ().

Defining differential accessibility regions

To capture more subtle differences in accessibility, rather than the more stringent non-overlapping peak definitions of subelement-specific peaks (e.g., ilium-specific peaks) used above, we implemented a differential accessibility methodology to define regions as being more/less accessible in a given pelvic subelement (relative to all others). All IDR-filtered, brain-filtered peaks from each pelvic subelement were pooled, with overlapping peaks merged using bedtools “merge.” We next defined 3-kb regions centered on each pooled peak (1.5 kb upstream/downstream); along these regions, a 1.5-kb sliding window (chosen based on the average called peak size across individual pelvic elements) was slid in 50-bp increments to generate a set of overlapping 1.5-kb regions. ATAC-seq read coverage within these regions was calculated using the “bedcov” function of samtools version 1.5 for each mapped .bam file corresponding to individual ATAC-seq samples. This sliding-window approach is used to identify the area of strongest cross-sample signal (i.e., ATAC-seq read coverage) around a pooled ATAC-seq peak, improving our confidence in defining differential accessibility (i.e., avoiding regions around an ATAC-seq peak with lower read coverage, which may make our analysis more sensitive to noise). For sets of overlapping windows corresponding to a single ATAC-seq peak, we applied a smoothening algorithm to eliminate extreme values occurring in overlapping windows for individual samples. Windows whose read coverage fell outside 1 SD of the set (of sliding windows for a single peak) were assigned the average read coverage of the two windows adjacent. This adjustment does not affect the later differential accessibility analysis but will affect the choice of sliding window assigned to represent this given peak and is done to avoid consistently using the most extreme windows to define differential accessibility (which may otherwise bias our results). Following this smoothening, for each window, we calculated the 75th percentile read coverage value across all pelvic subelements and samples (this was found to be more robust than mean read coverages, even after smoothening adjustment). The window with the greatest 75th percentile coverage was then selected as the representative region for that pooled ATAC-seq peak. Subsequently, raw read coverages for all optimized windows across all pelvic subelements/samples were imported as a matrix into DESeq2 version 1.26.0 using the “DESeqDataSetFromMatrix” function, with differential accessibility calculated using the “DESeq” function with tissue type as the main variable. For each individual pelvic subelement (i.e., tissue), we then calculated the log fold change (logFC) increase/decrease in accessibility for windows relative to the average of all other tissues (e.g., ilium-specific logFC relative to all other pelvic subelements) using contrast vectors and the “results” DESeq2 function.

Tissue-specific regulatory scores

We defined a “tissue-specific regulatory score” for all genes captured in our RNA-seq datasets, a concept inspired by methods for integrating ChIP-seq and RNA-seq datasets (). For each gene, we captured all chromatin accessibility regions (i.e., the optimized windows above) within 100 kb upstream/downstream of the TSS. For each accessibility region, we took the calculated tissue-specific logFC value (e.g., ilium-specific logFC accessibility for a given putative enhancer) and scaled it by distance to TSS using the following formula:where the second distance-scaling term is taken from Tang et al. (2011) (). These per-region regulatory scores were summed across all regions to give a single tissue-specific regulatory score. Thus, a single gene will have four such scores, one for each pelvic subelement.

Accessibility region conservation

phyloP20ways per–base pair conservation scores () were averaged over the lengths of all regions used in calculating the above tissue-specific regulatory scores using the “bigWigAverageOverBed” program from UCSC. Across all genes, we scaled these region conservation scores by distance to nearby gene TSS (as above), summing to a single score of conservation per gene.

Accessibility region sequence divergence

For all regions used in calculating the above tissue-specific regulatory and conservation scores, we calculated human-chimp %ID (as previously described). As a more stringent filter for human-derived modifications, we further calculated chimp-gorilla %ID for these regions [using aligned gorilla data from UCSC () and the same methods as for human-chimp]. Across all regions, we performed a linear regression between human-chimp and chimp-gorilla %ID, for which a significantly positive correlation was observed (Pearson’s correlation: 0.47, P < 1 × 10−16). We then calculated the residualized human-chimp %ID for each region using this regression—positive values being “more human-chimp similarity than expected based on chimp-gorilla” and negative values being “less human-chimp similarity than expected” (the latter suggesting possible human-derived modifications). Per-region residualized %ID values were subsequently scaled for distance to TSS for all genes (as described above) and summed for a single per-gene score of human-chimp regulatory divergence.

Rank aggregation

To isolate those genes that are of greatest interest in possible tissue-specific regulatory network evolution, we filtered for several key characteristics. Genes that are differentially expressed in a given pelvic subelement, and whose nearby putative enhancers meet the following criteria: (i) they are biased in their accessibility between pelvic subelements; (ii) they are largely conserved in their sequences across primates, suggesting functionality in determining pelvic morphology; and (iii) they show divergence between humans and chimpanzees, suggesting that natural selection has acted on them in shaping the human pelvis. Our three different per-gene scores, based on the presence of nearby accessibility regions (tissue-specific regulatory score, conservation score, and human-chimp divergence score), were combined with RNA-seq data for these genes, calculated as logFC expression in a given pelvic subelement relative to all others. For each individual pelvic subelement, we considered all genes with increased expression in the said tissue and subsequently ranked this set of genes on the basis of these four metrics to yield four sets of rankings. These rankings were aggregated using the RobustRankAggreg library in R version 1.1 (), which seeks to identify genes that are consistently ranked higher across different ranking sets, assigning a significance value per gene for their rankings relative to a null hypothesis of uncorrelated inputs. For an individual pelvic subelement, we retained all genes for which robust rank aggregation (RRA) significance values were below 0.05. We subsequently looked for gene-set enrichments of these significant genes using clusterProfiler version 3.16.1, testing for enriched GO terms using as a background all genes originally considered in the RRA analysis. To collapse semantically similar GO terms, we used the “simplify” function in clusterProfiler using default settings (table S4).

To prioritize gene sets for which aggregate ranks biased toward human-chimp divergence, for each enriched (Padj < 0.05) simplified GO term, we calculated the average human-chimp divergence rank of all genes associated with said term. Terms were subsequently sorted on the basis of these averaged ranks, such that terms with greater average human-chimpanzee divergence were prioritized. For simplicity, we then asked whether the top 10 terms had significantly greater divergence than expected. For each term, we generated a set of 10,000 randomized gene sets, matching the number of genes associated with that term. The average human-chimp divergence ranks of these randomized sets were then calculated to establish a background rank distribution (e.g., see Fig. 4G). Target and background values were standardized, and statistical significance was assessed using a CDF of the standard normal distribution as implemented in the “pnorm” function in R (version 4.0.3), with significance values subsequently adjusted using BH correction (table S4). Similar statistical results were obtained when using raw divergence score values in addition to divergence ranks described above.

Species diversity patterns

Variation data from the 1KG3 () (n = 2504 individuals) in .vcf.gz format were obtained and intersected with pelvic subelement-specific sets and with the pelvic common region set and paired subelement sets using tabix (version 1.9) () to obtain variants occurring within these putative regulatory regions. The per–base pair distributions of sequence conservation (phyloP20ways) noted above (fig. S9A) indicated that the central portion of regions tended to be the most conserved across species, suggesting functional constraint particularly in these areas. Therefore, for the analysis of species diversity patterns, we used a fixed peak size of 500 bp (that used in the above analyses of accelerated region intersections). Chimpanzee (n = 25) and gorilla (n = 31) sequence data were similarly obtained via the Great Ape Genome Diversity Project (GADP) (). Peak sets were lifted over from hg19 to hg18 for use with the GADP datasets with the UCSC LiftOver utility and relevant LiftOver chain file. Resulting subset Variant Call Format (VCF) files were converted to tab format with the following Unix command, using bcftools (version 1.8) ():

bcftools query -f ´%CHROM\t%POS\t%ID\t%REF\t%ALT[\t%SAMPLE=%TGT]\n´ -o out.vcf in.vcf. Variant data for all region sets were downsampled to n = 25 (with replacement, 5 resamples for gorilla and 200 resamples for the human set) to match sample size for all comparisons based on the least-sampled species (chimp), using a custom R script.

Common variants were defined using a minor allele frequency (MAF) threshold of ≥0.05 for all datasets, filtering tab-formatted files using a custom Python script. Counts data were defined as the number of variants intersecting a given regulatory region and were averaged over resampled variant sets (see below). Counts data across apes were then compared within a given subelement set (e.g., ilium-specific regions) to compare intraspecies diversity of putative regulatory regions. Hurdle modeling was used to test for significant differences in both total number of sequences containing variants (hurdle) and degree of variation between species (counts), implemented using the “hurdle” function from the pscl () package in R (version 1.5.5). A binomial model was applied for the initial hurdle/zero-counts step, with subsequent counts modeling done using a negative binomial regression model. Tukey post hoc testing was performed using the emmeans package in R (version 1.4.7) for both hurdle/zero-counts and counts models, with significance assessed at Padj < 0.05 (table S4). In addition, subelement sets were compared to one another (e.g., ilium-specific versus pubis-specific) within a given species using the above methods. Variant per-sequence counts were visualized as boxplots using ggplot2 (version 2.3.3.2) with logged values (Fig. 5C and fig. S12B).

To look at sequence constraint of putative regulatory regions within humans, the set of pelvic subelement-specific regions were pooled, and a set of 10,000 randomly generated region sets, consisting of 2000 regions of 500-bp size, was generated using the bedtools random function. These sets were subsequently pooled, sorted, and merged using bedtools, with the resulting bed file used to extract variants from the 1KG3 set with tabix (version 1.9). The pooled set of pelvic subelement-specific regions was also used to extract variants from the 1KG3 set. In addition, several genomic features were extracted from the HOMER (version 4.0.4) set of genomic annotations provided with the program, including the following sets: exon, intronic, promoter-TSS, and TTS. Regions from RepeatMasker were also obtained from the UCSC Table Browser. These additional sets were used to extract variants from the 1KG3 set. The resulting files were filtered for duplicate variants and, subsequently, MAF ≥ 0.05 with bcftools (version 1.8). Variants falling within particular genomic regions in the random background, target (i.e., subelement-specific regions), and genomic annotation sets were then extracted using tabix. Variants extracted for each set were counted using vcftools (version 0.1.15) “--counts2 --stdout” arguments (). Variant counts were then adjusted to account for the number of base pairs within a given set. The background distribution of these values was investigated using the “qqnorm” (R base) and “qqPlot” (car package) functions to look for visible deviations from normality, for which no obvious deviations were observed. Values were standardized, and statistical significance was assessed using a CDF of the standard normal distribution as implemented in the “pnorm” function in R (version 4.0.3). P values for significant deviations from the background distribution were corrected for the number of sets (n = 11) tested using BH correction. Significance was defined as Padj < 0.05 (table S4). As a robusticity check against possible non-normal deviations in our background set, we also performed enrichment testing using a fitted beta distribution (following 1000 bootstrapping samples for distribution parameters) using the methodology described above. Significance testing was again performed with “pbeta” (base R) for all target intersections/base pair values, with BH correction applied for the number of sets tested (n = 11) (table S4).

A similar sequence constraint analysis was also performed for chimpanzees. Pelvic subelement-specific regions were pooled and lifted over to hg18 using the LiftOver utility; a set of 10,000 randomly generated region sets, consisting of 2000 sequences of 500 bp, was generated using the bedtools random function. Randomized sequence sets were subsequently pooled, sorted, and merged using bedtools, with the resulting bed file used to extract variants from the GADP set with tabix. Several genomic features were extracted from chimpanzee HOMER genomic annotations, including the following sets: intronic, promoter-TSS, TTS, and exon. These additional sets were lifted over to hg18 (flags as indicated above) and used to extract variants from the GADP. The resulting files were filtered for duplicate variants and, subsequently, MAF ≥ 0.05 with bcftools (version 1.8). Variants falling within particular genomic regions in the random background, target, and genomic annotation sets were then extracted using tabix. Variants per set were counted using vcftools (version 0.1.15) --counts2 --stdout arguments. Variant counts were then adjusted to account for the number of base pairs within a given set. The background distribution of these values was investigated using the “qqnorm” (R base) and “qqPlot” (car package) functions to look for visible deviations from normality, for which an obvious left skew was observed. For comparison with the above human analysis, background values were standardized, and statistical significance was assessed using a CDF of the standard normal distribution with the “pnorm” function in base R. P values for significant deviations from the background distribution were corrected for the number of sets (n = 11) tested using BH correction. Significance was defined as Padj < 0.05. See table S4. Given the left-skew distribution, the “fitdistrplus” package (version 1.1.1) was used to determine an additional distribution to fit the data, as described for the accelerated region enrichment analysis above, with a beta distribution selected on the basis of goodness-of-fit statistics. Distribution parameters were fit using a bootstrap method (“bootdist” from “fitdistrplus”), with the median parameter estimates from 1000 samples used to define the distribution for significance testing of adjusted variant counts with “pbeta” (base R). Results differed negligibly from those obtained using a normal CDF (table S4).

Human variant properties analysis

ARGweaver ()–estimated allele ages, based on the European subset of the 1000 Genomes project, were obtained from http://compgen.cshl.edu/ARGweaver/CG_results/download/bigWigs/?C=S;O=A and assigned to individual variants using the “bigWigAverageOverBed” utility from UCSC. Variants for which estimated allele ages were not available were excluded from subsequent analyses. Precomputed LINSIGHT () scores were obtained for the hg19 genome from https://github.com/CshlSiepelLab/LINSIGHT. These were similarly assigned to individual variants using the “bigWigAverageOverBed” utility. As a standard set, genome-wide variants were subset to those described in available UK Biobank summary statistics. For a given ATAC-seq regulatory region, we considered the presence of SNPs both within and nearby these regions—this was done to capture the possible effects of local linkage disequilibrium, in which a strongly associated variant may not fall immediately within a region, but a nearby proxy SNP (which may be the causal variant for the association signal) does intersect. This was done using the “window” function in bedtools to consider all SNPs falling within 1000 bp of a given region for all region sets (specific pelvic subelement sets and embryonic brain). We also considered all SNPs not falling within/nearby ATAC-seq regions as an additional, genome-wide control set. Per-SNP metrics (for LINSIGHT and ARGweaver) were compared across sets using the “aov” function in base R. Tukey post hoc analysis (controlling for the number of pairwise comparisons) was performed with the “TukeyHSD” function in base R (see table S4). To visualize these results (Fig. 5B), we used the “plotmeans” function from gplots version 3.1.1.