the gibbon genome contains all previously described classes of transposable 
 elements that are mostly shared with the 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 .
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 ( fdr=0.031 , p - value=0.001 ) ( supplementary information s7 
 and file 6 ) ( edf 3 ) .
additional analyses with meta - pathway 
 database tools refined this enrichment to pathways related to 
 chromosome segregation , including 
establishment of sister chromatid 
 cohesion and  mitotic metaphase and anaphase ( table st7.3 ) .
genes with lava 
 insertions include proteins that function as checkpoints for cell division and for 
 spindle integrity / architecture ( e.g. , map4 , cep164 , bub1b ) , 
 participate in kinetochore assembly and attachment to the spindle ( e.g. , mad1l1 , 
 clasp2 ) , and play a role in chromosome segregation during 
 cell division ( e.g. , kifap3 , kif27 ) ( edt 1 ) .
intragenic lava insertions were skewed toward introns ( pearson chi - squared , 
 p=0.0001 ) and were less frequent than expected when within < 1 kbp of the 
 nearest exon junction ( edf 3 ) .
the majority 
 ( 74% ) of intronic lava elements were found in the antisense orientation .
we 
 hypothesized that intronic antisense lava insertions may cause early transcription 
 termination ( ett ) by providing a polyadenylation site in antisense orientation , as 
 previously described for l1 elements ( edf 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 .
luciferase 
 activity exceeding background level by ~50% was observed from the lava_f 
 reporter construct ( fig .
further , 3 rapid 
 amplification of cdna ends ( race ) experiments confirmed that the transcription 
 termination site had been supplied from the lava element ( edf 3 ) .
thus antisense intronic lava insertions can cause ett 
 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 ett in vivo by analyzing 
 rna - seq data generated for asia ( 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 presence of a poly(a ) tail was preventing full - length 
 alignment .
this analysis revealed that elements from a variety of sub - families have 
 the potential to cause ett , including those identified for lava elements inserted in 
 the microtubule cytoskeleton genes ( e.g. b2r2 , c4b , b1r2 ) ( edt 1 ) .
of note , we observed that ett occurred at relatively low 
 levels as we identified a significant number of read pairs indicative of normal 
 transcription and splicing for lava - terminated genes ( table st7.5 ) .
this is to be 
 expected , as full inactivation of many of these genes would be incompatible with 
 life . on the other hand , as alternative splicing and rna - pol ii transcript 
 termination/ polyadenylation are tightly coupled processes , lava - mediated ett could 
 also act by differently affecting distinct isoforms and/or influence the ratio 
 between isoforms . finally , lava insertions may also impact 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 s7 , table 
 st7 .
since 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 ( edf 4 ) 
 and found that most of the insertions ( 15/23 ) were shared by the four gibbon genera 
 ( supplementary file 6 ) . 

eleven of the genes match the structural requirements for ett and five of them are 
 also shared .
these genes include map4 , involved in spindle 
 architecture , and cep164 , a g2/m checkpoint whose inactivation 
 results in an aberrant spindle during cell division ( edt 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 presence / absence of 200 lava loci from among the evolutionarily youngest 
 elements in each subfamily ( edf 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 s7 ; table 
 st7.1 - 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 mya and ~17 mya , respectively ( 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 mya ( 95% ci : 15.9 - 17.6 
 mya ) assuming a split time with macaque of 29 mya ( supplementary information 
 s4 ) .
hence , the lava element likely 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 ~15x )
wgs short read data for two 
 individuals from each of the four genera , including two different 
 hylobates species ( h. moloch and h. 
 pileatus ) ( table 
 st2.1 - 2 ) . while 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 ( fig .
neighbor joining trees constructed from pairwise sequence divergence , 
 k , across ~11,000 genic ( 200 bp ) and ~12,000 
 non - genic ( 1 kb ) autosomal loci supported a supermatrix sequence topology of 
 ( ( ( siamang ( ssy ) , hoolock ( hle ) ) , 
 nomascus ( nle ) ) , ( h. pileatus ( hpl ) ) , 
 h. moloch ( hmo ) ) ( fig . 

4a ) , though 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 kbp sliding windows .
however all 15 possible rooted topologies 
 for the four genera were observed at considerable frequencies ( edf 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 employed a novel coalescent - based abc 
 methodology using the autosomal nongenic and genic loci ( veeramah et al . 

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 .
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 .
assuming an overall autosomal mutation rate of 1 x 
 10/site / year , we placed the beginning of the speciation 
 process at ~5 mya under both models , with the two hylobates 
 species diverging ~1.5 mya .
consistent with the abc analysis , ssy and hle share the largest number of 
 alleles across the whole genome ( 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 ( ne ) 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 ( fig st8.2 ) , suggests a larger 
 long - term ne 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 
 ne during the late pleistocene era ( 500 - 100 thousand years ago ( kya ) 
 followed by a subsequent decrease in ne 100 - 50 kya ( fig .
it is important to point out that 
 fluctuation in ne 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=81e-06 ; odds 
 ratio of 2.74 ( 1.794.07 , 95% ci ) ) .
consistent with this finding , a go 
 enrichment test for genes within + /100 kbp of each gibar ( in comparison with 
 background genes ) revealed enrichment for the  chromosome 
 organization category ( benjamini - hochberg fdr < 5% ) ( edf 6 ) . given evidence of functional roles 
 gathered for human accelerated regions , we speculate that the gibars may create functional elements 
 ( e.g. , enhancer , protein - binding domains ) to modulate the transcriptional effect of 
 local lava insertions ( supplementary information s12 and file 9 ) .
we assessed the potential presence of positive selection in 13,638 human 
 genes with one - to - one orthologs in gibbon using a branch - site likelihood ratio 
 test ( supplementary information 
 s10 ) .
one of the most striking features of gibbons is their use of 
 brachiation ( i.e. , 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.39e-11 ) , the fibril - forming 
 collagen 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 ( edf 
 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 , sds 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 . of note
, the link between erroneous 
 chromosome segregation and increased chromosomal rearrangement has been recently 
 demonstrated by others through in vitro experiments .
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 .
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 shattered gibbon genome helped gain insight into the 
 mechanisms of chromosome evolution and uncovered a novel 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 
 s1 ) .
the source dna for the sequencing was derived from a single female 
 ( asia ; studbook no .
short - read libraries were constructed at the oregon health & science 
 university ( ohsu ) following standard illumina protocols and sequenced on an illumina 
 hiseq 2000 .
a , the table compares the gibbon assembly statistics to 
 those of other primates sequenced with a similar strategy .
b , 
 the plot represents the percentage of the 10,734 single - copy gene hmms 
 ( hidden markov models ) for which just one gene ( blue ) is found in the 
 different mammalian genomes in ensembl 70 .
the missing hmms ( cyan ) either do not match any protein or the 
 score is within the range of what can be expected for unrelated proteins . 
 the remaining category ( green ) represents hmms for which the best matching 
 gene scores better than unrelated proteins but not as well as expected .
a , the image shows a representative gibbon - only wssd 
 ( whole - genome shotgun sequence detection ) call by sanger read depth .
the 
 duplication identified in this case overlaps with the gene 
 chad that codes for a cartilage matrix protein . 

b , examples of fish hybridizations on gibbon metaphases 
 using duplicated human fosmid clones that were identified by the ( wgs ) 
 detection strategy ( red signals ) .
c ) intrachromosomal tandem 
 duplication confirmed using cohybridization with a single control probe 
 ( blue signals ) .
c , megabases of lineage - specific and shared 
 duplications for primates based on grchr37 read depth analysis . 

moreover , insertions are significantly enriched in introns and depleted in 
 exons , most likely as a result of selection against insertions in exons . 

b , schematic representation of the mechanism through which 
 lava intronic insertions in anti - sense orientation might cause early 
 termination of transcription : a ) truncated transcript ; b ) normal transcript 
 ( pa = polyadenylation site ) .
c , we calculated the distance to the 
 nearest exon for each intronic lava and compared this to what would be 
 expected for random insertions ( i.e. , background ) .
we found fewer insertions 
 than expected by chance within 1 kbp of the nearest exon .
lava_f 3 tsd is highlighted by dark background ; the major 
 antisense lava_f polyadenylation signal ( maps ) is in red .
a , screenshots from the integrative genomics viewer 
 ( igv ) browser for loci map4 , rabgap1 , and 
 bbs9 .
each column shows portions of the igv 
 visualization of a lava insertion locus identified in nleu1.0 and its 
 flanking sequence .
read pairs are colored in red when their insert size is 
 larger - than - expected , indicating the presence of a lava insertion . 
 map4 is a shared lava insertion , while 
 rabgap1 and bbs9 are 
 nomascus - specific .
b , lava elements 
 containing at least 300 bp of the la section of lava elements were selected 
 and reanalyzed using repeatmasker to determine subfamily affiliation and 
 divergence from the consensus sequence .
lava elements are grouped based upon 
 their subfamily affiliations ( see legend top right ) .
the x - axis shows the 
 percent divergence from the respective consensus sequence , and the y - axis 
 shows the number of elements with a certain percent divergence from the 
 consensus sequence . a , neighbor - joining trees for gibbons using non - genic 
 loci .
b , upgma trees for 100 kbp nonoverlapping sliding windows 
 moving along the gibbon genome reporting the top 15 topologies ( see also 
 supplementary table 
 st8.3 ) .
different 
 sequence types are shown on the x - axis , and the y - axis displays the fraction 
 of gibars and candidate regions annotated to the respective class .
gibars 
 are significantly enriched in intergenic regions ( p = 4.7e-6 ) and 
 significantly depleted in exons ( p = 7.3e-6 ) .
p - values for each class were 
 calculated with the fisher 's exact test .
introns are comparably prevalent in 
 candidates and gibars , while in utr and flanking region counts are too low 
 to draw meaningful conclusions ( data not shown ) .
b , treemap 
 from revigo for goslim biological process terms with a benjamini - hochberg 
 fdr of 5% .
each rectangle is a cluster representative ; larger rectangles 
 represent  superclusters including loosely related terms .
microtubule cytoskeleton go category 
 with lava insertions genes highlighted in gray carry lava insertions that are 
 shared , antisense , and carry a perfect antisense polyadenylation 
 site .
 gibbons are small arboreal apes that display an accelerated rate of 
 evolutionary chromosomal rearrangement and occupy a key node in the primate 
 phylogeny between old world monkeys and great apes . here 
 we present the assembly 
 and analysis of a northern white - cheeked gibbon ( nomascus 
 leucogenys ) genome . 
 we describe the propensity for a 
 gibbon - specific retrotransposon ( lava ) to insert into chromosome segregation 
 genes and alter transcription by providing a premature termination site , 
 suggesting a possible molecular mechanism for the genome plasticity of the 
 gibbon lineage . 
 we further show that the gibbon genera 
 ( nomascus , hylobates , 
 hoolock and symphalangus ) experienced a 
 near - instantaneous radiation ~5 million years ago , coincident with major 
 geographical changes in southeast asia that caused cycles of habitat compression 
 and expansion . 
 finally , we identify signatures of positive selection in genes 
 important for forelimb development ( tbx5 ) and connective 
 tissues ( col1a1 ) that may have been involved in the adaptation 
 of gibbons to their arboreal habitat . 
LAVA insertions in the gibbon genome
The complex evolutionary history of gibbons
Functional sequence evolution
DISCUSSION
METHODS
Extended Data
Supplementary Material