Article Towards complete and error-free genome assemblies of all vertebrate species ­­­­ https://doi.org/10.1038/s41586-021-03451-0 A list of authors and their affiliations appears at the end of the paper. Received: 22 May 2020 Accepted: 12 March 2021 High-quality and complete reference genome assemblies are fundamental for the Published online: 28 April 2021 application of genomics to biology, disease, and biodiversity conservation. However, such assemblies are available for only a few non-microbial species1–4. To address this Open access issue, the international Genome 10K (G10K) consortium5,6 has worked over a five-year Check for updates period to evaluate and develop cost-effective methods for assembling highly accurate and nearly complete reference genomes. Here we present lessons learned from generating assemblies for 16 species that represent six major vertebrate lineages. We confirm that long-read sequencing technologies are essential for maximizing genome quality, and that unresolved complex repeats and haplotype heterozygosity are major sources of assembly error when not handled correctly. Our assemblies correct substantial errors, add missing sequence in some of the best historical reference genomes, and reveal biological discoveries. These include the identification of many false gene duplications, increases in gene sizes, chromosome rearrangements that are specific to lineages, a repeated independent chromosome breakpoint in bat genomes, and a canonical GC-rich pattern in protein-coding genes and their regulatory regions. Adopting these lessons, we have embarked on the Vertebrate Genomes Project (VGP), an international effort to generate high-quality, complete reference genomes for all of the roughly 70,000 extant vertebrate species and to help to enable a new era of discovery across the life sciences. Chromosome-level reference genomes underpin the study of func- Towards this end, having learned the lessons of having too many tional, comparative, and population genomics within and across spe- variables that make conclusions more difficult to reach in the G10K cies. The first high-quality genome assemblies of human1 and other from the G10K Assemblathon 2 effort15, we first evaluated multiple model species (for example, Caenorhabditis elegans2, mouse3, and genome sequencing and assembly approaches extensively on one zebrafish4) were put together using 500–1,000-base pair (bp) Sanger species, the Anna’s hummingbird (Calypte anna). We then deployed sequencing reads of thousands of hierarchically organized clones with the best-performing method across sixteen species representing six 200–300-kilobase (kb) inserts, and chromosome genetic maps. This major vertebrate classes, with a wide diversity of genomic character- approach required tremendous manual effort, software engineering, istics. Drawing on the principles learned, we improved these methods and cost, in decade-long projects. Whole-genome shotgun approaches further, discovered parameters and approaches that work better for simplified the logistics (for example, in human7 and Drosophila8), and species with different genomic characteristics, and made biological later next-generation sequencing with shorter (30–150-bp) sequenc- discoveries that had not been possible with the previous assemblies. ing reads and short insert sizes (for example, 1 kb) ushered in more affordable and scalable genome sequencing9. However, the shorter reads resulted in lower-quality assemblies, fragmented into thousands Complete, accurate assemblies require long reads of pieces, where many genes were missing, truncated, or incorrectly We chose a female Anna’s hummingbird because it has a relatively assembled, resulting in annotation and other errors10. Such errors small genome (about 1 Gb), is heterogametic (has both Z and W sex can require months of manual effort to correct individual genes and chromosomes), and has an annotated reference of the same individual years to correct an entire assembly. Genomic heterozygosity posed built from short reads16. We obtained 12 new sequencing data types, additional problems, because homologous haplotypes in a diploid including both short and long reads (80 bp to 100 kb), and long-range or polyploid genome are forced together into a single consensus by linking information (40 kb to more than 100 Mb), generated using eight standard assemblers, sometimes creating false gene duplications11–14. technologies (Supplementary Table 1). We benchmarked all technolo- To address these problems, the G10K consortium5,6 initiated the gies and assembly algorithms (Supplementary Table 2) in isolation and Vertebrate Genomes Project (VGP; https://vertebrategenomesproject. in many combinations (Supplementary Table 3). To our knowledge, org) with the ultimate aim of producing at least one high-quality, near this was the first systematic analysis of many sequence technologies, error-free and gapless, chromosome-level, haplotype-phased, and assembly algorithms, and assembly parameters applied on the same annotated reference genome assembly for each of the 71,657 extant individual. We found that primary contiguous sequences (contigs) named vertebrate species and using these genomes to address funda- (pseudo-haplotype; Supplementary Note 1) assembled from Pacific mental questions in biology, disease, and biodiversity conservation. Biosciences continuous long reads (CLR) or Oxford Nanopore long Nature | Vol 592 | 29 April 2021 | 737 Article reads (ONT) were approximately 30- to 300-fold longer than those Figs. 2a, 3b). We set initial minimum assembly metric goals of: 1 Mb assembled from Illumina short reads (SR), regardless of data type com- contig NG50; 10 Mb scaffold NG50; assigning 90% of the sequence bination or assembly algorithm used (Fig. 1a, Supplementary Table 3). to chromosomes, structurally validated by at least two independent The highest contig NG50s for short-read-only assemblies were about lines of evidence; Q40 average base quality; and haplotypes assem- 0.025 to 0.169 Mb, whereas for long reads they were about 4.6 to 7.66 Mb bled as completely and correctly as possible. When these metrics were (Fig. 1a); contig NG50 is an assembly metric based on a weighted median achieved, most genes were assembled with gapless exon and intron of the lengths of its gapless sequences relative to the estimated genome structures11, and fewer than 3% had frame-shift base errors identified size. After fixing a function in the PacBio FALCON software17 that caused in annotation. Q40 is the mathematical inflection point at which genes artificial breaks in contigs between stretches of highly homozygous and go from usually containing an error to usually not21. Of the curated heterozygous haplotype sequences (Supplementary Note 1, Supple- assemblies (Supplementary Table 10, Supplementary Note 2), 16 of mentary Table 2), contig NG50 nearly tripled to 12.77 Mb (Fig. 1a). These 17 achieved the desired continuity metrics (Extended Data Table 1). findings are consistent with theoretical predictions18 and demonstrate Scaffold NG50 was significantly correlated with genome size (Fig. 2a), that, given current sequencing technology and assembly algorithms, suggesting that larger genomes tend to have larger chromosomes. On it is not possible to achieve high contig continuity with short reads average, 98.3% of the assembled bases had reliable block NG50s ranging alone, as it is typically impossible to bridge through repeats that are from 2.3 to 40.2 Mb; collapsed repeat bases22 with abnormally high CLR longer than the read length. read coverage (more than 3 s.d.) ranged from 0.7 to 31.4 Mb per Gb; and the completeness of the genome assemblies ranged from 87.2 to 98.1%, with less than 4.9% falsely duplicated regions, consistent with the false Iterative assembly pipeline duplication rate we found for the conserved BUSCO vertebrate gene Scaffolds generated with all three scaffolding technologies (that is, set (Extended Data Table 1, Supplementary Tables 11, 12). 10X Genomics linked reads (10XG), Bionano optical maps (Opt.), and Arima Genomics, Dovetail Genomics, or Phase Genomics Hi-C) were approximately 50% to 150% longer than those generated using one Repeats markedly affect continuity or two technologies, regardless of whether we started with short- or For assemblies generated using our automated pipeline (Extended long-read-based contigs (Fig. 1b, Extended Data Fig. 1a, Supplemen- Data Fig. 3a) before manual curation, all but 2 (the thorny skate and tary Table 3). These findings include improvements we made to each channel bull blenny) of the 17 assemblies exceeded the desired continu- approach (Supplementary Note 1, Supplementary Tables 4, 5, Supple- ity metrics (Supplementary Table 13). In searching for an explanation mentary Fig. 1). Despite similar scaffold continuity, the short-read-only of these results, we found that contig NG50 decreased exponentially assemblies had from about 18,000 to about 70,000 gaps, whereas the with increasing repeat content, with the thorny skate having the high- long-read assemblies had substantially fewer (about 400 to about est repeat content (Fig. 2b, Supplementary Table 13). Consequently, 4,000) gaps (Fig. 1c). Many gaps in the short-read assemblies were after scaffolding and gap filling, we observed a significant positive in repeat or GC-rich regions. Considering the curated version of this correlation between repeat content and number of gaps (Fig. 2c). The assembly to be more accurate, we also identified roughly 5,000 to 8,000 kākāpō parrot, which had 15% repeat content, had about 325 gaps per mis-joins in short-read-based assemblies, whereas long-read-based Gb, including 2 of 26 chromosomes with no gaps (chromosomes 16 assemblies had only from 20 to around 700 mis-joins (Fig. 1d). These and 18) and no evidence of collapses or low support, suggesting that mis-joins included chimeric joins and inversions. After we curated the chromosomal contigs were complete (bStrHab1 in Supplemen- this assembly for contamination, assembly errors, and Hi-C-based tary Table 6). By contrast, the thorny skate, with 54% repeat content, chromosome assignments (Fig. 1e, f), the final hummingbird assembly had about 1,400 gaps per Gb (Extended Data Table 1); none of its 49 had 33 scaffolds that closely matched the chromosome karyotype in chromosomal-level scaffolds contained fewer than eight gaps, and all number (33 of 36 autosomes plus sex chromosomes) and estimated had some regions that contained collapses or low support (sAmbRad1 sizes (approximately 2 to 200 Mb; Fig. 1g, h), with only 1 to 30 gaps per in Supplementary Table 6). Even after curation and other modifications autosome (bCalAnn1 in Supplementary Table 6). Of the five autosomes to increase assembly quality (Supplementary Note 2), the number of with only one gap each, three (chromosomes 14, 15, and 19) had com- collapses, their total size, missing bases, and the number of genes in the plete spanning support by at least two technologies (reliable blocks, collapses all correlated with repeat content (Extended Data Fig. 4a–d). Extended Data Fig. 1c; bCalAnn1 in Supplementary Table 6), indicat- The average collapsed length, however, correlated with average CLR ing that the chromosome contigs were nearly complete. However, read lengths (10–35 kb; Extended Data Fig. 4e). There were no cor- they were missing long arrays of vertebrate telomere repeats within relations between the number of collapsed bases and heterozygosity 1 kb of their ends (Extended Data Fig. 1c; bCalAnn1 in Supplementary or genome size (Extended Data Fig. 4f, g). Depending on species, 77.4 Tables 6, 7). to 99.2% of the collapsed regions consisted of unresolved segmental duplications (Extended Data Fig. 4h). The remainder were high-copy repeats, mostly of previously unknown types (Extended Data Fig. 4i), Assembly pipeline across vertebrate diversity and of known types such as satellite arrays, simple repeats, long termi- Using the formula that gave the highest-quality hummingbird nal repeats (LTRs), and short and long interspersed nuclear elements genome, we built an iterative VGP assembly pipeline (v1.0) with (SINEs and LINEs), depending on species (Extended Data Fig. 4j). We haplotype-separated CLR contigs, followed by scaffolding with linked found that repeat masking before generating contigs prevented some reads, optical maps, and Hi-C, and then gap filling, base call polishing, repeats from making it into the final assembly (Supplementary Note 3). and finally manual curation (Extended Data Figs. 2a, 3a). We systemati- All of the above findings quantitatively demonstrate the effect that cally tested our pipeline on 15 additional species spanning all major repeat content has on the ability to produce highly continuous and vertebrate classes: mammals, birds, non-avian reptiles, amphibians, complete assemblies. teleost fishes, and a cartilaginous fish (Supplementary Tables 8, 9, Supplementary Note 2). For the zebra finch, we used DNA from the same male as was used to generate the previous reference genome19, Detection and removal of false duplications and included a female trio for benchmarking haplotype completeness, During curation, we discovered that one of the most common assem- where sequenced reads from the parents were used to bin parental bly errors was the introduction of false duplications, which can be haplotype reads from the offspring before assembly20 (Extended Data misinterpreted as exon, whole-gene, or large segmental duplications. 738 | Nature | Vol 592 | 29 April 2021 a v1.0* curated 12.77 b 44.7 c 429 d 0 v1.0* 12.77 72.7 434 20 CLR + linked + Opt. + Hi-C 4.53 35.2 3,083 730 CLR + Opt. + Hi-C 4.53 25.5 3,998 686 Long read CLR + linked + Hi-C 4.62 22 575 658 CLR + linked + Opt. 4.53 21.3 2,980 692 CLR + Opt. 4.53 15.5 3,798 576 CLR + Hi-C 4.54 15.8 432 680 CLR + linked 4.62 12.9 372 560 CLR 4.62 4.6 0 458 CLR + 1D and 2D 7.66 7.7 0 559 NRGene 0.06 21.9 37,501 5,193 Linked (Supernova2) 0.17 21.5 18,749 6,998 Short read SR + Opt. + Hi-C 0.025 31.4 70,676 7,982 SR + Opt. 0.025 14.8 70,422 7,922 SR + Hi-C 0.025 11.7 70,511 7,909 SR 0.025 4 70,084 7,861 0 5 10 15 0 25 50 75 100 0 25,000 50,000 75,000 100,000 0 3,000 6,000 9,000 Contig NG50 (Mb) Scaffold NG50 (Mb) No. gaps No. mis-joins e f g h 100 Karyotype size (Mb) Mis-join Z 1 2 3 4 5 6 7 8 30 W 9 10 11 12 13 14 15 16 10 17 18 19 20 21 22 23 24 3 25 26 27 28 29 30 31 32 3 10 30 100 33 34 35 36 Z W v1.0* v1.0* curated Assembled chr. size (Mb) Fig. 1 | Comparative analyses of Anna’s hummingbird genome assemblies 34 chromosomes. Grid lines indicate scaffold boundaries. Red arrow, example with various data types. a, Contig NG50 values of the primary mis-join that was corrected during curation. g, Karyotype of the identified pseudo-haplotype. b, Scaffold NG50 values. c, Number of joins (gaps). chromosomes (n = 36 + ZW), consistent with previous findings70. h, Correlation d, Number of mis-join errors compared with the curated assembly. The curated between estimated chromosome sizes (in Mb) based on karyotype images in g assembly has no remaining conflicts with the raw data and thus no known and assembled scaffolds in Supplementary Table 4 (bCalAna1) on a log–log mis-joins. *Same as CLR + linked + Opt. + Hi-C, but with contigs generated with scale. v1.0, VGP assembly v1.0 pipeline; linked, 10X Genomics linked reads; an updated FALCON17 version and earlier Hi-C Salsa version (v2.0 versus v2.2; Hi-C, Hi-C proximity ligation; 1D, 2D, Oxford Nanopore long reads; NRGene, Supplementary Table 2) for less aggressive contig joining. e, f, Hi-C NRGene paired-end Illumina reads; SR, paired-end Illumina short reads. interaction heat maps before and after manual curation, which identified We observed two types of false duplication: 1) heterotype duplications, we applied these tools, the primary assembly sizes and the k-mer and which occurred in regions of increased sequence divergence between BUSCO gene duplication rates were all reduced, and their correlations paternal and maternal haplotypes, where separate haplotype contigs with heterozygosity and repeat content were also reduced or eliminated were incorrectly placed in the primary assembly (Extended Data Fig. 5a); (Fig. 2d–h). These findings indicate that it is essential to properly phase and 2) homotype duplications, which occurred near contig boundaries or haplotypes and to obtain high consensus sequence accuracy in order to under-collapsed sequences caused by sequencing errors (Extended Data prevent false duplications and associated biologically false conclusions. Fig. 5b). False heterotype duplications appeared to occur with higher heterozygosity. For example, during curation of the female zebra finch genome, we found an approximately 1-Mb falsely duplicated heterozy- Curation is needed for a high-quality reference gous sequence (Extended Data Fig. 6a). This zebra finch individual had Each automated scaffolding method introduced tens to thousands the highest heterozygosity (1.6%) relative to all other genomes (0.1–1.1%). of unique joins and breaks in contigs or scaffolds (Supplementary Homotype duplications often occurred at contig boundaries, and were Table 14). Depending on species, the first scaffolding step with linked approximately the same length as the sequence reads (Extended Data reads introduced about 50–900 joins between CLR-generated contigs. Fig. 6b, c). We identified and removed false duplications during curation Optical maps introduced a further roughly 30–3,500 joins, followed by using read coverage, self-, transcript-, optical map- and Hi-C-alignments, Hi-C with about 30–700 more joins, and each identified up to several and k-mer profiles (Extended Data Fig. 6, Supplementary Fig. 2). dozen joins that were inconsistent with the previous scaffolding step. Before we purged false duplications, the primary assembly genome Manual curation resulted in an additional 7,262 total interventions for size correlated positively with estimated percentage heterozygosity; 19 genome assemblies or 236 interventions per Gb of sequence (Sup- more heterozygous genomes tended to have assembly sizes bigger plementary Table 15). When a genome assembly was available for the than the estimated haploid genome size (Fig. 2d). Similarly, the extra same or a closely related species, it was used to confirm putative chro- duplication rate in the primary assembly, measured using k-mers23 or mosomal breakpoints or rearrangements (Supplementary Table 15). conserved vertebrate BUSCO genes24, varied from 0.3% to 30% and These interventions indicate that even with current state-of-the-art trended towards correlation with heterozygosity (Fig. 2e, f, Supple- assembly algorithms, curation is essential for completing high-quality mentary Table 13). Apparent false gene duplication rates correlated reference assemblies and for providing iterative feedback to improve more strongly with the overall repeat rate in the assemblies (Fig. 2g, h). assembly algorithms. A further description of our curation approach To remove these false duplications automatically, we initially used and analyses of VGP genomes are presented elsewhere25. Purge_Haplotigs13, which removed retained falsely duplicated contigs that were not scaffolded (Extended Data Fig. 5; VGP v1.0–1.5). Later, we developed Purge_Dups14 to remove both falsely retained contigs Hi-C scaffolding and cytological mapping and end-to-end duplicated contigs within scaffolds (Extended Data Most large assembled scaffolds of each species spanned entire chromo- Fig. 5; VGP v1.6), which reduced the amount of manual curation. After somes, as shown by the relatively clean Hi-C heat map plots across each Nature | Vol 592 | 29 April 2021 | 739 Article a 600 b 30 c 1,500 r 2 = 0.4, d 110 * r 2 = 0.4, e** r 2 = 0.7, for Scaffold NG50 (Mb) Contig NG50 (Mb) 500 P = 1.9 x 10–5 P = 0.0044 P = 0.0025 Be Primary size (%) 100 Gaps per Gb r 2 = 0.4, 400 20 1,000 * ged P = 0.0053 ed Pur d *** r at 300 te Cu 90 a C ur 200 10 500 r 2 = 0.3, Curated** 80 P = 8.9 × 10–4 100 r 2 = 0.079, Before** P = 0.14 0 0 0 70 0 1 2 3 4 5 6 0 20 40 60 0 10 20 30 40 50 60 0 0.5 1.0 1.5 Genome size (Gb) Repeat (%) Repeat (%) Heterozygosity (%) e 30 f g h 20 20 30 BUSCO duplication (%) BUSCO duplication (%) k-mer duplication (%) r 2 = 0.048, r 2 = 0.15, r 2 = 0.45, r 2 = 0.30, k-mer duplication (%) P = 0.20 P = 0.07 P = 0.0019 P = 0.013 r 2 = 0.029, 15 r 2 = 0.058, r 2 = 0.47, 15 r 2 = 0.035, 20 P = 0.24 P = 0.18 20 P = 0.0014 P = 0.0069 10 ** 10 re re e* Bef o fo for 10 Be fore 10 Be Be 5 5 Purged ** * Purged Purged Purged* 0 0 0 0 0 0.5 1.0 1.5 0 0.5 1.0 1.5 0 10 20 30 40 50 60 0 10 20 30 40 50 60 Heterozygosity (%) Heterozygosity (%) Repeat (%) Repeat (%) Fig. 2 | Impact of repeats and heterozygosity on assembly quality. heterozygosity before and after purging of false duplications. g, h, As in e, f, but a, Correlation between scaffold NG50 and genome size of the curated with whole-genome repeat content before and after purging of false assemblies. b, Nonlinear correlation between contig NG50 and repeat content, duplications. Genome size, heterozygosity, and repeat content were estimated before and after curation. c, Correlation between number of gaps per Gb from 31-mer counts using GenomeScope71, except for the channel bull blenny, assembled and repeat content. d, Correlation between primary assembly size as the estimates were unreliable (see Methods). Repeat content was measured relative to estimated genome size (y axis) and genome heterozygosity (x axis), by modelling the k-mer multiplicity from sequencing reads. Sequence before and after purging of false duplications. Assembly sizes above 100% duplication rates were estimated with Merqury23 using 21-mers. *P < 0.05; indicate the presence of false duplications and those below 100% indicate **P < 0.01; ***P < 0.001, of the correlation coefficient: P values and adjusted r2 collapsed repeats. e, f, Correlations between genome duplication rate using from F-statistics. n = 17 assemblies of 16 species. k-mers23 (e) and conserved BUSCO vertebrate gene set (f), and genome scaffold after curation (Extended Data Fig. 7), near perfect correlation Data Fig. 9a) using our VGP trio pipeline (Extended Data Fig. 3b). Com- between chromosomal scaffold length and karyotypically determined pared to the non-trio assembly of the same individual, the trio version chromosome length (Fig. 1h), and the presence of telomeric repeat had seven- to eightfold fewer false duplications (k-mer and BUSCO dups motifs on some scaffold ends (Supplementary Table 7). In our VGP zebra in Supplementary Tables 11, 12), well-preserved haplotype-specific vari- finch assembly, all inferred chromosomes were consistent with previ- ants (k-mer precision/recall 99.99/97.08%), and higher base call accu- ously identified linkage groups in the Sanger-based reference, except racy, exceeding Q43 for both haplotypes (Extended Data Table 1). The for chromosomes 1 and 1B (Extended Data Fig. 8a). Their join in the VGP trio-based assembly was the only assembly with nearly perfect (99.99%) assembly was supported by both single CLR reads and optical maps separation of maternal and paternal haplotypes, determined using through the junction. We also corrected nine inversion errors and filled k-mers specific to each23. We identified haplotype-specific structural in large gaps at some chromosome ends. In the platypus, we identified variants, including inversions of 4.5 to 12.5 Mb on chromosomes 5, 11, 18 structural differences in 13 scaffolds between the VGP assembly and and 13 that were not readily identifiable in the non-trio version (Extended the previous Sanger-based reference anchored to chromosomes using Data Fig. 10a–e). Moving forward, the VGP is prioritising the collection fluorescence in situ hybridization (FISH) physical mapping (Extended of mother–father–offspring trios where possible, or single parent–off- Data Fig. 8b, Supplementary Table 16). Of these 18, all were supported spring duos, to assist with diploid assembly and phasing, as well as the with Hi-C, and seven were also supported by both CLR and optical maps development of improved methods for the assembly of diploid genomes in the VGP assembly. Our platypus assembly also filled in many large in the absence of parental genomic data, as described in another study29. (approximately 1–30 Mb) gaps and corrected many inversion errors (Extended Data Fig. 8b). Furthermore, we identified seven additional chromosomes (chromosomes 30–36) in the zebra finch, and eight Effects of polishing on accuracy (chromosomes 8, 9, 14, 15, 17, 19, 21, and X4; Extended Data Fig. 8a, b) Despite their increased continuity and structural accuracy, CLR-based in the platypus26,27. Relative to the VGP assembly, the earlier short-read assemblies required at least two rounds of short-read consensus polish- Anna’s hummingbird assembly was highly fragmented (Extended Data ing to reach 99.99% base-level accuracy (one error per 10 kb, Phred30 Fig. 8c), despite being scaffolded with seven different Illumina libraries Q40; Supplementary Table 5). Before polishing, the per-base accuracy spanning a wide range of insert sizes (0.2–20 kb). The previous climbing was Q30–35 (calculated using k-mers). The most common errors were perch assembled chromosomes were even more fragmented and also short indels from inaccurate consensus calling during CLR contig for- had large gaps of missing sequence (Extended Data Fig. 8d). On average, mation, which resulted in amino acid frameshift errors. Using our com- 97% ± 3% (s.d.) of the assembled bases were assigned to chromosomes bined approach of long-read and short-read polishing applied on both (Extended Data Table 1), compared with 76% and 32% in the prior zebra primary and alternate haplotype sequences together, we polished from finch and platypus references, respectively. We believe the comparable 82% to 99.7% of the primary and about 91.3% of the alternate assembly or higher accuracy of Hi-C relative to genetic linkage or FISH physical (Supplementary Table 17). Of the remaining unpolished sequence, mapping is due to the higher sampling rate of Hi-C pairs across the one haplotype was sometimes reconstructed at substantially lower genome. Nonetheless, visual karyotyping is useful for complementary quality, because most reads aligned to the higher quality haplotype validation of chromosome count and structure28. (Extended Data Fig. 11a). False duplications had similar effects, where the duplicated sequence acted as an attractor during the read mapping. Haplotypes in the more homozygous regions tended to be collapsed by Trios help to resolve haplotypes FALCON-Unzip17. All such cases recruited reads from both haplotypes We were able to assemble the trio-based female zebra finch contigs into and thereby caused switch errors, which we confirmed in the trio-based separate maternal and paternal chromosome-level scaffolds (Extended assembly and fixed when excluding read pairs from the other haplotype 740 | Nature | Vol 592 | 29 April 2021 a b c d Fig. 3 | Improvements to alignments and No. partial coding genes P < 2.0 × 10–16 Avg. reads mapped (%) 100 100 50,000 6,000 P< P = 1.27 × 10–6 annotations in VGP assemblies relative to prior 2.0 × 10–16 supported CDS 90 90 40,000 5,000 references. a, b, Average percentage of RNA-seq Total or fully P= 5.45 × 10–12 4,000 80 80 30,000 transcriptome samples (a; n = 44, mean ± s.e.m.) 3,000 Assembly 70 70 20,000 and ATAC–seq genome reads (b; n = 12) that align to 2,000 Previous 60 60 10,000 the previous and VGP zebra finch assemblies. Trio 1,000 VGP 50 50 0 0 Unique reads mapped to only one location in the Total Unique Total Unique Anna’s Zebra Platypus Anna’s Zebra Platypus RNA-seq ATAC−seq hummingbird finch hummingbird finch assembly. Total is the sum of unique and multi-mapped reads. P values are from paired e SPC25 (SPC25 component of NDC80 f GABRG2 (gamma-aminobutyric t-test. c, d, Total number of coding sequence (CDS) kinetochore complex) acid type A receptor gamma2 subunit) Taeniopygia_guttata-3.2.4 (previous) Taeniopygia_guttata-3.2.4 (previous) transcripts (full bar) and portion fully supported 3′ 5′3′ 5′ 5′ 3′ 5′ 3′ (inner bar) (c) and the number of RefSeq coding SPC25-like SPC25 (LOC100231268) 2–7 GABRG2-like (LOC101232861) 2–5 GABRG2 genes annotated as partial (d) in the previous and 3′ 5′ 5′ 3′ VGP assemblies using the same input data. e–h, SPC25 GABRG2 bTaeGut1_v1.p (VGP) bTaeGut1_v1.p (VGP) Examples of assembly and associated annotation errors in previous reference assemblies corrected g VTG2 (vitellogenin-2) h YIPF6 (Yip domain family, member 6) in the new VGP assemblies. See main text for Ornithorhynchus_anatinus_5.0.1 (previous) Anabas_testudineus ASM90030266v1 (previous) descriptions. i, Gene synteny around the VTR2C 5′ 3′ 5′ 3′ receptor in the platypus shows completely missing CAPN13 intron VTG2-like VTG2-like (No annotation available) (LOC100083241) (LOC100680960) genes (NUDT16), truncated and duplicated 6 12 25,26 32 1 5′ 3′ 3′ 3 2 5′ ARHGAP4, and many gaps in the earlier VTG2-like (LOC100083241) YIPF6 Sanger-based assembly compared with the filled in mOrnAna1.p.v1 (VGP) fAnaTes1.1 (VGP) and expanded gene lengths in the new VGP Scaffold Gap (N) Coding exon Exon with error Corrected in VGP Exon blast hit with no annotation UTR Intron assembly. Assembly accessions are in Supplementary Table 19. i Ornithorhynchus_anatinus_5.0.1 (previous) NW_001794330.1 NW_001704642.1 NW_001786488.1 NW_001665319.1 SRPK3 SSR4 PDZD4 ARHGAP4 L1CAM VTR2C ARHGAP4 NAA10 RENBP IDH3G CAV2 SRPK3 SSR4 PDZD4 L1CAM VTR2C ARHGAP4 NAA10 IDH3G CAV2 NUDT16 RENBP mOrnAna1.p.v1 (VGP) Chromosome 6 (19,648,353–19,890,187 bp) during polishing (Extended Data Fig. 11b). These findings indicate that Data Fig. 9e). These duplications were verified using single-molecule both sequence read accuracy and careful haplotype separation are CLR reads that spanned the duplication junctions or even the entire MT important for producing accurate assemblies. genome. Their absence in previous MT references34,35 is likely to result from the inability of Sanger or short reads to correctly resolve large duplications. More details on the MT-VGP pipeline and new biological Sex chromosomes and mitochondrial genomes discoveries are reported elsewhere33. Sex chromosomes have been notoriously difficult to assemble, owing to their greater divergence relative to autosomes and high repeat con- tent31. We successfully assembled both sex chromosomes (Z, W) for all Improvements to read alignment and annotation three avian species, the first W chromosome (to our knowledge) for Compared to previous Sanger (zebra finch and platypus) and Illumina vocal learning birds (Extended Data Figs. 7, 9b), the X and/or Y chromo- (Anna’s hummingbird and climbing perch) assemblies, we added about some in placental mammals (Canada lynx and two bat species), the X 42–176 Mb of missing sequence and placed 68.5 Mb (zebra finch) to chromosome in the thorny skate, and for the first time, to our knowl- 1.8 Gb (platypus) of previously unplaced sequence within chromo- edge, all ten sex chromosomes (5X and 5Y) in the platypus26 (Extended somes. We corrected about 7,800–64,000 mis-joins, and closed 55,177– Data Fig. 9c). The completeness and continuity of the zebra finch Z and 193,137 gaps per genome (Supplementary Table 19). Consistent with W chromosomes were further improved by the trio-based assembly these improvements, both transcriptome RNA sequencing (RNA-seq) (Extended Data Fig. 9b). However, the sex chromosome assemblies data (Fig. 3a) and genome assay for transposase-accessible chromatin were still more fragmented than the autosomes, probably owing to using sequencing (ATAC–seq) data (Fig. 3b) aligned with about 5 to their lower sequencing depth and high repeat content. 10% greater mapability to our new VGP assemblies compared with the Mitochondrial (MT) genomes, which are expected to be 11–28 kb previous assemblies. The NCBI RefSeq and EBI Ensembl annotations in size32, were initially found in only six assemblies (Supplementary revealed: 5,434 to 14,073 more protein-coding transcripts per species, Table 18). The MT-derived raw reads were present, but they failed to with 94.1 to 97.8% fully supported (Fig. 3c, Supplementary Table 20); assemble, in part because of minimum read-length cutoffs for the start- only about 100 to 300 partially assembled coding genes, compared ing contig assembly. Furthermore, if the MT genome was not present with about 1,600 to 5,600 (Fig. 3d); more orthologous coding genes during nuclear genome polishing, the raw MT reads were attracted to shared with human; and fewer transcripts that required corrections nuclear MT sequences (NuMTs), incorrectly converting them to the to compensate for premature stop codons or frame-shift indel errors full organelle MT sequence (Extended Data Fig. 11c). To address these (Extended Data Table 2). The total number of genes annotated went issues, we developed a reference-guided MT pipeline and included the down in the VGP assemblies (Extended Data Table 2), partly because MT genome during polishing33 (Extended Data Fig. 3c; VGP v1.6). With there were fewer false duplications (Supplementary Table 19). Sup- these improvements, we reliably assembled 16 of 17 MT genomes (Sup- porting these results, the VGP assemblies had 0 to 13% higher k-mer plementary Table 18) and discovered 2 kb of an 83-bp repeat expansion completeness (95% mean ± 3.5% s.d. versus 88 ± 4.3%; Extended Data within the control region in the kākāpō (Extended Data Fig. 9d), and Table 2, Supplementary Table 19; P = 0.0047, n = 4 prior and 17 VGP Nad1 and trnL2 gene duplications in the climbing perch (Extended assemblies, unpaired t-test). Nature | Vol 592 | 29 April 2021 | 741 Article Gene body Fig. 4 | VGP assemblies reveal GC content Exon patterns in protein-coding genes. a, Average GC Upstream (kb) 5′ UTR First Internal Last 3′ UTR Downstream (kb) a 70 content (n = 14,000–18,000 annotated coding Sequence 60 Found genes; Extended Data Table 2) in VGP assemblies GC content (%) Missing (black) and the percentage of genes with missing Missing (%) 50 40 Species sequence in the earlier references (red) based on a 30 Platypus Zebra finch Cactus alignment, in 100-bp blocks, 2 kb on either 20 Anna’s hummingbird side of all protein-coding genes (left and right), 10 Climbing perch and for UTRs, exons, and introns (middle). 0 2 1 0 0 1 2 Clade b, Average GC content (mean ± s.d. for lineages b 70 Bird with more than one species) of the six major Mammal Reptile vertebrate lineages sequenced, for 30 kb GC content (%) 60 Skate Amphibian upstream and downstream (in 100-bp blocks, log Fish 50 scale; left and right) and of the UTR, exons, and Type introns (middle). c, d, Left, specialized expression 40 Exon or up/downstream (arrows) shown by in situ hybridization of DRD1B 30 Intron in the zebra finch striatum (c) and ER81 in the 30 20 10 0 5′ UTR First Last 3′ UTR 0 10 20 30 arcopallium (d), from Jarvis et al.47; the cerebellum Intron c GC content (%) was removed from the ER81 image. Right, ATAC– DRD1B 0 5060 100 seq profiles in the GC-rich promoter regions of Chr. 4 13.580 13.581 13.582 13.583 13.584 Mb Striatum these genes, showing each gene’s GC content (red GC content (%) is high), the ATAC–seq peaks in striatum (purple) CpG islands Striatum ATAC or arcopallium (yellow) neurons, and portions of Arcopallium ATAC missing sequence (black) in the previous DRD1B VGP annotation 5′ 3′ reference assembly (grey). DRD1B Previous annotation 5′ 3′ Gaps in previous d ER81 Chr. 2 28.585 28.575 28.565 28.555 28.545 28.535 28.525 Mb Arcopallium GC content (%) CpG islands Striatum ATAC Arcopallium ATAC ER81 VGP annotation 5′ 3′ ER81 Previous annotation Gaps in previous An example of a whole-gene heterotype false duplication in the Ref- details in refs. 27,44), demonstrating that assembly quality has a critical Seq annotation of the previous zebra finch reference19 is the BUSCO effect on subsequent annotations and functional genomics. gene SPC2536, for which each haplotype was correctly placed in the VGP primary and alternate assemblies (Fig. 3e). The GABRG2 receptor, which shows specialized expression in vocal learning circuits37, had a partial GC-rich regulatory regions of coding genes tandem duplication of four of its ten exons, resulting in annotated We tested whether the higher-quality VGP assemblies enabled new partial false tandem gene duplications (GABRG2 and GABRG2-like; biological discoveries. Notably, beginning about 1.5 kb upstream of Fig. 3f). The vitellogenin-2 (VTG2) gene, a component of egg yolk in protein-coding genes, in 100-bp blocks, there was a steady increase all egg-laying species38, was distributed across 14 contigs in 3 different from about 6–20% to about 30–55% of genes having missing sequence scaffolds in the previous platypus assembly (Fig. 3g). Two of these scaf- in previous references (Fig. 4a); similarly high proportions of genes folds received two corresponding VTG2-like gene annotations, and the were missing their subsequent 5′ untranslated regions (UTRs) and first third was included as false duplicated intron in CAPN-13 (red), together exons. This fluctuation in missing sequence was directly proportional causing false amino acid sequences in five exons (blue). The BUSCO to GC content (Fig. 4a). We therefore studied the GC content pattern YIPF6 gene, which is associated with inflammatory bowel disease39, across all protein-coding genes in all 16 new VGP assemblies and found was split between two different scaffolds and is thus presumed to be a genome-wide signature: a rapid rise in GC content in the roughly 1.5 kb a gene loss in the earlier climbing perch assembly40 (Fig. 3h). Each of before the transcription start site, in the 5′ UTR, and in the first exon, these genes is now present on long VGP contigs, within validated blocks, followed by a steady decrease in subsequent exons and returning to with no gaps and no false gene gains or losses (Supplementary Table 21). near intergenic background levels in the 3′ UTR and about 1.5 kb after Going beyond individual genes, a ten-gene synteny window surround- the transcription termination site (Fig. 4b). The introns had lower GC ing the vasotocin receptor 2C gene (VTR2C; also known as AVPR2), which content, closer to the intergenic background. The intergenic GC content is involved in blood pressure homeostasis and brain function41,42, was was stable within 30 kb on either side of each gene (Fig. 4b). Mammals, split into 34 contigs on four scaffolds, one of which contained a false birds, and reptiles had the highest increase (around 20%) in GC content haplotype duplication of ARHGAP4 in the previous platypus assembly43 near the start site, followed by the amphibian and skate with medium (Fig. 3i). In our VGP assembly, all eleven genes were in one 37-Mb-long levels (around 10%). Teleost fishes showed an initial decrease, fol- contig within the approximately 50 Mb chromosome 6 scaffold. Fur- lowed by weaker increase (about 5%) from an already lower GC content thermore, eight of the eleven genes were remarkably increased in size (Fig. 4b). Given that the skate represents the sister branch to all other owing to the addition of previously unknown missing sequences. This vertebrate lineages sequenced, these findings suggest that teleosts chromosomal region was more GC-rich (54%) than the entire chromo- lost at least 5% GC content genome-wide, while maintaining most of the some 6 (46%). Thousands of such false gains and losses in previous GC content pattern in protein-coding genes. Although it is known that reference assemblies have been corrected in our VGP assemblies (more promoter regions can be CpG rich, and GC content can vary between 742 | Nature | Vol 592 | 29 April 2021 a Human (GRCh38) b Cretaceous–Palaeogene c 4 6 1 3 2 3 Chr. 6 p22.1 (28,238,097–28,930,763) boundary AN L N2 1 N3 AN N1 N2 9 27 Canada lynx SC AP CA BD CA SC CA CA X6 X5 ED IM 0.78 ZK NK ZS PG ZS ZK ZS ZS GP GP ZB TR Greater horseshoe bat Human (GRCh38) Boreoeutheria 0.31 Chr. B2 Pale spear-nosed bat 0.88 Canada lynx Platypus 0.31 Chr. 22 Chr. 3 Laurasiatheria Greater horseshoe bat Goode’s thornscrub tortoise 0.83 0.38 S17 S1 Anna’s hummingbird Yinpterochiroptera Egyptian fruit bat 1.00 1 Zebra finch (male) 0.26 Chr. 15 Chr. 4 Chiroptera Pale spear-nosed bat Kakapo 1.35 Yangochiroptera 0.78 S20 S8 Two-lined caecilian Velvety free-tailed bat 1.89 Blunt-snouted clingfish 0.45 S18 S56 S64 2.33 Greater mouse-eared bat Eastern happy Mammalia 1.89 0.39 S24 S96 Flier cichlid Kuhl’s pipistrelle Channel bull blenny 0.17 100 kb Platypus Zig-zag eel 0.34 Chromosome End of scaffold or Middle of scaffold or Climbing perch Chicken orientation chromosome chromosome Gene intact Pseudogene Gene completely lost Thorny skate 312 240 160 80 0 Divergence time (Myr) Fig. 5 | Chromosome evolution among bats and other vertebrates. (breaks per Myr) on different branches. Red, higher rates than average (0.84); a, Chromosome synteny maps across the species sequenced based on BUSCO blue, lower than average. c, Summary of alignment, gene organization, and gene alignments. Chromosome sizes (bar lengths) are normalized to genome functional gene status surrounding a bat interchromosomal EBR involving the size, to make visualization easier. Genes (lines) are coloured according to the homologue of human chromosome 6. End of scaffold (S) or chromosome (Chr.) human chromosome to which they belong; those on human chromosome 6 are means that the breakpoint is located at a chromosome arm end; middle means highlighted in blue and other chromosomes are in lighter shades. The that it is located within a scaffold or chromosome. Scale is relevant for human cladogram is from the TimeTree database72. b, Phylogenetic relationship of the Chr. 6 only. Actual gene sizes in the non-human species may differ and were mammalian species sequenced and their inferred chromosome EBR rates drawn to match the annotated human gene sizes for simplicity. exons and introns45,46, such a systematic pattern, the lineage-specific more massive chromosome rearrangements since divergence from differences within vertebrates, and the magnitude of these differences their most recent common ancestor with tetrapods, consistent with had not been previously described, to our knowledge. a proposed higher rearrangement rate in teleosts49. We tested whether the newly assembled GC-rich promoter regions To determine the precise locations of chromosome rearrangements contained novel regulatory sequences. Analysing the zebra finch brain, between species, we focused on a shorter evolutionary distance of we found that genes with upregulated expression specific to the stria- around 180 Myr among mammals, and added four additional bat spe- tum (for example, DRD1B, which encodes a dopamine receptor) had cies described in our Bat1K study50, the human genome reference51 ATAC-seq peaks in the GC-rich promoter and 5′ UTR region in striatal (GRCh38.p12), and a recently upgraded long-read chicken reference52 neurons, but not in arcopallium neurons (Fig. 4c); conversely, genes (galGal6a) as an outgroup. Pairwise whole-genome alignments to the (for example, the ER81 transcription factor) with upregulated expres- human reference defined homologous synteny blocks and evolu- sion in the arcopallium (mammalian cortex layer 5 equivalent47) had tionary breakpoint regions (EBRs) among the species. We found that ATAC–seq peaks in the GC-rich region in arcopallium neurons but not breakpoint rates (EBRs per Myr) tripled among bats soon after the striatal neurons (Fig. 4d). These GC-rich regions were missing in the last mass extinction event (about 66 million years ago (Mya)), a time earlier assembly. In addition, the missing region in DRD1B led to a false of rapid bat superfamily divergences53 (about 60 Mya; Fig. 5b). Some annotation as a two-exon gene48, whereas the VGP assembly revealed a rearrangements affected genes. For example, a 1.3-Mb inversion in single-exon gene (Fig. 4c). These GC-rich promoter regions are candi- greater horseshoe bat chromosome 28 (homologous to 29.5 Mb of dates for driving cell-type-specific expression. These findings demon- human chromosome 15; Extended Data Fig. 12a) disrupted STARD5, a strate the importance of using sequencing chemistry that reads through gene involved in cholesterol homeostasis in liver cells54. The rearrange- GC-rich regions, like the CLR method. The earlier hummingbird genome ment separated exons 1–5 from exon 6, and disrupted splicing of the assembly was generated using Illumina TruSeq3 chemistry16, which was transcripts (Extended Data Fig. 12b). Another example was an EBR that designed to read through GC-rich regions, and yet about 55% of the genes involved fission of an ancestral bat chromosome homologue of human were missing the 100-bp GC-rich region before the start site (Fig. 4a). chromosome 6 (boreoeutherian mammal chromosome 555) and was Another paper contains additional findings on missing regions27. later reused among the different bat lineages in rearrangements that involved the ancestral homologues of human chromosomes 1, 2 and 6 (Fig. 5c, Extended Data Fig. 12c). We also noted a fission in this region in Chromosomal evolution the mouse, rat, and dog genomes55. On the basis of the conserved gene We next investigated whether we could gain new insights into chromo- order in human and Canada lynx, we inferred that the boreoeutherian some evolution among vertebrates. Given the more than 430 million ancestral mammal locus corresponding to human 6p22.1 contained years (Myr) of evolutionary divergence among the species sampled 12 genes, including four ZSCAN and two ZKSCAN transcription fac- here, it was difficult to generate whole genome-to-genome align- tors, and two GPX enzyme genes, all associated with sequentially ments across all species. Thus, we focused our initial analyses on 1,147 increasing independent gene losses in bats (Fig. 5d). For example, the highly conserved BUSCO vertebrate genes that are shared among greater horseshoe bat lost only ZSCAN12 and GPX6 to pseudogeniza- our assemblies of all 16 species and the human reference (GRCh38). tion, whereas Kuhl’s pipistrelle lost all 12 genes. ZSCAN and ZKSCAN are Human chromosomes mapped with greater orthology to 3.7 ± 1.3 (s.d.) involved in cell differentiation, migration and invasion, proliferation, chromosomes on average in other mammals, compared to 5.6 ± 2.2 apoptosis, and innate immunity56. We speculate that loss of ZSCAN12 in in amphibians and 9.6 ± 3.3 in teleost fishes (Fig. 5a, Supplementary all six bats could contribute to their immune tolerance to pathogens50. Table 22). The skate chromosome arrangement was more conserved Other biological findings using these VGP assemblies are published with tetrapods, mapping to 2.9 ± 1.4 chromosomes on average, com- elsewhere, and include: 1) more accurate synteny across species, lead- pared to 4.8 ± 2.5 in teleost fishes. These findings indicate that, along ing to a better understanding of the evolution of and thus a univer- with a reduction in GC content, the teleost lineage has experienced sal nomenclature for the vasotocin (also known as vasopressin) and Nature | Vol 592 | 29 April 2021 | 743 Article Table 1 | Proposed standards and metrics for defining genome assembly quality Quality category Metric Finished VGP-2020 VGP-2016 B10k-2014 This study Notation x.y.P.Q.C c.c.Pc.Q60.C100 7.c.P6.Q50.C95 6.7.P5.Q40.C90 4.5.Q30 Continuity Contig NG50 (x) = Chr. NG50 >10 Mb >1 Mb >10 kb 1–25 Mb Scaffolds NG50 (y) = Chr. NG50 = Chr. NG50 >10 Mb >100 kb 23–480 Mb Gaps per Gb No gaps <200 <1,000 <10,000 75–1,500 Structural accuracy Reliable blocks = Chr. NG50 >10 Mb >1 Mb Not required 2.3–40.2 Mb False duplications 0% <1% <5% <10% 0.2–5.0% Curation Conflicts resolved Manual Manual Not required Manual Base accuracy Base pair QV (Q) >60 >50 >40 >30 39–43 k-mer completeness 100% complete >95% >90% >80% 87–98% Haplotype phasing Phase block NG50 (P) = Chr. NG50 >1 Mb >100 kb Not required 1.6 Mba Functional completeness Genes >98% complete >95% complete >90% >80% 82–98% Transcript mappability >98% >90% >80% >70% 96% Chromosome status Assigned (C) >100% >95% >90% Not required 94.4–99.9% Sex chromosomes Right order, no gaps Localized homo At least one shared Fragmented At least one pairs (for example, X or Z) shared Organelles (for example, One complete allele One complete Fragmented Not required One MT) allele complete allele The six broad quality categories in the first column are split into sub-metrics in the second column. The recommendations for draft to finished qualities (columns 3–6) are based on those achieved in past studies16,19,63, this study, and what we aspire to. In the x.y.P.Q.C notation, x = log10[contig NG50]; y = log10[scaffold NG50]; P = log10[haplotype phased NG50 block]; Q = Phred base accuracy QV; and C = percentage of the assembly assigned to chromosomes. c denotes ‘complete’ telomere-to-telomere continuity. The VGP assemblies (last column) satisfy the 6.7.6.Q40.C90 standard, but some come close to achieving a higher 7.c.7.Q50.C95 standard. These metrics apply to genomes about 1 Gb or bigger. a Phase blocks calculated for the zebra finch non-trio assembly using haplotype specific k-mers from parental data20; the trio assemblies had NG50 phase blocks of 17.3 Mb (maternal) and 56.6 Mb (paternal). oxytocin ligand and receptor gene families57; 2) greater understanding study chromosomal evolution, then the VGP-2016 minimum metric of the evolution of the carbohydrate 6-O sulfotransferase gene family, 6.7.P5.Q40.C95, with high structural and base accuracies and more which encodes enzymes that modify secreted carbohydrates58; 3) the than 95% assigned to chromosomes (Table 1), would be necessary. first Bat1K study50, which generated a genome-scale phylogeny that If having GC-rich promoter regions and complete 5′ exons in most better resolves the relationships between bats and other mammals, and genes is essential, then long-read approaches that sequence through which identified changes in bat genes that are involved in immunity and these regions are necessary. ‘Finished’ quality (Table 1) is obviously life span, including genes that are relevant to the COVID-19 pandemic59; the ideal assembly result, but this level of quality is currently routine 4) deleterious mutations that have been purged from the last surviving only for bacterial and non-vertebrate model organisms with smaller isolated and inbred population of the critically endangered kākāpō60; genome sizes that lack large centromeric satellite arrays65–67 and for and 5) more complete resolution of the evolution of the complex sex organelle genomes, as presented here33. The possibility of achieving chromosomes in platypus and echidna26. These discoveries were not complete, telomere-to-telomere assemblies of vertebrate and other possible with the previous reference assemblies, and we expect many eukaryotic species is foreseeable, given some assembled avian and bat future discoveries to follow. chromosomes with zero gaps in this study, and the recent complete assembly of two human chromosomes68,69. Proposed assembly quality metrics Drawing on the lessons learned from this work, we propose that assem- The Vertebrate Genomes Project bly quality should be summarized using 14 metrics under 6 categories Building on this initial set of assembled genomes and the lessons (Table 1; full details in Supplementary Note 4). We summarize the learned, we propose to expand the VGP to deeper taxonomic phases, most critical and commonly used metrics using the simple notation beginning with phase 1: representatives of approximately 260 verte- x.y.P.Q.C, where: x = log10[contig NG50], y = log10[scaffold NG50], brate orders, defined here as lineages separated by 50 million or more P = log10[haplotype phase block NG50], Q = QV base accuracy, and years of divergence from each other. Phase 2 will encompass species C = percentage of the assembly assigned to chromosomes (Table 1). that represent all approximately 1,000 vertebrate families; phase 3, Our current minimum VGP standard, for example, is 6.7.P5.Q40.C90. all roughly 10,000 genera; and phase 4, nearly all 71,657 extant named This revises our prior notation50,61,62, which reported log-scaled con- vertebrate species (Supplementary Note 5, Supplementary Fig. 3). To tinuity measured in ‘kilobases’ rather than ‘bases’. The thresholds accomplish such a project within 10 years, we will need to scale up to we chose were based on empirical and quantitative observations completing 125 genomes per week, without sacrificing quality. This between what is achievable currently and what is aspirational, and includes sample permitting, high molecular weight DNA extractions, the question the assemblies are meant to answer. For example, the sequencing, meta-data tracking, and computational infrastructure. We short-read paired-end library-based assemblies of the B10K Phase 1 will take advantage of continuing improvements in genome sequencing genomes in 201416 and the 10XG linked-read assembly of the Anna’s technology, assembly, and annotation, including advances in PacBio hummingbird presented here would be categorized as a 4.5.P7.Q50 HiFi reads, Oxford Nanopore reads, and replacements for 10XG reads assembly, with low continuity but high base accuracy (Table 1). Such (Supplementary Note 6), while addressing specific scientific questions a genome would be suitable for use in phylogenomics63 and for at increasing levels of phylogenetic refinement. Genomic technology population-scale SNP surveys64. If, instead, a genome is to be used to advances quickly, but we believe the principles of our pipeline and 744 | Nature | Vol 592 | 29 April 2021 the lessons learned will be applicable to future efforts. Areas in which 25. Howe, K. et al. Significantly improving the quality of genome assemblies through curation. Gigascience 10, giaa153 (2021). improvement is needed include more accurate and complete hap- 26. Zhou, Y. et al. Platypus and echidna genomes reveal mammalian biology and evolution. lotype phasing, base-call accuracy, and resolution of long repetitive Nature https://doi.org/10.1038/s41586-020-03039-0 (2021). regions such as telomeres, centromeres, and sex chromosomes. The 27. Kim, J. et al. False gene and chromosome losses affected by assembly and sequence errors. Preprint at https://doi.org/10.1101/2021.04.09.438906 (2021). VGP is working towards these goals and making all data, protocols, and 28. Lewin, H. A., Graves, J. A. M., Ryder, O. A., Graphodatsky, A. S. & O’Brien, S. J. Precision pipelines openly available (Supplementary Notes 5, 7). nomenclature for the new genomics. Gigascience 8, giz086 (2019). Despite remaining imperfections, our reference genomes are the 29. Kronenberg, Z. N. et al. Extended haplotype phasing of de novo genome assemblies with FALCON-Phase. Nat. Commun. https://doi.org/10.1038/s41467-020-20536-y (2021). most complete and highest quality to date for each species sequenced, 30. Ewing, B., Hillier, L., Wendl, M. C. & Green, P. Base-calling of automated sequencer traces to our knowledge. When we began to generate genomes beyond the using phred. I. Accuracy assessment. Genome Res. 8, 175–185 (1998). Anna’s hummingbird in 2017, only eight vertebrate species in Gen- 31. Tomaszkiewicz, M., Medvedev, P. & Makova, K. D. Y and W chromosome assemblies: approaches and discoveries. Trends Genet. 33, 266–282 (2017). Bank had genomes that met our target continuity metrics, and none 32. Kolesnikov, A. A. & Gerasimov, E. S. Diversity of mitochondrial genome organization. were haplotype phased (Supplementary Table 23). The VGP pipeline Biochem. (Mosc.) 77, 1424–1435 (2012). introduced here has now been used to complete assemblies of more 33. Formenti, G. et al. Complete vertebrate mitogenomes reveal widespread repeats and gene duplications. Genome Biol. (in the press). than 130 species of similar or higher quality (Supplementary Note 5; 34. Harrison, G. L. A. et al. Four new avian mitochondrial genomes help get to basic BioProject PRJNA489243). We encourage the scientific community to evolutionary questions in the late cretaceous. Mol. Biol. Evol. 21, 974–983 (2004). use and evaluate the assemblies and associated raw data, and to provide 35. Zhao, H. et al. The complete mitochondrial genome of the Anabas testudineus (Perciformes, Anabantidae). Mitochondrial DNA A DNA Mapp. Seq. Anal. 27, 1005–1007 feedback towards improving all processes for complete and error-free (2016). assembled genomes of all species. 36. Suzuki, A. et al. How the kinetochore couples microtubule force and centromere stretch to move chromosomes. Nat. Cell Biol. 18, 382–392 (2016). 37. Pfenning, A. R. et al. Convergent transcriptional specializations in the brains of humans and song-learning birds. Science 346, 1256846 (2014). Online content 38. Robinson, R. For mammals, loss of yolk and gain of milk went hand in hand. PLoS Biol. 6, e77 (2008). Any methods, additional references, Nature Research reporting sum- 39. Brandl, K. et al. Yip1 domain family, member 6 (Yipf6) mutation induces spontaneous maries, source data, extended data, supplementary information, intestinal inflammation in mice. Proc. Natl Acad. Sci. USA 109, 12650–12655 (2012). acknowledgements, peer review information; details of author con- 40. Malmstrøm, M. et al. Evolution of the immune system influences speciation rates in teleost fishes. Nat. Genet. 48, 1204–1210 (2016). tributions and competing interests; and statements of data and code 41. Japundžić-Žigon, N., Lozić, M., Šarenac, O. & Murphy, D. Vasopressin & oxytocin in control availability are available at https://doi.org/10.1038/s41586-021-03451-0. of the cardiovascular system: an updated review. Curr. Neuropharmacol. 18, 14–33 (2020). 42. Cataldo, I., Azhari, A. & Esposito, G. A review of oxytocin and arginine-vasopressin 1. International Human Genome Sequencing Consortium. Initial sequencing and analysis of receptors and their modulation of autism spectrum disorder. Front. Mol. Neurosci. 11, 27 the human genome. Nature 409, 860–921 (2001). (2018). 2. Sulston, J. et al. The C. elegans genome sequencing project: a beginning. Nature 356, 37– 43. Warren, W. C. et al. Genome analysis of the platypus reveals unique signatures of 41 (1992). evolution. Nature 453, 175–183 (2008). 3. Mouse Genome Sequencing Consortium. Initial sequencing and comparative analysis of 44. Ko, B. J. et al. Widespread false gene gains caused by duplication errors in genome the mouse genome. Nature 420, 520–562 (2002). assemblies. Preprint at https://doi.org/10.1101/2021.04.09.438957 (2021). 4. Howe, K. et al. The zebrafish reference genome sequence and its relationship to the 45. Lemaire, S. et al. Characterizing the interplay between gene nucleotide composition bias human genome. Nature 496, 498–503 (2013). and splicing. Genome Biol. 20, 259 (2019). 5. Genome 10K Community of Scientists. Genome 10K: a proposal to obtain whole-genome 46. Zhang, L., Kasif, S., Cantor, C. R. & Broude, N. E. GC/AT-content spikes as genomic sequence for 10,000 vertebrate species. J. Hered. 100, 659–674 (2009). punctuation marks. Proc. Natl Acad. Sci. USA 101, 16855–16860 (2004). 6. Koepfli, K.-P., Paten, B., the Genome 10K Community of Scientists & O’Brien, S. J. The 47. Jarvis, E. D. et al. Global view of the functional molecular organization of the avian Genome 10K Project: a way forward. Annu. Rev. Anim. Biosci. 3, 57–111 (2015). cerebrum: mirror images and functional columns. J. Comp. Neurol. 521, 3614–3665 (2013). 7. Venter, J. C. et al. The sequence of the human genome. Science 291, 1304–1351 (2001). 48. Kubikova, L., Wada, K. & Jarvis, E. D. Dopamine receptors in a songbird brain. J. Comp. 8. Adams, M. D. et al. The genome sequence of Drosophila melanogaster. Science 287, Neurol. 518, 741–769 (2010). 2185–2195 (2000). 49. Sémon, M. & Wolfe, K. H. Rearrangement rate following the whole-genome duplication in 9. Shendure, J. & Ji, H. Next-generation DNA sequencing. Nat. Biotechnol. 26, 1135–1145 teleosts. Mol. Biol. Evol. 24, 860–867 (2007). (2008). 50. Jebb, D. et al. Six reference-quality genomes reveal evolution of bat adaptations. Nature 10. Yin, Z.-T. et al. Revisiting avian ‘missing’ genes from de novo assembled transcripts. BMC 583, 578–584 (2020). Genomics 20, 4 (2019). 51. Schneider, V. A. et al. Evaluation of GRCh38 and de novo haploid genome assemblies 11. Korlach, J. et al. De novo PacBio long-read and phased avian genome assemblies correct demonstrates the enduring quality of the reference assembly. Genome Res. 27, 849–864 and add to reference genes generated with intermediate and short reads. Gigascience 6, (2017). 1–16 (2017). 52. Warren, W. C. et al. A new chicken genome assembly provides insight into avian genome 12. Kelley, D. R. & Salzberg, S. L. Detection and correction of false segmental duplications structure. G3 (Bethesda) 7, 109–117 (2017). caused by genome mis-assembly. Genome Biol. 11, R28 (2010). 53. Meredith, R. W. et al. Impacts of the Cretaceous Terrestrial Revolution and KPg extinction 13. Roach, M. J., Schmidt, S. A. & Borneman, A. R. Purge Haplotigs: allelic contig on mammal diversification. Science 334, 521–524 (2011). reassignment for third-gen diploid genome assemblies. BMC Bioinformatics 19, 460 54. Rodriguez-Agudo, D. et al. StarD5: an ER stress protein regulates plasma membrane and (2018). intracellular cholesterol homeostasis. J. Lipid Res. 60, 1087–1098 (2019). 14. Guan, D. et al. Identifying and removing haplotypic duplication in primary genome 55. Kim, J. et al. Reconstruction and evolutionary history of eutherian chromosomes. Proc. assemblies. Bioinformatics 36, 2896–2898 (2020). Natl Acad. Sci. USA 114, E5379–E5388 (2017). 15. Bradnam, K. R. et al. Assemblathon 2: evaluating de novo methods of genome assembly 56. Lin, B., Dutta, B. & Fraser, I. D. C. Systematic investigation of multi-TLR sensing identifies in three vertebrate species. Gigascience 2, 10 (2013). regulators of sustained gene activation in macrophages. Cell Syst. 5, 25–37.e3 (2017). 16. Zhang, G. et al. Comparative genomics reveals insights into avian genome evolution and 57. Theofanopoulou, C., Gedman, G. L., Cahill, J. A., Boeckx, C. & Jarvis, E. D. Universal adaptation. Science 346, 1311–1320 (2014). nomenclature for oxytocin-vasotocin ligand and receptor families. Nature https://doi. 17. Chin, C.-S. et al. Phased diploid genome assembly with single-molecule real-time org/10.1038/s41586-020-03040-7 (2021). sequencing. Nat. Methods 13, 1050–1054 (2016). 58. Ocampo Daza, D. & Haitina, T. Reconstruction of the carbohydrate 6-O sulfotransferase 18. Bresler, G., Bresler, M. & Tse, D. Optimal assembly for high throughput shotgun gene family evolution in vertebrates reveals novel member, CHST16, lost in amniotes. sequencing. BMC Bioinformatics 14 (Suppl. 5), S18 (2013). Genome Biol. Evol. 12, 993–1012 (2020). 19. Warren, W. C. et al. The genome of a songbird. Nature 464, 757–762 (2010). 59. Damas, J. et al. Broad host range of SARS-CoV-2 predicted by comparative and structural 20. Koren, S. et al. De novo assembly of haplotype-resolved genomes with trio binning. analysis of ACE2 in vertebrates. Proc. Natl Acad. Sci. USA 117, 22311–22322 (2020). Nat. Biotechnol. (2018). 60. Dussex, N. et al. Population genomics reveals the impact of long-term small population 21. Koren, S., Phillippy, A. M., Simpson, J. T., Loman, N. J. & Loose, M. Reply to ‘Errors in size in the critically endangered kākāpō. Cell Genom. (in the press). long-read assemblies can critically affect protein prediction’. Nat. Biotechnol. 37, 127–128 61. Teeling, E. C. et al. Bat biology, genomes, and the Bat1K project: to generate (2019). chromosome-level genomes for all living bat species. Annu. Rev. Anim. Biosci. 6, 23–46 22. Vollger, M. R. et al. Long-read sequence and assembly of segmental duplications. (2018). Nat. Methods 16, 88–94 (2019). 62. Lewin, H. A. et al. Earth BioGenome Project: sequencing life for the future of life. Proc. 23. Rhie, A., Walenz, B. P., Koren, S. & Phillippy, A. M. Merqury: reference-free quality, Natl Acad. Sci. USA 115, 4325–4333 (2018). completeness, and phasing assessment for genome assemblies. Genome Biol. 21, 245 63. Jarvis, E. D. et al. Whole-genome analyses resolve early branches in the tree of life of (2020). modern birds. Science 346, 1320–1331 (2014). 24. Waterhouse, R. M. et al. BUSCO applications from quality assessments to gene prediction 64. Li, S. et al. Genomic signatures of near-extinction and rebirth of the crested ibis and other and phylogenomics. Mol. Biol. Evol. 35, 543–548 (2018). endangered bird species. Genome Biol. 15, 557 (2014). Nature | Vol 592 | 29 April 2021 | 745 Article 65. Koren, S. & Phillippy, A. M. One chromosome, one contig: complete microbial genomes Massachusetts Cooperative Fish and Wildlife Research Unit, Amherst, MA, USA. 22School of from long-read sequencing and assembly. Curr. Opin. Microbiol. 23, 110–120 (2015). Biological Science, The Environment Institute, University of Adelaide, Adelaide, South 66. Jenjaroenpun, P. et al. Complete genomic and transcriptional landscape analysis using Australia, Australia. 23Bond Life Sciences Center, University of Missouri, Columbia, MO, USA. third-generation sequencing: a case study of Saccharomyces cerevisiae CEN.PK113-7D. 24 Department of Biology, East Carolina University, Greenville, NC, USA. 25UQ Genomics, Nucleic Acids Res. 46, e38 (2018). University of Queensland, Brisbane, Queensland, Australia. 26Department of Biological 67. Tyson, J. R. et al. MinION-based long-read sequencing and assembly extends the Sciences, Clemson University, Clemson, SC, USA. 27The Genetic Rescue Foundation, Caenorhabditis elegans reference genome. Genome Res. 28, 266–274 (2018). 68. Miga, K. H. et al. Telomere-to-telomere assembly of a complete human X chromosome. Wellington, New Zealand. 28Kākāpō Recovery, Department of Conservation, Invercargill, New Nature 585, 79–84 (2020). Zealand. 29Department of Zoology, University of Otago, Dunedin, New Zealand. 30University of 69. Logsdon, G. A. et al. The structure, function and evolution of a complete human Arizona Genetics Core, Tucson, AZ, USA. 31Department of Life Sciences, Natural History chromosome 8. Nature https://doi.org/10.1038/s41586-021-03420-7 (2021). Museum, London, UK. 32School of Natural Sciences, Bangor University, Gwynedd, UK. 70. Beçak, M. L., Beçak, W., Roberts, F. L., Shoffner, R. N. & Volpe, P. (eds.) Chromosome Atlas: 33 Department of Biology, University of Konstanz, Konstanz, Germany. 34Department of Fish, Amphibians, Reptiles, and Birds Vol. 2 (Springer, 1973). Organismic and Evolutionary Biology, Harvard University, Cambridge, MA, USA. 35Department 71. Vurture, G. W. et al. GenomeScope: fast reference-free genome profiling from short of Marine and Environmental Sciences, Northeastern University Marine Science Center, reads. Bioinformatics 33, 2202–2204 (2017). Nahant, MA, USA. 36Department of Biology, University of Antwerp, Antwerp, Belgium. 72. Kumar, S., Stecher, G., Suleski, M. & Hedges, S. B. TimeTree: a resource for timelines, timetrees, and divergence times. Mol. Biol. Evol. 34, 1812–1819 (2017). 37 Naturalis Biodiversity Center, Leiden, The Netherlands. 38Institute of Biology, Karl-Franzens University of Graz, Graz, Austria. 39Florida Museum of Natural History, University of Florida, Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Gainesville, FL, USA. 40Center for Systems Biology, Dresden, Germany. 41Zoological Institute, published maps and institutional affiliations. University of Basel, Basel, Switzerland. 42Tag.bio, San Francisco, CA, USA. 43UC Santa Cruz Genomics Institute, University of California, Santa Cruz, CA, USA. 44San Diego Zoo Global, Open Access This article is licensed under a Creative Commons Attribution Escondido, CA, USA. 45Pacific Biosciences, Menlo Park, CA, USA. 46Digital BioLogic, 4.0 International License, which permits use, sharing, adaptation, distribution Ivanić-Grad, Croatia. 47Bionano Genomics, San Diego, CA, USA. 48Arima Genomics, San Diego, and reproduction in any medium or format, as long as you give appropriate CA, USA. 49Dovetail Genomics, Santa Cruz, CA, USA. 50Independent Researcher, Santa Cruz, credit to the original author(s) and the source, provide a link to the Creative Commons license, CA, USA. 51CNAG-CRG, Centre for Genomic Regulation, Barcelona Institute of Science and and indicate if changes were made. The images or other third party material in this article are Technology, Barcelona, Spain. 52Universitat Pompeu Fabra, Barcelona, Spain. 53Department of included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your Computer Science, University of Maryland College Park, College Park, MD, USA. 54School of intended use is not permitted by statutory regulation or exceeds the permitted use, you will Computer Science and Technology, Center for Bioinformatics, Harbin Institute of Technology, need to obtain permission directly from the copyright holder. To view a copy of this license, Harbin, China. 55Department of Psychology, Institute for Mind and Biology, University of visit http://creativecommons.org/licenses/by/4.0/. Chicago, Chicago, IL, USA. 56Department of Genetics and Biochemistry, Clemson University, Clemson, SC, USA. 57Department of Behavioral Neuroscience, Oregon Health and Science © The Author(s) 2021 University, Portland, OR, USA. 58Max Planck Institute for the Physics of Complex Systems, Dresden, Germany. 59Monash University Malaysia Genomics Facility, School of Science, Arang Rhie1,103, Shane A. McCarthy2,3,103, Olivier Fedrigo4,103, Joana Damas5, Selangor Darul Ehsan, Malaysia. 60Tropical Medicine and Biology Multidisciplinary Platform, Giulio Formenti4,6, Sergey Koren1, Marcela Uliano-Silva7,8, William Chow3, Monash University Malaysia, Selangor Darul Ehsan, Malaysia. 61Qatar Falcon Genome Project, Arkarachai Fungtammasan9, Juwan Kim10, Chul Lee10, Byung June Ko11, Mark Chaisson12, Doha, Qatar. 62Department of Biosciences, University of Milan, Milan, Italy. 63eGnome, Inc., Gregory L. Gedman6, Lindsey J. Cantin6, Francoise Thibaud-Nissen13, Leanne Haggerty14, Seoul, Republic of Korea. 64LOEWE Centre for Translational Biodiversity Genomics, Frankfurt, Iliana Bista2,3, Michelle Smith3, Bettina Haase4, Jacquelyn Mountcastle4, Sylke Winkler15,16, Sadye Paez4,6, Jason Howard17, Sonja C. Vernes18,19,20, Tanya M. Lama21, Frank Grutzner22, Germany. 65Senckenberg Research Institute, Frankfurt, Germany. 66Goethe-University, Faculty Wesley C. Warren23, Christopher N. Balakrishnan24, Dave Burt25, Julia M. George26, of Biosciences, Frankfurt, Germany. 67BGI-Shenzhen, Shenzhen, China. 68Department of Matthew T. Biegler6, David Iorns27, Andrew Digby28, Daryl Eason28, Bruce Robertson29, Biology, Pennsylvania State University, University Park, PA, USA. 69Center for Medical Taylor Edwards30, Mark Wilkinson31, George Turner32, Axel Meyer33, Andreas F. Kautt33,34, Genomics, Pennsylvania State University, University Park, PA, USA. 70Center for Paolo Franchini33, H. William Detrich III35, Hannes Svardal36,37, Maximilian Wagner38, Computational Biology and Bioinformatics, Pennsylvania State University, University Park, PA, Gavin J. P. Naylor39, Martin Pippel15,40, Milan Malinsky3,41, Mark Mooney42, Maria Simbirsky9, USA. 71Department of Computer Science and Engineering, Pennsylvania State University, Brett T. Hannigan9, Trevor Pesout43, Marlys Houck44, Ann Misuraca44, Sarah B. Kingan45, University Park, PA, USA. 72Department of Biochemistry and Molecular Biology, Pennsylvania Richard Hall45, Zev Kronenberg45, Ivan Sović45,46, Christopher Dunn45, Zemin Ning3, State University, University Park, PA, USA. 73Hoonygen, Seoul, Korea. 74Department of Alex Hastie47, Joyce Lee47, Siddarth Selvaraj48, Richard E. Green43,49, Nicholas H. Putnam50, Migration, Max Planck Institute of Animal Behavior, Radolfzell, Germany. 75Department of Ivo Gut51,52, Jay Ghurye49,53, Erik Garrison43, Ying Sims3, Joanna Collins3, Sarah Pelan3, James Torrance3, Alan Tracey3, Jonathan Wood3, Robel E. Dagnew12, Dengfeng Guan2,54, Biological Sciences, Universidad de los Andes, Bogotá, Colombia. 76Center for Evolutionary Sarah E. London55, David F. Clayton56, Claudio V. Mello57, Samantha R. Friedrich57, Hologenomics, The GLOBE Institute, University of Copenhagen, Copenhagen, Denmark. Peter V. Lovell57, Ekaterina Osipova15,40,58, Farooq O. Al-Ajli59,60,61, Simona Secomandi62, 77 University Museum, NTNU, Trondheim, Norway. 78China National Genebank, BGI-Shenzhen, Heebal Kim10,11,63, Constantina Theofanopoulou6, Michael Hiller64,65,66, Yang Zhou67, Shenzhen, China. 79Villum Center for Biodiversity Genomics, Section for Ecology and Robert S. Harris68, Kateryna D. Makova68,69,70, Paul Medvedev69,70,71,72, Jinna Hoffman13, Evolution, Department of Biology, University of Copenhagen, Copenhagen, Denmark. 80State Patrick Masterson13, Karen Clark13, Fergal Martin14, Kevin Howe14, Paul Flicek14, Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Brian P. Walenz1, Woori Kwak63,73, Hiram Clawson43, Mark Diekhans43, Luis Nassar43, Academy of Sciences, Kunming, China. 81Center for Excellence in Animal Evolution and Benedict Paten43, Robert H. S. Kraus33,74, Andrew J. Crawford75, M. Thomas P. Gilbert76,77, Genetics, Chinese Academy of Sciences, Kunming, China. 82Institute of Molecular and Cell Guojie Zhang78,79,80,81, Byrappa Venkatesh82, Robert W. Murphy83, Klaus-Peter Koepfli84, Biology, A*STAR, Biopolis, Singapore, Singapore. 83Centre for Biodiversity, Royal Ontario Beth Shapiro85,86, Warren E. Johnson84,87,88, Federica Di Palma89, Tomas Marques-Bon et90,91,92,93, Emma C. Teeling94, Tandy Warnow95, Jennifer Marshall Graves96, Museum, Toronto, Ontario, Canada. 84Smithsonian Conservation Biology Institute, Center for Oliver A. Ryder44,97, David Haussler43,85, Stephen J. O’Brien98,99, Jonas Korlach45, Species Survival, National Zoological Park, Washington, DC, USA. 85Department of Ecology Harris A. Lewin5,100,101, Kerstin Howe3,104 ✉, Eugene W. Myers15,40,102,104 ✉, and Evolutionary Biology, University of California Santa Cruz, Santa Cruz, CA, USA. 86Howard Richard Durbin2,3,104 ✉, Adam M. Phillippy1,104 ✉ & Erich D. Jarvis4,6,86,104 ✉ Hughes Medical Institute, Chevy Chase, MD, USA. 87The Walter Reed Biosystematics Unit, Museum Support Center MRC-534, Smithsonian Institution, Suitland, MD, USA. 88Walter Reed Army Institute of Research, Silver Spring, MD, USA. 89Department of Biological Sciences, 1 Genome Informatics Section, Computational and Statistical Genomics Branch, National Earlham Institute, University of East Anglia, Norwich, UK. 90Institute of Evolutionary Biology Human Genome Research Institute, National Institutes of Health, Bethesda, MD, USA. (UPF-CSIC), PRBB, Barcelona, Spain. 91Catalan Institution of Research and Advanced Studies 2 Department of Genetics, University of Cambridge, Cambridge, UK. 3Wellcome Sanger (ICREA), Barcelona, Spain. 92Centre for Genomic Regulation (CRG), Barcelona Institute of Institute, Cambridge, UK. 4Vertebrate Genome Lab, The Rockefeller University, New York, NY, Science and Technology (BIST), Barcelona, Spain. 93Institut Català de Paleontologia Miquel USA. 5The Genome Center, University of California Davis, Davis, CA, USA. 6Laboratory of Crusafont, Universitat Autònoma de Barcelona, Barcelona, Spain. 94School of Biology and Neurogenetics of Language, The Rockefeller University, New York, NY, USA. 7Leibniz Institute Environmental Science, University College Dublin, Dublin, Ireland. 95Department of Computer for Zoo and Wildlife Research, Department of Evolutionary Genetics, Berlin, Germany. 8Berlin Science, The University of Illinois at Urbana-Champaign, Urbana, IL, USA. 96School of Life Center for Genomics in Biodiversity Research, Berlin, Germany. 9DNAnexus Inc., Mountain Science, La Trobe University, Melbourne, Victoria, Australia. 97Department of Evolution, View, CA, USA. 10Interdisciplinary Program in Bioinformatics, Seoul National University, Seoul, Behavior, and Ecology, University of California San Diego, La Jolla, CA, USA. 98Laboratory of Republic of Korea. 11Department of Agricultural Biotechnology and Research Institute of Genomics Diversity-Center for Computer Technologies, ITMO University, St. Petersburg, Agriculture and Life Sciences, Seoul National University, Seoul, Republic of Korea. 12University Russian Federation. 99Guy Harvey Oceanographic Center, Halmos College of Natural of Southern California, Los Angeles, CA, USA. 13National Center for Biotechnology Sciences and Oceanography, Nova Southeastern University, Fort Lauderdale, FL, USA. Information, National Library of Medicine, NIH, Bethesda, MD, USA. 14European Molecular 100 Department of Evolution and Ecology, University of California Davis, Davis, CA, USA. 101John Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, Muir Institute for the Environment, University of California Davis, Davis, CA, USA. 102Faculty of UK. 15Max Planck Institute of Molecular Cell Biology and Genetics, Dresden, Germany. Computer Science, Technical University Dresden, Dresden, Germany. 103These authors 16 DRESDEN-concept Genome Center, Dresden, Germany. 17Novogene, Durham, NC, USA. contributed equally: Arang Rhie, Shane A. McCarthy, Olivier Fedrigo. 104These authors jointly 18 Neurogenetics of Vocal Communication Group, Max Planck Institute for Psycholinguistics, supervised this work: Kerstin Howe, Eugene W. Myers, Richard Durbin, Adam M. Phillippy, Nijmegen, The Netherlands. 19Donders Institute for Brain, Cognition and Behaviour, Nijmegen, Erich D. Jarvis. ✉e-mail: kj2@sanger.ac.uk; gene@mpi-cbg.de; rd109@cam.ac.uk; adam. The Netherlands. 20School of Biology, University of St Andrews, St Andrews, UK. 21University of phillippy@nih.gov; ejarvis@rockefeller.edu 746 | Nature | Vol 592 | 29 April 2021 Methods A for 1 h at 37 °C. HMW gDNA was purified with two washes of phe- nol–chloroform-IAA equilibrated to pH 8.0, followed by two washes of Genome assembly naming chloroform-IAA, and precipitated in ice-cold 100% ethanol. Filamen- For each completed assembly of an individual, we gave that assembly tous HMW gDNA was either spooled with shepherds hooks or collected an abbreviated name with the following rules: Lineage/GenusSpecies/ by centrifugation. HMW gDNA was washed twice with 70% ethanol, Individual#.Assembly#. The first letter, in lowercase, identifies the dried for 20 min at room temperature and eluted in TE. For the flier particular lineage: m, mammals; b, birds; r, reptiles; a, amphibians; f, cichlid muscle gDNA sample used for PacBio CLR and 10XG libraries, teleost fish; and s, sharks and other cartilaginous fishes. The next three glycogen was precipitated by adding 1/10 (v/v) 0.3 M sodium acetate, letters (first in caps) identify the species scientific genus name; the pH 6.0 to the extracted genomic DNA, mixing carefully and spinning at next three letters (first in caps) identifies the specific species name. room temperature at 10,000g. PFGE revealed thatDNA molecule length In the last position is the genome identifier, where integers (1, 2, 3, …) was between 50 and 300 kb—often lower in size than that obtained with represent different individuals of the same species, and decimals (1.1, the agarose plug but sufficient for long-range sequencing of CLR and 1.2, 1.3, …) represent different assemblies of the same individual. For linked read data types. example, the first submission of the curated Anna’s hummingbird (Calypte anna) assembly is bCalAnn1.1, and an updated assembly for Others. We also used the Qiagen MagAttract HMW DNA kit (cat no. the same individual is bCalAnn1.2. When the abbreviated lineage or 67563) and the KingFisher Cell and Tissue DNA kit (Thermo Scientific; genus and species names for two or more species were identical, we cat no. 97030196), following the manufacturers’ guidelines. These replaced the subsequent letters (fourth, fifth and so on) of the genus protocols yielded HMW DNA ranging from 30 to 50 kb. The Genomic or species name until they could be differentiated. We have created Tip (Qiagen) kit was also used for tissue-based extraction of HMW DNA. abbreviated names for all 71,657 vertebrate species (http://vgpdb.snu. ac.kr/splist/; https://id.tol.sanger.ac.uk/). Libraries and sequencing PacBio libraries and sequencing. DNA obtained from agarose plugs Sample collection was sheared down to ~40 kb fragment size with a MegaRuptor device The production of high-quality genome assemblies required us to obtain (Diagenode, Belgium) and fragmented using Covaris g-tubes (520079) high-quality cells or tissue that would yield high-molecular-weight or by needle shearing. PacBio large insert libraries were prepared with (HMW) DNA for long-read sequencing technologies (CLR and ONT) either the SMRTbell Template Prep Kit 1.0‐SPv3 (no.100‐991‐900) or and optical mapping (Bionano). Therefore, we obtained fresh-frozen the SMRTbell Express Template Prep Kit v1 (no. 101‐357‐000). Libraries samples of various tissues (Supplementary Table 8). All samples were were size-selected between 12 and 25 kb using Sage BluePippin (Sage obtained according to approved protocols of the respective animal Science, USA), depending on the DNA quality and extraction method. care and use committees or permits obtained by the respective persons These libraries were sequenced on either RSII or Sequel I instruments, at and institutions listed in Supplementary Table 8. Additional details of least 60× coverage per species using Sequel Binding Kit and Sequencing the samples are on their respective BioSample pages (https://www. Plate versions 2.0 and 2.1 with 10-h movie time (Supplementary Table 9). ncbi.nlm.nih.gov/biosample; accession numbers in Supplementary Table 8). All tissue types tested yielded a sufficient quantity and quality 10X Chromium libraries and sequencing. Unfragmented HMW DNA of DNA for sequencing and assembly, but we found that blood worked from the agarose plugs was used to generate linked read libraries on best for species that have nucleated red blood cells (that is, bird and the 10X Genomics Chromium platform (Genome Library Kit & Gel Bead reptiles), and spleen or cultured cells worked best for mammals, as of Kit v2 PN-120258, Genome Chip Kit v2 PN-120257, i7 Multiplex Kit PN- to date. Analysis of different tissue types will be presented elsewhere 120262) following the manufacturer’s guidelines. We sequenced the (in preparation). 10X libraries at ~60× coverage per species on an Illumina NovaSeq S4 150-bp PE lane. Isolation of high-molecular-weight DNA Agarose plug DNA isolation. For tissue, HMW DNA was extracted using Bionano libraries and optical map imaging. Unfragmented the Bionano animal tissue DNA isolation fibrous tissue protocol (cat ultra-HMW DNA from the agarose plugs was labelled using either two no. RE-013-10; document number 30071), according to the manufac- different nicking enzymes (BspQI and BssSI) or a direct labelling en- turer’s guidelines. A total of 25–30 mg was fixed in 2% formaldehyde zyme (DLE1) following the Bionano Prep Labelling NLRS (document and homogenized using the Qiagen TissueRuptor or manual tissue number 30024) and DLS protocols, respectively (document number disruption. For nucleated blood, 27–54 μl was used with an adapted 30206). Labelled samples were then imaged on a Bionano Irys or on a protocol (Bionano, personal communication) of the Bionano Prep Bionano Saphyr instrument. For all species, we aimed for at least 100× Blood and Cell Culture DNA Isolation Kit (cat no. RE-130-10). Lysates coverage per label (Supplementary Table 9). were embedded into agarose plugs and treated with Proteinase K and RNase A. Plugs were then purified by drop dialysis with 1× TE. DNA Hi-C libraries and sequencing. Chromatin interaction (Hi-C) librar- quality was assessed using pulse field gel electrophoresis (PFGE) (Pip- ies were generated using either Arima Genomics, Dovetail Genom- pin Pulse, SAGE Science, Beverly, MA) or the Femto Pulse instrument ics, or Phase libraries on muscle, blood, or other tissue with in vivo (Agilent). PFGE revealed that we isolated ultra-high-molecular-weight cross-linking (Supplementary Table 9) and sequenced on Illumina in- DNA between ~100 and ~500 kb long. struments. Arima-HiC preparations were performed by Arima Genom- ics (https://arimagenomics.com/) using the Arima-HiC kit that uses two Phenol–chloroform gDNA extraction. For some samples, we per- enzymes (P/N: A510008). The resulting Arima-HiC proximally ligated formed phenol–chloroform extractions for HMW gDNA. Snap-frozen DNA was then sheared, size-selected around 200–600 bp using SPRI tissue was pulverized into a fine powder with a mortar and pestle in beads, and enriched for biotin-labelled proximity-ligated DNA using liquid nitrogen. The powdered tissue was lysed overnight at 55 °C in streptavidin beads. From these fragments, Illumina-compatible librar- high-salt tissue lysis buffer (400 mM NaCl, 20 mM Tris base (pH 8.0), ies were generated using the KAPA Hyper Prep kit (P/N: KK8504). The 30 mM EDTA (pH 8.0), 0.5% SDS, 100 μg/ml Proteinase K), and pow- resulting libraries were PCR amplified and purified with SPRI beads. The dered lung tissue was lysed overnight in Qiagen G2 lysis buffer (cat quality of the final libraries was checked with qPCR and Bioanalyzer, no. 1014636, Qiagen, Hilden, Germany) containing 100 μg/ml Pro- and then sequenced on Illumina HiSeq X at ~60× coverage following the teinase K at 55 °C. RNA was removed by incubation in 50 μg/ml RNase manufacturer’s protocols. Dovetail-HiC preparations were performed Article by Dovetail using a single-enzyme (DpnII) proximity ligation approach. support. For conservatively filling sequences, we compared different Phase-HiC libraries were made by Phase Genomics using a Proximo Hi-C parameters in output stage with --minreads 1 and --minreads 4 in ad- Library single-enzyme reaction. dition to no restrictions. We found that the number of gaps closed was similar to the gaps filled with Arrow76 (Supplementary Table 4) and Quality control chose not to run PBJelly82 for future assemblies. Before we performed any assembly, all genomic data of all data types from each sample were used to screen potential outlier libraries, outlier Short-read polishing. Illumina polishing benchmarking was per- sequencing runs, or accidental species contamination with Mash73 by formed using Longranger83 2.1.3 and Pilon84 1.21 with --fix bases, local measuring sequence similarity (Supplementary Fig. 4). When run- option (Supplementary Table 5). Later, for the VGP pipeline, we used ning Mash, we used 21-mers to generate sketches with sketch size of FreeBayes85 as Pilon84 was not computationally scalable for large ge- 10,000 and compared among each sequencing run, and then differ- nomes with the updated Longranger 2.2.2. ences assessed between sequencing sets. Base-level accuracy estimate. Base-level accuracy was measured Genome size, repeat content, and heterozygosity estimations using a mapping-based approach and later using the k-mer-based ap- These estimations were made with k-mer-based methods applied to the proach23. To determine the number of rounds to polish, we used Illu- Illumina short reads obtained from 10XG linked sequencing libraries. mina paired-end reads from the hummingbird16. After trimming off barcodes during scaff10x74 preprocessing, canonical 31-mer counts were collected using Meryl23. With the resulting 31-mer Mis-joins and missed-joins. The curated hummingbird assembly was histogram, GenomeScope71 was used to estimate the haploid genome mapped to the target assemblies with MashMap286 with --filter_mode length, repeat content, and heterozygosity. The thorny skate linked one-to-one --pi 95 using 5 kb segments (-s 5000) for CLR assemblies read data failed quality control, which we suspect was due to low com- and 1 kb (-s 1000) for SR assemblies to compensate for the shorter con- plexity sequences from the high repeat content (54.1%) of the genome; tig sizes, as contigs smaller than a segment size will be excluded from so k-mers were collected later from Illumina whole-genome sequenc- the alignment. The number of mis-joins and missed joins were identified ing reads instead. The genome size and repeat content of the channel using the assembly_comparison.pl used in the ‘Curation’ section below bull blenny were estimated from an alternative method that looks at (Supplementary Methods, Supplementary Fig. 5). the mode of long read overlap coverage and WindowMasker75, as the estimated genome size from GenomeScope was almost doubling the VGP standard genome assembly pipeline 1.0 to 1.6 known haploid genome size (1.29 Gb versus 0.6 Gb) and repeat content All 17 genomes were assembled with the VGP pipeline (Extended Data (28.0% versus 58.0%), for reasons related to either the quality of the Fig. 2a) for benchmark purposes, with some uncurated. The pale 10X data or species differences. spear-nosed bat, greater horseshoe bat, Canada lynx, platypus, male and female zebra finch, kākāpō, Anna’s hummingbird, Goode’s thorn- Benchmarking assembly steps with the Anna’s hummingbird scrub tortoise, flier cichlid, and blunt-snouted clingfish assemblies were To develop the VGP standard pipeline, we compared various scaffold- generated using the VGP pipeline 1.0 to 1.6 and curated for submission ing, gap filling, and polishing tools. Default options were used unless to NCBI and EBI public archives. The curated and submitted two-lined otherwise noted. Detailed software versions are listed in Supplemen- caecilian, zig-zag eel, climbing perch, channel bull blenny, eastern tary Table 2. happy, and thorny skate assemblies were generated using a similar process developed in parallel (Supplementary Note 2). Two submitted Contigging and scaffolding. FALCON76 and FALCON-Unzip17 (smr- curated versions of the female zebra finch were made, one using the tanalysis 3.0.0) were used to generate contigs that used CLR. Canu77 standard VGP pipeline and the other using the VGP trio pipeline, so that 1.5+67 was used to generate the combined PacBio CLR and Oxford comparative analyses could be performed by others. Nanopore ONT assembly. To benchmark scaffolding with linked reads, we used scaff10x74 2.0. For the linked read-only assembly, Supernova Contigging. For PacBio data, contigs were generated from subreads 278 was used. For the optical maps, two-enzyme hybrid scaffolding using FALCON76 and FALCON-Unzip17, with one round of Arrow pol- was used in the Bionano Solve v3.2.1 software, using BspQI and BssSI ishing (smrtanalysis 5.1.0.26412). A minimum read length of 2 kb or a initially, as well as DLE1 later when the technology was developed. For cutoff at which reads longer than the cutoff include 50× coverage was benchmarking Hi-C in scaffolding, Salsa 2.279 was used for scaffolding used, whichever was longer. For calculating read coverage, we used results in Fig. 1a, with Hi-C reads generated from Arima Genomics. estimated genome size from http://www.genomesize.com/ when avail- Additional comparisons for the Hi-C libraries were performed using able, or from the literature (Supplementary Table 11) while waiting for assemblies provided by Dovetail Genomics and Phase Genomics (Sup- 10XG sequencing to estimate genome size using k-mers. FALCON and plementary Table 3). We used Hi-C from Arima Genomics as it had the FALCON-Unzip were run with default parameters, except for computing smallest number of PCR duplicates and better coverage for short and the overlaps. Raw read overlaps were computed with DALIGNER param- long interactions at the time of comparison (Supplementary Fig. 1). eters -k14 -e0.75 -s100 -l2500 -h240 -w8 to better reflect the higher Assembly statistics from HiRise, Proximo HiC, 3D-DNA80 and Arima error rate in early PacBio sequel I and II. Pread (preassembled read) Hi-C are available in Supplementary Table 3. We concluded that all Hi-C overlaps were computed with DALIGNER parameters -k24 -e.90 -s100 scaffolding algorithms had similar performance. We decided to use -l1000 -h600 intending to collapse haplotypes for the FALCON step to Salsa, as HiRise and Proximo HiC were not open access, and 3D-DNA better unzip genomes with high heterozygosity rate. FALCON-Unzip was computationally expensive on the DNAnexus platform. For short outputs both a pseudo-haplotype and a set of alternate haplotigs that read assemblies, other than Supernova and the NRGene assembly, the represent the secondary alleles. We refer to these outputs as the primary assembly GCA_000699085.116 was used for benchmarking, which was contig set (c1) and alternate contig set (c2). generated with Illumina paired-end, multiple mate-pair libraries and the SoapDeNovo81 assembler. The NRGene assembly was provided by Purging false duplications. Heterotype false duplications occurred the company with DeNovo Magic. despite setting FALCON76 parameters to resolve up to 10% haplotype di- vergence. FALCON-Unzip17 also incorrectly retained some secondary al- Gap filling. We ran PBJelly with support --capturedOnly --spanOnly leles in the primary contig set, which appeared as false duplications. To parameters, to avoid greedy gap closures with no spanning read reduce these false duplications, we ran Purge_Haplotigs13, first during curation (VGP v1.0 pipeline) and then later after contig formation (VGP QV for all genomes herein considered, we noticed that it was particu- v1.5 pipeline). To do the former, Purge_Haplotigs was run on the primary larly sensitive to the coverage cutoff parameter (-x). This is because contigs (c1), and identified haplotigs were mapped to the scaffolded Arrow generates a de novo consensus from the mapped reads without primary assembly with MashMap286 for removal. In the latter, identified explicitly considering the reference sequence. Later, we found that the haplotigs were moved from the primary contigs (c1) to the alternate second round of Arrow polishing sometimes reduced the QV accuracy haplotig set (p2). The remaining primary contigs were referred to as p1; for some species. Upon investigation, this issue was traced back to p2 combined with c2 was referred to as q2. Later, in the VGP v1.6 pipe- option -x 5, which requires at least 5 reads to call consensus. Such low line, we replaced Purge_Haplotigs with Purge_Dups14, a new program minimum requirements can lead to uneven polishing in low coverage developed by several of the authors in response to Purge_Haplotigs not regions. To avoid this behaviour, we suggest to increase the -x close to removing partial false duplication at contig boundaries. Purging also the half sequence coverage (for example, 30× when 60× was used for removes excessive low-coverage ( junk) and high-coverage (repeats) assembly) and check QV before moving forward. contigs. To calculate the presence and overall success of purging false For genomes with a combined assembly size larger than 4 Gb, we duplications, we used a k-mer approach (Supplementary Methods, used Minimap290 with parameters -ax map-pb instead of Blasr91 to Supplementary Fig. 6). overcome reference index size limitations. Two more rounds of base-pair polishing were performed with linked Scaffolding with 10XG linked reads. The 10X Genomics linked reads reads. The reads were aligned with Longranger align 2.2.2, which incor- were aligned to the primary contigs (p1), and an adjacency matrix was porates the Lauriat for barcode-aware alignment83. From the align- computed from the barcodes using scaff10x74 v2.0–2.1. Two rounds of ments, homozygous mismatches (variants) were called with FreeBayes83 scaffolding were performed. The first round was run with parameters v1.2.0 using default options. Consensus was called with bcftools con- -matrix 2000 -reads 12 -link 10, and the second round with parameters sensus92 with -i’QUAL>1 && (GT=’’AA’’ || GT = ‘’Aa’’)’ -Hla. -matrix 2000 -reads 8 -link 10. A gap of 100 bp (represented with ‘N’s) was inserted between joined contigs. The resulting primary scaffold VGP Trio Pipeline v1.0–v1.6. The trio pipeline is similarly designed to set was named s1. the standard pipeline, except for the use of parental data (Extended Data Fig. 3b). When parental genomes are available, the child’s CLR Scaffolding with Bionano optical maps. Bionano cmaps were gen- reads are binned to maternal and paternal haplotypes, and assembled erated using the Bionano Pipeline in non-haplotype assembly mode separately as haplotype-specific contigs (haplotigs) using TrioCanu20. and used to further scaffold the s1 assembly with Bionano Solve In brief, parental specific marker k-mers were collected using Meryl23 v3.2.187. We began with a one-enzyme nick map (BspQI), followed from the parental Illumina WGS reads of the parents. These markers by a two-enzyme nick map (BspQI and BssSI), and then with a DLE-1 were filtered and used to bin the child’s CLR read. A haplotype was as- one-enzyme non-nicking approach when the later data type became signed given the markers observed, normalized by the total markers in available (Supplementary Table 9). Scaffold gaps were sized accord- each haplotype. The subsequent purging, scaffolding, and polishing ing to the software estimate. The resulting scaffold set was named s2. steps were similarly updated with the use of Purge_Dups14 (v1.6). We extended binning to linked reads and Hi-C reads, by excluding read Scaffolding with Hi-C reads. Hi-C reads were aligned to the s2 scaf- pairs that had any parental-specific marker. The binned Hi-C reads were folds using the Arima Genomics mapping pipeline88. In brief, both used to scaffold its haplotype assembly, and polished with the binned ends of a read pair were mapped independently using BWA-MEM89 linked reads from the observation of haplotype switching using the with the parameter -B8, and filtered when mapping quality was <10. standard polishing approach. During curation, one of the haplotype Chimeric reads containing a restriction enzyme site were trimmed assemblies with the higher QV and/or contiguity was chosen as the from the restriction site onward, leaving only the 5′ end. The filtered representative haplotype. The heterogametic sex chromosome from single-read alignments were then rejoined as paired read alignments. the unchosen haplotype was added to the representative assembly. The processed alignments were then used for scaffolding with Salsa279, However, while curating several trios, we found that in regions of low which analyses the normalized frequency of Hi-C interactions between divergence between shared parental homogametic sex chromosomes all pairs of contig ends to determine a likely ordering and orienta- (that is, X or Z), a small fraction of offspring CLR data was mis-assigned tion of each. We used parameters -m yes -i 5 -p yes to allow Salsa2 to to the wrong haplotype. This mis-alignment resulted in a duplicate, break potentially mis-assembled contigs and perform five iterations low-coverage offspring X or Z assembly in the paternal (for mammals) of scaffolding. After feedback from curation, later versions of Salsa or maternal (for birds) haplotype, respectively, which required removal were developed, which more conservatively determine the number during curation. We are working on methods to improve the binning of iterations (v2.1) and actively break at mis-assemblies (v2.2), and accuracy for resolution of this issue going forward. run for the Canada lynx, Goode’s thornscrub tortoise, and two-lined For the female zebra finch in particular, contigs were generated caecilian. The restriction enzyme(s) used to generate each library were before the binning was automated in the Canu assembler as Trio- specified using parameters -e GATC,GANTC for Arima and -e GATC for Canu1.7, and therefore a manual binning process was applied as Dovetail and Phase Genomics Hi-C data. The resulting Hi-C scaffolded described in the original Trio-binning paper20 (Supplementary Meth- assembly was named s3. ods). Contigs were assembled for each haplotype using the binned reads, excluding unclassified reads. The contigs were polished with Consensus polishing. To polish bases in both haplotypes with two rounds of Arrow polishing using the binned reads, and scaffolded minimal alignment bias, we concatenated the alternate haplotig set following the v1.0 pipeline with no purging. Additional scaffolding (c2 in v1.0 or q2 in v1.5–1.6) to the scaffolded primary set (s3) and the as- rounds with Bionano (s4) and Hi-C were applied. Scaffolds were sembled mitochondrial genome (mitoVGP in v1.6). We then performed renamed according to the primary scaffold assembly of the same another round of polishing with Arrow (smrtanalysis 5.1.0.26412) individual (s5), with sex chromosomes grouped as Z in the paternal using PacBio CLR reads, aligning with pbalign --minAccuracy=0.75 assembly and W in the maternal assembly following synteny to the Z --minLength=50 --minAnchorSize=12 --maxDivergence=30 –con- chromosome from the curated male zebra finch VGP assembly. Two cordant --algorithm=blasr --algorithmOptions=--useQuality --max- rounds of SR polishing were applied using linked reads, by mapping Hits=1 --hitPolicy=random --seed=1 and consensus polishing with on both haplotypes. After haplotype switches were discovered, addi- variantCaller --skipUnrecognizedContigs haploid -x 5 -q 20 -X120 tional rounds of polishing were applied using binned linked reads –v --algorithm=arrow. While this round of polishing resulted in higher (Supplementary Methods). Article organelle assembly pipeline. Organelle matches embedded in nuclear Mitochondrial genome assembly. Similar to other recent methods93,94, sequences that were found to be NuMTs were kept. we developed a reference-guided MT assembly pipeline. MT reads in the raw CLR data were identified by mapping the whole read set to an False duplication removal. Retained false duplications were identi- existing reference sequence of the specific species or of closely re- fied using Purge_Haplotigs13 run either after scaffolding and polishing lated species using Blasr. Filtered mtDNA CLRs were assembled into a (Anna’s hummingbird, kākāpō, male zebra finch, female zebra finch, single contig using Canu v1.8, polished with Arrow using CLR and then platypus, pale spear-nosed bat, and greater horseshoe bat) or on the FreeBayes v1.0.2 together with bcftools v1.9 using short reads from the c1 before scaffolding (two-lined caecilian, flier cichlid, Canada lynx, 10XG data (Extended Data Fig. 3c). The overlapping sequences at the and Goode’s thornscrub tortoise). Subsequent manual curation identi- ends of the contig were trimmed, and the remaining contig sequence fied additional haplotypic duplications for the listed assemblies and circularized. The mitoVGP pipeline is made available at https://github. also those that were not treated with Purge_Haplotigs (Eastern happy, com/VGP/vgp-assembly/tree/master/mitoVGP. A more detailed pro- climbing perch, zig-zag eel). The evidence used included read cover- tocol description of the assembly pipeline and new discoveries from age, sequence self-comparison, transcript alignments, Bionano map the MT assemblies are published elsewhere33. alignments and Hi-C 2D maps, all confirming the superfluous nature of one allele. The identified haplotype duplications were moved from Curation the primary to the alternate assembly. The VGP genome assembly pipeline produces high quality assemblies, yet no automated method to date is free from the production of errors, Chromosome assignment. For a scaffold to be annotated as a chro- especially during the scaffolding stages. To minimize the impact of mosome, we used evidence from Hi-C as well as genetic linkage or FISH the remaining algorithmic shortcomings, we subjected all assemblies karyotype mapping when available. For Hi-C evidence, we considered to rigorous manual curation. All data generated for a species in this a scaffold as a complete chromosome (albeit with gaps) when there study and other publicly available data (for example, genetic maps, was a clear unbroken diagonal in the Juicebox or HiGlass plots for that gene sets and genome assemblies of the same or closely related spe- scaffold and no other large scaffolds that could be joined to that same cies) were collated, aligned to the primary assembly and analysed scaffold; if present and no unambiguous join was possible, we named it in gEVAL95 (https://vgp-geval.sanger.ac.uk/index.html), visualizing as an unlocalized scaffold for that chromosome. When we could not find discordances in a feature browser and issue lists. In parallel, Hi-C data evidence of a complete chromosome, we kept the scaffold number for were mapped to the primary assembly and visualized using Juice- its name. We named all evidence-validated scaffolds as chromosomes box96 and/or HiGlass97. With these data, genome curators identified down to the smallest Hi-C box unit resolution allowed with these char- mis-joins, missed joins and other anomalies, and corrected the primary acteristics. When there was an established chromosome terminology assembly accordingly. No change was made without unambiguous for a given species or set of species, we use the established terminology evidence from available data types; for example, a Hi-C suggested except when our new assemblies revealed errors in the older assembly, join would not be made unless supported by BioNano maps, long-read such as scaffold/chromosome fusions, fissions, rearrangements, and data, or gene alignments. When sequencing the heterogametic sex, non-chromosome names. For species without an established chromo- we identified sex chromosomes based on half coverage, homology some terminology, we named the scaffolds as chromosomes numbers alignments to sex chromosomes in other species, and the presence 1, 2, 3 …, in descending order of scaffold size. For the sex chromosomes, of sex chromosome-specific genes. we used the letters X and Y for mammals and Z and W for birds. Contamination removal. A succession of searches was used to identify Using comparative genomics to assess assembly structure. In cases potential contaminants in the generated assemblies. where a high-quality chromosome-level genome was available for a 1) A megaBLAST98 search against a database of common contami- closely related species, comparative genome analysis was performed. nants (ftp://ftp.ncbi.nlm.nih.gov/pub/kitts/contam_in_euks.fa.gz) The polished primary assembly (t3.p) was mapped to the related ge- requiring e ≤ 1 × 10−4, reporting matches with ≥98% sequence iden- nome using MashMap286 with --pi 75 -s 300000. The number of chro- tity and match length 50–99 bp, ≥94% and match length 100–199 bp, mosomal differences was identified using a custom script available at or ≥90% and match length 200 bp or above. https://github.com/jdamas13/assembly_comparison. This resulted 2) A vecscreen (https://www.ncbi.nlm.nih.gov/tools/vecscreen/) in the identification of ~60 to ~450 regions for each genome assem- search against a database of adaptor sequences (ftp://ftp.ncbi.nlm. bly flanking putative misassemblies or lineage-specific genome rear- nih.gov/pub/kitts/adaptors_for_screening_euks.fa) rangements. To identify which were real misassemblies, the identified 3) After soft-masking repeats using Windowmasker75, a megaBLAST discrepancies were communicated to the curation team for manual search against chromosome-level assemblies from RefSeq requiring verification (see above). e ≤ 1 × 10−4, match score ≥100, and sequence identity ≥98%; regions To identify any possible remaining mis-joins, each curated avian matching highly conserved rDNAs were ignored. and mammalian assembly was compared with the zebra finch (tae- Manual inspection of the results was necessary to differentiate Gut2) or human (hg38) genomes, respectively. Pairwise alignments contamination from conservation and/or horizontal gene transfer. between each of the VGP assemblies and the clade reference were Adaptor sequences were masked; other contaminant sequences were generated with LastZ99 (version 1.04) using the following parameters: removed. Assemblies were also checked for runs of Ns at the ends of C = 0 E = 30 H = 2000 K = 3000 L = 2200 O = 400. The pairwise scaffolds, created as artefacts of the iterative scaffolding process, and alignments were converted into the UCSC ‘chain’ and ‘net’ for- when found they were trimmed. mats with axtChain (parameters: -minScore = 1000 -verbose = 0 -linearGap = medium) followed by chainAntiRepeat, chainSort, chain- Organelle genomes. These were detected by a megaBLAST search PreNet, chainNet and netSyntenic, all with default parameters100. Pair- against a database of known organelle genomes requiring e ≤ 1 × 10−4, wise synteny blocks were defined using maf2synteny101 at 100-, 300-, sequence identity ≥90%, and match length ≥500; the databases are and 500-kb resolutions. Evolutionary breakpoint regions were detected available at ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/mito.nt.gz and and classified using an ad hoc statistical approach102. This analysis ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plastid/*genomic.fna.gz. identified 2 to 90 genomic regions per assembly that could be flank- Only scaffolds consisting entirely of organelle sequences were assumed ing misassemblies, lineage-specific chromosome rearrangements, or to be organelle genomes, and replaced by the genome from the separate reference-specific chromosome rearrangements (116 in the human and 26 in the zebra finch). Determining the underlying cause for each 88. Arima Genomics, Inc. Arima Genomics Mapping Pipeline. https://github.com/ ArimaGenomics/mapping_pipeline. of the flagged regions will need further verification. All alignments 89. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler are available for visualization at the Evolution Highway comparative transform. Bioinformatics 25, 1754–1760 (2009). chromosome browser (http://eh-demo.ncsa.illinois.edu/vgp/). 90. Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100 (2018). 91. Chaisson, M. J. & Tesler, G. Mapping single molecule sequencing reads using basic local Annotation alignment with successive refinement (BLASR): application and theory. BMC NCBI and Ensembl annotation pipeline used in this study are described Bioinformatics 13, 238 (2012). 92. Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, in the Supplementary Methods. 2078–2079 (2009). 93. Dierckxsens, N., Mardulyn, P. & Smits, G. NOVOPlasty: de novo assembly of organelle Evaluation genomes from whole genome data. Nucleic Acids Res. 45, e18 (2017). 94. Soorni, A., Haak, D., Zaitlin, D. & Bombarely, A. Organelle_PBA, a pipeline for assembling Detailed methods for other types of evaluation, including BUSCO chloroplast and mitochondrial genomes from PacBio DNA sequencing data. BMC runs, mis-join and missed-join identification, reliable blocks, collapsed Genomics 18, 49 (2017). repeats, telomeres, RNA-seq and ATAC–seq mapping, and false gene 95. Chow, W. et al. gEVAL — a web-based browser for evaluating genome assemblies. Bioinformatics 32, 2508–2510 (2016). duplications are in the Supplementary Methods. No statistical meth- 96. Durand, N. C. et al. Juicebox provides a visualization system for Hi-C contact maps with ods were used to predetermine sample size, the experiments were not unlimited zoom. Cell Syst. 3, 99–101 (2016). randomized, and the investigators were not blinded to group during 97. Kerpedjiev, P. et al. HiGlass: web-based visual exploration and analysis of genome interaction maps. Genome Biol. 19, 125 (2018). experiments and outcome assessment. 98. Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinformatics 10, 421 (2009). Reporting summary 99. Harris, R. S. Improved Pairwise Alignment of Genomic DNA. Thesis, Pennsylvania State Univ. (2007). Further information on research design is available in the Nature 100. Kent, W. J., Baertsch, R., Hinrichs, A., Miller, W. & Haussler, D. Evolution’s cauldron: Research Reporting Summary linked to this paper. duplication, deletion, and rearrangement in the mouse and human genomes. Proc. Natl Acad. Sci. USA 100, 11484–11489 (2003). 101. Kolmogorov, M., Raney, B., Paten, B. & Pham, S. Ragout—a reference-assisted assembly tool for bacterial genomes. Bioinformatics 30, i302–i309 (2014). Data availability 102. Farré, M. et al. Novel insights into chromosome evolution in birds, archosaurs, and reptiles. Genome Biol. Evol. 8, 2442–2451 (2016). All raw data, intermediate and final assemblies are publicly available via 103. Guan, D. Asset. https://github.com/dfguan/asset. GenomeArk (https://vgp.github.io/genomeark), archived on NCBI/EBI 104. Tarailo-Graovac, M. & Chen, N. Using RepeatMasker to identify repetitive elements in BioProject under accession PRJNA489243 with annotations, and brows- genomic sequences. Curr. Protoc. Bioinformatics 25, 4.10.1–4.10.14 (2009). 105. Krumsiek, J., Arnold, R. & Rattei, T. Gepard: a rapid and sensitive tool for creating dotplots able on the UCSC Genome Browser (https://hgdownload.soe.ucsc.edu/ on genome scale. Bioinformatics 23, 1026–1028 (2007). hubs/VGP/). The final primary assembly from the automated pipeline 106. Harry, E. PretextView. https://github.com/wtsi-hpag/PretextView. before curation is browsable on gEVAL (https://vgp-geval.sanger.ac.uk) 107. Kurtz, S. et al. Versatile and open software for comparing large genomes. Genome Biol. 5, R12 (2004). with all four raw data mappings. The VGP assembly pipeline is available 108. Nattestad, M. Dot. https://github.com/MariaNattestad/dot. as a stand-alone pipeline (https://github.com/VGP/vgp-assembly) as well as a workflow on DNAnexus (https://platform.dnanexus.com/). Acknowledgements We thank the following persons for feedback and support: R. Johnson, A VGP-specific assembly hub portal in the U.C. Santa Cruz browser is E. Karlsson, K. Lindblad Toh, W. Jun, I. Korf, W. Haerty, G. Etherington, B. Clavijo, and A. available as a gateway to access all VGP genome assemblies and annota- Komissarov for discussions in the early stages of the project; R. Fuller for help with the G10K tions (https://hgdownload.soe.ucsc.edu/hubs/VGP). website maintenance, and H. Segal for help with with VGP website development; M. Linh Pham for help with initial grant writing; L. Shalmiyev for administrative help; D. Church, G. Kol, K. Baruch, O. Barad, I. Liachko, E. Muzychenko, S. Garg, and M. Kolmogorov for preliminary analyses performed on one or more genomes; K. Oliver, C. Corton and J. Skelton for data Code availability generation; E. Harry for technical support in scaff10x and Pretext; C. Mazzoni for coordinating students and training at Leibniz Institute for Zoo and Wildlife Research and Berlin Center for All codes used in the VGP Assembly Pipeline and the VGP Trio Pipeline Genomics in Biodiversity Research; and M. Driller, C. Caswara, M. Vafadar, N. Hill, D. De Panis, are publicly available at https://github.com/VGP/vgp-assembly/tree/ A. Whibley, B. Maloney, C. Mitchell, G. Gallo, J. Gaige, K. Amoako-Boadu, M. Jose Gomez, M. master/pipeline. Montero, D. Ratnikov, S. Brown, S. Zylka, S. Marcus, and T. Carrasco for completing training and testing the VGP pipeline by producing ordinal representative genome assemblies not described in this manuscript. We thank our company partners (listed below), NCBI, EBI, and 73. Ondov, B. D. et al. Mash: fast genome and metagenome distance estimation using Amazon AWS, including AWS for sponsoring sequence storage. J. Fekecs and D. Leja created MinHash. Genome Biol. 17, 132 (2016). the animal images, and J. Kim modified them to silhouettes. We thank them for their 74. Ning, Z. & Harry, E. Scaff10X https://github.com/wtsi-hpag/Scaff10X. permission to publish. A.R., S.K., B.P.W. and A.M.P. were supported by the Intramural Research 75. Morgulis, A., Gertz, E. M., Schäffer, A. A. & Agarwala, R. WindowMasker: window-based Program of the NHGRI, NIH (1ZIAHG200398). A.R. was also supported by the Korea Health masker for sequenced genomes. Bioinformatics 22, 134–141 (2006). Technology R&D Project through KHIDI, funded by the Ministry of Health & Welfare, Republic 76. Chin, C.-S. et al. Nonhybrid, finished microbial genome assemblies from long-read SMRT of Korea (HI17C2098). S.A.M., I.B. and R.D. were supported by Wellcome Trust grant WT207492; sequencing data. Nat. Methods 10, 563–569 (2013). W.C., M. Smith, Z.N., Y.S., J.C., S. Pelan, J.T., A.T., J.W. and Kerstin Howe by WT206194; L.H., F.M., 77. Koren, S. et al. Canu: scalable and accurate long-read assembly via adaptive k-mer Kevin Howe and P. Flicek by WT108749/Z/15/Z, WT218328/B/19/Z and the European Molecular weighting and repeat separation. Genome Res. 27, 722–736 (2017). Biology Laboratory. O.F. and E.D.J. were supported by Howard Hughes Medical Institute and 78. Weisenfeld, N. I., Kumar, V., Shah, P., Church, D. M. & Jaffe, D. B. Direct determination of Rockefeller University start-up funds for this project. J.D. and H.A.L. were supported by the diploid genome sequences. Genome Res. 27, 757–767 (2017). Robert and Rosabel Osborne Endowment. M.U.-S. received funding from the European Union’s 79. Ghurye, J. et al. Integrating Hi-C links with assembly graphs for chromosome-scale Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant assembly. PLoS Comput. Biol. 15, e1007273 (2019). agreement (750747). F.T.-N., J. Hoffman, P. Masterson and K.C. were supported by the 80. Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions reveals Intramural Research Program of the NLM, NIH. C.L., B.J.K., J. Kim and H.K. were supported by folding principles of the human genome. Science 326, 289–293 (2009). the Marine Biotechnology Program of KIMST, funded by the Ministry of Ocean and Fisheries, 81. Luo, R. et al. SOAPdenovo2: an empirically improved memory-efficient short-read Republic of Korea (20180430). M.C. was supported by Sloan Research Fellowship (FG-2020- de novo assembler. Gigascience 1, 18 (2012). 12932). S.C.V. was funded by a Max Planck Research Group award from the Max Planck Society, 82. English, A. C. et al. Mind the gap: upgrading genomes with Pacific Biosciences RS and a Human Frontiers Science Program (HFSP) Research grant (RGP0058/2016). T.M.L., W.E.J. long-read sequencing technology. PLoS ONE 7, e47768 (2012). and the Canada lynx genome were funded by the Maine Department of Inland Fisheries & 83. Bishara, A. et al. Read clouds uncover variation in complex regions of the human Wildlife (F11AF01099), including when W.E.J. held a National Research Council Research genome. Genome Res. 25, 1570–1580 (2015). Associateship Award at the Walter Reed Army Institute of Research (WRAIR). C.B. was 84. Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial variant detection supported by the NSF (1457541 and 1456612). D.B. was funded by The University of Queensland and genome assembly improvement. PLoS ONE 9, e112963 (2014). (HFSP - RGP0030/2015). D.I. was supported by Science Exchange Inc. (Palo Alto, CA). H.W.D. 85. Garrison, E. & Marth, G. Haplotype-based variant detection from short-read sequencing. was supported by NSF grants (OPP-0132032 ICEFISH 2004 Cruise, PLR-1444167 and OPP- Preprint at http://arxiv.org/abs/1207.3907 (2012). 1955368) and the Marine Science Center at Northeastern University (416). G.J.P.N. and the 86. Jain, C., Koren, S., Dilthey, A., Phillippy, A. M. & Aluru, S. A fast adaptive algorithm for thorny skate genome were funded by Lenfest Ocean Program (30884). M.P. was funded by the computing whole-genome homology maps. Bioinformatics 34, i748–i756 (2018). German Federal Ministry of Education and Research (01IS18026C). M. Malinsky was supported 87. Bionano Genomics, Inc. Bionano Software Downloads. https://bionanogenomics.com/ by an EMBO fellowship (ALTF 456-2016). The following authors’ contributions were supported support/software-downloads/. by the NIH: S. Selvaraj (R44HG008118); C.V.M., S.R.F., P.V.L. (R21 DC014432/DC/NIDCD); K.D.M. Article (R01GM130691); H.C. (5U41HG002371-19); M.D. (U41HG007234); and B.P. (R01HG010485). D.G. J.L., B.P.W., M. Malinsky. Generated/modified software: S.K., A.R., S.B.K., R.H., Z.K., J. Korlach, was supported by the National Key Research and Development Program of China I.S., C.D., Z.N., A.H., J.L., J.G., E.G., C.V.M., S.R.F., N.H.P. Pipeline development: A.R., S.A.M., G.F., (2017YFC1201201, 2018YFC0910504 and 2017YFC0907503). F.O.A. was supported by S.K., M.U.-S., A.F., M. Simbirsky, B.T.H., T.P., M.P., E.W.M., R.D., A.M.P. Generated MT assemblies: Al-Gannas Qatari Society and The Cultural Village Foundation-Katara, Doha, State of Qatar and G.F., J. Korlach. Curation: Kerstin Howe, W.C., Y.S., J.C., S. Pelan, J.T., A.T., J.W., Y.Z., J.D., H.A.L. Monash University Malaysia. C.T. was supported by The Rockefeller University. M. Hiller was Sex chromosomes: Y.Z., R.S.H., K.D.M., P. Medvedev, J. M. Graves. Hummingbird karyotype supported by the LOEWE-Centre for Translational Biodiversity Genomics (TBG) funded by the analyses: M. Houck, A. Misuraca, M.P., E.W.M., E.D.J. Annotation: F.T.-N., L.H., J. Hoffman, P. Hessen State Ministry of Higher Education, Research and the Arts (HMWK). H.C. was supported Masterson, K.C., F.M., Kevin Howe, P. Flicek, D.B. Evaluation analysis: A.R., J.D., M.U.-S., J. Kim, by the NHGRI (5U41HG002371-19). R.H.S.K. was funded by the Max Planck Society with C.L., B.J.K., M.C., G.L.G., L.J.C., F.T.-N., L.H., J. M. George, J.G., R.E.D., D.G., S.E.L., D.F.C., C.V.M., computational resources at the bwUniCluster and BinAC funded by the Ministry of Science, S.R.F., P.V.L., E.O., F.O.A.-A., S. Secomandi, C.T., M. Hiller, H.K., Kerstin Howe, E.W.M., R.D., A.M.P., Research and the Arts Baden-Württemberg and the Universities of the State of E.D.J. Biological findings: J.D., J. Kim, C.L., B.J.K., G.L.G., L.J.C., H.A.L., A.R., E.D.J. Data Baden-Württemberg, Germany (bwHPC-C5). B.V. was supported by the Biomedical Research availability: A.R., S.A.M., W.C., A.F., S. Paez, M. Simbirsky, B.T.H., B.P.W., W.K., H.C., M.D., L.N., Council of A*STAR, Singapore. T.M.-B. was funded by the European Research Council under B.P., A.M.P., E.D.J. G10K council, founders, and coordination of VGP: T.M.-B., A.J.C., F.D.P., R.D., the European Union’s Horizon 2020 research and innovation programme (864203), MINECO/ M.T.P.G., E.D.J., K.-P.K., H.A.L., R.W.M., E.W.M., E.C.T., B.V., G.Z., A.M.P., S. Paez, J. M. Graves, FEDER, UE (BFU2017-86471-P), Unidad de Excelencia María de Maeztu, AEI (CEX2018- O.A.R., D.H., S.J.O., T.W. and B.S. All authors reviewed the manuscript. 000792-M), a Howard Hughes International Early Career award, Obra Social “La Caixa” and Secretaria d’Universitats i Recerca and CERCA Programme del Departament d’Economia i Competing interests During the contributing period, B.T.H., M. Simbirsky, A.F. and M. Mooney Coneixement de la Generalitat de Catalunya (GRC 2017 SGR 880). E.C.T. was supported by the were employees of DNAnexus Inc. S.B.K., R.H., Z.K., J. Korlach, I.S. and C.D. were full-time European Research Council (ERC-2012-StG311000) and an Irish Research Council Laureate employees at Pacific Biosciences, a company developing single-molecule long read Award. M.T.P.G. was supported by an ERC Consolidator Award 681396-Extinction Genomics, sequencing technologies. R.E.G., N.H.P., and J.G. were affiliated with Dovetail Genomics, a and a Danish National Research Foundation Center Grant (DNRF143). T.W. was supported by company developing genome assembly tools, including Hi-C. I.G. was affiliated with Oxford the NSF (1458652). J. M. Graves was supported by the Australian Research Council Nanopore Technologies, a company generating long read sequencing technologies. A.H. and (CEO561477). E.W.M. was partially supported by the German Federal Ministry of Education and J.L were employees of Bionano Genomics, a company developing optical maps for genome Research (01IS18026C). Complementary sequencing support for the Anna’s hummingbird and assembly. S. Selvaraj was an employee of Arima Genomics, a company developing Hi-C data several genomes was provided by Pacific Biosciences, Bionano Genomics, Dovetail Genomics, for genome assemblies. R.D. is a scientific advisory board member of Dovetail Inc. P. Flicek is a Arima Genomics, Phase Genomics, 10X Genomics, NRGene, Oxford Nanopore Technologies, member of the Scientific Advisory Boards of Fabric Genomics, Inc., and Eagle Genomics, Ltd. Illumina, and DNAnexus. All other sequencing and assembly were conducted at the H.C. receives royalties from the sale of UCSC Genome Browser source code, LiftOver, GBiB, Rockefeller University, Sanger Institute, and Max Planck Institute Dresden genome labs. Part of and GBiC licenses to commercial entities. S.K. has received travel funds to speak at symposia this work used the computational resources of the NIH HPC Biowulf cluster (https://hpc.nih. organized by Oxford Nanopore. M.D. and L.N. receive royalties from licensing of UCSC gov). We acknowledge funding from the Wellcome Trust (108749/Z/15/Z) and the European Genome Browser. For W.E.J., the content here is not to be construed as the views of the DA or Molecular Biology Laboratory. We thank Le Comité Scientifique Régional du Patrimoine DOD. All other authors declare no competing interests. Naturel and Direction de l’Environnement, de l’Aménagement et du Logement, Guyanne for research approvals and export permits. Additional information Supplementary information The online version contains supplementary material available at Author contributions Wrote the paper and co-coordinated the study: A.R., E.D.J., A.M.P., R.D., https://doi.org/10.1038/s41586-021-03451-0. E.W.M., Kerstin Howe, S.A.M., O.F. Coordination with vendors: J. Korlach, S. Selvaraj, R.E.G., Correspondence and requests for materials should be addressed to Kerstin Howe., E.W.M., A.H., M. Mooney. Collected samples: M.T.P.G., W.E.J., R.W.M., G.Z., B.V., M.T.B., J. Howard, S.C.V., R.D., A.M.P. or E.D.J. T.M.L., F.G., W.C.W., D.B., J. M. George, M.T.B., D.I., A.D., D.E., B.R., T.E., M. Wilkinson, G.T., A. Peer review information Nature thanks Michael Schatz, Justin Zook and the other, anonymous, Meyer, A.F.K., P. Franchini, H.W.D., H.S., M. Wagner, G.J.P.N., R.D., E.D.J., E.C.T., R.H.S.K. reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are Generated genome data: O.F., I.B., M. Smith, B.H., J.M., S.W., C.B., A. Meyer, A.F.K., P. Franchini, available. I.G., D.F.C., C.V.M. Generated genome assemblies: A.R., S.A.M., S.K., M.P., S.B.K., R.H., J.G., Z.N., Reprints and permissions information is available at http://www.nature.com/reprints. Extended Data Fig. 1 | Assessment of completeness of the Anna’s shows the number of counts of low coverage or signs of collapsed repeats. hummingbird assembly. a, b, Steps and NG50 continuity values of the VGP Larger chromosomal scaffolds (1–19) have fewer gaps and low coverage or assembly pipeline that gave the highest quality assembly for Anna’s collapsed regions compared with the micro chromosomes (20–33). hummingbird (a) and Canada lynx (b) in this study. The specific steps are Chromosomes 14, 15 and 19 of the Anna’s hummingbird were the most outlined further in Extended Data Fig. 2a, and Methods. c, Whole-genome structurally reliable scaffolds, having only one gap each with no low-support alignment of CLR (red), linked reads (green), optical maps (blue), and Hi-C regions. We defined reliable blocks as those supported by at least two reads (purple) of the Anna’s hummingbird, along with telomere motif (TTAGGG technologies. Reliable blocks excluded regions with structural assembly and its reverse complement, yellow) and gaps (grey) using Asset software103. errors, such as collapsed repeats or unresolved segmental duplications. For each data type, the first row shows the mapped coverage, and the second Low-support regions are those where the reliable blocks row has a peak. Article Extended Data Fig. 2 | VGP assembly pipeline applied across multiple the assemblies of 16 vertebrate species, weighted by the fraction of bases in species. a, Iterative assembly pipeline of sequence data types (coloured as in each length bin (log scaled). Molecule length above 1 kb was measured from b) with increasing chromosomal distance. Thin bars, sequence reads; thick read length for CLR, estimated molecule coverage for linked reads, raw black bars, assembled contigs; black bars with space and arcing links, scaffolds; molecule length for optical maps, and interaction distance for Hi-C reads. For grey bars, gaps placed by previous steps; thick red border, tracking of an each species, the fragment length distribution of each data type was similar to example contig in the pipeline. The curation step shows an example of a those for the Anna’s hummingbird, with differences primarily influenced by mis-assembly break identified by sequence coverage (grey, left) and an tissue type, preservation method, and collection or storage conditions example of an inversion error (right) detected by the optical map. (unpublished data). b, Intra-molecule length distribution of the four data types used to generate Extended Data Fig. 3 | Flow charts of assembly pipelines used to generate assembly errors using gEVAL 25,95 b, Standard VGP trio assembly pipeline when high-quality assemblies in this study. a, Standard VGP assembly pipeline DNA is available for a child and parents20. Dashed line indicates that the other when sequencing data of one individual, that generated the highest quality haplotype went through the same steps before curation. In addition to the assemblies: generate primary pseudo-haplotype and alternate haplotype curated assemblies of both haplotypes, a representative haplotype with both contigs with CLR using FALCON-Unzip17; generate scaffolds with linked reads sex chromosomes is submitted. c, Mitochondrial assembly pipeline. Figure key using Scaff10x74; break mis-joins and further scaffold with optical maps using applies to a–c. Steps newly introduced in v1.5–v1.6 are highlighted in light blue. Solve87; generate chromosome-scale scaffolds with Hi-C reads using Salsa279; c, contigs; p, purged false duplications from primary contigs; q, purged fill in gaps and polish base-errors with CLR using Arrow (Pacific BioSciences); alternate contigs; s, scaffolds; t, polished scaffolds. Further details and perform two or more rounds of short-read polishing with linked reads using instructions are available elsewhere33 and at https://github.com/VGP/ FreeBayes85; and perform expert manual curation to correct potential vgp-assembly. Article Extended Data Fig. 4 | Relationship between collapses and genomic coloured by species with its size relative to the length of the collapse as it characteristics. a, Correlation between the total number of collapses and appears in the assembly. Collapses above the horizontal bar (>90%) are further percentage repeat content estimated in the submitted curated versions of classified as collapsed high-copy repeats, and those below the horizontal bar n = 17 genomes from 16 species. b, Correlation between total number of bases are classified as segmental duplications (low-copy repeats). i, Major repeat in collapsed regions per Gb and repeat content. c, Correlation between total types in collapsed high-copy repeats. Most of the repeats were masked only missing bases collapsed per Gb and repeat content. d, Correlation between with WindowMasker75, with no annotation available by RepeatMasker104. total number of genes (coding and non-coding) in the collapsed regions and j, Minor repeat types in collapsed repeats. This is a breakdown of the repeats repeat content. e, Lack of correlation between the average collapsed size and categorized as ‘Others’ in i, owing to the smaller scale. Bar colours in i and j are repeat content. f, Lack of correlation between the total number of collapses as in h. Note smaller scale in j compared with i. Collapsed satellite arrays were and percentage heterozygosity. g, Lack of correlation between the total almost exclusively found in the platypus, comprising ~2.5 Mb. Collapsed simple number of collapses per Gb and genome size. Genome size, heterozygosity, repeats were the major source in the thorny skate (~400 kb). There was a higher and repeat content were estimated from 31-mer counts using GenomeScope71. proportion of LTRs in birds, LINEs and SINEs in mammals, and DNA repeats in Reported are adjusted r2 and P values from F-statistics. h, Cumulative collapsed the amphibian. Among the genes in the collapses, many were repetitive short bases per Gb in each collapse and percentage repeat masked. Each circle is non-coding RNAs. P values from F-statistics. Extended Data Fig. 5 | False duplication mechanisms in genome assembly. homotype duplications occur where a sequence from the same genomic locus a, False heterotype (haplotype) duplications occurs when more divergent is duplicated, and are of two types: 1. Overlapping sequences at contig sequence reads from each haplotype A (blue) and B (red) (maternal and boundaries: in current overlap-layout-consensus assemblers, branching paternal) form greater divergent paths in the assembly graph (bubbles), while sequences in assembly graphs that are not selected as the primary path have a nearly identical homozygous sequences (black) become collapsed. When the small overlapping sequence (purple), dovetailing to the primary path where it assembly graph is properly formed and correctly resolved (green arrow), one originated a branch. The size of the duplicated sequence is often the length of a of the haplotype-specific paths (red or blue) is chosen for building a ‘primary’ corrected read. Subsequent scaffolding results in tandem duplicated pseudo-haplotype assembly and the other is set apart as an ‘alternate’ sequences with a gap between. 2. Under-collapsed sequences: sequencing assembly. When the graph is not correctly resolved (purple arrow), one of four errors in reads (red x) randomly or systematically pile up, forming types of pattern are formed in the contigs and subsequent scaffolds. under-collapsed sequences. Subsequent duplication errors in the scaffolding Depending on the supporting evidence, the scaffolder either keeps these are similar to the heterotype duplications. Purge_haplotigs13 align sequences haplotype contigs on separate scaffolds or brings them together on the same to themselves to find a smaller sequence that aligns fully to a larger contig or scaffold, often separated by gaps: 1. Separate contigs: both contigs are scaffold, and removes heterotype duplication types 1, 3, and 4. Purge_dups14 retained in the primary contig set, an error often observed when additionally uses coverage information to detect heterotype duplication type haplotype-specific sequences are highly diverged. 2. Flanking contigs: the 2 and homotype duplications. We distinguished the two types of duplications assembly graph is partially formed, connecting the homozygous sequence of by: 1) haplotype-specific variants in reads aligning at half coverage to each the 5′ side to one haplotype (blue) and the 3′ side to the other haplotype (red). 3. heterotype duplication; 2) differing consensus quality that resulted from read Partial flanking contigs: only one haplotype (blue) flanks one side of the coverage fluctuations when aligning reads to homotype duplications; and 3) homozygous sequence. 4. Failed connecting of contigs: all haplotype k-mer copy number anomalies in which homotype duplications were observed sequences fail to properly connect to flanking homozygous sequences. b, False in the assembly with more than the expected number of copies. Article Extended Data Fig. 6 | False duplication examples fixed during manual gEVAL95. b, Example of a homotype duplication found in the hummingbird curation. a, An example of a heterotype duplication in the female zebra finch, assembly. These were caused by an algorithm bug in FALCON, which was later non-trio assembly. Left, a self-dot plot of this region generated with Gepard105, fixed. c, Example of a combined duplication involving both heterotype (green) with sequences coloured by haplotypes. Gaps, duplicated sequences (green and homotype (orange) duplications. Assembly graph structure is shown on and purple), and haplotype-specific marker densities are indicated at the top. the left for clarity, highlighting the overlapping sites at the contig boundary Right, a detailed alignment view of the green haplotype duplication with shaded following the duplication type. Assembly errors including the above paternal and maternal markers, self-alignment components, transcripts false duplications were detected and fixed during the curation process. annotated, contigs, bionano maps, and repeat components displayed in Extended Data Fig. 7 | Evidence of near-complete chromosome scaffolds in matches to smaller scaffolds are not linked because of ambiguous order or the VGP assemblies. Shown are Hi-C interaction heat maps for each species orientation, and are instead submitted as ‘unlocalized’ belonging to the after curation, visualized with PretextView106. A scaffold is considered a relevant chromosome. Bands at the top of each heat map show scaffolds putative arm-to-arm chromosome when all Hi-C read pairs in a row and column identified as X, Z (blue) or Y, W (red) sex chromosomes. The Hi-C map of map to a square (that is, an assembled chromosome) on the diagonal without fAstCal1 is not included as we had no remaining tissue left of the animal used to any other interactions off the diagonal. Those with remaining off-diagonal generate Hi-C reads. Article Extended Data Fig. 8 | Comparison of chromosomal organization between reference. The lines within each block depict orientation relative to the previous and new VGP assemblies. a, Zebra finch male compared to a reference; a positive slope is the same orientation as the reference, whereas a previous reference assembly of the same animal. b, Platypus male compared negative slope is the inverse orientation. Gaps are white boxes with no lines, in with a previous reference female assembly (so the Y chromosomes are absent the reference relative to the VGP assembly. A white box for the entire in the previous reference). c, Hummingbird female compared to a previous chromosome means a newly identified chromosome in the VGP assembly. Top reference of the same animal. d, Climbing perch compared to a previous 20 is the longest 20 scaffolds of the hummingbird and climbing perch reference. Each row represents a VGP-generated chromosome for the target assemblies. Accession numbers of the assemblies compared are listed in species. Colours depict identity with the reference (see key to the right); more Supplementary Table 19. than one colour indicates reorganization in the VGP assembly relative to the Extended Data Fig. 9 | Haplotype-resolved sex chromosomes and mammals (platypus, Canada lynx, pale spear-nosed bat, and greater horseshoe mitochondrial genomes. a, Alignment scatterplot, generated with MUMmer bat) showing a higher density of repeats within the mammalian X chromosome NUCmer107, visualized with dot108, of maternal and paternal chromosomes from than the avian Z chromosome. d, VGP kākāpō mitochondrial genome assembly the female zebra finch trio-based assembly. Blue, same orientation; red, reveals previously missing repetitive sequences (adding 2,232 bp) in the origin inversion; orange, repeats between haplotypes. The paternal Z chromosome is of replication region, containing an 83-bp repeat unit. e, VGP climbing perch highly divergent from the maternal W, and thus mostly unaligned. b, Alignment mitochondrial genome assembly showing a duplication of trnL2 and partial scatterplot of assembled Z and W chromosomes across the three bird species, duplication of Nad1, which were absent from the prior reference. Orange approximated with MashMap286. Segments of 300 kb (green), 500 kb (blue), arrows and red lines, tRNA genes and their alignments; dark grey arrows and and 1 Mb (purple) are shaded darker with higher sequence identity, with a grey shading, all other genes and their alignments; black, non-coding regions; minimum of 85%. The smaller size and higher repeat content of the W green line, conventional starting point of the circular sequence. chromosome are clearly visible. c, X and Y chromosome segments of the Article Extended Data Fig. 10 | Large haplotype inversions with direct evidence in inversion found on chromosome 11 and a complicated 8.1-Mb rearrangement the zebra finch trio assembly. a, Two inversions (green and red) in on chromosome 13 between maternal and paternal haplotypes. d, No chromosome 5 found from the MUMmer NUCmer107 alignment of the maternal mis-assembly signals were detected from the binned Hi-C interaction plots, and paternal haplotype assemblies, visualized with dot108. b, Hi-C interaction indicating that the haplotype-specific inversions are real. e, Half the PacBio plot showing that the trio-binned Hi-C data remove most of the interactions CLR span and Bionano optical maps agree with the inversion breakpoints in from the other haplotype (red arrows), which could be erroneously classified as chromosome 11, supporting the haplotype-specific inversion. a mis-assembly if only one haplotype was used as a reference. c, An 8.5-Mb Extended Data Fig. 11 | Polishing artefacts. a, An example of uneven mapping reads to both haplotype assemblies, and polishing by mapping paternally coverage in the primary and alternate sequence pair of the Anna’s binned reads to the paternal assembly. Polishing improves QV, but introduces hummingbird assembly. In this example, the alternate (alt) sequence was built haplotype switch errors when using reads from both haplotypes as shown in at higher quality, attracting all linked-reads for polishing. The matching locus row 2. This can be avoided when using haplotype binned reads for polishing. in the primary (pri) assembly was left unpolished, resulting in frameshift errors c, Example of over-polishing. The nuclear mitochondria (NuMT) sequence was in the TLK1 gene. b, Haplotype-specific markers (red for maternal, blue for transformed as a full mitochondria (MT) sequence during long-read polishing paternal) and error markers found in the assembly on the Z chromosome owing to the absence of the MT contig, where the NuMT attracted all long reads (inherited from the paternal side) of the trio-binned female zebra finch from the MT. In comparison, the trio-binned assembly had the MT sequence assembly. Each row shows markers before short-read polishing, mapping all assembled in place, preventing mis-placing of MT reads during read mapping. Article Extended Data Fig. 12 | Chromosome evolution among the bat species fission of the boreoeutherian chromosome 5 in the bat ancestor; blue star, the sequenced. a, Genes surrounding an inversion in the greater horseshoe bat, region upstream of the breakpoint in the bats; green star, the region relative to human chromosome 15 (red highlight). The STARD5 gene is directly downstream of the breakpoint in the bats. The red starred breakpoint was disrupted by this inversion, which separates exons 1–5 from exon 6 in the confirmed as reused, as opposed to assembly errors, in chromosomal greater horseshoe bat. b, RNA-seq tracks showing the lack of RNA splicing rearrangements of the pale spear-nosed bat, Egyptian fruit bat, and greater evidence of STARD5 transcripts in the greater horseshoe bat (bottom) in horseshoe bat. There is no evidence of reuse for the velvety free-tailed bat. We comparison to the pale spear-nosed bat where the STARD5 gene is not disrupted could not confirm breakpoint reuse in the greater mouse-eared bat or Kuhl’s (top). c, Circos plots of chromosome organization relationships between the pipistrelle at the chromosomal scale because they were on small scaffolds that each of the analysed bats and segments of the human chromosomes 1, 2, 6 and may not be completely assembled. 10. Red star, breakpoint location in human chromosome 6, depicting the Extended Data Table 1 | Summary metrics of the curated and submitted vertebrate species assemblies Colour shading indicates degree of heterozygosity or repeats (red), primary assembly sizes and relative size of alternate haplotypes (orange), continuity measures (green), gaps and collapses (blue), and base call accuracy (purple). A dash indicates that the sex chromosomes were not found or are not known. Accessions are available in Supplementary Table 10. Article Extended Data Table 2 | Annotation summary statistics in previous and newly assembled VGP reference genomes Annotation results of VGP assemblies and previous reference assemblies with the NCBI Eukaryotic Genome Annotation Pipeline, using the same RNA-seq data and nearly identical sets of transcripts and proteins on input. Highlighted rows are plotted in Fig. 5c, d. nature research | reporting summary Erich D. Jarvis, Arang Rhie, Kerstin Howe, Corresponding author(s): Eugene Myers, Richard Durbin, Adam Phillippy Last updated by author(s): Jan 8, 2021 Reporting Summary Nature Research wishes to improve the reproducibility of the work that we publish. This form provides structure for consistency and transparency in reporting. For further information on Nature Research policies, see our Editorial Policies and the Editorial Policy Checklist. Statistics For all statistical analyses, confirm that the following items are present in the figure legend, table legend, main text, or Methods section. n/a Confirmed The exact sample size (n) for each experimental group/condition, given as a discrete number and unit of measurement A statement on whether measurements were taken from distinct samples or whether the same sample was measured repeatedly The statistical test(s) used AND whether they are one- or two-sided Only common tests should be described solely by name; describe more complex techniques in the Methods section. A description of all covariates tested A description of any assumptions or corrections, such as tests of normality and adjustment for multiple comparisons A full description of the statistical parameters including central tendency (e.g. means) or other basic estimates (e.g. regression coefficient) AND variation (e.g. standard deviation) or associated estimates of uncertainty (e.g. confidence intervals) For null hypothesis testing, the test statistic (e.g. F, t, r) with confidence intervals, effect sizes, degrees of freedom and P value noted Give P values as exact values whenever suitable. For Bayesian analysis, information on the choice of priors and Markov chain Monte Carlo settings For hierarchical and complex designs, identification of the appropriate level for tests and full reporting of outcomes Estimates of effect sizes (e.g. Cohen's d, Pearson's r), indicating how they were calculated Our web collection on statistics for biologists contains articles on many of the points above. Software and code Policy information about availability of computer code Data collection We collected genome sequence data of various technologies from individuals representing 16 vertebrate species. Data analysis All software used is at at https://github.com/VGP/vgp-assembly, as of November 2020. The assembly software and versions used are listed in Supplementary Table The specific VGP assembly pipeline used is at https://github.com/VGP/vgp-assembly/tree/master/pipeline. For manuscripts utilizing custom algorithms or software that are central to the research but not yet described in published literature, software must be made available to editors and reviewers. We strongly encourage code deposition in a community repository (e.g. GitHub). See the Nature Research guidelines for submitting code & software for further information. Data Policy information about availability of data All manuscripts must include a data availability statement. This statement should provide the following information, where applicable: - Accession codes, unique identifiers, or web links for publicly available datasets - A list of figures that have associated raw data April 2020 - A description of any restrictions on data availability All raw sequence and assembly data are available at Genome Ark https://vgp.github.io/genomeark/ and the NCBI BioProject page PRJNA489243 https:// www.ncbi.nlm.nih.gov/bioproject/489243 1 nature research | reporting summary Field-specific reporting Please select the one below that is the best fit for your research. If you are not sure, read the appropriate sections before making your selection. Life sciences Behavioural & social sciences Ecological, evolutionary & environmental sciences For a reference copy of the document with all sections, see nature.com/documents/nr-reporting-summary-flat.pdf Life sciences study design All studies must disclose on these points even when the disclosure is negative. Sample size For most analyses, we analyzed data from a sample size of n = 16 vertebrate species, and n = 17 individuals; we had two individuals for one species (the zebra finch). For the chromosomal evolution analyses, we added an additional n = 4 bat species. Data exclusions Only sequence data that failed quality control were excluded or repeated. Replication We confirmed the ability to replicate all code using multiple rounds of assembly or analyses. Randomization No analyses required generating randomize data sets. Blinding No analyses required being blind to groups Reporting for specific materials, systems and methods We require information from authors about some types of materials, experimental systems and methods used in many studies. Here, indicate whether each material, system or method listed is relevant to your study. If you are not sure if a list item applies to your research, read the appropriate section before selecting a response. Materials & experimental systems Methods n/a Involved in the study n/a Involved in the study Antibodies ChIP-seq Eukaryotic cell lines Flow cytometry Palaeontology and archaeology MRI-based neuroimaging Animals and other organisms Human research participants Clinical data Dual use research of concern Animals and other organisms Policy information about studies involving animals; ARRIVE guidelines recommended for reporting animal research Laboratory animals The care and collection of a laboratory female zebra finch sample was done under an approved IACUC protocol at the Rockefeller University. The male laboratory sample was obtained from the approval mention of the previous reference genome in Warren et al 2011 Nature. Wild animals Samples of the 15 other species were collected from wild animals, with approved permits of the source institutions and local governments involved. These sources, persons with the permits, geographic location, sex and relative age are listed in Supplementary Table 8 and the BioSample submissions in NCBI and ENA. Field-collected samples No laboratory work was conducted on field-collected animals Ethics oversight Rockefeller University for the zebra finch species. Note that full information on the approval of the study protocol must also be provided in the manuscript. April 2020 2