AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication Baoxing Songa,1 , Santiago Marco-Solab,c , Miquel Moretob,d , Lynn Johnsona, Edward S. Bucklera,e,f,1 , and Michelle C. Stitzera,g,1 a Institute for Genomic Diversity, Cornell University, Ithaca, NY 14853; bDepartment of Computer Sciences, Barcelona Supercomputing Center, Barcelona 08034, Spain; cDepartament d’Arquitectura de Computadors i Sistemes Operatius, Universitat Auto  noma de Barcelona, Barcelona 08193, Spain; dDepartament d’Arquitectura de Computadors, Universitat Polite cnica de Catalunya, Barcelona 08034, Spain; eSection of Plant Breeding and Genetics, Cornell University, Ithaca, NY 14853; fAgricultural Research Service, US Department of Agriculture, Ithaca, NY 14853; and gDepartment of Molecular Biology and Genetics, Cornell University, Ithaca, NY 14853 Edited by Douglas Soltis, University of Florida, Gainesville, FL; received July 23, 2021; accepted November 15, 2021 Millions of species are currently being sequenced, and their duplicated copies by deletion or pseudogenization (5) and chro- genomes are being compared. Many of them have more complex mosomal rearrangement, further complicating genome genomes than model systems and raise novel challenges for alignments. genome alignment. Widely used local alignment strategies often Seed-and-extend local alignment strategies (6, 7) have been produce limited or incongruous results when applied to genomes widely successful for comparing model genomes. Such strate- with dispersed repeats, long indels, and highly diverse sequences. gies generally trigger alignments with pairs of highly similar Moreover, alignment using many-to-many or reciprocal best hit k-mers (seeds) from two genomes. This strategy can trigger approaches conflicts with well-studied patterns between species false alignments when aligning genomes with many dispersed repeats or even fail to generate an alignment when repetitive EVOLUTION with different rounds of whole-genome duplication. Here, we introduce Anchored Wavefront alignment (AnchorWave), which genome sequences are masked. The seed-and-extend strategy performs whole-genome duplication–informed collinear anchor often fails when aligning regulatory elements with essential identification between genomes and performs base pair–resolved functions, even though those sequences are expected to be con- global alignment for collinear blocks using a two-piece affine gap served between species. For example, the core motifs of tran- cost strategy. This strategy enables AnchorWave to precisely iden- scription factor–binding sites (TFBSs) (8–10) are 6.8 base pairs tify multikilobase indels generated by transposable element (TE) (bp) on average— much shorter than the seed size used in presence/absence variants (PAVs). When aligning two maize genome alignment. Furthermore, the presence of highly diverse genomes, AnchorWave successfully recalled 87% of previously fragments can limit alignment extension (6) and confound local reported TE PAVs. By contrast, other genome alignment tools showed low power for TE PAV recall. AnchorWave precisely aligns up to three times more of the genome as position matches or Significance indels than the closest competitive approach when comparing diverse genomes. Moreover, AnchorWave recalls transcription One fundamental analysis needed to interpret genome factor–binding sites at a rate of 1.05- to 74.85-fold higher than assemblies is genome alignment. Yet, accurately aligning other tools with significantly lower false-positive alignments. regulatory and transposon regions outside of genes remains AnchorWave complements available genome alignment tools by challenging. We introduce Anchored Wavefront alignment showing obvious improvement when applied to genomes with (AnchorWave), which implements a genome duplication dispersed repeats, active TEs, high sequence diversity, and whole- informed longest path algorithm to identify collinear regions genome duplication variation. and performs base pair–resolved, end-to-end alignment for collinear blocks using an efficient two-piece affine gap cost sensitive genome alignment j whole-genome duplication j genome strategy. AnchorWave improves the alignment under a num- comparison j transposable element variation j regulatory element ber of scenarios: genomes with high similarity, large alignment genomes with high transposable element activity, genomes with many inversions, and alignments between species with deeper evolutionary divergence and different whole- G enome alignment tools are fundamental for comparative evolutionary analysis. Unlike initial genome sequencing efforts, which concentrated on cost-effective sequencing of genome duplication histories. Potential use cases include genome comparison for evolutionary analysis of nongenic sequences and population genetics of taxa with large, model species, fulfilling the goal of sequencing a million repeat-rich genomes. eukaryotic reference genomes (1) adds many species with larger or repeat-rich genomes (2). Aligning those genomes pro- Author contributions: B.S., E.S.B., and M.C.S. designed research; B.S. and M.C.S. vides a revolutionary opportunity to understand the evolution performed research; B.S., S.M.-S., M.M., L.J., and E.S.B. contributed new reagents/ of eukaryotic genomes. Pipelines that have successfully aligned analytic tools; B.S. and M.C.S. analyzed data; and B.S., E.S.B., and M.C.S. wrote the genomes of model species often do not work well among paper. genomes with many complex variants, especially plant species. The authors declare no competing interest. Although the alignment of genic regions is bolstered by their This article is a PNAS Direct Submission. This article is distributed under Creative Commons Attribution-NonCommercial- modest length and conservation of amino acid residues, genes NoDerivatives License 4.0 (CC BY-NC-ND). comprise only a minority of the nucleotides in a genome. Distal 1 To whom correspondence may be addressed. Email: esb33@cornell.edu, mcs368@ regulatory regions can be recalcitrant to alignment, often cornell.edu, or bs674@cornell.edu. pushed far away from the genes they regulate by transposable This article contains supporting information online at http://www.pnas.org/lookup/ element (TE) insertions (3, 4). In addition, recursive whole- suppl/doi:10.1073/pnas.2113075119/-/DCSupplemental. genome duplications (WGD) result in fractionation of Published December 21, 2021. PNAS 2022 Vol. 119 No. 1 e2113075119 https://doi.org/10.1073/pnas.2113075119 j 1 of 9 alignment. In addition, alignment using the affine gap cost step 4). AnchorWave concatenates base pair sequence alignment strategy does not model mechanisms that generate indels of dif- for each anchor and interanchor region and outputs the alignment ferent length distributions (11), and alignment extension termi- in Multiple Alignment Format (MAF) (Fig. 1, step 5). nates in front of long indels (i.e., TE presence/absence variants To benchmark AnchorWave, we used partially synthetic and [PAVs]). Finally, most genome alignment approaches generate real genomes under a number of scenarios: small genomes with many-to-many alignments or are limited to one-to-one align- high similarity, large genomes with high TE activity, genomes ments (12), which may not reflect the true evolutionary history with many inversions, and alignments between species with when comparing taxa with unshared WGD. varying evolutionary divergence and WGD histories. We Here, we developed AnchorWave (Anchored Wavefront focused on the alignment performance in terms of TE align- alignment), a whole-genome alignment method that utilizes ment, sensitivity in putative regulatory sequence, and computa- genome collinearity and can be guided by differing levels of tional resources. To test AnchorWave for the alignment of WGD to perform sensitive sequence alignment with high accu- highly similar genomes, we synthesized benchmark alignments racy and recall long indels. These features provide a significant by introducing variant calls of 18 Arabidopsis (Arabidopsis thali- improvement compared to current methods when aligning ana) accessions (23) to the TAIR10 reference genome. In these genomes with enriched dispersed repeats, high sequence diver- benchmark alignments, variants account for 2.13 to 3.60% of sity, high transposon activity, or WGD variation. Some of these genome sites for each accession, including 1.22 to 2.38% of complex genomic variations that are challenging for alignment genome sites caused by variants longer than 50 bp (SI can be found in vertebrates and insects (13–15) but are wide- Appendix, Fig. S1). We synthesized genomes with these variant spread among plant genomes (16). calls and aligned them against the TAIR10 reference genome using AnchorWave, minimap2 (19), LAST (12), MUMmer4 Results (24), and GSAlign (25) and compared the newly generated alignments with benchmark alignments. AnchorWave was the AnchorWave leverages collinear regions to improve genome only algorithm that aligned chromosomes end to end and alignment (Fig. 1). Syntenic or collinear arrangements among aligned highly diverse fragments that were not aligned in the taxa have been investigated by aligning protein-coding gene benchmark (SI Appendix, Fig. S2), leading to a slight decrease sequences (17). These collinear blocks are parsimoniously in precision with only minimap2 ranking higher. AnchorWave interpreted as being derived from a shared ancestor (17, 18). had the highest F-score and recall (Fig. 2A and SI Appendix, AnchorWave takes the reference genome sequence and gene Fig. S3 and Supporting Note 1). annotation as input and extracts the reference full-length cod- To evaluate the performance for detecting long indels in ing sequence (CDS) to use as anchors. We use a splice-aware repeat-rich genomes, we developed a benchmark by removing sequence alignment program [minimap2 (19) or GMAP (20)] ∼60% of long terminal repeat (LTR) retrotransposons from to lift over the start and end positions of the reference full- the maize (Zea mays L.) B73 v4 assembly (26). This synthetic length CDS to the query genome (Fig. 1, step 1). AnchorWave TE-removed genome had 84,271 deletions with lengths ranging then identifies collinear anchors using one of three user- from 1,144 to 33,730 bp (SI Appendix, Fig. S4). We counted the specified algorithm options (Fig. 1, step 2) and uses a newly number of TE deletions that could be correctly recalled by implemented two-piece affine gap cost strategy in a Wavefront aligning the TE-removed genome against the reference genome Algorithm (WFA) library (21) to perform alignment for each (SI Appendix, Supporting Note 2). The genome alignment anchor and interanchor interval (Fig. 1, step 4). Some interan- results from GSAlign, LAST, and MUMmer4 did not generate chor regions cannot be aligned via WFA due to high computa- any variant longer than 1 kilobase pair (Kbp). Minimap2 (19) tional costs. For these situations, AnchorWave either identifies recalled ∼21% of these long deletions correctly, likely benefit- novel anchors within long interanchor regions (Fig. 1, step 3) ing from the usage of global alignment between adjacent or, for those that cannot be split by novel anchors, aligns them anchors in a chain. AnchorWave recalled ∼95% of these dele- using either the dynamic programming global sequence align- tions correctly (Fig. 2B). ment function ksw_extd2 implemented in the minimap2 library To evaluate the performance of the alignment approaches (19) or a reimplemented sliding window approach (22) (Fig. 1, using a realistic polymorphism landscape where long indels are || || | - Fig. 1. Principle of the AnchorWave process. AnchorWave identifies collinear regions via conserved anchors (here, full-length CDS) and breaks collinear regions into shorter fragments (i.e., anchor and interanchor intervals). By merging shorter intervals together after performing sensitive sequence align- ment via a two-piece affine gap cost global sequence alignment strategy, AnchorWave generates a whole-genome alignment. AnchorWave implements commands to guide collinear block identification with or without chromosomal rearrangements and provides options to use known WGD to inform the alignment. 2 of 9 j PNAS Song et al. https://doi.org/10.1073/pnas.2113075119 AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication A 1.00 B 0.75 F score 0.50 precision recall 0.25 0.00 An mi mi mi LA MU GS ni m ni m ni m c ST Ali Mm ho gn ap ap ap on rW er4 2_ 2_ 2_ e av as as as to e m5 m1 m2 on 0 0 e Test of Arabidopsis genome alignment C D TE_PAV non TE_PAV gap unaligned 1.00 0.75 0.50 Propor 0.25 0.00 An mi mi mi MU GS nim nim nim ch Ali Mm orW gn ap a ap p2 e r4 2_ 2_ av _a as as e sm m5 m2 10 0 EVOLUTION Fig. 2. Comparison of genome alignment tools using genomes of different individuals in the same species. (A) Comparison of the performance of genome alignment tools at variant sites of 18 Arabidopsis accession alignment benchmarks. Genome alignments were performed using minimap2 with preset options asm5, asm10, and asm20, which terminate extension in regions with 5, 10, and 20% sequence divergence, respectively. GSAlign and MUM- mer4 alignments were performed with default parameters. LAST genome alignments with default parameters were termed as LAST many-to-many. LAST many-to-many alignments were processed with a chain and net procedure to generate LAST many-to-one alignments (each query genome nucleotide may be aligned multiple times, while each reference nucleotide can be aligned up to one time). LAST many-to-one alignments were filtered to generate LAST one-to-one alignments (SI Appendix, Supporting Note 1). (B) Recall and precision of TE deletions by aligning the TE-removed maize B73 genome against the reference genome. MUMmer4, GSAlign, and LAST one-to-one had zero recall ratio. (C) Overview of the maize B73 genome sites aligned to the maize Mo17 genome. In those TE regions which were previously reported as present in B73 and absent in Mo17 (TE PAV on the legend), no position match alignments are expected. A higher number of position matches in these regions (striped orange) indicates a higher false-positive ratio. (D) Two inversions were located using AnchorWave between the maize B73-Ab10 assembly and the B73 v4 reference genome. mixed together with single nucleotide polymorphisms (SNPs) cases of previously reported inversions using AnchorWave are and short indels, we aligned the genomes of two maize individ- included in SI Appendix, Supporting Note 4 and Dataset S2. uals [B73 v4 and Mo17 (27), SI Appendix, Supporting Note 3]. Genome collinearity breaks down as species diverge and as AnchorWave aligned 61% of the B73 genome as position match genomes fractionate after polyploidy events (30). We evaluated (defined as an ungapped alignment, either matched or mis- AnchorWave by aligning a number of representative genome- matched nucleotides, SI Appendix, Fig. S5), which is compara- s—autotetraploids, paleopolyploids, and genomes separated by ble to other tools (Fig. 2C). However, AnchorWave produced multiple rounds of WGD. the lowest number of position matches in TE insertions previ- Historically, genome assemblies have been produced as hap- ously reported to be specific to B73 (28). This suggests that loid chromosomes after collapsing heterozygous regions. AnchorWave generated the fewest false positives among com- Recent sequencing advances allow the haplotype-resolved pared tools, as these sites should not be matched (orange- assembly of polyploid genomes as shown for the tetraploid striped region of Fig. 2C). Moreover, AnchorWave aligned potato (Solanum tuberosum) (31). We aligned the tomato (Sola- much (37%) of the genome as deletions (gaps). Such gaps are num lycopersicum) genome to the potato genome, allowing expected for alignments between maize individuals in whom each tomato anchor to define four potato anchors. Anchor- indel variation arises from TE PAVs. Anderson et al. (28) Wave aligned 63.5% of the tomato genome as position matches resolved a subset (15,182) of TE PAVs between these two or indels in the potato genome, ranking third highest, lower genomes to base pair resolution using an alignment approach than that generated via LAST many-to-many (76.2%) and anchored on syntenic genes (28). AnchorWave increased this LAST many-to-one (66.5%) (SI Appendix, Fig. S7 and number of TE PAVs to 28,321, recovering 87% (13,181) (28). Supporting Note 5). Polyploidy is rare in vertebrates, but gold- Other tools had almost zero recall ratios of TE PAVs (SI fish (Carassius auratus) (32) have undergone a WGD ∼14 Mya. Appendix, Table S1), as they generated few gapped alignments Goldfish have twice as many chromosomes as zebrafish (Danio (Fig. 2C). rerio) (33), but the size of their genome assemblies are similar The frequent presence of inversions in eukaryotic genomes (∼1.3 gigabase pairs [Gbp]), suggesting that extensive DNA poses an obstacle to the end-to-end alignment of chromosomes. loss has occurred along the goldfish lineage. Using parameters By incorporating anchor strand information into the collinear that allowed each zebrafish anchor to define up to two collinear identification approach, AnchorWave efficiently identifies inver- blocks in goldfish, AnchorWave aligned ∼82.7% of the zebra- sions. As an example, we show the gene-level resolution of two fish genome sequence as position matches or indels (SI neighboring inversions of abnormal chromosome 10 of maize, a Appendix, Supporting Note 6), over twice as much as the cytologically known inversion that carries a meiotic drive locus second-highest, generated via LAST many-to-many. Often, the (29) (Fig. 2D, SI Appendix, Fig. S6, and Dataset S1). Other recall reduction to diploidy results in fractionation of subgenomes, a Song et al. PNAS j 3 of 9 AnchorWave: Sensitive alignment of genomes with high sequence diversity, https://doi.org/10.1073/pnas.2113075119 extensive structural polymorphism, and whole-genome duplication signal observed in older WGDs. The maize lineage has under- often reflect indels between genomes, as most have inserted gone a WGD since its divergence with sorghum (34), but subse- more recently than species divergence. While 70 to 85% (38, 39) quent chromosomal fusions resulted in these species having the of the maize genome is composed of structurally recognizable same chromosome number (n = 10). Allowing up to two collin- TEs, almost all have estimated insertion times more recent ear paths for each sorghum anchor (Fig. 3A and SI Appendix, than the divergence from sorghum (Fig. 3B). AnchorWave Fig. S8 and Supporting Note 7), AnchorWave aligned a signifi- aligns less of the maize genome as position matches in sor- cantly larger proportion of the maize genome as position match ghum than any LAST alignment and minimizes position or indels, 3.4 times that of the second-highest, generated via matches located in maize TE regions (orange hatched region; LAST many-to-many (Fig. 3C). Multiple rounds of WGDs can Fig. 3C). Other investigated species lack information about further complicate alignment. Soybeans (Glycine max) and the ages of individual TEs, but the turnover of LTR retro- common beans (Phaseolus vulgaris) share an ancient WGD ∼56 transposons suggests that many TEs in tomato and soybean Mya (35) and then diverged ∼19 Mya (35). Subsequently, the are younger than the divergence from the species to which we soybean lineage had an allotetraploidy event ∼13 Mya (36). We compare them (40). In the tomato–potato comparison, aligned the common bean genome to the soybean genome, although LAST approaches aligned more of the genome than allowing each common bean anchor to define two soybean AnchorWave, much of that was position matches in TE anchors to account for the unshared WGD. AnchorWave sequence (SI Appendix, Fig. S7). In the soybean comparison, aligned 78.06% of the soybean genome in collinear blocks with AnchorWave not only aligns more of the genome but aligned the common bean, which is ranked as the highest (SI Appendix, a smaller ratio in TE regions as position matches compared to Fig. S9 and Supporting Note 8). For the above alignments, non-TE regions (SI Appendix, Fig. S9). These TE comparisons AnchorWave aligned more base pairs as indels than all the show AnchorWave reduces false-positive alignments in repeti- other alignment approaches and fewer base pairs as position tive regions with low expected conservation. matches than LAST many-to-many and LAST many-to-one. Conservation in genic regions is expected across species, but An alignment that covers a large proportion of the genome since AnchorWave uses genes to guide genome alignments, we may not be optimal, as false-positive alignment can always did not use coding sequences to assess alignment quality. increase this proportion. Here, we assessed alignment quality Instead, we turn to putative regulatory sequences, which are based on biologically informed expectations about lack of expected to be more conserved and less affected by absence sequence conservation in TEs and sequence conservation in variants than the whole–genome background (41). The regula- putative regulatory sequences. TEs evolve more rapidly than tory sequence of tomato and soybean has been investigated by their host genomes, with independent TE movement and identifying accessible chromatin regions through ATAC-seq sequence divergence among species (37). Thus, TE regions (4, 42). AnchorWave aligned a higher proportion of accessible A chr4 chr5 B 1500 250M Cumulative TE length (Mbp) 200M 150M chr4 1000 100M maize 50M 0M 500 200M chr5 150M 100M 0 0M 20 40 60 0M 20 40 60 0 5 10 15 20 M M M M M M sorghum TE age (My) C TE non TE D 12.5 alignment to position match ratio gap Position match ratio in TFBS in non TFBS alignment ze 10.0 1.00 AnchorWave 0.75 minimap2_asm5 7.5 minimap2_asm10 0.50 minimap2_asm20 LAST many to many 0.25 5.0 LAST many to one Propor LAST one to one 0.00 MUMmer4 GSAlign 2.5 An mi mi mi LA LA MU GS nim nim nim ch ST ST Ali Mm orW gn ap ap ap ma ma e 2_ 2_ 2_ r4 0.0 av n n as as as e m5 m1 m2 0.0 0.2 0. 5 0.7 1.0 0 0 0 5 0 5 0 ny TFBS recall Fig. 3. A comparison of different genome alignment tools for aligning the maize B73 v4 and sorghum genomes. (A) The identified collinear anchors between the maize B73 v4 assembly and sorghum genome on chr4 and chr5. Each dot is plotted based on the start coordinate of the reference genome and query genome of each anchor. Collinear anchors on the same strand between the reference genome and query genome are shown in blue, otherwise red. (B) Cumulative distribution of TE length versus TE age in the B73 v4 maize genome with age measured in millions of years (My). The dashed line indi- cates 12 My, the estimated divergence time of maize and sorghum (34). TE age data are from Stitzer et al. (39); 371 TEs older than 20 My were not plot- ted, and the total length of these 371 TEs is 531 Kbp. (C) Sequence alignment between the maize B73 v4 genome and the sorghum genome. Minimap2, MUMmer4, and GSAlign generated many-to-many alignments. Since most maize TEs are not shared with sorghum, a higher number of position matches in maize TE regions (striped orange) indicates a higher false-positive ratio. AnchorWave aligns 88.7% of the maize genome to the sorghum genome, while the second highest is 28.0% generated by LAST many-to-many. (D) Comparison of the proportion of sites in maize TFBS that were aligned as a posi- tion match (recall) and the position match ratios (number of position match sites to number of aligned sites) in TFBS versus non-TFBS regions. 4 of 9 j PNAS Song et al. https://doi.org/10.1073/pnas.2113075119 AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication chromatin regions as position match than any other aligner. does not need highly similar seeds to trigger alignment in Moreover, AnchorWave alignments showed higher sequence diverse regions and increases the sequence alignment sensitiv- identity across these regions, where tomato accessible chroma- ity. AnchorWave uses collinear anchors to guide the alignment tin regions show a 2.21-fold higher match ratio compared to of repeat elements; thus, local duplications or translocations the whole-genome background (SI Appendix, Fig. S7), and without enough anchors to be identified as a collinear block are soybean accessible chromatin regions show 1.82-fold higher (SI expected to be aligned as indels. Chromosomal fusions after Appendix, Fig. S9). All other aligners show nearly equal match WGD can complicate the separation of subgenomes, but the ratios in ATAC-seq regions with the whole genome back- WGD informed collinear blocks identification function in ground. TFBSs have been identified directly for maize via chro- AnchorWave improves the alignment for genomes with WGD matin immunoprecipitation sequencing (ChIP-seq). (8). variation. The collinear anchor guided alignment, and the two- AnchorWave aligned the highest proportion of maize TFBS piece affine gap cost global sequence alignment strategy regions (8) as position matches with sorghum and generated a enabled the recall of long indels using AnchorWave. Based on higher match ratio in TFBS regions compared to the rest of the the observation that long indels are generally derived from TE genome (Fig. 3D), suggesting that AnchorWave aligns these movement, while different mechanisms introduce shorter indel conserved functional sequences more accurately and mutations (11), the two-piece affine gap cost strategy enables comprehensively. the alignment of short and long indels. The performance of We further explored the performance of AnchorWave on other genome alignment implementations may improve by fine- various species at differing levels of divergence and polyploidy tuning gap parameters (47), but, with the exception of mini- (43) and investigated the proportion of genomes that could be map2, all other tested tools align long indels and short indels identified in collinear blocks (SI Appendix, Supporting Note 9 using the same gap penalty profile. This precludes parameter and Table S2). Maize and rice are two grasses with a greater tuning as a solution to simultaneously optimize the alignment (50 My) phylogenetic distance than that between maize and of long indels and short indels. sorghum (34). AnchorWave identified 85.5% of the maize We highlight AnchorWave’s ability to align several polyploid genome as collinear with the rice genome (44) (Oryza sativa, species to diploid relatives (Fig. 3 and SI Appendix, Figs. S8 IRGSP-1.0), 3.5% less than that between maize and sorghum. and S9). Many of these species are ancient allopolyploids, and extensive fractionationation of subgenomes and chromosome EVOLUTION Arabidopsis and chocolate diverged more recently than did Ara- bidopsis and grape (45) (Vitis vinifera, 12X), and while 61.8% of fusions have occurred to generate the extant genomes (48). For the Arabidopsis genome was identified as collinear with choco- example, thousands of maize genes have remained duplicated late, only 50.0% was collinear with grape. Fractionation after relative to sorghum (48), but, without an explicit model of WGD breaks collinear blocks, so we examined extreme paleo- WGD, most genome alignment strategies fail to align these polyploid comparisons. Two rounds of WGD and one round of regions. By generating alignment paths in each subgenome, whole-genome triplication occurred between Arabidopsis and AnchorWave allows interpretation of the conservation and evo- tomato, and AnchorWave detected only 36.3% of the Arabidop- lution of these genes and their local regulatory regions. As sis genome as collinear with tomato, despite a similar diver- much of this putative regulatory sequence is embedded gence time as between Arabidopsis and grape. Six rounds of between repetitive TEs, AnchorWave reduces false-positive WGDs separate maize and banana (46) (Musa schizocarp), and alignments in this repetitive space. although the divergence time between maize and banana is While AnchorWave provides improved alignments for many smaller than that between Arabidopsis and tomato, only 11.6% complicated but real issues in genomics, relying on the exis- of the maize genome was classified as collinear with banana. tence of collinear regions can be a limitation. Collinearity is The joint action of sequence divergence and WGD resolution more limited in vertebrate systems than plants (49), yet, in com- parisons between the human and mouse genomes, AnchorWave affects the existence of collinear blocks between species. provided alignment over putative regulatory sequences (SI AnchorWave exploits the latest advances in sequence align- Appendix, Supporting Note 10). There are limits to aligning ment using the WFA (21), which reduces the overall memory noncollinear sequences—although the autosomes of humans and computational requirements for global pairwise sequence and chimpanzees share largely collinear genes, the gene- alignment. To limit memory usage for long sequence compari- depauperate Y chromosome could not be aligned (SI Appendix, sons, AnchorWave implements two strategies: identifying novel Supporting Note 11). Additionally, technical limitations such as anchors within long interanchor regions and, for interanchor fragmented genome assemblies can prevent the identification regions that lack sufficient homology, approximating alignment of collinearity (50), although this will likely pose less of a prob- using either a banded approach or sliding windows. Alignments lem with advances in long-read sequencing. between Arabidopsis genomes typically used less than 10 Gb When comparing genomes with different rounds of WGD, memory and only minutes of computing time, while the most AnchorWave significantly increased the proportion of the resource intensive maize–sorghum alignment required 85 Gb genome that was aligned compared to one-to-one alignments memory and 130 h (SI Appendix, Tables S3–S5), and alignment and reduced false-positive alignments compared to many-to- can be parallelized when more memory is available. Although many alignments (Fig. 3C). Compared to alignment approaches the execution time of AnchorWave is high for some experi- using a seed-and-extension strategy, AnchorWave increased ments, it remains comparable when considering that most other sensitivity for putative regulatory sequences and could recall methods fail to align many bases in the genome. long indel variants. AnchorWave’s collinear approach further reduces false-positive alignments from dispersed repeats. We Discussion showed that AnchorWave can generate whole-genome align- Genome evolution across the tree of life has resulted in species ments, facilitating studies of the evolution of regulatory ele- with vastly differentiated ploidy, chromosome number, and ments, TE polymorphisms, and chromosomal rearrangements genome organization. Beyond genome structure, species differ such as inversions. AnchorWave even allows duplicated collin- in the complexity and magnitude of nongenic repetitive sequen- ear blocks to be aligned, making it particularly relevant to ces. Despite this widespread variation, genomes are always plants, in which an estimated 35% of species are polyploids punctuated by evolutionarily conserved regulatory sequences (16). AnchorWave complements available genome alignment and genes. AnchorWave makes use of this conservation, utiliz- tools by improving the genome alignment of many plant ing gene collinearity to guide genome alignment. This approach species. Song et al. PNAS j 5 of 9 AnchorWave: Sensitive alignment of genomes with high sequence diversity, https://doi.org/10.1073/pnas.2113075119 extensive structural polymorphism, and whole-genome duplication Materials and Methods 4. For each me and mf pair, update mf’s chainScore based on the following: Collinear Anchor Identification. AnchorWave takes the reference genome in If matchScoref þ chainScoree > chainScoref : FASTA format, the reference genome gene annotation in GFF3 format, and chainScoref ¼ matchScoref þ chainScoree the query genome in FASTA format as input. AnchorWave extracts the full- preIndexf ¼ e: length CDS from the reference genome using the reference genome and Starting from the match that has the maximum chainScore, AnchorWave will annotation. The start and end positions of the reference full-length CDS to use the preIndex values to track back and produce a list of matches that give the query genome are lifted over using a splice-aware sequence alignment the highest score. program [minimap2 (19) was used in this manuscript]. AnchorWave then Longest path considering inversions. After sorting the matches in ascending implements a longest-path dynamic programming algorithm to identify col- order based on the reference anchor start positions, to consider inversions, we linear anchors. Base pair resolution sequence alignments within each anchor create currentScore, which is a cumulative score, and maxScore, which is the and interanchor region are conducted using the two-piece affine gap cost maximum value that currentScore has ever reached for each round of iteration. global sequence alignment strategy, and these alignments are concatenated When a match is encountered on the negative strand, we assign its match- together to generate the alignment for each collinear block. Score to currentScore. If the next match is on the negative strand and the A longest-path dynamic programming approach is applied to a pair of query start position is smaller than the query start position of the current one, chromosomes. On the reference chromosome is a list of n anchors: we add this matchScore to the currentScore. Otherwise, we subtract this matchScore from the currentScore. When the currentScore drops below 0, the Q ¼ q0 , q1 , …, qn1 : iteration is terminated, and the next iteration starts. If the maximum cumula- tive score (maxScore) is larger than a preset threshold for all matches in a kept On the query chromosome is a list of m lift over hits: group, we reverse their order in the list of matches. Those reversed matches have increasing query start positions and decreasing reference positions and T ¼ t0 , t1 , …, tm1 : thus remain on the diagonal. A similar longest-path dynamic programming algorithm is applied to find For each anchor and its hit, q and t, the start and end positions, respectively, an end-to-end chain as described for sequences without inversions or rear- are identified from the SAM file output from the splice-aware setting of mini- rangements, except that the end position of anchor qif is smaller than the start map2. We set up a list of o anchor matches: position of anchor qie on the reference chromosome when strande and strandf are negative. With 0 ≤ e < f < o, for a previous element, me = (qie, tje, M ¼ m0, m1 , …, mo 1 : matchScoree, strande, chainScoree, preIndexe) and a current element, mf = (qif, tjf, matchScoref, strandf, chainScoref, preIndexf), we iterate the list of matches, Individual matches are defined as the following: incrementing e and f, while the following conditions are true:  1. If the current match and previous match are on the negative strand, the m0 ¼ qi 0 , tj 0 , matchScore0 , strand0 , chainScore0 , preIndex0  m1 ¼ qi 1 , tj 1 , matchScore1 , strand1 , chainScore1 , preIndex1 end position of reference anchor qif is smaller than the start position of ref- … erence anchor qie. Otherwise, the end position of reference anchor qie is  mo 1 ¼ qio1 , tjo1 , matchScoreo1 , strando 1 , chainScoreo 1 , preIndexo 1 : smaller than the start position of reference anchor qif. 2. The end position of the query anchor tie is smaller than the start position of the query anchor tif. tjk is the lift over hit of qik identified via minimap2, 0 ≤ k < o, 0 ≤ i < n, 0 ≤ j < 3. For each me and mf pair, update the mf’s chainScore based on the following: m. matchScorek is the sequence similarity (ratio of the number of identical nucleotides to the length of the reference full-length CDS). If qik and tjk are on If ðmatchScoref þ chainScoree > chainScoref Þ : chainScoref ¼ matchScoref þ chainScoree the same strand, strandk is set as positive; otherwise, it is set as negative. chain- preIndexf ¼ e: Scorek is initialized with matchScorek. preIndexk is the index of the previous anchor on the chain and is initialized as 1, meaning that the current anchor Starting from the match that has the maximum chainScore, AnchorWave will is the first one on a chain. use the preIndex values to track back and output the list of anchor matches Three longest-path approaches have been developed for genomes with that give the highest score. different types of chromosomal rearrangements. The first two approaches For each pair of homologous chromosomes, at the base pair sequence align- expect chromosome-level assemblies and the homologous chromosomes in ment step, sequences between two anchors on opposite strands were skipped. the input reference and query files to be named identically, as is common for If the first collinear anchor is on the negative strand, the sequence upstream intraspecific comparisons. Those two approaches try to align each pair of the first anchor would not be aligned. If the last collinear anchor is on the neg- homologous chromosomes from the first base pair to the last base pair. The ative strand, the sequence downstream the last anchor would not be aligned. Longest path considering inversions, rearrangements, and WGDs. Anchor- third approach is the most flexible and identifies collinear blocks using a local Wave implements a function to constrain the alignment depths for both the longest-path algorithm and then performs base pair resolution global align- reference and query genome, which is useful when there may be multiple col- ment for each collinear block. The third method gains this flexibility by not linear paths (i.e., genomes with chromosomal translocations, chromosome requiring prior information about homologous chromosomes and may gener- fusions, and varying numbers of rounds of WGDs). This function does not ate alignments for fewer base pairs compared to the first two approaches. assume there are homologous chromosome pairs. Instead, it identifies homol- The first approach should be used when no chromosomal rearrangements ogous collinear blocks by applying a local longest-path algorithm on anchors. or inversions occur between homologous chromosomes and the second when The matches are sorted in ascending order based on the reference anchor inversions disrupt collinearity. When a user has limited background knowl- start positions. Then, with 0 ≤ e < f < o, for a previous element, me = (qie, tje, edge on chromosomal contiguity of the genomes being aligned, the third matchScoree, strande, chainScoree) and a current element, mf = (qif, tjf, approach can be used as the default choice. matchScoref, strandf, chainScoref), iterate the list of matches, incrementing e Longest-path approach for genome sequences without inversions or rear- and f, while the following conditions are true: rangements. The target is to select a subset of positive-strand, nonoverlap- ping anchor matches from M that give a maximum value of chainScore in 1. The current match and previous match are on the same strand. which the positions of anchors on the reference and query chromosomes 2. The end position of anchor qie is smaller than the start position of anchor increase. The matches in M are first sorted in ascending order by reference qif on the reference chromosome. anchor start positions. Then, with 0 ≤ e < f < o, for a previous element, me = 3. If strande is positive, the end position of anchor tie is smaller than the start (qie, tje, matchScoree, strande, chainScoree, preIndexe) and for a current ele- position of anchor tif on the query chromosome. If strande is negative, the ment, mf = (qif, tjf, matchScoref, strandf, chainScoref, preIndexf), the list of end position of anchor tie is smaller than the start position of anchor tif on matches is iterated, incrementing e and f, while the following conditions the query chromosome. are true: 4. For each me and mf pair, update the mf’s chainScore based on the following: 1. The strands of mf and me are positive. 2. The end position of reference anchor qie is smaller than the start position if ðmatchScoref þ maxð0, chainScoree þ O þ E  NumberOfGapsðe, f ÞÞ of reference anchor qif. > chainScoref Þ: 3. The end position of the query anchor tie is smaller than the start position chainScoref ¼ matchScoref þ maxð0, chainScoree þ O þ E  of the query anchor tif. NumberOfGapsðe, fÞÞ 6 of 9 j PNAS Song et al. https://doi.org/10.1073/pnas.2113075119 AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication preIndexf ¼ e; Identifying Additional Anchors to Reduce the Size of Interanchor Intervals. in which O is a gap opening penalty, and E is a gap To improve sequence alignment quality and computational efficiency, addi- tional anchors are needed in long interanchor intervals. AnchorWave used the extension penalty. Let a be the number of anchors “mm_map” function of the minimap2 library with a single-piece affine gap between anchor qie and qif, and let b be the number cost setting to perform local alignment in collinear interanchor regions. In of anchors between tje and tjf. each interanchor region, AnchorWave selects the mm_map primary alignment and specifies this as a novel anchor. This step is iterated until either all interan- NumberOfGapsðe, fÞ ¼ ða þ b þ ja  bjÞ=d: chor intervals are shorter than a settable threshold, the sequence similarity of NumberOfGaps reflects the number of anchor mismatches and PAVs of the new primary alignment is lower than a threshold (by default, we do not anchors. ja  bj is used to penalize differences in the number of anchors filter novel anchors using similarity, and this is set as 0), or no new mm_map between two sequences (in contrast to a similar number of mismatched matches can be found. anchors). d is a settable parameter (by default: 3) to calculate a normalized value for the three parts (a, b, ja  bj); higher values of d allow more continu- Base Pair Resolution Sequence Alignment Using the Two-Piece Affine Gap ous chaining across gaps, which may introduce false-positive chains. To avoid Cost Strategy. Based on the assumption that pairs of sequences in each extremely long indels or mismatches, if a or b is larger than a settable thresh- anchor or interanchor region are passed down from a common ancestor, old D (by default: 30), the current chaining stops. The values of d, D, O, and E AnchorWave performs base pair resolution global sequence alignment for were set empirically by manually comparing the dot plots of minimap2 splice- each anchor and interanchor interval using the two-piece affine gap cost aware mapping results and collinear anchors for all genomes used in this man- strategy (19). We implemented the two-piece affine gap cost strategy in the uscript, and default values were used throughout the manuscript. WFA (21) library, and sequence alignments are conducted using WFA by AnchorWave selects the chain with the maximum chainScore. If the maxi- default. If the WFA library requires more memory than a preset threshold, the mum chainScore is larger than a settable threshold, then output the chain. To “ksw_extd2” functions implemented in minimap2 are called using a calcu- align genomes with WGD variants, all the reference and query anchors that lated bandwidth. The memory cost threshold of the WFA library is calculated fall into the chain range will be counted. If there are multiple hits for a refer- using the “-w” parameter of AnchorWave. The bandwidth of ksw_extd2 is cal- ence anchor in the chain range, the reference anchor is counted once. Anchor culated using the “-w” parameter of AnchorWave and anchor/interanchor matches that fall into chain ranges that are counted as larger than a settable sequence length. For longer sequences with calculated ksw_extd2 bandwidth alignment depth threshold are marked and will not be used for the next smaller than “-w,” we implement the two-piece affine gap cost strategy with round of iteration. Then, the next iteration starts using all the unmarked a sliding window approach (22), which generates approximate sequence matches until the maximum chainScore is smaller than a settable threshold. alignments. The sliding window size (-w) was set as 38,000 in this study, which EVOLUTION The user-settable alignment depth thresholds are the parameters “-R” and was also used as the minimum bandwidth for ksw_extd2. “-Q.” They are used to control the alignment depth for the reference genome Long indels are generally derived from TE movements, while different and query genome, especially when the reference and query genomes have mechanisms introduce shorter indel mutations (11). Since most available WGD variants. “-R” is the maximum alignment depth for the reference genome alignment tools align long indels and short indels using the same gap genome, and “-Q” is the maximum alignment depth for the query genome. penalty profile, long indels generally fail to be aligned. Here, we implemented For example, the maize genome has an additional round of WGD compared the two-piece gap cost strategy to align long indels following equation 4 to the sorghum genome, so we set the maximum alignment depth of the described by Li (19). Let the reference r = r0, r1 … rn1 and the query q = q0, maize genome as 1 and the maximum alignment depth of the sorghum q1 … qm1 be a pair of sequence fragments in an anchor or interanchor inter- genome as 2. We use parameters “-R 1 -Q 2” to align the sorghum genome val from the reference genome and the query genome, respectively, with against the maize genome. length jrj = n and jqj = m. O1 is the first piece affine gap open penalty, E1 is the first piece affine extend open penalty, O2 is the second piece affine gap open penalty, and E2 is the second piece affine extend open penalty. Let l be Filtering Anchors to Improve Alignment Quality. A correct lift over of full- the gap length; minfO1 + jlj × E1, O2 + jlj × E2g was used as the indel penalty length CDSs from the reference genome to the query genome is central to the for the dynamic programming sequence alignment approach. We always AnchorWave pipeline. If the full-length CDSs of two or more genes are identi- assume O1 + E1 < O2 + E2. On the condition that E1 > E2, it applies cost O1 + jlj × cal, AnchorWave ignores all of them in the subsequent analysis because of E1 to gaps shorter than d(O2  O1)/(E1  E2)e and applies O2 + jlj × E2 to longer their ambiguous mappings. We used the splice-aware function of minimap2 gaps. to lift over full-length CDSs across this manuscript. Minimap2 misaligns small The values of the first piece affine gap cost and second piece affine gap exons (https://github.com/lh3/minimap2#limitations). To reduce this side cost were selected based on the finding that TE copies are longer than 50 bp effect, when extracting full-length CDSs, AnchorWave ignores CDS exon (39) (SI Appendix, Fig. S10). records <20 bp (although this limit is a user-settable parameter). All exons can be used when more accurate splice-aware aligners [e.g., GMAP (20)] are used Data Sources and Methods for Validation. All the AnchorWave results shown for lift over. in this manuscript are repeatable using AnchorWave version 1.0.0. The We first identify full-length CDSs that might produce incorrect lift overs. In parameter “-IV” of the AnchorWave “genoAli” command was used to iden- this study, we used minimap2 to map extracted full-length CDSs to the refer- tify inversions between maize de novo assemblies against the maize B73 v4 ence genome and then rank the hits of each full-length CDS based on reference genome. To align the sorghum genome against the maize similarity (ratio of the number of identical DNA base pairs to the length of the genome, we used the “proali” command of AnchorWave with parameters full-length reference CDS). To minimize the effect of tandemly duplicated “-R 1 -Q 2” to utilize the knowledge that the maize lineage has been anchors, AnchorWave uses two thresholds: a number of mapping hits through a WGD since its divergence with sorghum (34). The value of “-Q,” threshold (e, 1 by default) and a similarity ratio threshold (y, 0.6 by default). “-R,” and “-e” were also adjusted for other polyploid and paleopolyploid Any full-length CDS with >e mapped hits on a single reference chromosome genome alignments using AnchorWave (see SI Appendix, Supporting Notes sequence is further investigated using the similarity ratio threshold. From a 1–11 for more details). similarity-sorted list of all hits of this full-length CDS, we calculate the similarity The values of asm5, asm10, and asm20 of the -x setting of minimap2 v2.16- ratio by dividing the e þ 1 hit similarity by the highest hit similarity. If this ratio r922 (19) were used separately. When aligning the TE-removed maize B73 is above the similarity ratio threshold (y), we drop this full-length CDS and its genome sequence using a single thread, the setting asm20 gave an hits on any chromosome. We then use the longest-path approach for genome “insufficient memory” error on a computer with 2 terabytes of available sequences without inversions or rearrangements to further filter the hits. We memory, and we did not obtain the corresponding result. compare the coordinates of remaining hits with the original GFF3 file; any full- The “lastal” function from the LAST toolkit v932 (12) was used to perform length CDS with different coordinates between the original GFF3 file and lift genome alignment with default parameters; the results were termed LAST over reference hits is placed on an unwanted list. many-to-many. Following previously described methods (51), the LAST many- We then use minimap2 to lift over extracted full-length CDSs to the query to-many results were transformed into PSL format using the “maf-convert” genome. Full-length CDSs in the unwanted list are not used. Anchors are fur- command of LAST, and the psl files were fed into the chain-net pipeline, ther filtered to reduce the impact of tandem duplications using the same “axtChain -linearGap=loose,” “chainMergeSort,” “chainPreNet,” “chainNet,” approach and parameters as described for the reference genome. After filter- “netToAxt,” “axtSort,” and “axtToMaf” in sequential order to generate the ing, the remaining anchors are fed into a user-specified longest-path algo- LAST many_to_one results. The LAST many_to_one results were further proc- rithm to identify collinear blocks. essed via the “last-split j maf-swap j last-split j maf-swap” command to Song et al. PNAS j 7 of 9 AnchorWave: Sensitive alignment of genomes with high sequence diversity, https://doi.org/10.1073/pnas.2113075119 extensive structural polymorphism, and whole-genome duplication generate LAST one_to_one results. “last-split” and “maf-swap” are compo- Data Availability. The source code of AnchorWave is hosted by GitHub at nents of the LAST genome alignment toolkits. https://github.com/baoxingsong/anchorwave (53). All the commands, plotting The parameters –sam-short of MUMmer4 (24) were used to produce scripts, and plot source data for this manuscript can be found at: https:// genome alignments in SAM format. We used the parameter “-fmt 1” of GSA- github.com/baoxingsong/genomeAlignment (54). lign (25) to perform genome alignments with default settings and output in MAF. ACKNOWLEDGMENTS. This project is supported by the United States Depart- To calculate the proportion of the reference genome that was aligned and ment of Agriculture Agricultural Research Service, NSF No. 1822330, NSF No. matched, all the alignments in MAF were reformatted into bam files using the 1854828, the European Union's Horizon 2020 Framework Programme under the DeepHealth project [825111], the European Union Regional Development “maf-convert sam” command of LAST and SAMtools v1.10 (52). We used the Fund within the framework of The European Regional Development Fund “depth” command of SAMtools to calculate how many base pairs of the refer- Operational Program of Catalonia 2014 to 2020 with a grant of 50% of total ence genome were aligned. We used the “samtools depth j awk '$3>0fprint cost eligible under the DRAC project [001-P-001723], and National Natural Sci- $0g' j wc -l” command to calculate how many base pairs of a reference ence Foundation of China No. 31900486. M.C.S. was supported by NSF Post- genome has a matched position in the query genome. doctoral Research Fellowship in Biology No. 1907343. M.M. was partially sup- Because of the limit of available computational resources (maximum of 2 ported by the Spanish Ministry of Economy, Industry, and Competitiveness terabytes of memory), we split the maize B73 v4 genome and maize Mo17 under Ramo  n y Cajal (RYC) fellowship number RYC-2016-21104. We thank the genome into individual chromosomes and performed alignments using LAST members of the E.S.B. laboratory (Cornell University) for helpful discussions. We thank Travis Wrightsman (Cornell University) for suggesting the name and minimap2 for each pair of homologous chromosomes independently. AnchorWave and submitting it to bioconda and Hequan Sun (Max Planck When aligning chromosome 1 using minimap2 asm20, we set “-w 19” to Institute for Plant Breeding Research) for sharing the tetraploid potato assem- reduce memory usage to less than 2 terabytes. The outputs of minimap2, bly prerelease. We thank Merritt Khaipho-Burch, Qi Sun, Cheng Zou, Minghui MUMmer4, and GSAlign were filtered as one-to-one alignment for subse- Wang, Zack Miller, Sara Miller (Cornell University), Jeffrey Ross-Ibarra (Univer- quent analysis using the “last-split j maf-swap j last-split j maf-swap” com- sity of California, Davis), and Shi Huang (University of California, San Diego) mands. More detailed description about the data and methods could be for helpful comments. found in SI Appendix, Supporting Notes 1–11. 1. H. A. Lewin et al., Earth BioGenome project: Sequencing life for the future of life. 27. S. Sun et al., Extensive intraspecific gene order and gene structural variations Proc. Natl. Acad. Sci. U.S.A. 115, 4325–4333 (2018). between Mo17 and other maize genomes. Nat. Genet. 50, 1289–1295 (2018). 2. M. Exposito-Alonso, H.-G. Drost, H. A. Burbano, D. Weigel, The Earth BioGenome 28. S. N. Anderson et al., Transposable elements contribute to dynamic genome content project: Opportunities and challenges for plant genomics and conservation. Plant J. in maize. Plant J. 100, 1052–1065 (2019). 102, 222–229 (2020). 29. R. J. Mroczek, J. R. Melo, A. C. Luce, E. N. Hiatt, R. K. Dawe, The maize Ab10 meiotic 3. B. Wei et al., Genome-wide characterization of non-reference transposons in crops drive system maps to supernumerary sequences in a large complex haplotype. Genet- suggests non-random insertion. BMC Genomics 17, 536 (2016). ics 174, 145–154 (2006). 4. Z. Lu et al., The prevalence, evolution and chromatin signatures of plant regulatory 30. H. Tang et al., Synteny and collinearity in plant genomes. Science 320, 486–488 elements. Nat. Plants 5, 1250–1259 (2019). (2008). 5. M. Freeling, M. J. Scanlon, J. E. Fowler, Fractionation and subfunctionalization fol- 31. H. Sun et al., Chromosome-scale and haplotype-resolved genome assembly of a tetra- lowing genome duplications: Mechanisms that drive gene content and their conse- ploid potato cultivar. bioRxiv [Preprint] https://doi.org/10.1101/2021.05.15.444292. quences. Curr. Opin. Genet. Dev. 35, 110–118 (2015). Accessed 17 May 2021. 6. S. F. Altschul, W. Gish, W. Miller, E. W. Myers, D. J. Lipman, Basic local alignment 32. Z. Chen et al., NISC Comparative Sequencing Program, De novo assembly of the gold- search tool. J. Mol. Biol. 215, 403–410 (1990). fish (Carassius auratus) genome and the evolution of genes after whole-genome 7. H. Li, N. Homer, A survey of sequence alignment algorithms for next-generation duplication. Sci. Adv. 5, eaav0547 (2019). sequencing. Brief. Bioinform. 11, 473–483 (2010). 33. K. Howe et al., The zebrafish reference genome sequence and its relationship to the 8. X. Tu et al., Reconstructing the maize leaf regulatory network using ChIP-seq data of human genome. Nature 496, 498–503 (2013). 104 transcription factors. Nat. Commun. 11, 5089 (2020). 34. Z. Swigonova  et al., Close split of sorghum and maize genome progenitors. Genome 9. R. C. O’Malley et al., Cistrome and epicistrome features shape the regulatory DNA Res. 14, 1916–1923 (2004). landscape. Cell 165, 1280–1292 (2016). 35. J. Schmutz et al., A reference genome for common bean and genome-wide analysis 10. B. Song et al., Conserved noncoding sequences provide insights into regulatory of dual domestications. Nat. Genet. 46, 707–713 (2014). sequence and loss of gene expression in maize. Genome Res. 31, 1245–1257 (2021). 36. J. Schmutz et al., Genome sequence of the palaeopolyploid soybean. Nature 463, 11. J. L. Bennetzen, J. Ma, K. M. Devos, Mechanisms of recent genome size variation in 178–183 (2010). flowering plants. Ann. Bot. 95, 127–132 (2005). 37. C. Vitte, O. Panaud, LTR retrotransposons and flowering plant genome size: Emer- 12. S. M. Kiełbasa, R. Wan, K. Sato, P. Horton, M. C. Frith, Adaptive seeds tame genomic gence of the increase/decrease model. Cytogenet. Genome Res. 110, 91–107 sequence comparison. Genome Res. 21, 487–493 (2011). (2005). 13. Z. Li et al., Multiple large-scale gene and genome duplications during the evolution 38. M. B. Hufford et al., De novo assembly, annotation, and comparative analysis of 26 diverse of hexapods. Proc. Natl. Acad. Sci. U.S.A. 115, 4713–4718 (2018). maize genomes. bioRxiv [Preprint] https://doi.org/10.1101/2021.01.14.426684. Accessed 16 14. G. Tsagkogeorga, J. Parker, E. Stupka, J. A. Cotton, S. J. Rossiter, Phylogenomic analy- January 2021. ses elucidate the evolutionary relationships of bats. Curr. Biol. 23, 2262–2267 (2013). 39. M. C. Stitzer, S. N. Anderson, N. M. Springer, J. Ross-Ibarra, The genomic ecosystem of 15. C. Sun et al., LTR retrotransposons contribute to genomic gigantism in plethodontid transposable elements in maize. PLoS Genet. 17, e1009768 (2021). salamanders. Genome Biol. Evol. 4, 168–183 (2012). 40. P. Jedlicka, M. Lexa, E. Kejnovsky, What can long terminal repeats tell us about the 16. T. E. Wood et al., The frequency of polyploid speciation in vascular plants. Proc. Natl. age of LTR retrotransposons, gene conversion and ectopic recombination? Front Acad. Sci. U.S.A. 106, 13875–13879 (2009). Plant Sci. 11, 644 (2020). 17. H. Tang et al., Screening synteny blocks in pairwise genome comparisons through 41. E. T. Dermitzakis, A. G. Clark, Evolution of transcription factor binding sites in mam- integer programming. BMC Bioinformatics 12, 102 (2011). malian gene regulatory regions: Conservation and turnover. Mol. Biol. Evol. 19, 18. D. Liu, M. Hunt, I. J. Tsai, Inferring synteny between genome assemblies: A systematic 1114–1121 (2002). evaluation. BMC Bioinformatics 19, 26 (2018). 42. K. A. Maher et al., Profiling of accessible chromatin regions across multiple plant spe- 19. H. Li, Minimap2: Pairwise alignment for nucleotide sequences. Bioinformatics 34, cies and cell types reveals common gene regulatory principles and new control mod- 3094–3100 (2018). ules. Plant Cell 30, 15–36 (2018). 20. T. D. Wu, C. K. Watanabe, GMAP: A genomic mapping and alignment program for 43. J. Schnable, E. Lyons, Plant paleopolyploidy (2015) https://doi.org/10.6084/m9. mRNA and EST sequences. Bioinformatics 21, 1859–1875 (2005). figshare.1538627.v1. 21. S. Marco-Sola, J. C. Moure, M. Moreto, A. Espinosa, Fast gap-affine pairwise align- 44. Y. Kawahara et al., Improvement of the Oryza sativa Nipponbare reference ment using the wavefront algorithm. Bioinformatics 37, 456–463 (2020). genome using next generation sequence and optical map data. Rice (N. Y.) 6, 4 22. B. Song et al., Complement genome annotation lift over using a weighted sequence (2013). alignment strategy. Front. Genet. 10, 1046 (2019). 45. O. Jaillon et al.; French-Italian Public Consortium for Grapevine Genome Characteri- 23. X. Gan et al., Multiple reference genomes and transcriptomes for Arabidopsis thali- zation, The grapevine genome sequence suggests ancestral hexaploidization in ana. Nature 477, 419–423 (2011). major angiosperm phyla. Nature 449, 463–467 (2007). 24. G. Marçais et al., MUMmer4: A fast and versatile genome alignment system. PLoS 46. C. Belser et al., Chromosome-scale assemblies of plant genomes using nanopore long Comput. Biol. 14, e1005944 (2018). reads and optical maps. Nat. Plants 4, 879–887 (2018). 25. H.-N. Lin, W.-L. Hsu, GSAlign: An efficient sequence alignment tool for intra-species 47. Y. Wu et al., A multiple genome alignment workflow shows the impact of repeat genomes. BMC Genomics 21, 182 (2020). masking and parameter tuning on alignment of functional regions in plants. 26. Y. Jiao et al., Improved maize reference genome with single-molecule technologies. bioRxiv [Preprint] https://doi.org/10.1101/2021.06.01.446647. Accessed 2 June Nature 546, 524–527 (2017). 2021. 8 of 9 j PNAS Song et al. https://doi.org/10.1073/pnas.2113075119 AnchorWave: Sensitive alignment of genomes with high sequence diversity, extensive structural polymorphism, and whole-genome duplication 48. J. C. Schnable, N. M. Springer, M. Freeling, Differentiation of the maize subgenomes 52. H. Li et al.; 1000 Genome Project Data Processing Subgroup, The sequence align- by genome dominance and both ancient and ongoing gene loss. Proc. Natl. Acad. Sci. ment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009). U.S.A. 108, 4069–4074 (2011). 53. B. Song et al., Data for "AnchorWave: Sensitive alignment of genomes with high 49. A. Coghlan, E. E. Eichler, S. G. Oliver, A. H. Paterson, L. Stein, Chromosome evolution sequence diversity, extensive structural polymorphism, and whole-genome in eukaryotes: A multi-kingdom perspective. Trends Genet. 21, 673–682 (2005). duplication." GitHub. https://github.com/baoxingsong/anchorwave. Deposited 50. T. Zhao, M. E. Schranz, Network-based microsynteny analysis identifies major differ- 15 October 2021. ences and genomic outliers in mammalian and angiosperm genomes. Proc. Natl. 54. B. Song et al., Data for AnchorWave: Sensitive alignment of genomes with high Acad. Sci. U.S.A. 116, 2165–2174 (2019). sequence diversity, extensive structural polymorphism, and whole-genome dupli- 51. L. Kistler et al., Multiproxy evidence highlights a complex evolutionary legacy of cation. GitHub. https://github.com/baoxingsong/genomeAlignment. Deposited maize in South America. Science 362, 1309–1313 (2018). 1 November 2021. EVOLUTION Song et al. PNAS j 9 of 9 AnchorWave: Sensitive alignment of genomes with high sequence diversity, https://doi.org/10.1073/pnas.2113075119 extensive structural polymorphism, and whole-genome duplication