Gibbons (Hylobatidae) are critically endangered¹ small apes that inhabit the tropical forests of southeast Asia (Fig. 1) and belong to the superfamily Hominoidea along with great apes and humans. In the primate phylogeny, gibbons diverged between Old World monkeys and great apes, providing a unique perspective from which to study the origins of hominoid characteristics.

Gibbons have several distinctive traits, the most striking of which is the unusually high number of large-scale chromosomal rearrangements in comparison to the inferred ancestral ape karyotype². The four gibbon genera (Nomascus, Hylobates, Hoolock and Symphalangus) occupy different regions of southeast Asia and bear distinctive karyotypes, with diploid chromosome numbers ranging from 38 to 52 (Fig. 1). Given the relatively recent differentiation of these genera (4–6 million years ago (Myr ago), this constitutes an extraordinarily fast rate of karyotype change.

In order to investigate the mechanisms behind the plasticity of the gibbon genome, understand the evolutionary relationships among the four extant gibbon genera and study the evolution of putatively functional sequences related to gibbon-specific adaptations, we sequenced and assembled the genome of a female northern white-cheeked gibbon (Nomascus leucogenys) named ‘Asia’. The reference assembly (Nleu1.0) provides on average 5.7-fold Sanger read coverage over 2.9 gigabase pairs (Gb) (Table 1 and Supplementary Table ST1.1). Our quality assessment (Extended Data Fig. 1) confirmed its equivalence to other Sanger sequence-based non-human primate draft assemblies (such as the orangutan or rhesus macaque^{3, 4}) (Supplementary Information section S1, Supplementary Data Files 1 and 2). We also obtained ~15× whole-genome shotgun (WGS) short-read data (Illumina) for two individuals of each gibbon genus and high-coverage exome data (>60×) for two of the same individuals in order to derive error models for single nucleotide polymorphism (SNP) calls (Supplementary Information section S2; Supplementary Tables ST2.1–2.3).

Nleu1.0 scaffolds were aligned against the human reference (GRCh37) to be ordered and oriented into 26 chromosomes (Nleu3.0) under extensive guidance by cytogenetic data. The reshuffled nature of the gibbon genome was especially evident when human–gibbon chromosome alignments were compared with those between human and great apes, rhesus macaque (Old World monkey) and marmoset (New World monkey) (Fig. 2a). This higher rate of reshuffling applied only to large-scale chromosomal rearrangements (>10 megabases (Mb)), whereas smaller-scale rearrangements (10–100 kilobases (kb)) were comparable with other species (Fig. 2b) (Supplementary Information section S1).

We identified 96 gibbon–human synteny breakpoints in Nleu1.0 and classified them as to whether they could be defined at the base-pair level (class I, n = 42) or only narrowed to an interval due to greater complexity (class II, n = 54). As previously reported⁵, breakpoints were significantly depleted of genes (Supplementary Fig. SF5.2 and Supplementary Data File 3) and breakpoint intervals contained a mixture of repetitive sequences that inserted exclusively into the gibbon genome^{2, 5, 6} (Fig. 2c). To assess breakpoint segmental duplication content, we identified gibbon-specific segmental duplication using in silico methods followed by experimental validation (Extended Data Fig. 2, Supplementary Fig. SF3.1, Supplementary Information section S3 and Supplementary Data File 4). Of note, both gibbon-specific segmental duplication and gene family expansion analyses suggested the gibbon genome has not undergone a greater rate of duplication than other hominoids, further supporting a model in which accelerated evolution has been limited to gross chromosomal rearrangements (Supplementary Information section S6, Supplementary Fig. SF6.1).

Segmental duplication enrichment was the best predictor of gibbon–human synteny breakpoints, as shown through permutation analyses (P value < 0.0001); however, breakpoints were also enriched for Alu elements (Supplementary Table ST5.1; Supplementary Information section S5; Supplementary Fig. SF5.2). Although non-allelic homologous recombination between highly similar sequences can mediate large-scale rearrangements⁷, the majority of gibbon chromosomal breakpoints bore signatures of non-homology based mechanisms (Fig. 2c). These included the insertion of non-templated sequences (2–51 nucleotides (nt)) and/or the absence of identity, suggesting non-homologous end joining. The presence of micro-homologies (2–26 nt) in a small portion of the breakpoints (13/42) pointed to additional alternative mechanisms such as microhomology-mediated end joining⁸ or microhomology-mediated break-induced replication⁹. The origin of the complex structure of breakpoint intervals (class II) was less obvious and reinforced the observation that repeats have the tendency to accumulate at the breakpoints.

To explore the possibility that chromatin conformation, rather than sequence, might predispose regions to breakage, we investigated the relationship between gibbon breakpoints and CCCTC-binding factor (CTCF), an evolutionarily conserved protein with multiple functions, including mediating intra- and interchromosomal interactions¹⁰. We performed chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) of CTCF-bound DNA using lymphoblast cell lines established from eight gibbon individuals (Supplementary Information section S5). We observed an enrichment of gibbon–human breakpoints in CTCF-binding events (P value = 0.0028), which increased when we considered a ~20 kb window centred around each breakpoint (P value of < 0.0001). Notably, this enrichment was maintained only for CTCF-binding events shared with other primates (human, orangutan and rhesus macaque)¹¹ but not those specific to gibbon (P value = 0.0019) (Supplementary Fig. SF5.4).

Thus, gibbon–human breakpoints co-localized with distinct genomic features and epigenetic marks; however, as many of these features were shared with other primates, other factors unique to the gibbon lineage must have been present to trigger the increased frequency of chromosomal rearrangements.

The gibbon genome contains all previously described classes of transposable elements that are mostly also present in other primates. One exceptional addition is the LAVA element, a novel retrotransposon that emerged exclusively in gibbons¹² and has a composite structure comprised of portions of other repeats (3′-L1-AluS-VNTR-Alu_-like-5′) (Fig. 3a). Searches of Nleu1.0 retrieved 1,797 LAVA insertions, 1,256 of which were 3′ intact elements, many carrying signs of target-primed reverse transcription (TPRT)¹³. The distribution of 3′ intact LAVA elements uncovered a significant overlap with genes (Pearson chi-squared, P = 0.017) and Gene Ontology (GO) analyses using the database for annotation, visualization, and integrated discovery (DAVID)¹⁴ showed a significant functional enrichment exclusive to the ‘microtubule cytoskeleton’ category (false discovery rate = 0.031, P value = 0.001) (Supplementary Information section S7 and Supplementary Data File 6) (Extended Data Fig. 3). Additional analyses with meta-pathway database tools^{15, 16} refined this enrichment to pathways related to chromosome segregation, including ‘establishment of sister chromatid cohesion’ and ‘mitotic metaphase and anaphase’ (Supplementary Table ST7.3). Genes with LAVA insertions include proteins that function as checkpoints for cell division and for spindle integrity/architecture (such as MAP4, CEP164 and BUB1B)^{17, 18, 19}, participate in kinetochore assembly and attachment to the spindle (for example, MAD1L1 and CLASP2)^{20, 21}, and have a role in chromosome segregation during cell division (for example, KIFAP3 and KIF27)²² (Extended Data Table 1).

Intragenic LAVA insertions were skewed toward introns (Pearson chi-squared, P = 0.0001) and were less frequent than expected when within <1 kb of the nearest exon junction (Extended Data Fig. 3). The majority (74%) of intronic LAVA elements were found in the antisense orientation. We speculated that intronic antisense LAVA insertions may cause early transcription termination by providing a polyadenylation site in the antisense orientation, as previously described for L1 elements^{23, 24} (Extended Data Fig. 3). Indeed, we found 84.1% of the 3′-intact LAVA elements encoded a perfect polyadenylation signal at their 3′ end in antisense orientation.

To obtain experimental evidence that LAVA elements disrupt transcription, we performed a reporter assay in which the 3' end of a luciferase gene construct lacking a transcriptional termination site was fused to the 3'-terminal fragments of LAVA_E and LAVA_F elements, mimicking the arrangement observed in gibbon genes (Fig. 3b, left). Luciferase activity exceeding background level by ~50% was observed from the LAVA_F reporter construct (Fig. 3b, right), indicating faithful termination of luciferase transcription. Furthermore, 3′ rapid amplification of cDNA ends (RACE) experiments confirmed that the transcription termination site had been supplied from the LAVA element (Extended Data Fig. 3). Thus antisense intronic LAVA insertions can cause early transcription termination with some variability possibly due to the genomic context of the polyadenylation site, which explained the difference between the two reporter constructs.

We also investigated LAVA induced early transcription termination in vivo by analyzing RNA-seq data generated for the gibbon named Asia (Supplementary Table ST2.4). Specifically, we looked for paired-end reads only partially aligning to an antisense LAVA element due to untemplated residues and then identified cases for which the presence of a poly(A) tail was preventing full-length alignment. This analysis revealed that elements from a variety of subfamilies have the potential to cause early transcription termination, including those identified for LAVA elements inserted in the microtubule cytoskeleton genes (for example, LAVA_B2R2, LAVA_C4B, LAVA_B1R2) (Extended Data Table 1). We observed that early transcription termination occurred at relatively low levels as we identified a significant number of read pairs indicative of normal transcription and splicing for LAVA-terminated genes (Supplementary Table ST7.5). This is to be expected, as full inactivation of many of these genes would be lethal. On the other hand, as alternative splicing and RNA pol II transcript termination/polyadenylation are tightly coupled processes, LAVA-mediated early transcription termination could also act by differently affecting distinct isoforms and/or influencing the ratio between isoforms. Finally, LAVA insertions may also affect gene expression by functioning as exon traps, as shown for SVA elements²⁵. One putative example of an exon trapping event was identified for HORMAD2, a gene that monitors the formation of synapsis during crossover²⁶ (Supplementary Information section S7, Supplementary Table ST7.6, Supplementary Fig. SF7.1–7.2).

As genome reshuffling began in the common ancestor of all extant gibbon species, LAVA insertions must have occurred in key genes before the four genera diverged. We experimentally confirmed the mode and tempo of all 23 LAVA insertions in genes from the microtubule cytoskeleton category using both site-specific PCR and in silico methods (Extended Data Figure 4) and found that most of the insertions (15/23) were shared by the four gibbon genera (Supplementary Data File 6). Eleven of the genes match the structural requirements for early transcription termination and five of them are also shared. These genes include MAP4, involved in spindle architecture and CEP164, a G2/M checkpoint gene whose inactivation results in an aberrant spindle during cell division^{18, 19} (Extended Data Table 1).

We explored the relationship between LAVA family expansion and evolution of the gibbon lineage and, through analyses of diagnostic mutations, identified 22 LAVA subfamilies (Fig. 3c). In addition, we tested for the presence or absence of 200 LAVA loci from among the evolutionarily youngest elements in each subfamily (Extended Data Fig. 4) across 17 unrelated gibbon individuals and found that 52% of loci were shared among all four genera, whereas 27% were Nomascus specific. The remaining LAVA insertions showed a variety of confounding phylogenetic relationships consistent with incomplete lineage sorting (ILS) of ancestral polymorphisms, perhaps as a result of a rapid radiation of gibbon genera (Supplementary Information section S7; Supplementary Table ST7.1–7.2). We used a maximum likelihood method²⁷ to obtain age estimates for the 22 LAVA subfamilies. In the case of the two oldest subfamilies, LAVA_A1 and LAVA_A2, we obtained estimates of ~18 Myr ago and ~17 Myr ago, respectively (Supplementary Table ST7.3). A coalescent-based methodology implemented in the software G-PhosCS²⁸ using Nleu1.0 estimated a gibbon–great ape population divergence time of ~16.8 Myr ago (95% confidence intervals (CI): 15.9–17.6 Myr ago) assuming a split time with macaque of 29 Myr ago (Supplementary Information section S4). Hence, the LAVA element probably originated around the time of the divergence of gibbons from the ancestral great ape/human lineage.

The evolutionary history of the gibbon lineage and, in particular, the timing and order of splitting among the four genera, is still a subject of debate²⁹. To address this issue, we generated medium coverage (mean ~15×) WGS short read data for two individuals from each of the four genera, including two different Hylobates species (H. moloch and H. pileatus) (Supplementary Table ST2.1–2.2). Although phylogenetic analysis of assembled whole mitochondrial DNA genomes using BEAST³⁰ strongly supported monophyletic groupings for each gibbon genus, the branching order of the four genera remained unresolved (Supplementary Fig. SF9.1–9.2; Supplementary Information S9).

Neighbour-joining trees constructed from pairwise sequence divergence, k, across ~11,000 genic (200 base pairs (bp)) and ~12,000 non-genic (1 kilobase (kb)) autosomal loci supported a supermatrix sequence topology of (((Siamang (SSY), Hoolock (HLE)), Nomascus (NLE)), (H. pileatus (HPL)), H. moloch (HMO)) (Fig. 4a); nevertheless, bootstrap confidence for the node separating NLE and Hylobates was low (~52%). This topology was also the most frequently observed when constructing k-based unweighted pair group method with arithmetic mean (UPGMA) trees along the genome using non-overlapping 100-kb sliding windows. However, all 15 possible rooted topologies for the four genera were observed at considerable frequencies (Extended Data Fig. 5), consistent with the extensive ILS observed in the LAVA element analysis.

In order to infer the most likely bifurcating species topology amongst the four genera while taking into account ILS, we used a novel coalescent-based ABC methodology using the autosomal non-genic and genic loci (Veeramah et al., in the press) (Supplementary Information section S8). The topology described above had the highest combined posterior probability, though support was relatively low (P (model) = 17%) and other topologies, including one with NLE and Hylobates interchanged as the most external taxa, had comparable probabilities (Fig. 4a).

The estimated internal branch lengths under the best species topology using our ABC framework and G-PhoCS were very short, supporting a rapid speciation process for the four gibbon genera (Fig 4b, right). Given this observation and uncertainty in the best topology, we also estimated parameters under an instantaneous speciation model (Fig. 4b, left). Assuming an overall autosomal mutation rate of 1 × 10⁻⁹ per site per year, we placed the beginning of the speciation process at ~5 Myr ago under both models, with the two Hylobates species diverging ~1.5 Myr ago.

Consistent with the ABC analysis, SSY and HLE share the largest number of alleles across the whole genome (Supplementary Table ST8.5). However, NLE and the two Hylobates samples are both significantly closer to SSY than HLE as assessed by the D-statistic³¹. This result could be explained by two independent gene flow events between SSY and both NLE and Hylobates. However, fertile intergenic hybrids have yet to be observed either in the wild or captivity³²; an alternative explanation would be long-term population structure in the gibbon ancestral population. Both the ABC and G-PhoCS analyses suggest that the ancestral gibbon effective population size (N_e) was large (80,000–130,000), but neither of these frameworks can distinguish this from a structured ancestral population.

The coalescent-based analysis (Fig. 4a), along with estimates of genome-wide heterozygosity (Supplementary Fig. ST8.2), suggests a larger long-term N_e for both N. leucogenys and H. moloch compared to the other species. Analysis using the pairwise sequentially Markovian coalescent (PSMC) model³³ indicates that these two species underwent an increase in N_e during the Late Pleistocene era (500–100 thousand years ago (kyr ago) followed by a subsequent decrease in N_e 100–50 kyr ago (Fig. 4c) (Supplementary Information section S8). Fluctuation in N_e could result from changes in the actual number of individuals in the population, changes in population structure, and/or variable gene flow.

Accelerated substitution rates are a hallmark of adaptive evolution, and genomic regions with excess lineage-specific substitutions have been found to have functional roles³⁴. We identified 240 short (153 bp) median length) regions with accelerated substitution rates in the gibbon lineage (gibARs). We observed that gibARs were primarily intergenic (66%) and tended to co-localize near the same genes as LAVA elements (P value = 81 × 10⁻⁶; odds ratio of 2.74 (95% CI: 1.79–4.07)). Consistent with this finding, a GO enrichment test for genes within ± 100 kb of each gibAR (in comparison with background genes) revealed enrichment for the ‘chromosome organization’ category (Benjamini–Hochberg false discovery rate <5%) (Extended Data Fig. 6). Given evidence of functional roles gathered for human accelerated regions³⁵, we speculate that the gibARs may create functional elements (for example, enhancers or protein-binding domains) to modulate the transcriptional effect of local LAVA insertions (Supplementary Information section S12 and Supplementary Data File 9).

We assessed the potential presence of positive selection in 13,638 human genes with one-to-one orthologues in gibbon using a branch-site likelihood ratio test³⁶ (Supplementary Information section S10). One of the most striking features of gibbons is their use of brachiation (arboreal locomotion using only the arms). We uncovered evidence related to traits possibly associated with this adaptation such as the gibbon’s longer arms, more powerful shoulder flexors, rotator muscles and elbow flexors³⁷. First, some genes whose functions relate to these anatomical specializations appear to have undergone positive selection in gibbons. They include TBX5 (P value = 0.00015), required for the development of all forelimb elements³⁸; COL1A1 (pro-alpha1 chains of type I collagen) (P value = 3.39 × 10⁻¹¹), the fibril-forming collagen that is the main protein of bones, tendons and teeth³⁹; and CHRNA1 (acetylcholine receptor subunit alpha precursor) (P value = 0.00039), involved in skeletal muscle contraction⁴⁰. These genes have not been identified as positively selected in other primates to date. We also observed that some genes involved in chondrogenesis (SNX19, ID2 and EXT1) were associated with gibARs. Finally, the chondroadherin gene (CHAD)⁴¹ coding for a cartilage matrix protein is specifically duplicated in all gibbon genera (Extended Data Fig. 2).

Our sequencing, assembling and analysis of the gibbon genome has provided numerous insights into the accelerated evolution of the gibbon karyotype and identified genetic signatures related to gibbon biology. First, segmental duplications and repetitive sequences were the best predictors of gibbon–human breakpoints, although we excluded a causal role given the predominance of non-homology-based repair signatures. Furthermore, accelerated rearrangement was confined to large-scale chromosomal events, pointing to a mechanism responsible for causing gross chromosomal changes, rather than global genomic instability. This is in line with our hypothesis that the high rate of chromosomal rearrangements may have been due to LAVA-induced premature transcription termination of chromosome segregation genes. This effect may have occurred at a low enough level to be compatible with life but sufficient to increase the frequency of chromosome segregation errors. The link between erroneous chromosome segregation and increased chromosomal rearrangement has been recently demonstrated by others through in vitro experiments^{25, 26}.

The question remains how such a high number of chromosomal rearrangements could become fixed in such a relatively short time. One possibility is that a combination of geographic isolation and post-mating reproductive barriers accelerated the radiation of the four gibbon genera. Our estimates dated the lineage-splitting event to the Miocene–Pliocene transition, when major changes in the distribution of tropical and subtropical forests were caused by the elevation of the Yunnan plateau and rise in sea levels^{42, 43}. Furthermore, fluctuation in sea levels beginning in the Early Pliocene appears to have brought about cycles of forest fragmentation and amalgamation, leading to alternating range compression and expansion for many mammalian groups⁴⁴.

Together, these results advance our knowledge of the unique traits of the small apes and highlight the complex evolutionary history of these species. Moreover, our analyses of the rearranged gibbon genome help to provide insight into the mechanisms of chromosome evolution as well as uncovering a new source of genome plasticity.

Sanger-based whole-genome sequencing was performed as described for other species. The genome assembly was generated using the ARACHNE genome assembler assisted with alignment data from the human genome (Supplementary Information section S1). The source DNA for the sequencing was derived from a single female (Asia; studbook no. 0098, ISIS no. NLL605) housed at the Virginia Zoo in Norfolk, Virginia. Short-read libraries were constructed at the Oregon Health & Science University (OHSU) following standard Illumina protocols and sequenced on an Illumina HiSeq 2000. Analyses were performed with custom analysis pipelines. See Supplementary Information for additional information about the methods.

1. Mittermeier, R. A., Rylands, A. B. & Wilson, D. E. Handbook of the Mammals of the World Vol. 3 (Lynx Edicions,. (2013)
2. Carbone, L. et al. A high-resolution map of synteny disruptions in gibbon and human genomes. PLoS Genet. 2, e223 (2006)
   • CAS
   • PubMed
   • Article
3. Locke, D. P. et al. Comparative and demographic analysis of orang-utan genomes. Nature 469, 529–533 (2011)
   • CAS
   • ISI
   • PubMed
   • Article
4. Gibbs, R. A. et al. Evolutionary and biomedical insights from the rhesus macaque genome. Science 316, 222–234 (2007)
   • CAS
   • ISI
   • PubMed
   • Article
5. Girirajan, S. et al. Sequencing human–gibbon breakpoints of synteny reveals mosaic new insertions at rearrangement sites. Genome Res. 19, 178–190 (2009)
   • CAS
   • PubMed
   • Article
6. Carbone, L. et al. Evolutionary breakpoints in the gibbon suggest association between cytosine methylation and karyotype evolution. PLoS Genet. 5, e1000538 (2009)
   • CAS
   • PubMed
   • Article
7. Bailey, J. A. & Eichler, E. E. Primate segmental duplications: crucibles of evolution, diversity and disease. Nature Rev. Genet. 7, 552–564 (2006)
   • Article
8. Yan, C. T. et al. IgH class switching and translocations use a robust non-classical end-joining pathway. Nature 449, 478–482 (2007)
   • CAS
   • ISI
   • PubMed
   • Article
9. Hastings, P. J., Ira, G. & Lupski, J. R. A microhomology-mediated break-induced replication model for the origin of human copy number variation. PLoS Genet. 5, e1000327 (2009)
   • CAS
   • PubMed
   • Article
10. Merkenschlager, M. & Odom, D. T. CTCF and cohesin: linking gene regulatory elements with their targets. Cell 152, 1285–1297 (2013)
    • CAS
    • PubMed
    • Article
11. Schwalie, P. C. et al. Co-binding by YY1 identifies the transcriptionally active, highly conserved set of CTCF-bound regions in primate genomes. Genome Biol. 14, R148 (2013)
    • PubMed
    • Article
12. Carbone, L. et al. Centromere remodeling in Hoolock leuconedys (Hylobatidae) by a new transposable element unique to the gibbons. Genome Biol. Evol. 4, 648–658 (2012)
    • CAS
    • PubMed
    • Article
13. Luan, D. D., Korman, M. H., Jakubczak, J. L. & Eickbush, T. H. Reverse transcription of R2Bm RNA is primed by a nick at the chromosomal target site: a mechanism for non-LTR retrotransposition. Cell 72, 595–605 (1993)
    • CAS
    • ISI
    • PubMed
    • Article
14. Huang da W, Sherman, B. T. & Lempicki, R. A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protocols 4, 44–57 (2009)
    • CAS
    • PubMed
    • Article
15. Mostafavi, S., Ray, D., Warde-Farley, D., Grouios, C. & Morris, Q. GeneMANIA: a real-time multiple association network integration algorithm for predicting gene function. Genome Biol. 9 (Suppl. 1). S4 (2008)
    • CAS
    • PubMed
    • Article
16. Kamburov, A., Wierling, C., Lehrach, H. & Herwig, R. ConsensusPathDB—a database for integrating human functional interaction networks. Nucleic Acids Res. 37, D623–D628 (2009)
    • CAS
    • PubMed
    • Article
17. Baker, D. J., Jin, F., Jeganathan, K. B. & van Deursen, J. M. Whole chromosome instability caused by Bub1 insufficiency drives tumorigenesis through tumor suppressor gene loss of heterozygosity. Cancer Cell 16, 475–486 (2009)
    • CAS
    • ISI
    • PubMed
    • Article
18. Samora, C. P. et al. MAP4 and CLASP1 operate as a safety mechanism to maintain a stable spindle position in mitosis. Nature Cell Biol. 13, 1040–1050 (2011)
    • Article
19. Leber, B. et al. Proteins required for centrosome clustering in cancer cells. Sci. Transl. Med. 2, 33ra38 (2010)
    • CAS
    • PubMed
    • Article
20. Schuyler, S. C., Wu, Y. F. & Kuan, V. J. The Mad1–Mad2 balancing act—a damaged spindle checkpoint in chromosome instability and cancer. J. Cell Sci. 125, 4197–4206 (2012)
    • CAS
    • PubMed
    • Article
21. Maia, A. R. et al. Cdk1 and Plk1 mediate a CLASP2 phospho-switch that stabilizes kinetochore-microtubule attachments. J. Cell Biol. 199, 285–301 (2012)
    • CAS
    • PubMed
    • Article
22. Haraguchi, K., Hayashi, T., Jimbo, T., Yamamoto, T. & Akiyama, T. Role of the kinesin-2 family protein, KIF3, during mitosis. J. Biol. Chem. 281, 4094–4099 (2006)
    • CAS
    • PubMed
    • Article
23. Han, J. S., Szak, S. T. & Boeke, J. D. Transcriptional disruption by the L1 retrotransposon and implications for mammalian transcriptomes. Nature 429, 268–274 (2004)
    • CAS
    • ISI
    • PubMed
    • Article
24. Wheelan, S. J., Aizawa, Y., Han, J. S. & Boeke, J. D. Gene-breaking: a new paradigm for human retrotransposon-mediated gene evolution. Genome Res. 15, 1073–1078 (2005)
    • CAS
    • PubMed
    • Article
25. Damert, A. et al. 5′-Transducing SVA retrotransposon groups spread efficiently throughout the human genome. Genome Res. 19, 1992–2008 (2009)
    • CAS
    • ISI
    • PubMed
    • Article
26. Wojtasz, L. et al. Meiotic DNA double-strand breaks and chromosome asynapsis in mice are monitored by distinct HORMAD2-independent and -dependent mechanisms. Genes Dev. 26, 958–973 (2012)
    • CAS
    • PubMed
    • Article
27. Marchani, E. E., Xing, J., Witherspoon, D. J., Jorde, L. B. & Rogers, A. R. Estimating the age of retrotransposon subfamilies using maximum likelihood. Genomics 94, 78–82 (2009)
    • CAS
    • PubMed
    • Article
28. Gronau, I., Hubisz, M. J., Gulko, B., Danko, C. G. & Siepel, A. Bayesian inference of ancient human demography from individual genome sequences. Nature Genet. 43, 1031–1034 (2011)
    • Article
29. Wall, J. D. et al. Incomplete lineage sorting is common in extant gibbon genera. PLoS ONE 8, e53682 (2013)
    • CAS
    • PubMed
    • Article
30. Drummond, A. J. & Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7, 214 (2007)
    • CAS
    • PubMed
    • Article
31. Durand, E. Y., Patterson, N., Reich, D. & Slatkin, M. Testing for ancient admixture between closely related populations. Mol. Biol. Evol. 28, 2239–2252 (2011)
    • CAS
    • PubMed
    • Article
32. Hirai, H., Hirai, Y., Domae, H. & Kirihara, Y. A most distant intergeneric hybrid offspring (Larcon) of lesser apes, Nomascus leucogenys and Hylobates lar. Hum. Genet. 122, 477–483 (2007)
    • PubMed
    • Article
33. Li, H. & Durbin, R. Inference of human population history from individual whole-genome sequences. Nature 475, 493–496 (2011)
    • CAS
    • ISI
    • PubMed
    • Article
34. Prabhakar, S. et al. Human-specific gain of function in a developmental enhancer. Science 321, 1346–1350 (2008)
    • CAS
    • ISI
    • PubMed
    • Article
35. Pollard, K. S. et al. An RNA gene expressed during cortical development evolved rapidly in humans. Nature 443, 167–172 (2006)
    • CAS
    • ISI
    • PubMed
    • Article
36. Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591 (2007)
    • CAS
    • ISI
    • PubMed
    • Article
37. Michilsens, F., Vereecke, E. E., D’Août, K. & Aerts, P. Functional anatomy of the gibbon forelimb: adaptations to a brachiating lifestyle. J. Anat. 215, 335–354 (2009)
    • PubMed
    • Article
38. Browne, M. L. et al. Evaluation of genes involved in limb development, angiogenesis, and coagulation as risk factors for congenital limb deficiencies. Am. J. Med. Genet. A. 158A, 2463–2472 (2012)
39. Marini, J. C. et al. Consortium for osteogenesis imperfecta mutations in the helical domain of type I collagen: regions rich in lethal mutations align with collagen binding sites for integrins and proteoglycans. Hum. Mutat. 28, 209–221 (2007)
    • CAS
    • ISI
    • PubMed
    • Article
40. Masuda, A. et al. hnRNP H enhances skipping of a nonfunctional exon P3A in CHRNA1 and a mutation disrupting its binding causes congenital myasthenic syndrome. Hum. Mol. Genet. 17, 4022–4035 (2008)
    • CAS
    • PubMed
    • Article
41. Hessle, L. et al. The skeletal phenotype of chondroadherin deficient mice. PLoS ONE 8, e63080 (2013)
    • CAS
    • PubMed
    • Article
42. Cane, M. A. & Molnar, P. Closing of the Indonesian seaway as a precursor to east African aridification around 3–4 million years ago. Nature 411, 157–162 (2001)
    • CAS
    • ISI
    • PubMed
    • Article
43. Xu J.-X, Ferguson D. K, Li, C.-S. & Wang Y.-F Late Miocene vegetation and climate of the Lühe region in Yunnan, southwestern China. Rev. Palaeobot. Palynol. 148, 36–59 (2008)
    • Article
44. Woodruff D. S & Turner L. M The Indochinese–Sundaic zoogeographic transition: a description and analysis of terrestrial mammal species distributions. J. Biogeogr. 36, 803–821 (2009)
    • Article
45. Harvey, P. H., Martin, R. D. & Clutton-Brock, T. H. in Primate Societies (eds Smuts B. B., et al.) Life histories in comparative perspective. 181–196 (Chicago Univ. Press, 1987)
46. Kim, S. K. et al. Patterns of genetic variation within and between Gibbon species. Mol. Biol. Evol. 28, 2211–2218 (2011)
    • CAS
    • PubMed
    • Article

Download references

The gibbon genome project was funded by the National Human Genome Research Institute (NHGRI) including grants U54 HG003273 (R.A.G.) and U54 HG003079 (R.K.W.) with further support from National Institutes of Health NIH/NIAAA P30 AA019355 and NIH/NCRR P51 RR000163 (L.C.), R01_HG005226 (J.D.W., M.F.H.), NIH P30CA006973 (S.J.W.), a fellowship from the National Library of Medicine Biomedical Informatics Research Training Program (N.H.L.), R01 GM59290 (M.A.B.) and U41 HG007497-01 (M.A.B, M.K.K.), R01 MH081203 (J.M.S.), HG002385 (E.E.E.), National Science Foundation (NSF) CNS-1126739 (B.U., M.A.B., M.K.K.) and DBI-0845494 (M.W.H.), PRIN 2012 (M.R.), Futuro in ricerca 2010 RBFR103CE3 (M.V.), ERC Starting Grant (260372) and MICINN (Spain) BFU2011-28549 (T.M.-B.), grant of the Ministry of National Education, CNCS – UEFISCDI, project number PN-II-ID-PCE-2012-4-0090 (A.D.), grant of the Deutsche Forschungsgemeinschaft SCHU1014/8-1 (G.G.S.), ERC Starting and Advanced Grant and EMBO Young Investigator Award (Z.I., N.V.F.), ERC Starting Grant and EMBO Young Investigator Award (D.T.O.), Commonwealth Scholarship Commission (M.C.W.). E.E.E. is an investigator of the Howard Hughes Medical Institute. We acknowledge the contributions of the staff of the HGSC, including the operations team: H. Dinh, S. Jhangiani V. Korchina, C. Kovar; the library team: K. Blankenburg, L. Pu, S. Vattathil; the assembly team: D. Rio-Deiros, H. Jiang; the submissions team: M. Batterton, D. Kalra, K. Wilczek-Boney, W. Hale, G. Fowler, J. Zhang; the quality control team: P. Aqrawi, S. Gross, V. Joshi, J. Santibanez; and the sequence production team: U. Anosike, C. Babu, D. Bandaranaike, B. Beltran, D. Berhane-Mersha, C. Bickham, T. Bolden, M. Dao, M. Davila, L. Davy-Carroll, S. Denson, P. Fernando, C. Francis, R. Garcia III, B. Hollins, B. Johnson, J. Jones, J. Kalu, N. Khan, B. Leal, F. Legall III, Y. Liu, J. Lopez, R. Mata, M. Obregon, C. Onwere, A. Parra, Y. Perez, A. Perez, C. Pham, J. Quiroz, S. Ruiz, M. Scheel, D. Simmons, I. Sisson, J. Tisius, G. Toledanes, R. Varghese, V. Vee, D. Walker, C. White, A. Williams, R. Wright, T. Attaway, T. Garrett, C. Mercado, N. Ngyen, H. Paul and Z. Trejos. We thank Z. Ivics for providing some of the reagents. We additionally acknowledge the Production Sequencing Group at The Genome Institute. Wellcome Trust (grant numbers WT095908 and WT098051), NHGRI (U41HG007234) and European Molecular Biology Laboratory. For the production of next-generation sequences, we acknowledge the Massively Parallel Sequencing Shared Resources (MPSSR) at OHSU, the National Center of Genomic Analyses (CNAG) (Barcelona, Spain), the University of Arizona Genetics Core (UAGC), and the UCSF sequencing core. We also acknowledge the Louisiana Optical Network Institute (LONI). We thank the Gibbon Conservation Center and the Fort Wayne Children’s Zoo for providing the gibbon samples. The MAKER annotation pipeline is supported by NSF IOS-1126998.We thank T. Brown for proofreading and editing the manuscript.

The N. leucogenys WGS project has been deposited in GenBank under the project accession ADFV00000000.1. All short-read data have been deposited into the Short Read Archive (http://www.ncbi.nlm.nih.gov/sra) under the accession number SRP043117. Resources for exploring the gibbon genome are available at UCSC (http://genome.ucsc.edu), Ensembl (http://ensembl.org), NCBI (http://ncbi.nlm.nih.gov), and the Baylor College of Medicine Human Genome Sequencing Center (https://www.hgsc.bcm.edu/non-human-primates/gibbon-genome-project). This paper is dedicated to the memory of Alan R. Mootnick (1951–2011).