{'fraction_non_alphanumeric': 0.02804756568297846, 'fraction_numerical': 0.026615614493836384, 'mean_word_length': 5.2524328532502915, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The evolution and population diversity of human-specific segmental duplications

SUMMARY Segmental duplications contribute to human evolution, adaptation and genomic instability but are often poorly characterized. We investigate the evolution, genetic variation and coding potential of human-specific segmental duplications (HSDs). We identify 218 HSDs based on analysis of 322 deeply sequenced archaic and contemporary hominid genomes. We sequence 550 human and nonhuman primate genomic clones to reconstruct the evolution of the largest, most complex regions with protein-coding potential (n=80 genes/33 gene families). We show that HSDs are non-randomly organized, associate preferentially with ancestral ape duplications termed “core duplicons”, and evolved primarily in an interspersed inverted orientation. In addition to Homo sapiens-specific gene expansions (e.g., TCAF1/2), we highlight ten gene families (e.g., ARHGAP11B and SRGAP2C) where copy number never returns to the ancestral state, there is evidence of mRNA splicing, and no common gene-disruptive mutations are observed in the general population. Such duplicates are candidates for the evolution of human-specific adaptive traits.

; a = ancestral and b = derived), annotated SDs, and CN heatmaps from Illumina sequence mapping of chimpanzee, gorilla, and orangutan, and 10 African individuals (see Figure 1 for color coding index). We note timing estimates shown here are based on an evolutionary human-chimpanzee divergence time of 6 mya, though recent estimates suggest as much as 12 mya. Muliple-species sequence alignments were generated using the human reference genome (GRCh38) or CH17 BAC contigs for human paralogs and orthologs from NHPs, including chimpanzee (panTro2 or CH251 BACs), gorilla (CH277 BACs), and orangutan (ponAbe2 or CH276 BACs). Trees were generated using MEGA6 2 using the neighbor-joining method 3 . See the confidence probability (multiplied by 100) that the interior branch length is greater than 0, as estimated using the bootstrap test (500 replicates shown next to the branches 4,5 ). The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Kimura 2-parameter method 6 and are in the units of the number of base substitutions per site. All positions containing gaps and missing data were eliminated. The alignments used to generate pictured trees were inputs into Tajima's relative rate test 7 . Subsequent timing estimates were performed using pairwise comparisons of human and chimpanzee orthologs using the complete deletion option (with all outgroups removed from the alignment). Alignment lengths and timing estimates are shown in Supplementary Table 8.
Supplementary Figure 9. Human copy number diversity. Overall average CN was calculated per individual from read depth produced from Illumina mappings across a set region defining each duplication (Supplementary  Figure 10). The ancestral ROCK1A resides on human chromosome 18q11 while the human-specific ROCK1B shows polymorphism at 18p11.32, with one individual homozygously deleted for the paralog. The red signal adjacent to the blue signal indicates human-specific ARHGAP11B and the red signal adjacent to the green signal indicates ancestral ARHGAP11A. This same procedure was performed for individuals with suggestive rare homozygous HSD deletions for GPRIN2/NPYR4 (not pictured), MIR4267 (not pictured) and ROCK1 ( Figure 3B). See the Supplementary Supplementary Figure 13. Gene expression analysis of duplicated genes using a k-mer approach. Expression analysis using the GTEx RNA-seq dataset 8 . The heatmap shows a log 2 scale of the median normalized read counts per subtissue, based on unique 30-mers that differentiate between the paralogs (Methods; Supplementary Table 16 Red hexagons represent locations of identified nonsense and frameshift mutations with the total number of cases and controls with variants shown above each site. Colored boxes represent predicted functional domains (using SMART (http://smart.embl-heidelberg.de/)) with labels beneath. Schematic created using IBS 10 . Details of mutations are described in Supplementary Table 20.
Supplementary Figure 16. Pairwise comparison of genomic organization of human and chimpanzee at the 7q11.23 distal breakpoint using Miropeats. A Miropeats comparison of the human and chimpanzee contigs shows the pairwise differences between the orthologous regions. Lines connect stretches of homologous regions (threshold s = 1000) and match the arrow colors when they connect SD blocks (defined in Figure 6A). Additional annotations include: wholegenome shotgun sequence detection (WSSD) in human and chimpanzee, indicating duplicated regions identified by sequence read depth 11 ; DupMasker 12 ; and exons of genes. Gaps in the chimpanzee contig are shown as red boxes.
Supplementary Figure 17. SD configuration of chromosome 7q11.23 in great apes. Tiling paths of BACs were selected and sequenced using Illumina Nextera, capillary, and SMRT methods. Finished BAC sequences were compared and assembled into contigs (shown as gray lines beneath each haplotype). In some cases, BACs could not be assembled due to complex structural configurations contained within the clone, and Illumina Nextera sequencing was used to deduce the structural configuration and connections of larger contigs (shown as dashed gray line). SDs were identified within contigs from each species, as previously defined 13 , and shown as colored arrows: A is orange, B is red, and C is green. Additional SDs from a distal region are also depicted as blue (D) and purple (E) arrows in gorilla, as they were transplanted to BP2 via large-scale inversion.  Figure 18. Pairwise comparison of human reference (GRCh38) and newly constructed contig using CH17 BAC clones. A Miropeats comparison of the human and chimpanzee contigs shows the pairwise differences between the orthologous regions. Lines connect stretches of homologous regions (threshold s = 1000) and match the arrow colors when they connect SD blocks (defined in Figure 3C). Additional annotations include: WSSD in human and chimpanzee, indicating duplicated regions identified by sequence read depth 11 ; DupMasker 12 ; and exons of genes. In the newly created contig, errors in the human reference genome were corrected, including removal of a gap (pictured as a red box) and ~29 kbp of extra sequence (chr7:143336323-143347897 and chr7:143397898-143415538; hg38).

Supplementary Figure 19. An example of HSD gene family (CFC1) showing evidence of interlocus gene conversion from paralog-specific duplication analysis.
We leveraged singly unique nucleotide k-mers (SUNKs) calculated per individual from read depth produced from Illumina mappings across a set region defining each duplication in human populations, including the HGDP (N = 236; GRCh38), NHPs, archaic humans, a Denisovan and a Neanderthal. SUNK analysis of CFC1/TISP43, which we originally determined was "fixed" by our aggregate diploid copy number analysis, revealed that nearly all non-African individuals carry only the B duplicate paralog and almost no copies of the A paralog. Africans show a more diverse distribution with the B locus ranging from 1 to 4 copies while the A locus ranges from 0 to 3. In contrast, Neanderthal and Denisova both carry two copies of each paralog indicating that gene conversion has likely occurred within modern humans making the two paralogs nearly identical among the out-of-Africa populations (Supplementary Tables 22 and 23).
Orangutan (17) Gorilla (32) Chimpanzee (23) Bonobo (14) Denisovan & Neanderthal (2) Archaics ( Orangutan (17) Gorilla (32) Chimpanzee (23) Bonobo (14) Denisovan & Neanderthal (2) Archaics (3) Western European (59) Siberian (22) South Asian (23) Oceanic (21) East Asian (45) American (18) African ( (2) Archaics (3) Western European (59) Siberian (22) South Asian (23) Oceanic (21) East Asian (45) American (18) African ( (2) Archaics (3) Western European (59) Siberian (22) South Asian (23) Oceanic (21) East Asian (45) American (18) African ( Table  19), which would normally impact the RhoGAP domain of the full-length ARHGAP11A protein. Figure 25. Complex models of chromosome 10q11.23 HSD evolutionary history. BACs tiling across human chromosome 10q11.23 region were sequenced and assembled (representing human and additional great apes) and supercontigs were created. Estimates of sizes and evolutionary timing (human-chimpanzee distance; Supplementary Table 8) of events are denoted between each predicted intermediate genomic structure. SD organization is depicted as colored arrows and the orientations of intervening regions are shown with arrows. Models of the predicted evolutionary histories of the HSDs at all loci are depicted starting with the predicted human-chimpanzee common ancestor to the most common haplotype present in modern humans. A Miropeats comparison of the human and chimpanzee contigs shows the pairwise differences between the orthologous regions. Lines connect stretches of homologous regions based on a chosen threshold (s), defined as the number of matching bases minus the number of mismatching bases (s = 500) and match the arrow colors when they connect SD blocks. Additional annotations include WSSD in human and chimpanzee, indicating duplicated regions identified by sequence read depth 11 , DupMasker 12 , and genes. Figure 26. SD configuration of chromosome 10q11 in great apes. Tiling paths of BACs were selected and sequenced using Illumina Nextera, capillary, and SMRT methods. Finished BAC sequences were compared and assembled into contigs. SDs were identified within contigs from each species, as previously defined 13 , and shown as colored arrows: PTPN20/FRMPD2 as purple and GPRIN2/NPY4R as dark green. Additional SDs from a distal region are also depicted as light blue and light green. Gaps in the contigs are shown as red boxes.  Figure  3 and Supplementary Figure 11. See accompanying Supplementary Tables Excel file.

SUPPLEMENTARY TABLES Supplementary
Supplementary

SUPPLEMENTARY NOTE
Proposed mechanism and evolution of SDs on human chromosome 1. Previous directed sequencing efforts using the haploid BAC library have completely corrected this locus in the most recent build (GRCh38) 15,16 . This hotspot of duplication has been the target of at least six independent HSD events, including four primary duplications and two secondary duplications (Supplementary Figure 7). The majority of events occurred within bands 1p11.2 to 1q21.2 and contain seven HSD gene families, many of which have been previously implicated with neurological diseases and phenotypes. HYDIN, a gene in which homozygous mutations lead to congenital hydrocephalus in mouse 17 and primary ciliary dyskinesia in humans 18 , duplicated from chromosome 16q22.2 (357 kbp) to the 1q21.1 locus ~3.1 million years ago (mya) 19 creating partial paralog HYDIN2 (Supplementary Table 8, assuming a human-chimpanzee split of 6 mya). An 87 kbp duplication occurred ~4.7 mya, including GPR89 and a partial PDZK1, flank the 1q21.1 microdeletion locus in direct orientation, creating a potential substrate for non-allelic homologous recombination (NAHR). Our analysis also provided additional insight into the evolutionary history of SRGAP2, a gene previously shown to have duplicated in the human lineage 15 15 . We note timing estimates are based on an evolutionary human-chimpanzee divergence time of 6 mya, though recent estimates suggest as much as 12 mya.

Human copy number diversity
A notable exception to the general finding that higher copy number gene families have greater variance are duplication blocks containing GTF2I, GTF2IRD2, and NCF1 (diploid CN = 6; HGDP standard deviation (s.d.) = 0.27), located within the Williams-Beuren syndrome region on chromosome 7q11.23. This may be due to functional constraints of the genes.

Homo sapiens-specific duplications of BOLA2 and TCAF1/TCAF2
Comparisons with the genomes of Neanderthal and Denisova allowed us to identify the most recent duplications to emerge specifically on the lineage of Homo sapiens since divergence from archaic hominins. This includes duplications of the TRP channel-associated factors, TCAF1 and TCAF2, as well as the recently characterized BOLA2 gene family 21,22 , which shows some of the most significant expression differences between human and chimpanzee induced pluripotent stem cells 23 and whose SD mediates large-scale copy number variation associated with 1% of autism cases 24,25 . It is interesting that both gene families are involved with pathways related to "environmental sensing". For example, BOLA2 regulates intracellular iron levels while TCAF1/TCAF2 are associated with posttranslational regulation of the primary cold sensor, TRPM8 26 . Our sequence analysis suggests that both loci have undergone extensive expansion and restructuring creating not only additional copies but the potential for novel fusion genes and truncated copies to have emerged specifically Homo sapiens 21 .
For the novel TCAF-duplications finding, we characterized copy number of the TCAF locus and show that among modern humans copy number is uniform over a ~131 kbp cassette (involving A, B, and C SDs) and ranges from 2 to 7 ( Figure 4). The only exception to this pattern of copy number polymorphism was a Western European individual (HGDP00798) (A and B segments were discordant with C copy number), which we predict arose as a result of a NAHR-mediated deletion between directly oriented B1 and B2. Further, we identified 17 individuals (n = 2,367) where the copy number was consistent with the non-duplicated ancestral state observed in Denisova and Neanderthal. We remapped Illumina data of HGDP individuals to our newly constructed human contig of this locus in order to determine the paralogs experiencing polymorphism but were unable to perform SUNK analysis due to the extremely high similarity of the duplicate paralogs.

The complexity of HSD evolutionary history
To highlight the complex evolutionary history associated with HDS regions, we selected three loci (TCAF1/2 locus at chromosome 7q35 (main text); Williams-Beuren disease region at chromosome 7q11.23 (main text); and genomic hotspot-associated chromosome 10q11.23) for further investigation taking advantage of high-quality NHP BAC sequence generated as part of this study (N = 196; Supplementary Table 6) or taken from GenBank (N = 35). At the chromosome 10q11.21 locus, large-scale deletions and duplications have been identified in children with developmental delay with variable expressivity and penetrance 27,28 . HSD genes FRMPD2, PTPN20, GPRIN2, and NPY4R reside within two separate SDs proximal to the disease-associated region ( Supplementary Figures 25 and 26). Our data predict an initial inversion of 589 kbp, which resulted in a duplication of a 122 kbp segment containing a partial paralog of FRMPD2B and a full-length paralog of PTPN20B (0.389 ± 0.019 human-chimpanzee distance). The inversion breakpoint maps within PTPN20A and truncates the likely ancestral version of this gene by removing the last two exons, leaving the human duplicate PTPN20B as the only functional paralog. A 489 kbp duplication containing full-length GPRIN2 and NPY4R along with additional great ape-duplicated genes occurred potentially concurrently with the previous event (0.380 ± 0.014 human-chimpanzee distance). Comparing this region in human and chimpanzee identified an additional 550 kbp inversion that included the GPRIN2/NPY4R SD and adjoining proximal region.
Our BAC-based targeted sequencing of the 7q35 locus (containing TCAF1 and TCAF2) not only eliminated the gap annotated in the reference sequence (GRCh38) but also dramatically reorganized the structure of the region removing incorrectly assigned paralogous sequence, including ~29 kbp of extra sequence (Supplementary Figure 18). Errors also existed within the sequence itself as evidenced by dramatically different sequence identities of paralogs between the reference and the corrected contig, which would have incorrectly estimated the timing of the final TCAF1 and TCAF2 duplications before the split with Denisova and Neanderthal.

Additional important HSD features emerging from targeted sequencing of loci
We distinguish primary duplications (as those ancestrally derived from the unique ortholog in NHP species) from secondary duplications (those that have duplicated from an HSD). Secondary duplications (average 437 kbp) are twice as large when compared to primary duplications (average 136 to 196 kbp). The size discrepancy between older and younger duplications may be explained, in part, by subsequent internal deletions accruing over time within older duplications, as was the case for the SRGAP2B paralog 15 . Nevertheless, unlike primary duplications, secondary duplications are always distributed intrachromosomally in inverted orientation with the majority mapping within 5 Mbp of their progenitor.
We also found an enrichment of HSDs near primate ancestral "core duplicons." It is possible that the complex genomic architecture associated with core duplicons residing in proximity to HSDs acted to perturb the replisome leading to forkstalling and template switching (FoSTeS) or represented sites of fragility predisposing to microhomology-mediated breakinduced repair (MMBIR). The enrichment of inverted SDs also emphasizes the intimate association between inversions and the dispersal of SDs [29][30][31] . The association of inversion breakpoints with SDs has contributed to this type of variation being underestimated in studies of human genomes [32][33][34] .

Additional information on putatively functional HSD genes
While the functional relevance of most transcripts mapping to HSDs remains to be determined, several recent studies have suggested that these regions may encode genes relevant to human neurocognitive and neuroanatomical adaptation. The human-specific duplicate SRGAP2C, for example, has been shown in vivo to alter dendrite formation and potentially spine density in developing neurons 15,20,35 , while the HSD gene ARHGAP11B appears to promote apical basal radial progenitor amplification in the subventricular zone 36 . Microinjection of the ARHGAP11B RNA into the developing mouse brain increases the number of basal radial glial divisions leading to cortical expansion as well as gyrification.
There are also a handful of HSD-associated genes implicated in immune response, including CD8B (CD8 antigen, beta polypeptide), which encodes the beta chain of the heterodimeric CD8 glycoprotein responsible for recognizing antigenic peptides on the surface of immune T cells. This gene was previously thought to have duplicated in a common ancestor of human, gorilla and chimpanzee 14 . Our analysis shows a gorilla-specific 40 kbp duplication of CD8B exon 1 not shared with human, chimpanzee, bonobo or orangutan (Supplementary Figure 20). The gorilla duplication has no overlap with our HSD, which encompasses CD8B exons 2 to 5; our timing estimate of CD8B-B, which places the duplication at just after the human-chimpanzee divergence from a common ancestor, supports this finding ( Figure 2).

Genomic characterization of SMN1 locus
When we compared our new chimpanzee reference for chromosome 5q13.3 containing SMN1 37,38 , we found that human and chimpanzee differ by 1.3 Mbp (Supplementary Figure 21). Despite having complete sequences of the two species, we were unable to delineate the precise mechanism of duplications leading from the assumed ancestral structure in the chimpanzee genome to the configuration in the human reference today. Our efforts were hindered by the presence of high copy (CN > 8), polymorphic, nearly identical, palindromic duplications peppered throughout this locus that likely arose in the last two million years. Sequencing additional human haplotypes will be necessary to delineate the precise mechanisms leading to the rapid expansion of this complex region.

Evidence of ancestral paralog deletions
We leveraged SUNKs to investigate copy number variation diversity at the level of individual paralogs 39 (Supplementary  Tables 10, 23 and 24). Among the 72 gene paralogs assayed here, we identified apparent homozygous deletions for 24 paralogs in at least one human individual. Genotyping implicated three ancestral paralogs as showing evidence of heterozygous deletions (Supplementary Figure 22); closer inspection of copy number heatmaps found that FRMPD2A and OCLN-A appeared to be truly deleted, though in the latter case gene conversion may be at play based on the known disease relevance of this gene 40 (Supplementary Figure 23). Meanwhile, though deletions appear to exist within portions of the HSD-containing CHRNA7, the gene itself appeared undeleted. These results highlight certain limitations of this SUNK approach, especially when paralogs share high-sequence similarity 15,41 .

Alternative derived gene model of ARHGAP11B places common LGD variant in intron of gene
A common LGD variant was identified in ARHGAP11B, an HSD previously implicated in neuronal migration 36 , suggesting the gene may not be functionally relevant (Supplementary Tables 10 and 19). Closer inspection shows the LGD variant falling within an ARHGAP11B-specific intron due to an alternative gene model of ARHGAP11B compared to its ancestral counterpart (Supplementary Figure 24). We propose that systematic characterization of gene models for other HSD genes is necessary to fully understand the impact of LGD variants on function.

Identification of HSD regions from Illumina data
Whole-genome Illumina short-read mappings against the human reference genome (GRCh37/hg19) of a diverse, highcoverage set of 236 human genomes from the HGDP 42 and a set of 86 NHP genomes 43 were used to estimate aggregate copy number in 500 bp windows using previously described methods 39 . Windows with >90% of human genomes at copy number >2.5 and >90% of NHP genomes at copy number <2.8 were identified as HSDs. The copy number cutoff for NHP genomes was empirically determined using previously identified HSD genes to reduce false negatives 39,44 . These regions were merged if within 1 kbp of each other with the final HSD regions, including merged windows 5 kbp in size or greater.
To connect adjacent HSD regions punctuated by higher ancient "core" duplicons 45,46 (found duplicated across all great apes), HSDs within 20 kbp of each other were additionally merged if, in all the windows between the regions, >90% of human genomes had copy number >2.8.

Validation of HSD regions
We intersected HSD regions with SDs identified using whole-genome shotgun sequence detection (WSSD) using Sanger sequence read depth 11 and whole-genome analysis comparison (WGAC) 1 methods (Supplementary Table 1). We compared the 37 regions that did not intersect previously identified SDs to a genome assembly of the CHM1 haploid hydatidiform mole (NCBI Assembly PacBioCHM1_r2_GenBank_08312015) using BLASR to look for orthogonal evidence that these regions might be resolved or collapsed duplications 47,48 . We counted a query region as resolved if it had multiple mappings with match length minus edit distance greater than 90% of the query length. To find collapsed duplications, we calculated coverage across the CHM1 assembly in 100 bp windows and identified regions >5 kbp of elevated coverage, where the threshold for elevated coverage was set at the third quartile plus two times the interquartile range (or 72.7X).

Sequencing of BAC clones
DNA from CH17, CH251, CH276 and CH277 BAC clone libraries was isolated, prepped into barcoded genomic libraries, and sequenced (150 bp paired end) on an Illumina MiSeq using a Nextera protocol 31 . Sequencing data (~300-fold coverage) were mapped with mrFAST 49 to the human reference genome (GRCh37) and SUNKs were used to discriminate between highly identical SDs 39 . Some clones were subject to a hierarchical clone-based strategy with highquality capillary fluorescent-based sequencing. This entailed the construction of genomic libraries, sequencing of paired-end shotgun libraries, and assembly of inserts into a finished sequencing contig. The majority of clones were subject to PacBiobased sequencing entailing the preparation of SMRTbell libraries sequenced using PacBio RSII C2P4 or C4P6 chemistry (one SMRT cell/BAC sample with two 45-minute movies) for a subset of clones spanning HSD regions based on Illumina sequence analyses. Inserts were assembled using Quiver and HGAP (Hierarchical Genome Assembly Process) as described 50 . BACs were assembled into contigs with PacBio-and capillary-sequenced clones using Sequencher and compared to the human reference genome using Miropeats 51 and BLAST 52 . Sequences for contigs not included in the human reference can be found in Supplementary Dataset 1. Additional BAC clones mapping to HSD regions in human and NHPs were identified within GenBank, previously sequenced by the Wellcome Trust Sanger Center and McDonnell Genome Institute at Washington University as part of separate projects, which we included in breakpoint and evolutionary analyses (Supplementary Tables 3 and 6). These clones were not included in counts of BACs sequenced for this study.

Breakpoint identification
A combination of BLAST 52 , BLAT 53 , and WGAC 1 methods were used to identify HSD paralogous regions. Sequences (±500 bp) were subsequently extracted from the human reference genome (GRCh38) or from BAC-assembled contigs and pairwise comparisons performed and visualized using BLAST and Miropeats 51 to identify the maximal duplicon breakpoints. When precise breakpoints could not be defined due to flanking core duplicons, minimal and maximal breakpoints were reported for both ancestral and duplicate paralogs. When comparing sizes of HSDs, we used the greater of the two sizes between ancestral and duplicate paralogs. In size estimates and comparisons of HSDs, we excluded the NAIP-C HSD due to uncertainty of its status as a primary or secondary duplication. We found that some sets of duplication sizes were drawn from a non-normal distribution using the Shapiro-Wilk test (e.g., maximum size of HSDs in period III); hence, we compared all duplications sizes using nonparametric tests. Specifically, we applied the Wilcoxon-Mann-Whitney test-using the wilcox.test() function in R (version 2.15.0)-when comparing primary versus secondary duplication sizes and the Kruskal-Wallis rank sum test to assess size differences across all three evolutionary periods. After identifying significant differences between time periods of minimum primary duplication sizes, we applied a Wilcoxon-Mann-Whitney test post hoc to identify the duplication waves that were significantly different and adjusted for multiple comparisons using the Holm method. All cumulative sizes were reported as medians.

HSD clustering simulations
We simulated a null distribution by shuffling 218 HSD regions (defined in Supplementary Table 1; GRCh37) and an unbiased set of 18 primary HSDs (of our total 24) containing previously identified HSD genes 39 (derived genomic coordinates only, Supplementary Table 5; GRCh38) within the same chromosome 10 million times using BEDTools shuffle (v2.23.0) and, when multiple duplications occurred on a single chromosome, calculated the distance to the nearest duplication using midpoint coordinates. In this and all subsequent simulations, shuffled intervals were not allowed to intersect non-scaffold gap regions. We calculated the median distance for each iteration of the simulation and compared this distribution to the empirical value. We also determined midpoint distances of the 18 HSD primary duplication to the nearest non-primary HSD defined by WGAC 1 except shuffling was performed one million times. All analyses were repeated allowing shuffling genome-wide.
We recalculated the sequence and location of core duplicons using coordinates of all duplicon clades previously defined 45 . As the original definition of a core duplicon required all duplicons to occur on the same chromosome, we limited our analysis to clades with intrachromosomal duplicons. We extracted the GRCh35/hg17 sequence corresponding to each duplicon and aligned all duplicon sequences against themselves with BLASR using parameters tuned for high-quality queries (-affineAlign -affineOpen 8 -affineExtend 0 -bestn 30 -maxMatch 30 -sdpTupleSize 13) 54 . We retained all alignments >100 bp between duplicons from the same clade except for alignments of each duplicon to itself. For each clade, we clustered all alignments that reciprocally overlapped by 50% or more and selected the cluster with the most components (i.e., alignments from other duplicons in the same clade) as the representative core for the clade. We aligned the core sequences to GRCh38/hg38 with BLASR using the same parameters described above and filtering matches <75% of the query length. A null distribution of median distance to the nearest core duplicon was created using simulations performed as described above with shuffling of the 18 HSD primary duplications 10 million times.

Evolutionary analysis
Sequences from HSD orthologs were identified and extracted from genome reference or BAC assemblies for chimpanzee (panTro4 and CH251), gorilla (gorGor4.1 and CH277), and orangutan (ponAbe2 and CH276). Multiple sequence alignments (MSAs) were generated using MAFFT and included the maximal shared genomic regions of human paralogs and nonhuman orthologs excluding any flanking core duplicons 55 . Alignments were visualized for manual editing using Jalview 56 . MSAs are available in Supplementary Dataset 1. Phylogenetic analyses were performed using MEGA6 2 (Supplementary Dataset 2). If a full-length alignment did not pass the Tajima's D relative rate test 7 , using orangutan as the outgroup, pairwise sequence identities were calculated across 500 bp sliding windows with 100 bp increments using the PopGenome statistical package 57 and visualized using ggplot in R. Portions of the alignment that exhibited aberrant spikes in sequence identity were excluded and phylogenetic analyses were repeated on a refined region (Supplementary Table 8). In the case of HSD regions containing SRGAP2, corrections were made to distance estimates to account for differences in substitution rates of paralogous regions as previously described 15 . Duplication mechanisms were predicted using a combined approach of defining ancestral paralogs/configurations using genomic synteny taken from chimpanzee and/or orangutan and evolutionary timing estimates to predict the order of rearrangements.

Detection of interlocus gene conversion
We implemented two approaches to detect signatures consistent with recent interlocus gene conversion among HSD paralogs. First, we created pairwise sequence alignments between all annotated HSDs containing the ancestral and duplicate paralogs. We next calculated the identity of two aligned sequences over 2 kbp sliding windows across the alignment with a stepwise increment of 100 bp (Supplementary Dataset 1). Pairwise sequence identity was plotted against the length of the alignment and potential interlocus gene conversion events were identified by the presence of sharp sequence identity transitions from low <99% to high >1.0. All sequence identity plots are included in Supplementary Dataset 1. We next used the program GENECONV to scan each MSA used in the phylogenetic analysis of HSD regions. GENECONV identifies pairs of sequences with longer than expected tracks of 100% sequence identity conditioning on the overall pattern of variable sites in the alignment 58,59 . The program was run using default parameters and tracks with a global P < 0.05 were considered significant for follow-up analysis 58 . Tracks of perfect sequence identity were then hardmasked from the original MSAs using BEDTools and the resulting branch lengths were then recalculated from the phylogenies (Supplementary Table 8).