Published online 16 November 2020 NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 1 doi: 10.1093/nargab/lqaa089 Characterization of nuclear mitochondrial insertions in the whole genomes of primates Gargi Dayama 1 , Weichen Zhou 1 , Javier Prado-Martinez 2 , Tomas Marques-Bonet3,4,5,6 and Ryan E. Mills 1,7,* 1 Department of Computational Medicine and Bioinformatics, University of Michigan Medical School, Ann Arbor, MI 48109, USA, 2 Wellcome Sanger Institute, Hinxton, CB10 ISA, Cambridge, UK, 3 Institute of Evolutionary Biology (UPF-CSIC), PRBB, Dr. Aiguader 88, 08003 Barcelona, Spain, 4 Catalan Institution of Research and Advanced Studies (ICREA), Passeig de Lluı́s Companys, 23, 08010, Barcelona, Spain, 5 CNAG-CRG, Centre for Genomic Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, 08028 Barcelona, Spain, 6 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/ Columnes s/n, 08193 Cerdanyola del Vallès, Barcelona, Spain and 7 Department of Human Genetics, University of Michigan Medical School, Ann Arbor, MI 48109, USA Received February 21, 2020; Revised May 04, 2020; Editorial Decision October 12, 2020; Accepted October 15, 2020 ABSTRACT INTRODUCTION The transfer and integration of whole and partial mi- Polymorphic Nuclear Mitochondrial Insertions (NumtS) tochondrial genomes into the nuclear genomes of are fragments of mitochondrial DNA (mtDNA) that have eukaryotes is an ongoing process that has facili- been transferred and integrated into the nuclear genome of tated the transfer of genes and contributed to the an organism. NumtS vary in size from smaller fragments to evolution of various cellular pathways. Many previ- full-length mitochondrial integrations and can share vary- ing degrees of homology with the genome of their parent ous studies have explored the impact of these inser- mitochondria, depending on the age and time of insertion. tions, referred to as NumtS, but have focused pri- Once these pieces of mtDNA have inserted themselves, their marily on older events that have become fixed and mutation rate decreases by an order of magnitude to the are therefore present in all individual genomes for background rate of the organism’s nuclear genome, essen- a given species. We previously developed an ap- tially fossilizing the fragment and providing a snapshot of proach to identify novel Numt polymorphisms from the ancestral mitochondria from when it was inserted (1– next-generation sequence data and applied it to thou- 3). This information can be very beneficial in phylogenetic sands of human genomes. Here, we extend this anal- analysis (4,5), and indeed such inserted sequences have been ysis to 79 individuals of other great ape species in- identified in previous studies exploring their prevalence in cluding chimpanzee, bonobo, gorilla, orang-utan and humans (5–7), chimpanzees (5,8–9) and other organisms (5,10). As with other forms of genetic variation, NumtS can also an old world monkey, macaque. We show that be both fixed within all genomes of a species as well as poly- recent Numt insertions are prevalent in each species morphic between individuals. While the majority of prior though at different apparent rates, with chimpanzees studies have focused on older, fixed insertions, recent efforts exhibiting a significant increase in both polymorphic have begun to explore the impact of segregating Numt alle- and fixed Numt sequences as compared to other les within human populations (2,7). great apes. We further assessed positional effects It has been shown that all parts of the mtDNA are able to in each species in terms of evolutionary time and transfer into the nuclear genome, though relative rates for rate of insertion and identified putative hotspots on individual regions are still unknown (7). Prior studies have chromosome 5 for Numt integration, providing in- reported a deficit of NumtS from the hypervariable (HV) sight into both recent polymorphic and older fixed regions of mtDNA, specifically HV2 within the mitochon- reference NumtS in great apes in comparison to hu- drial control region (MCRF ) (1). In particular, a negative correlation was reported between the prevalence of NumtS man events. and the proportion of variability of the site, thus increas- ing the difficulty of detecting events in the nuclear genome, which exhibits a slower rate of mutation (1). More recent * To whom correspondence should be addressed. Tel: +1 734 647 9628; Fax: +1 734 615 6553; Email: remills@umich.edu  C The Author(s) 2020. Published by Oxford University Press on behalf of NAR Genomics and Bioinformatics. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. 2 NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 work focusing specifically on Numt polymorphisms, how- The reference versions used were as follows; for chimpanzee ever, reported a slight enrichment of the D-loop for the ‘panTro4’, for bonobo ‘panPan1’ (assembled by (24)), for polymorphic events (7). Likewise, the patterns of insertional gorilla ‘gorGor4’, for orang-utan ‘ponAbe2’ and for rhesus preferences within chromosomal DNA have not been well ‘rheMac3’, and comparisons were made to previously iden- defined. Earlier work reported that younger insertions are tified NumtS in humans using the GRCh37 version of the more prevalent within intronic regions, while older NumtS human reference genome. tend to be intergenic (11). More recent analyses have sug- gested that they may insert non-randomly, with a tendency Numt identification and assembly to insert in regions significantly deficit of transposable el- ements (12,13), while also seeming to have a local affinity To reduce misalignments from variability in mtDNA se- toward low GC regions and A+T dinucleotides in humans quences within individuals of each species, modified ref- (7). An abundance of human specific NumtS on chromo- erences were generated for each sample by concatenating some Y was reported in one study, though these appear to individual mtDNA to its species-specific reference. The Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 have arisen through the duplication of existing NumtS and sample-specific mtDNA were assembled using mtArchitect, were likely not driven through new insertions (1). However, a method previously developed (25) which uses iterative others have reported no pattern or specific insertion sites mapping whereby sequence reads initially aligning to the (14). reference mtDNA for each individual were isolated and as- Previous studies have also found discrepancies in the rate sembled into a single mtDNA contig. These sample-specific of insertion for NumtS across different species. For exam- mtDNA sequences were then inserted into each respective ple, some groups have observed a consistent rate of insertion reference genomes as an additional contig and sequence throughout the evolution of great apes (1,6,15), whereas reads were realigned and processed as described above. The others reported variability between species (9,16). Indeed, updated BAM files were used for polymorphic Numt dis- earlier observations reported a higher rate of insertions in covery using dinumt (7), which utilizes aberrant/discordant both gorillas (17) and chimpanzees (12) as compared to hu- reads (where one read maps to the mtDNA and the other mans. Additional research has reported a burst of inser- elsewhere on the nuclear genome), read orientation and tion at the divergence of old-world and new-world mon- various other filters to define insertion breakpoints. Fixed keys (5,15,18–19). Further, an interesting stratification in NumtS in each species-specific reference genome were re- the evolution of fixed Numt insertions was also reported be- analyzed using the method previously described by (26) to tween primate and non-primate species (5). Overall, there is generate a consistent and updated set across organisms. much ambiguity regarding the interspecies rate of insertion, A three step process was utilized for the de novo assem- copy number and length and these discrepancies are likely bly of non-reference NumtS detected by our approach: (i) due in part to the predominant use of older, species-specific Clustering: the aberrant and soft-clipped read pairs from NumtS as well as limited sample size in the prior analysis. the flanking region of non-reference Numt insertion sites Here, we present a holistic analysis of NumtS using both were extracted and collated from all the samples for which older, fixed insertions as well as new polymorphic events with the insertion was detected. (ii) Assembly: where pos- detected by our previously developed methodology (7). We sible, these polymorphic NumtS were then assembled into derive the prevalence and insertion rates of NumtS across contigs using CAP3 (27). (iii) Annotation: the resulting con- species and explore positional effects and sequence con- tigs were subsequently aligned against their respective ref- text. We further de novo assemble a subset of Numt in- erence mtDNA genome to annotate the specific Numt frag- sertions and analyze their evolutionary distance from con- ment. structed ancestral mtDNA sequences for the various spe- ciation events. This represents the first exploration and as- Orthogonal validation of NumtS using long-read sequences sessment of polymorphic events and their significance in the evolution of the primates. To assess the accuracy of putative Numt insertion sites ob- tained from the short-read data sets, we made use of ex- MATERIALS AND METHODS isting long-read whole genome sequences that were avail- able for a subset of the species we interrogated. Long- Data acquisition and processing read sequencing data from Pacific Biosciences (PacBio) Whole genome, paired-end sequencing data were ob- Single Molecule, Real Time (SMRT) technology were tained in SRA format from the short read archive (https: obtained from NCBI under the project accession num- //www.ncbi.nlm.nih.gov/sra) for various non-human pri- bers PRJNA369439 (chimpanzee, orang-utan) and PR- mates: chimpanzee––Pan troglodytes; bonobo––Pan pani- JEB10880 (gorilla) (28,29). We used pbmm2 (https://github. cus; gorilla–– Gorilla gorilla and orang-utan––Pongo abelii com/PacificBiosciences/pbmm2) to align the PacBio raw and Pongo pygmaeus, were obtained from the Great Ape reads of three primate genomes to their relative references. Genomes Project (20) and rhesus––Macaca mulatta data We customized an approach for identifying mobile ele- were obtained elsewhere as an outgroup (21). The SRA ment insertions, PALMER (30) (https://github.com/mills- files (SRP018689, ERP002376) were converted to fastq lab/PALMER), to allow the discovery of NumtS in PacBio format using the SRA Toolkit (version 2.3.5-2 released aligned data. We required at least four sub-reads at a April 2014) and then subsequently aligned to its respec- given locus with NumtS-related sequences and performed tive species-specific reference and processed using BWA (22) additional error correction (sequencing error <1%) using and GATK (23) to generate alignment files in BAM format. CANU (31) to correct PacBio sub-reads that support the NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 3 insertions. These were subsequently re-aligned locally to the Enrichment analysis predicted insertion sites using pbmm2, followed by another We tested for any potential preferential bias of MT genes of PALMER to better define the sequences and integration being inserted as NumtS. To assess this, we performed a sites of the detected NumtS. We then examined the inter- permutation analysis by generating 1000 sets of random in- section between calls from PacBio and standard short reads, sertion coordinates based on the list of derived polymor- allowing an extension of ±500 bp for insertion site position. phic and reference Numt fragment sizes. The frequency of complete gene insertions was calculated for all MT genes Frequency and characterization of insertion site across all the for each species. As some of the primate species were lack- species ing the MT gene annotation, we defined the genes us- To investigate whether insert-size or coverage had an ef- ing human mtDNA reference (GRCh37) for each primate fect on the number of insertions we observed, we derived species through the lift-over tool and chain files from UCSC and compared the frequency of Numt insertions relative Genome Browser (32). A Z-score was calculated for each Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 to their respective reference for every non-human primate gene to determine trends of MT gene enrichment or deple- sample. The defined Numt breakpoints (site of insertion in tion across all human and non-human primates as follows: the genome relative to the reference) were then unified to Z − score = (x − μ) /σ the human reference (GRCh37) using the lift-over tool and chain files from UCSC Genome Browser (32) to compare where x is the true values of total gene insertions, ␮ is the the sites of insertion across all primates including humans mean of permutation count and ␴ is the standard deviation (7) for both reference and polymorphic events. They were of the permutation count. then visualized using a circos (33) plot. A similar strategy was used to assess the subset of assem- bled NumtS. Their aligned coordinates were unified to the Analysis of putative hotspots of Numt insertion human reference (GRCh37) as described above and a cir- The entire nuclear genome was divided into 25 Mbp bins cos plot was generated to display integration preferences and the frequency for Numt insertions was calculated and for mtDNA fragments into each species. The human ref- compared across the homologous positions in each species erence mtDNA (GRCh37) was binned into genes and the to determine any hotspots above a permutation threshold average rate of insertion for each bin was calculated across that defines 5% FDR (false discovery rate) as follows: all species.   n We also utilized these liftover coordinates as unified to the   2 human reference (GRCh37) from all other species to define  xi − xi ancestral NumtS positions. Polymorphic and fixed NumtS i =1 from each species were combined to conduct a pairwise assessment of shared events across all species. Ancestral Where x− is the mean frequency of NumtS across all the bins NumtS at each node of divergence were determined by com- for a given species i and x is the frequency of NumtS in each bining all species that diverged at the given node and screen- bin, equated across n number of species. A total of 25 Mbp ing for shared events relative to humans. Furthermore, we was determined to be a lower bound for window size given also examined shared events between all species relative to the sparsity of reference and polymorphic NumtS in each each other in a pairwise fashion, omitting bonobo to avoid species. any bias introduced due to their unusually low number of shared events. Software availability Characterization and rate of insertions for NumtS The software and analysis pipelines used in this analysis can be found at https://github.com/mills-lab/primate numt and To investigate species-specific differences in insertional https://github.com/mills-lab/dinumt. characteristics, we compared the fragment length and GC content for all the assembled polymorphic and reference NumtS (extracted from the reference genome of each RESULTS species). In addition, the rate of insertion was estimated for Discovery and frequency of NumtS across species assembled polymorphic NumtS using the method described previously (7). The ancestral mitochondrial sequences for We identified a total of 409 polymorphic NumtS in 79 sam- each primate genus were obtained from ENSEMBL Com- ples across 5 different non-human primate species (Addi- para release 77 based on 8 primates EPO (Enredo-Pecan- tional file 1). Chimpanzee was observed to encompass the Ortheus) pipeline (34). A profile of nucleotide changes was largest number of polymorphic Numt insertions with a to- generated between these ancestral mtDNA and the consen- tal of 187 events discovered in 24 samples across four dif- sus mtDNA for each species. The assembled non-reference ferent subgroups at an observed rate of 20.8 average per NumtS were then annotated for these specific positions to sample. Pan troglodytes troglodytes, P. troglodytes ellioti, calculate an allele ratio. Finally, after post-correction of any Pan troglodytes schweinfurthii exhibited an average of 22 mismatches, these allele ratios were normalized by the evo- events per sample, though P. troglodytes verus showed a lutionary distance since their specific node divergence, thus much lower rate of 10.2. Macaque had the highest over- calculating the approximate age of each assembled inser- all average per sample rate of 21.2 insertions per sample tion. and a total of 56 unique events called in five samples, while 4 NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 orang-utans showed the lowest at 12.9 average per sample, to construct assembled contigs for 54% of our predictions with 43 total unique events in 10 samples from two differ- across each species. We assessed the accuracy of our as- ent species, though with differing rates between Pongo pyg- sembled sequences through comparison with high-quality maeus (17.2 avg/sample) and Pongo abelii (8.6 avg/sample). NumtS we had previously sequenced and observed a very The bonobo individuals in aggregate had 67 total inser- high accuracy at concordant sites (99.9% average iden- tions with an average of 15.46 per sample and 27 samples tity). Using this approach, we were able to recapitulate the of gorilla had a total of 55 Numt called in two different underlying insertion sequences for 222 NumtS across all species and three different subspecies with an overall 14.6 species (Additional Files 4–8). The majority of NumtS that avg/sample. Both Gorilla gorilla gorilla and Gorilla gorilla we were able to assemble were found in multiple samples. diehli had approximately 15 events per sample, but the sub- Longer/full-length NumtS were more difficult to assemble species Gorilla beringei graueri only had 2.5 avg/sample, due to the circular nature of mtDNA as well as the lim- similar to human polymorphic events (Table 1). We reex- itations of using short reads, though we were able to de- amined the reference genomes for each species as many rive the estimated length and mitochondrial region of origin Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 had been updated since previous studies on NumtS had for some longer insertions based on the partial assembly of been conducted and cataloged a total of 4274 reference each breakpoint. (fixed) NumtS in great apes and macaque (chimpanzee: 974; bonobo: 809; gorilla: 657, orang-utan: 963 and macaque: Validation of Numt calls using long-read sequencing 871) (Additional File 2). We next assessed the frequency of sample-specific NumtS PacBio long-read sequencing can generate longer (10 Kbp and observed a common trend in chimpanzee, bonobo and or higher on average) contiguous DNA sequences than macaque for a higher rate of sample-specific, relatively re- standard short-read sequencing technologies and can iden- cent insertions compared to other species. In contrast, most tify full-length insertions that fall completely within an indi- insertions in gorilla and orang-utan were shared among vidually sequenced molecule (30). We made use of three pri- other members of the species (Figure 1). We also discov- mate genomes (one chimpanzee, one gorilla and one orang- ered several events that were found in almost all samples utan) sequenced to >65× by PacBio to provide orthogonal suggesting they may represent older insertions. Overall, a validation for our call sets of NumtS (28). Although these higher frequency of polymorphic NumtS was observed in sequences were obtained from different individual samples chimpanzee subspecies, with the exception of P. troglodytes and were not present in the short-read cohort we analyzed, verus exhibiting a frequency similar to the other primates. we postulated that it would still be likely that we would iden- Furthermore, we did not observe any effect of insert-size tify some species-specific polymorphisms that were shared or coverage on the number of insertions (Additional File 3: between the sets that could be used as a basis for compar- Supplementary Figure S1). ison. Indeed, we observed orthogonal evidence for seven NumtS in chimpanzee, four NumtS in orang-utan and one in gorilla (Additional File 9) from PacBio data. The PacBio Spatial organization of Numt insertions into nuclear primate long reads also provided augmented information for these genomes calls, including the specific sequence and segment coordi- We examined the chromosomal locations of Numt inte- nates within the mtDNA as well as insertion orientation gration across homologous regions of each species to ex- (Additional File 9). In addition to verifying the insertion plore shared insertional preferences. We found no indica- loci in the nuclear genome, we further compared our as- tion of preferential bias toward specific chromosomes or re- sembled NumtS to the sequences present in the respective gions for either NumtS present in the respective references PacBio long reads and found a 99.5% average identity con- of each species or polymorphic events (Additional File 3: cordance between the two technologies. Together, the re- Supplementary Figure S2A and B), with the exception of sults showed a high consistency between calls by PacBio chromosome Y. No polymorphic insertions were detected data and our assembled Numt insertions by standard short- even for species in which the Y chromosome was available read technology with regard to insertion loci, segment co- for analysis (humans, chimpanzees and orang-utans). There ordinates within the mtDNA, and the insertion sequences were also several gaps on some chromosomes in various suggesting we obtained a high confidence set of polymor- species that may be due to variability in assembly, partic- phic assembled Numt insertions. ularly around the centromeric regions, potentially resulting in a dearth of Numt detection. Characterization and enrichment of NumtS We previously reported that polymorphic NumtS in hu- De novo assembly of Numt insertions mans exhibited a higher GC% than their parent mito- We next sought to derive the underlying sequences of the chondrial genome (7), and thus we sought to determine Numt insertions. In contrast to our earlier work with hu- whether this observation extended to non-human primates man cell lines, the underlying primate samples from which as well. We calculated the GC% for the subset of NumtS we the whole genome sequences were derived were either un- were able to assemble for each species and, unlike humans, available or very sparse and thus prohibit the direct molec- observed no significant change in any of the polymorphic ular characterization of these events. We thus applied an as- insertions compared to their parent mtDNA with the ex- sembly strategy (see ‘Materials and Methods’ section and ception of those present in chimpanzee (Figure 3). How- Figure 2) around each putative insertion site and were able ever, the GC% for all polymorphic events differed between NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 5 Table 1. Frequency of NumtS Frequency of Numt insertions Species N Coverage Reference Polymorphic Ave/sample Spp Sub Spp Homo sapiens* 999 4-12 767 141 1.5 Pan troglodytes 24 4.4 918 187 20.8 Pan troglodytes troglodytes 4 6.2 67 25.7 Pan troglodytes ellioti 9 3.5 98 21.8 Pan troglodytes verus 5 5.7 41 10.2 Pan troglodytes schweinfurthii 6 5.2 78 23.6 Pan panicus 13 5.9 736 67 15.4 Gorilla 27 6.4 657 55 13.6 Gorilla gorilla gorilla 24 6.8 48 14.4 Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 Gorilla gorilla diehli 1 6.2 16 16 Gorilla beringei graueri 2 6.3 5 2.5 Pongo 10 6.4 923 44 12.9 Pongo abelii 5 6.1 34 8.6 Pongo pygmaeus 5 6.6 40 17.2 Macaca mulata 5 9 872 57 19.4 Total 1078 4873 432 The frequency of Numt insertions in reference genomes and Numt insertions relative to the respective reference, discovered in several (N = number of samples) samples of six different groups (human, chimpanzee, bonobo, gorilla, orang-utan and macaque). Chimpanzee and gorilla are further divided into subspecies. The total polymorphic NumtS for the species are listed under ‘spp’ whereas the frequency for subspecies are listed under ‘sub spp’. Average coverage for each group is noted in the third column. The last column has the average Numt insertions in each group. *Reported previously in (7) Figure 1. Shared and sample-specific Numt frequency. Each bar is representing the frequency of Numt insertion in that sample. The red portion of the bar indicates the number of events that are shared between samples within that species, whereas the blue bar shows the Numt insertions unique to that sample. Ass three-letter code is used for grouping the subspecies (in case of chimpanzee and gorilla) and population (for humans), whereas a two-letter code is used for denoting the species (for bonobo, orang-utan and macaque). species, with those within macaque, gorilla, and bonobo ments for recently inserted polymorphic NumtS. We further exhibiting a lower GC% in contrast to the higher GC% generated a phylogenetic tree using the reference mtDNA observed in human, chimpanzee and orang-utans (Addi- for each species and their respective ancestral mtDNA to tional File 3: Supplementary Figure S3). A similar ob- explore ancestral relationships and found species-specific servation was also noted in average GC% for reference mtDNA to cluster perfectly with their ancestral mtDNA mtDNA in all species. An increase in GC% of polymor- (Additional file 3: Supplementary Figure S4). Additionally, phic events compared to their respective reference mtDNA we also found that the average GC% for macaque, orang- was observed for chimpanzees and humans. This suggests utan and human ancestral mtDNA was lower compared that while polymorphic events mimic the GC% of their to GC% of polymorphic NumtS and reference mtDNA, parent mtDNA in all other species, chimpanzees and hu- whereas it was similar or higher for the rest of the species mans might have a selective bias for higher GC% frag- (Figure 3). 6 NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 and matched our predicted mtDNA coordinates, enabling us to infer that ∼5% of the events were likely full length. In contrast, reference insertions across all species were up to 8 Kbp in size, with 3 (>8 Kbp) outliers only in humans. To get an overview of fixed Numt evolution spanning over millions of years, we characterized the ancestral NumtS at each node of species divergence. A total of 530 shared events between human, chimpanzee and bonobo were ob- served, of which only two human polymorphic NumtS were shared with chimpanzee reference NumtS suggesting they may actually represent human-specific deletion polymor- phisms. At the next node of divergence, we found a total of 426 events shared between human, chimpanzee, bonobo Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 and gorilla, further decreasing to 355 where orang-utan diverged with the rest of primates. Finally, a total of 132 shared NumtS were characterized by the common ancestor between macaque and primates (Additional File 3: Supple- mentary Figure S5A). The human coordinates of all the 530 shared NumtS with each of the five other species, including the two polymorphic shared events (denoted by ‘HP’) can be found in Additional File 10. As expected, we observed the least number of shared NumtS at this node, given the 25MY’s of divergence. Due to the unusually low number of shared events between bonobo and the rest of the species, we decided to omit this group from further analysis to ex- amine the shared events between groups. We observed the largest number of potentially fixed novel insertions occur- ring in chimpanzee and orang-utan (Additional File 3: Sup- plementary Figure S5B). Gorilla exhibited the fewest inser- tions, though this could also be due to its less-well charac- terized reference genome that still consists of several gaps. Finally, we conducted a pair-wise comparison between each species to look at the shared polymorphic and fixed NumtS. We observed that humans and chimpanzees have the most shared NumtS (Additional File 3: Supplementary Figure S5C). In contrast, bonobo showed a very low number of shared events with all other species, which was quite un- expected given the relatively short divergence period from chimpanzees. We also observed several additional shared polymorphic events between various species, suggesting ei- ther a deletion of reference events in the species or mischar- acterization. Characterization of Numt insertional preferences in the nu- Figure 2. Pipeline for Numt assembly. The graphical representation is di- vided into three main parts of the Numt assembly pipeline. (A) Utilizing clear genome BAM files to extract the flanking reads and split reads around the break- We examined insertional preferences for any particular re- points of Numt insertion from all the samples that they were discovered in. (B) These reads are then combined into one file to assemble into con- gion of mtDNA, with a focus on Numt fragments contain- tigs using CAP3. (C) The contigs are then annotated using UCSC BLAT ing complete genes or the D-loop (Figure 4A and B). We feature. compared detected NumtS in the reference sequences of each species to random background models of insertions by converting gene annotations for mitochondrial genes for We next assessed whether there were length differences each species relative to human mtDNA using the liftOver between our polymorphic NumtS and those present in tool (see ‘Materials and Methods’ section) and did not ob- their respective reference sequences (Figure 3). The major- serve any obvious preferential bias toward any region or ity (∼95%) of the assembled polymorphic insertions were genes from the mtDNA relative to the human reference smaller than 1 Kb in length, which is expected given the lim- mtDNA sequence across all the species. We next exam- itations of short-read assembly and is consistent with our ined polymorphic and reference NumtS identified individ- previous results (7). Even though we were not able to fully ually at the species level and observed several regions of assemble these long insertions, we were able to partially as- mtDNA that exhibited enrichment or depletion across mul- semble contigs that crossed over the breakpoints precisely tiple genes (Additional File 3: Supplementary Figures S6– NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 7 Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 Figure 3. Rate of insertion, length of NumtS and GC content. (A) Rate of Numt insertion represented by the red bars divided per MYs for all 6 groups. (B) %GC content in assembled polymorphic Numts compared (represented by blue notched box plot) with %GC content of reference mtDNA (orange bar) and ancestral mtDNA (green bar). (C) Length of assembled polymorphic NumtS represented by the blue bars compared to the length of reference NumtS represented by overlain purple density plot. 8 NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 A B Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 Figure 4. Spatial representation of species-specific insertions relative to the human reference genome. The breakpoints for all the species are lifted over to human reference (GRCh37/hg19). (A) The normalized frequency of each mitochondrial gene across all the reference Numt fragments. A unique color is assigned for each genus:purple for humans, green for chimpanzee, blue for bonobo, red for gorilla, black for orang-utan and orange for macaque. Using the same color scheme figure (B) was generated to illustrate the frequency of polymorphic NumtS. 11). We assessed these regions across all species to identify high abundance of fixed reference events in all six species, any selective biases toward particular MT genes and did not with orang-utans exhibiting the largest number of inser- observe any pan-species MT gene bias for reference NumtS tions and humans the least. Additionally, to understand the (‘Materials and Methods’ section, Figure 5A). However, for evolutionary timeline of these insertions, we looked at how polymorphic events, we found an enrichment of genes lo- many of the shared (ancestral) events were a part of the cated toward the center of mtDNA and depletion closer hotspot on chromosome 5. We observed the same number to the D-loop (Figure 5B). This was consistent across all (n = 6) of shared events between ‘humans, chimpanzees or species except bonobo, which interestingly exhibited the op- bonobo’ and ‘humans, chimpanzees or bonobos and goril- posite, showing a depletion of genes toward the center of las’. This decreased to four shared events when including mtDNA and enrichment of D-loop containing fragments. orang-utan with these species. Only one ancestral shared We next investigated whether there was any enrichment event was observed in the hotspot region of chromosome 5 for polymorphic and reference NumtS within and across between all six species, implying this region began exhibit- the different species The human genome (GRCh37) was seg- ing multiple insertions after the divergence of macaque from mented into bins of 25 Mbps to compare the frequency of the great apes. Numt insertion (lifted over to the human reference) for each Once a piece of mtDNA inserts itself into the nuclear bin across all the groups (human, non-human primates, genome, it typically follows the mutational rate of its host and Macaca mulatta) against a random background. Eight DNA that is an order of magnitude lower than the mtDNA, hotspots with a frequency greater than the threshold value giving a snapshot of when the insertion occurred (1–3). (permutation threshold that defines 5% FDR) were ob- Based on this, we calculated the rate of insertion per mil- served for reference NumtS, one of which on chromosome lion years normalized against the evolutionary time of di- 5 was also identified as the sole hotspot for polymorphic vergence. We observed that most polymorphic insertions insertions (Additional file 3: Supplementary Figure S12). occurred in the past million years for gorilla, bonobo, chim- We further explored this region on chromosome 5 (chr5:75 panzee and macaque, which was similar to what has been 000 000–100 000 000) to determine whether this was spe- observed for humans. orang-utan was the only group that cific to an individual species or common across primates. showed a different rate of insertion, where the insertions We observed that the majority of these insertions were pri- were spread evenly over the past 4 to 5 million years (Figure marily from chimpanzee, bonobo and gorilla, while human, 3). Interestingly, we also found a Numt that inserted ∼15 orang-utan and macaque were lacking any polymorphic in- MYA in the macaque. It is possible however that these in- sertions in this region. However, we did observe an overall sertions are fixed in the population but are not a part of the NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 9 A though we can posit no specific cause for the change in the rates of insertions, we can proffer several speculations. First, this might be due to Mitochondria’s rate of disintegration caused by factors such as oxidative phosphorylation (35) unique to each species. Alternatively, the prevalent mech- anism for Numt insertions is understood to be implicated with the repair of double-strand breaks (DSB) (8), however, there may be other mechanisms yet to be discovered that are specific to a given species. We also observed an absence of polymorphic Numt in- sertions on the Y chromosomes across all species. This is in part due to the bonobo, gorilla, and macaque reference genomes omitting the Y chromosome, however, we also Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 B did not observe any polymorphic insertions in orang-utan, chimpanzee or even humans whose references are more complete. This suggests some selective pressure against these insertions on Y chromosome; alternatively, if the in- sertions were to occur in the gamete stage when only the ovum contains mtDNA and not the sperm, the absence of any insertions on the Y chromosome might be explained. This hypothesis is further supported by the findings shown by (18) that all the recent reference NumtS on the Y chro- mosome were due to the result of duplication rather than novel insertions. We observed a higher GC% of polymorphic NumtS in chimpanzee and humans compared to their respective par- ent mtDNA, which might be due to two reasons; i) If there is an insertional bias toward mtDNA fragments with higher Figure 5. Mitochondrial gene enrichment or depletion across species. GC%, to balance the low GC content of the region where (A) Linear mtDNA representing the enriched and depleted (Z-score ≥ 1 the insertion occurs on nuDNA or ii) the mtDNA fragment and ≤ −1) major genes (genes coding for tRNA excluded) inserted as poly- with lower GC% disintegrating faster than the fragments morphic NumtS in all six species. Unique color and alphabet is assigned that have a higher GC, thus increasing the chances of these for each genus, purple ‘H’ for humans, green ‘C’ for chimpanzee, blue ‘B’ fragments being selected for repairing the non-homologous for bonobo, red ‘G’ for gorilla, black ‘O’ for orang-utan and orange ‘M’ is for macaque. Using the same color scheme figure (B) was generated to breaks. illustrate the significant Z-scores for polymorphic NumtS. We further report an enrichment of both reference and polymorphic Numt insertions on chromosome 5. While for reference events the hotspot was derived by a substantial reference due to weak characterization and gaps in some of number of insertions across all the species (especially orang- these primate species. utan with 50 insertions), polymorphic events were mainly derived by gorilla and bonobo. Combined with the low number of shared NumtS observed at each node of diver- DISCUSSION gence, it suggests that these hotspots were created by in- Advances in sequencing technologies have made it possible sertions that are an ongoing and independent occurrence to conduct studies at a finer scale, allowing the exploration across all of the species. We did observe that human refer- of previously understudied genomic variation in popula- ence (GRCh37) has a slightly lower average GC% at 34% tions across multiple species. Here, we have described one of for 100 bp flanks on either side of Numt breakpoints in this the first forays into the investigation of polymorphic nuclear region. Further, (36) examined structural variations and mitochondrial insertions (NumtS) in multiple samples from transposable elements by dividing the genome into 10 Kbp primates and an old world monkey. We report a shift in the bins and reported several hotspots in this 25 Mb region for frequency of polymorphic Numt insertions with a succes- human (9), chimpanzee (6) and orang-utan (6). Interest- sive increase from orang-utans to chimpanzees, followed by ingly, they did not find any hotspots in rhesus. So far the a decrease in humans by an order of magnitude. The high- only known mechanism for Numt insertions is through the est rate of insertion was observed in Macaca mulatta which repair of DSBs, but exploring this phenomenon could po- diverged from humans approximately 25 MYA. This sug- tentially lead to other mechanistic explanations for Numt gests that at some point after macaque-primate divergence insertions and their continued occurrence. It should be there was a considerable drop in the rate of insertion, as seen noted that the above mentioned observations might also in orang-utans. Indeed, a burst of Numt insertions in the possibly be because some of the primate references are not reference genome (and thus, primarily fixed in the species) well defined and we might be missing either newer NumtS has been previously reported (15) around the same period relative to the reference or calling events that are fixed in the as the divergence of old world monkeys, which is consis- species but missing in the reference. We also note that there tent with their presence as an outgroup in our analysis. Al- is likely an ascertainment bias however due to the preferen- 10 NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 tial assembly of smaller NumtS events from the short-read 5. Calabrese,F.M., Balacco,D.L., Preste,R., Diroma,M.A., Forino,R., sequence data available. Ventura,M. and Attimonelli,M. (2017) NumtS colonization in mammalian genomes. Sci. Rep., 7, 16357. Our study does have some limitations, such as the lack of 6. Mourier,T., Hansen,A.J., Willerslev,E. and Arctander,P. (2001) The direct molecular validation due to the unavailability of the human genome project reveals a continuous transfer of large samples we analyzed. However, we previously established mitochondrial fragments to the nucleus. Mol. Biol. Evol., 18, the high fidelity of our NumtS detection method in hu- 1833–1837. mans (7) and expect a similar accuracy for the non-human 7. Dayama,G., Emery,S.B., Kidd,J.M. and Mills,R.E. (2014) The genomic landscape of polymorphic human nuclear mitochondrial primates assessed here. We also made use of the long-read insertions. Nucleic Acids Res., 42, 12640–12649. PacBio sequences which, while limited, had very good con- 8. Hazkani-Covo,E., Einat,H.-C. and Shay,C. (2008) Numt-mediated cordance with both our predicted insertions as well as their double-strand break repair mitigates deletions during primate assembled sequences. This indicates that our de novo as- genome evolution. PLos Genet., 4, e1000237. 9. Hazkani-Covo,E. (2009) Mitochondrial insertions into primate sembly approach can serve as a valuable tool for when nuclear genomes suggest the use of numts as a tool for phylogeny. exploring such insertions in nuclear genomes for datasets Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 Mol. Biol. Evol., 26, 2175–2179. where the further molecular characterization of samples is 10. Hazkani-Covo,E., Zeller,R.M. and Martin,W. (2010) Molecular not possible. We also acknowledge the potential bias cre- poltergeists: mitochondrial DNA copies (numts) in sequenced nuclear ated due to poorly assembled references or gaps for certain genomes. PLos Genet., 6, e1000834. 11. Ricchetti,M., Tekaia,F. and Dujon,B. (2004) Continued colonization species. However, by characterizing shared events between of the human genome by mitochondrial DNA. PLoS Biol., 2, species, we have been able to derive the most comprehensive 1313–1324. overview till date of the evolution of these NumtS across 12. Jensen-Seaman,M.I., Wildschutte,J.H., Soto-Calderón,I.D. and primates and macaque spanning over several million years Anthony,N.M. (2009) A comparative approach shows differences in patterns of numt insertion during hominoid evolution. J. Mol. Evol., of divergence. 68, 688–699. 13. Tsuji,J., Frith,M.C., Tomii,K. and Horton,P. (2012) Mammalian SUPPLEMENTARY DATA NUMT insertion is non-random. Nucleic Acids Res., 40, 9073–9088. 14. Hazkani-Covo,E. and Graur,D. (2007) A comparative analysis of Supplementary Data are available at NARGAB Online. numt evolution in human and chimpanzee. Mol. Biol. Evol., 24, 13–18. 15. Bensasson,D., Feldman,M.W. and Petrov,D.A. (2003) Rates of DNA ACKNOWLEDGEMENTS duplication and mitochondrial DNA insertion in the human genome. We would like to thank Jeffrey M. Kidd for his insightful J. Mol. Evol., 57, 343–354. 16. Gherman,A., Adrian,G., Chen,P.E., Teslovich,T.M., Pawel,S., advice on the analysis and also Christian Brion for help with Marjorie,W., Kashuk,C.S., Aravinda,C., Lupski,J.R., Cutler,D.J. figure construction. et al. (2007) Population bottlenecks as a potential major shaping force of human genome architecture. PLos Genet., 3, e119. 17. Jensen-Seaman,M.I., Sarmiento,E.E., Deinard,A.S. and Kidd,K.K. FUNDING (2004) Nuclear integrations of mitochondrial DNA in gorillas. Am. J. National Institutes of Health (1R01HG007068-01A1 to Primatol., 63, 139–147. 18. Hazkani-Covo,E., Sorek,R. and Graur,D. (2003) Evolutionary R.E.M). T.M.B is supported by funding from the Eu- dynamics of large numts in the human genome: rarity of independent ropean Research Council (ERC) under the European insertions and abundance of post-insertion duplications. J. Mol. Union’s Horizon 2020 research and innovation pro- Evol., 56, 169–174. gramme (grant agreement No. 864203), BFU2017-86471-P 19. Tourmen,Y., Baris,O., Dessen,P., Jacques,C., Malthièry,Y. and (MINECO/FEDER, UE), “Unidad de Excelencia Marı́a Reynier,P. (2002) Structure and chromosomal distribution of human mitochondrial pseudogenes. Genomics, 80, 71–77. de Maeztu”, funded by the AEI (CEX2018-000792-M), 20. Prado-Martinez,J., Sudmant,P.H., Kidd,J.M., Li,H., Kelley,J.L., Howard Hughes International Early Career, Obra So- Lorente-Galdos,B., Veeramah,K.R., Woerner,A.E., O’Connor,T.D., cial ”La Caixa” and Secretaria d’Universitats i Recerca Santpere,G. et al. (2013) Great ape genetic diversity and population and CERCA Programme del Departament d’Economia i history. Nature, 499, 471–475. 21. Gokcumen,O., Tischler,V., Tica,J., Zhu,Q., Iskow,R.C., Lee,E., Coneixement de la Generalitat de Catalunya (GRC 2017 Fritz,M.H.-Y., Langdon,A., Stütz,A.M., Pavlidis,P. et al. (2013) SGR 880). Primate genome architecture influences structural variation Conflict of interest statement. None declared. mechanisms and functional consequences. Proc. Natl. Acad. Sci. U.S.A., 110, 15764–15769. 22. Li,H. and Durbin,R. (2009) Fast and accurate short read alignment REFERENCES with Burrows-Wheeler transform. Bioinformatics, 25, 1754–1760. 1. Soto-Calderón,I.D., Lee,E.J., Jensen-Seaman,M.I. and 23. DePristo,M.A., Banks,E., Poplin,R., Garimella,K.V., Maguire,J.R., Anthony,N.M. (2012) Factors affecting the relative abundance of Hartl,C., Philippakis,A.A., del Angel,G., Rivas,M.A., Hanna,M. nuclear copies of mitochondrial DNA (numts) in hominoids. J. Mol. et al. (2011) A framework for variation discovery and genotyping Evol., 75, 102–111. using next-generation DNA sequencing data. Nat. Genet., 43, 2. Lang,M., Sazzini,M., Calabrese,F.M., Simone,D., Boattini,A., 491–498. Romeo,G., Luiselli,D., Attimonelli,M. and Gasparre,G. (2012) 24. Prüfer,K., Munch,K., Hellmann,I., Akagi,K., Miller,J.R., Walenz,B., Polymorphic NumtS trace human population relationships. Hum. Koren,S., Sutton,G., Kodira,C., Winer,R. et al. (2012) The bonobo Genet., 131, 757–771. genome compared with the chimpanzee and human genomes. Nature, 3. Zischler,H., Hans,Z., Helga,G., von Haeseler,A. and Svante,P. (1995) 486, 527–531. A nuclear ‘fossil’ of the mitochondrial D-loop and the origin of 25. Lobon,I., Tucci,S., de Manuel,M., Ghirotto,S., Benazzo,A., modern humans. Nature, 378, 489–492. Prado-Martinez,J., Lorente-Galdos,B., Nam,K., Dabad,M., 4. Baldo,L., de Queiroz,A., Hedin,M., Hayashi,C.Y. and Gatesy,J. Hernandez-Rodriguez,J. et al. (2016) Demographic history of the (2011) Nuclear-mitochondrial sequences as witnesses of past genus pan inferred from whole mitochondrial genome interbreeding and population diversity in the jumping bristletail reconstructions. Genome Biol. Evol., 8, 2020–2030. Mesomachilis. Mol. Biol. Evol., 28, 195–210. NAR Genomics and Bioinformatics, 2020, Vol. 2, No. 4 11 26. Calabrese,F.M., Simone,D. and Attimonelli,M. (2012) Primates and 32. Speir,M.L., Zweig,A.S., Rosenbloom,K.R., Raney,B.J., Paten,B., mouse NumtS in the UCSC Genome Browser. BMC Bioinformatics, Nejad,P., Lee,B.T., Learned,K., Karolchik,D., Hinrichs,A.S. et al. 13(Suppl. 4), S15. (2016) The UCSC Genome Browser database: 2016 update. Nucleic 27. Huang,X. and Madan,A. (1999) CAP3: A DNA sequence assembly Acids Res., 44, D717–D725. program. Genome Res., 9, 868–877. 33. Krzywinski,M., Schein,J., Birol,I., Connors,J., Gascoyne,R., 28. Kronenberg,Z.N., Fiddes,I.T., Gordon,D., Murali,S., Cantsilieris,S., Horsman,D., Jones,S.J. and Marra,M.A. (2009) Circos: an Meyerson,O.S., Underwood,J.G., Nelson,B.J., Chaisson,M.J.P., information aesthetic for comparative genomics. Genome Res., 19, Dougherty,M.L. et al. (2018) High-resolution comparative analysis 1639–1645. of great ape genomes. Science, 360, 1–11. 34. Herrero,J., Muffato,M., Beal,K., Fitzgerald,S., Gordon,L., 29. Gordon,D., Huddleston,J., Chaisson,M.J.P., Hill,C.M., Pignatelli,M., Vilella,A.J., Searle,S.M.J., Amode,R., Brent,S. et al. Kronenberg,Z.N., Munson,K.M., Malig,M., Raja,A., Fiddes,I., (2016) Ensembl comparative genomics resources. Database, 2016, Hillier,L.W. et al. (2016) Long-read sequence assembly of the gorilla 1–17. genome. Science, 352, aae0344. 35. Shigenaga,M.K., Hagen,T.M. and Ames,B.N. (1994) Oxidative 30. Zhou,W., Emery,S.B., Flasch,D.A., Wang,Y., Kwan,K.Y., Kidd,J.M., damage and mitochondrial decay in aging. Proc. Natl. Acad. Sci. Moran,J.V. and Mills,R.E. (2019) Identification and characterization U.S.A., 91, 10771–10778. Downloaded from https://academic.oup.com/nargab/article/2/4/lqaa089/5983420 by guest on 10 March 2022 of occult human-specific LINE-1 insertions using long-read 36. Lin,Y.-L. and Gokcumen,O. (2019) Fine-scale characterization of sequencing technology. Nucleic Acids Res., 48, 1146–1163. genomic structural variation in the human genome reveals adaptive 31. Koren,S., Walenz,B.P., Berlin,K., Miller,J.R., Bergman,N.H. and and biomedically relevant hotspots. Genome Biol. Evol., 11, Phillippy,A.M. (2017) Canu: scalable and accurate long-read 1136–1151. assembly via adaptive k-mer weighting and repeat separation. Genome Res., 27, 722–736.