Genomic Repeat Abundances Contain Phylogenetic Signal

1Understanding aspects of comparative evolution, including at its simplest relationships between taxa at varying levels of classification, is being revolutionized by the advent of next-generation sequencing (NGS) technologies. All recent methods in this area are based on multiplexing samples from diverse taxa, thereby maximizing the number of taxa that can be sequenced in one lane or one plate of an NGS run (e.g., Illumina). NGS approaches have enabled a quantum leap in the amount of data available while becoming increasingly cost-effective. Approaches include amplicon sequencing (sequencing of specific genes or regions of interest) using barcoded primers, full mitochondrial and plastid genome sequencing, and phylogenomics based on the full complement of protein-coding genes. Many recent approaches are based on reduced-representation libraries (i.e., reducing genomic complexity/increasing recovery of homologous regions across taxa); in this arena RAD-sequencing based on restriction-site associated DNA fragments scattered across the genome and hybridization methods of targeted capture, so-called “pull-down” approaches, are two of the most common methodologies.

2However, in such phylogenetic/phylogenomics studies, and indeed in broader studies of comparative evolution, the repetitive portion of the genome is often discarded without consideration of any potential use. Repetitive elements in genomes consist of both tandem repeats and interspersed mobile elements (e.g., DNA transposons and retrotransposons). In angiosperms (flowering plants), such repeats are diverse and numerous, contributing up to 70%–80% of nuclear genomic DNA (gDNA), thus making flowering plants an excellent group in which to study the dynamics of repetitive element evolution. Genome sizes vary 2400-fold in angiosperms alone; aside from cases involving wholegenome duplication, much of this variability can be explained by differing amounts of repetitive DNA.

3NGS of a small, random sample of the genome (0.5–5% genome proportion [GP], i.e., genome coverage as a percentage) results in data consisting mainly of repetitive sequences; genic regions will not be adequately covered in such a dataset, but repeats present in thousands of copies will be well represented. Previous analyses have shown that low-coverage sequencing of gDNA, followed by graph-based clustering of sequence reads, is sufficient to provide characterization of many hundreds or thousands of well-represented repeats; these studies also provide detailed insights into patterns of genome evolution. Low-coverage gDNA sequencing (i.e., “genome skimming”; ) and repeat clustering are now both cost-effective and easy to implement. The proportion of sequence reads representing a particular repetitive element cluster has also been shown to accurately reflect genomic abundance. Repetitive elements are scattered across the genome and provide much of the characteristic differences between chromosomes and chromosomal subregions, including those in which the majority of genes are embedded. Thus, relative abundance of well-represented repeats is reflective of broad-scale genome composition. Localization of repeats on chromosomes and use of repeats as markers for fluorescence in-situ hybridization (FISH) in some groups has shown that often the most-parsimonious explanation for these localizations and rearrangements reflect hypotheses of the species tree derived from other data types, usually DNA sequence data (e.g.,; ). In structure and chromosomal position, repeats in closely related species are nearly identical, whereas more distantly related species diverge in repeat structure and location as genetic similarity decreases.

4Here we test the usefulness of a novel phylogenetic methodology based on the abundance of different repetitive elements, estimated through bioinformatic analysis of NGS reads from a small proportion of the genome. Previously similar studies have found that genomic signatures present in the frequency of short sequence repeats can be used to reconstruct phylogenetic relationships, i.e., tetranucleotide frequencies in microbial genomes and 2- to 5-nt repeats in birds. Here we use different criteria and a clustering method in order to identify homologous repeat classes. This method can essentially be viewed as a hybrid between molecular systematics and morphometric cladistics, as abundances of repetitive DNAs are used as continuously varying characters for phylogenetic inference. We utilize in combination graph-based clustering estimation of repeats and the computational methodology of in particular, which allows for analysis of continuous characters without assignment (coding) of arbitrarily circumscribed characters, implemented in the software “tree analysis using new technology,” TNT. Such a combined approach has been utilized successfully with eigenshape-based geometric morphometrics and continuous character phylogenetics in TNT. We investigate this method in six diverse groups – five orders of angiosperms and one insect group, with differing genome sizes and amounts of repetitive DNA, and show a high (but not always identical) level of congruence with previously hypothesized species trees (i.e., current knowledge from gene trees and morphological circumscriptions) at a variety of taxonomic levels.

5Plant materials (accession numbers), DNA extraction and Illumina sequencing details (including NCBI Short Read Archive [SRA] accession numbers) can be found in and.

6Plant materials (including voucher specimen details), DNA extraction and 454 sequencing details (including SRA accession numbers) for this dataset can be found in. Orobanchaceae is the largest family of parasitic flowering plants. Four genera were included in our dataset, representing a variety of life history strategies: Lindenbergia, autotrophic, nonparasitic; Schwalbea, parasitic but still photosynthetic; four species of Orobanche, nonphotosynthetic, parasitic, including one tetraploid species (O. gracilis); and three species of Phelipanche, nonphotosynthetic, parasitic.

7Seeds of Vicia tetrasperma (VIC726), V. hirsuta (VIC728), V. sylvatica (VIC63), and V. ervilia (ERV52) were obtained from the seed bank of the Leibniz Institute of Plant Genetics and Crop Plant Research (IPK), Germany. Seeds of Lathyrus sativus and L. latifolius were purchased from Fratelli Ingegnoli S.p.A., Milano, Italy (cat.no. 455) and SEMO Smrzice, Czech Republic (acc.no. 1-0040-68867-01), respectively. Lathyrus vernus was collected from a wild population at Vidov, Czech Republic (GPS 48°55'17.401”N, 14°29'44.158″E). Pisum fulvum (accession ICARDA IG64207) was provided by Petr Smykal, Palacky University, Olomouc, Czech Republic. In all species, genomic DNA was extracted from isolated leaf nuclei and sequenced on the Illumina platform (paired-end 100 nt reads) at Elim Biopharmaceuticals, Hayward, USA (P. fulvum) or GATC Biotech, Konstanz, Germany (all other species). Illumina sequencing of P. sativum was described in. Voucher specimens are available for all material sequenced at IPMB, CZ. All read data are available at the SRA with the following accession numbers: V. hirsuta—ERR413114; V. ervilia—ERR413112; V. sylvatica—ERR413113; V. tetrasperma—ERR413111; Lathyrus sativus—ERR413118 & ERR413119; L. vernus—ERR413116 & ERR413117; L. latifolius—ERR413120; Pisum sativum—ERR063464; P. fulvum—ERR413083. Tribe Fabeae Rchb. is a group of five genera and ∼380 species, containing several important crop species including pea (Pisum sativum). In our analysis we included species from three genera, although this includes species proposed to be members of a further two new genera.

8DNA extractions were sourced from the Royal Botanic Gardens, Kew DNA Bank (http://apps.kew.org/dnabank.homepage.html). At the University of Liverpool, 454 sequencing was performed by the Centre for Genomic Research; reads were trimmed to 100 bp prior to clustering, and any reads of < 100 bp were discarded. All read data are available at the SRA with the following accession numbers: F. affinis—ERR571997; F. alfredae subsp. glaucoviridis—ERR571998; F. davidii—ERR571999; F. imperialis—ERR572000; F. koidzumiana—ERR572001; F. maximowiczii—ERR572002; F. pluriflora—ERR572003; F. sewerzowii—ERR572004; F. tubiformis—ERR572005; Lilium pyrenaicum—ERR572006.

9Fritillaria is a genus of bulb-bearing petaloid monocots with species possessing some of the largest recorded genome sizes. It comprises approximately 140 species (Rix 2001) and is closely related to the genus Lilium. In our analysis, we had nine representatives of the genus from each of the two main clades—the North American and the Eurasian clades.

10Illumina reads for the following species were downloaded from the SRA: Drosophila bipectinata—SRR345542; Drosophila suzukii—SRR1002946; Drosophila biarmipes—SRR345536; Drosophila ananassae—SRR491410; Drosophila melanogaster—SRR1005465; Drosophila sechellia—SRR869587; Drosophila simulans—SRR580369.

11Illumina reads for species from the Sonoran Desert clade of Asclepias were downloaded from the SRA: A. macrotis 149—SRX384308; A. albicans × subulata 282—SRX384307; A. cutleri 382—SRX384306; A. subulata 423—SRX384305; A. macrotis 150—SRX384304; A. albicans 422—SRX384303; A. subulata 411—SRX384302; A. masonii 154—SRX384301; A. leptopus 137—SRX384300; A. cutleri 421—SRX384299; A. coulteri 45—SRX384298; A. subaphylla 272—SRX384297; A. subaphylla 271—SRX384296; A. albicans 003—SRX384295; A. syriaca 4885—SRX040889.

12Further details of raw data, quality filtering and resultant read datasets for all examples are provided in Online Appendix 1 (http://dx.doi.org/10.5061/dryad.vn0gc).

13In order to calculate the number of reads for each comparative clustering an accurate genome size should ideally be available for each species. For most datasets genome sizes were available from the Plant DNA C-Values database (http://data.kew.org/cvalues/) or were estimated using flow cytometry. For Drosophila, genome sizes were taken from the Animal Genome Size database (http://www.genomesize.com). For Asclepias genome sizes were assumed to be equal (423 Mb, as for A. syriaca; ), as data were unavailable for each species. Without accounting for genome size (e.g., taking the same number of reads), the abundance of each cluster is more likely to reflect genome size rather than the proportion of that repeat in the genome. Genome size ranges (1C) in each dataset are as follows: Nicotiana (1.51–5.32 Gb); Orobanchaceae (0.454.32 Gb); Fabeae (3.05–9.98 Gb); Fritillaria (30.1–75.7 Gb); Asclepias (n/a); Drosophila (0.16–0.20 Gb).

14Graph-based clustering of NGS reads was performed as described in using the latest Galaxy-based web server implementation of the pipeline, RepeatExplorer. In brief, all sequence reads (sequence data) are subjected to pair-wise (BLAST) comparison, and similarities are represented by a graph structure in which nodes represent sequence reads and overlapping reads are connected by edges. Edge weights represent the amount of sequence similarity (similarity scores). Clusters of nodes more frequently connected to one another than to outside nodes in the graph represent families of genomic repeats or their parts. Families would be likely to include sequences of the same length (or portions thereof) in which sequence variation is low, 90% similarity over at least 55% of their length.

15Combined datasets of reads (sequence data) were compiled as follows: (1) 5% genome proportion (GP) each of four diploid species of Nicotiana L. (N. sylvestris Spreng., N. tomentosiformis Goodsp., N. attenuata Torr., N. obtusifolia M.Martens & Galeotti); (2) 5% GP each of the four diploid species of Nicotiana in (1) and two species of allopolyploid section Repandae (N. repanda Sims and N. nudicaulis S.Watson); (3) 5% GP each of the four diploid species of Nicotiana in (1) and two types of N. tabacum L. (N. tabacum SR1A and N. tabacum TR1A synthetic); (4) 2.08% GP each of nine species of Orobanchaceae (Lindenbergia philippensis (Cham. and Schltd.) Benth., Schwalbea americana L., Phelipanche ramosa (L.) Pomel, Phelipanche purpurea (Jacq.) Soják, Phelipanche lavandulacea Pomel, Orobanche pancicii Beck, Orobanche gracilis Beck, Orobanche cumana Wallr., Orobanche crenata Forssk.); (5) 1% GP each of nine species of tribe Fabeae (Vicia sylvatica L., Vicia ervilia Willd., Vicia hirsuta (L.) Gray, Vicia tetrasperma (L.) Schreb., Pisum sativum L., Pisum fulvum Sibth. & Sm., Lathyrus sativus L., Lathyrus vernus (L.) Bernh., Lathyrus latifolius L.); (6) 0.01% GP each of nine species of Fritillaria L. and one of Lilium L. (Fritillaria affinis (Schult. & Schult.f.) Sealy, F. alfredae subsp. glaucoviridis (Turrill) Rix, F. davidii Franch., F. imperialis L., F. koidzumiana Ohwi, F. maximowiczii Freyn, F. pluriflora Torr. ex Benth., F. sewerzowii Regel, F. tubiformis Gren. & Godr., Lilium pyrenaicum Gouan); (7) 2% GP each of 15 Asclepias L. (Asclepias syriaca L. 4885, A. albicans S. Watson 003, A. albicans 422, A. coulteri A. Gray 45, A. cutleri Woodson 382, A. cutleri 421, A. leptopus I. M. Johnst. 137, A. macrotis Torr. 149, A. macrotis 150, A. masonii Woodson 154, A. subaphylla Woodson 271, A. subaphylla 272, A. subulata Decne. 411, A. subulata 423, A. albicans × subulata 282; (8) 5% GP each of 7 Drosophila species (D. ananassae, D. bipectinata, D. suzukii, D. biarmipes, D. melanogaster, D. sechellia, D. simulans). Different GP values were used across datasets due to genome size differences and the amount of sequencing data available or that could be clustered with the available computing power.

16Separate comparative analyses (i.e., simultaneous clustering of reads from all species in the dataset) were run for each dataset on RepeatExplorer, using default settings (i.e., similarity threshold of 90% over 55% of the read length). Reads were prefixed with codes specific to the taxon in question, enabling comparative analysis of repetitive element abundances in different taxa. Comparative counts of the number of reads in each cluster (which is proportional to their genomic abundance) were used for phylogenetic analyses. Plastid and mitochondrial reads were either filtered out prior to clustering (using BLAST and custom scripts) or were identified after clustering (BLAST to most closely related plastome currently available) and plastid clusters removed prior to phylogenetic inference.

17The extent to which clusters represent homologous entities is dictated by the similarity parameters specified. We use the default settings of RepeatExplorer, with a threshold similarity of 90% over 55% of the read length to be exceeded in order for a hit to be recorded. Clusters are then produced using a graph-based algorithm and a principle of maximum modularity, which results in clusters where most reads have a high similarity to one another within clusters and a low similarity between clusters (see for further details on the clustering process). Different repeats will form different clusters in the output of RepeatExplorer and are treated here as separate evolutionary entities (characters). The abundance of a repeat in a species, its genome proportion, depends on repeat copy number and genome size. Tandem repeats have variable monomer sizes up to 180 bp; those with monomer sizes shorter than the read length (typically 100 bp) will form a spherical graph, and those with monomer sizes greater than the read length will form a ring graph structure.

18Plant genomes in particular contain a large abundance of LTR retroelements (LTR-REs), which are typically several kb, up to 5kb. These repetitive elements are complex, often dispersed across the genome, and there may be a spectrum of related (or degraded products) of similar LTR retroelements. Based on the RepeatExplorer threshold and graphical algorithm, LTR-REs are often split into different clusters, as parts of these elements are less conserved (e.g., around the LTR) than others (e.g., the protein-coding domains). Sequence divergence within the protein-coding domains is insufficient for phylogenetic analysis, although the number of elements is variable and putatively indicative of evolutionary history. Although LTR-REs may be split over several clusters, they will be split in the same way for every species included in the same clustering run, thereby preserving phylogenetic signal, and each piece of LTR-RE would be expected to contain a uniform phylogenetic pattern.

19High-copy DNA sequences were assembled directly from short read data using the program MIRA (http://www.chevreux.org/projects_mira.html). The general settings used in the manifest file are provided in Online Appendix 1. The following assemblies were performed: (2) Nicotiana—large subunit rDNA and whole plastomes were assembled by mapping to Nicotiana tabacum sequences as a reference, using raw Illumina reads; (2) Fritillaria—whole plastome sequences were assembled directly from plastid 454 reads only (filtered using a custom perl script and BLAST), using the Lilium longiflorum plastid genome as a reference; (3) Orobanchaceae—whole plastomes for O. cumana, O. pancicii, O. crenata were assembled from raw 454 reads using the O. gracilis plastome as a reference and P. lavandulacea was assembled using P. ramosa as a reference.

20Data matrices consisting of the 1000 most abundant clusters, each representing a repetitive element family, were converted to legal TNT format (modified Hennig86). All abundances were transformed by a constant factor dependent upon the largest cluster abundance in the matrix. Each abundance was divided by this factor ( abundance/65) in order to make all numbers in the matrices ≤ 65, the maximal value for continuous character implementation in TNT tree searches. This factorial transformation does not affect the normal distribution of abundance for each cluster and is only necessary for efficient implementation in the TNT program, as described below.

21Trees were inferred using maximum parsimony (MP), utilizing the implementation of Farris' algorithm for the down-pass and Goloboff's algorithm for the up-pass, as described in. In such an approach, continuous characters are not arbitrarily recoded but are simply used as additive characters (i.e., count changes can be of noninteger differences). Implicit enumeration (branch and bound) tree searches were used for datasets in this study owing to the small number of taxa in each dataset. Resampling was performed using 100 000 replicates and symmetrical resampling, a modification of the standard bootstrap. Sequence trees were inferred using the same method for comparison, with gaps coded as missing data. The same phylogenetic reconstruction methodology was employed to enable direct comparison to the repeat trees.

22The following datasets were used: (1) full plastomes and 18S-5.8S-26S rDNA for Nicotiana diploids and section Repandae (assembled); (2) full plastomes for Fritillaria (assembled); (3) combined matrix of 17 mitochondrial and nuclear genes for Drosophila (28S, adh, amy, amr, cdc6, COI, COII, ddc, esc, gpd, h2s, hb, ITS, ND1, ND4, nup and ptc – see for GenBank accession numbers); (4) whole plastomes and complete 26S to 18S rDNA cistron for Asclepias—alignments taken from; (5) whole plastomes for Orobanchaceae (assembled); and (6) nuclear ITS rDNA and plastid trnL from for tribe Fabeae.

23Maximum likelihood (ML) trees were computed using gene frequency/continuous character implementation in Contml, part of the Phylip package. This method assumes that each character evolves independently and only in accordance with random genetic drift, using a Brownian motion model of likelihood. Matrices were transformed such that cluster abundances represented allele frequencies (0–1) by dividing all clusters by the largest cluster size; this is required for resampling prior to ML tree computation. Resampling from the matrix with replacement (bootstrapping) was first carried out using Seqboot for 1000 replicate datasets. ML analyses were then performed on all 1000 datasets using Contml, and bootstrap percentages mapped onto the strict consensus tree for each dataset computed using Consense.

24All trees were viewed in FigTree (http://tree.bio.ed.ac.uk/software/figtree/) and further edited in iDraw (Indeeo, Inc.). All ML trees are shown in Online Appendix 2. Reticulation in the Nicotiana tabacum dataset was explored using SplitsTree4, using 10 000 bootstrap trees from the MP analysis as input for filtered supernetwork analysis (filtering performed at 10% of all input trees). Nicotiana tabacum is a relatively recently formed allotetraploid; its two parents have been determined to be N. tomentosiformis and N. sylvestris.

25To test performance of the method several parameters were analyzed with the smallest clustering dataset (four diploid species of Nicotiana) and in TNT as above, but with modifications described below. In each case, the resultant tree was compared with the expected tree topology (Fig. 1a) and the symmetric bootstrap percentage recorded.

26To find the minimal GP necessary to resolve relationships, several sequence datasets were produced at 11 levels of GP from 0.005% to 5.120% (doubling of GP at each step). Three replicate clustering runs were computed for each GP. The mean support and standard error were calculated and used, in addition to tree support and topology, to observe how reproducibility varies with GP.

27To evaluate the number of characters (clusters) sufficient to resolve the tree, trees were built with different numbers of clusters, varying from 5 to 1000 for 25 datasets, and the tree inferred each time, comparing the topology and support percentages.

28Variance in phylogenetic signal across the cluster abundance matrix was tested by partitioning it into sets of 150 cluster abundances (this number chosen from the cluster number analysis above). Trees were then inferred from each set, and the resulting resolution and symmetric bootstrap support of the unrooted tree were then estimated. Three GPs were tested (2.00%, 0.32%, and 0.07%) in order to show how different partitions respond at different GPs (chosen to represent difference of three orders of magnitude).

29To evaluate the effect of sampling sequence data on tree building, trees were inferred from clustering of three random samplings of read data. The mean and its standard error were calculated for abundance of each cluster. A phylogenetic analysis based on the range of the mean ±1 standard error of the mean was conducted, as ranges may more accurately reflect the phylogenetic signal in continuous characters, thereby reducing the artefact of two taxa appearing as distinct when they are not. Clusters (repetitive element abundances) that have overlapping normal distributions result in a step count of 0, i.e., no change. Range analysis was tested with a GP of 0.32%.

30The relative informativeness of different repeat types was analyzed by creating subsets of the original matrix based on different repeat annotations. Annotations were assigned to the following categories based on BLAST hits to repbase, custom-protein domain database in the RepeatExplorer pipeline and graph structure: DNA transposon, Ty1/Copia LTR retrotransposon, Ty3/Gypsy LTR retrotransposon, rDNA, satellite, and other (e.g., nonLTR retrotransposon) including unclassified repeats. Matrices were created based on each repeat type, for each example taxon dataset, and trees were inferred as above. The mean bootstrap (as a proxy for tree resolution) was computed for each analysis (Fig. 5).

31The unrooted tree presented in Figure 1a contains four diploid species of Nicotiana and mirrors gene trees based on other nuclear DNA regions, including rDNA sequences reconstructed from the NGS data (Fig. 1). There is a different relationship for these four using the plastome tree constructed from these NGS data (Fig. 1a), in line with plastid gene trees published previously.

32In a further analysis of Nicotiana, tree resolution and topology were investigated using data from diploid species and allotetraploid species from two sections of Nicotiana (Nicotiana sections Repandae and Nicotiana). Nicotiana section Repandae is a group of four allopolyploids derived from a single allopolyploid formation event approximately 5 Ma. Genomes of sect. Repandae have experienced extensive genome turnover subsequent to their formation, and the genomes retain more similarity to the extant relative of their maternal progenitor, N. sylvestris, rather than the extant relative of their paternal progenitor, N. obtusifolia. This striking bias is supported in our analysis here (Fig. 1b), in which N. repanda and N. nudicaulis together (100 bpP) are strongly supported as sister (89 bp) to N. sylvestris. Previous analyses of repetitive DNA and genomic in situ hybridization (GISH) have shown that the genomes of sect. Repandae have diverged extensively since their formation, despite being each other's closest relatives, through loss of middle and lower-abundance repetitive elements. This is particularly evident in N. repanda (Fig. 1b).

33To contrast this example, we examined how the method performs with a different allopolyploid section of much more recent origin—Nicotiana sect. Nicotiana, which contains the familiar allotetraploid Nicotiana tabacum, the most common tobacco species in commerce. Nicotiana tabacum is estimated to have originated 200 000 years ago or less, and its formation involved entities closely related to extant N. tomentosiformis and N. sylvestris, its paternal and maternal progenitors, respectively. GISH is able to distinguish the progenitor genomes, the T-genome from N. tomentosiformis and the S-genome from N. sylvestris. Analyses of repetitive DNA show the genome of N. tabacum has preferentially lost paternal repeats and is much more similar to N. sylvestris, although its nrITS and IGS sequences of ribosomal DNA are identical to its paternal progenitor, N. tomentosiformis. In repeat phylogenetic analyses the tree with all four diploid species and N. tabacum shows that N. tabacum is more closely related to N. sylvestris than to N. tomentosiformis (Fig. 2a), reflecting that abundances of repetitive DNA in N. tabacum are in general more similar to those in the maternal parent, N. sylvestris (Fig. 2a). This result was previously found based on analyses of the GP in different clusters of repetitive DNA, which showed a preferential loss of paternal (i.e., N. tomentosiformis) repeats in N. tabacum. Nevertheless, the supernetwork (Fig. 2b) illustrates the presence of splits that group N. tabacum with N. tomentosiformis in addition to placing it with N. sylvestris, indicating that some repeats inherited from the paternal progenitor are still present. Analysis of 10 000 bootstrap trees reveals that N. tabacum groups with N. tomentosiformis in 17% of the trees; in the remaining 83% of trees it is sister to N. sylvestris. Additionally, the relatively long branch length separating the N. tabacum samples from that of N. sylvestris highlights the retention of characters conflicting with its position as sister to N. sylvestris (i.e., presence of paternal-type repeats).

34Lilium is the designated outgroup, following, and our analysis generally places species into their expected clades (Fig. 3a). Other than the placement of F. maximowiczii, the repeat tree is in agreement with trees based on plastid/plastome data (Fig. 3a;; ). The ML tree is partially resolved (Online Appendix 2). Owing to the huge genome sizes in this genus our analysis was based on a very low GP of 0.01%, and this may have had an impact on the tree building method—this level of GP should be sufficient, although it is on the cusp of being too low to be representative of repeat diversity and composition (Fig. 4a; see discussion below). However, it still reproduces species relationships in a similar manner to previous results, indicating that this number of reads still contains phylogenetic signal despite their low GP.

35Our analyses for the diverse fly genus Drosophila focused on seven species from the melanogaster subgroup. Drosophila simulans and D. sechellia are strongly supported as sister species, to which D. melanogaster is then sister (Fig. 3b). Drosophila suzukii and D. biarmipes form a clade, which is sister to the D. melanogaster clade. Drosophila ananassae is sister to the rest, with rooting on D. bipectinata. These results (Fig. 3b) mirror those found in many recent phylogenetic studies based on large amounts of sequence data including mtDNA and nuclear markers. The ML analysis mirrors these results with reasonably high levels of support (Online Appendix 2).

36To test our method on a difficult phylogenetic problem we investigated the Sonoran Desert clade (SDC) of Asclepias and present for comparison the rDNA and plastome results of, with A. syriaca as the outgroup. In the repeat tree, A. macrotis, A. coulteri and A. leptopus are supported as separate from the core SDC, which includes A. albicans, A. subulata, A. subaphylla, A. masonii and A. cutleri; in plastid and nrDNA analyses a strongly supported core SDC excludes A. cutleri (Fig. 3c). Otherwise the results are different from both rDNA and plastomes, which are in turn different from one another and that from mtDNA. Note that species are not monophyletic (as with the plastome tree), and a putative homoploid hybrid was included (A. albicans, A. subulata). Thus our method provides yet another novel hypothesis of relationships between species for this difficult phylogenetic problem, although these relationships are weakly supported. It should be noted however that our results may be compromised by a lack of genome size data with which to calibrate the number of input reads. If there are large GS differences between species, then the GP analyzed will be significantly different between species, which may influence the topology of the trees produced.

37The tree is rooted with Lindenbergia philippensis based on the analyses of and. Orobanche and Phelipanche are each monophyletic (Fig. 3d), and Schwalbea americana is then sister to them, generally exhibiting similar to greater levels of divergence than analyses based on DNA sequence data (Fig. 3d). Such high levels of divergence are only apparent in this dataset.

38Internal relationships in Orobanche and Phelipanche are entirely congruent with previous analyses based on nrITS and plastid rps2 sequence data, in all cases with similar or better support in our analysis (Fig. 3d;;; ). The position of O. gracilis/O. cumana is switched, however, relative to the tree generated from full plastome sequences (Fig. 3d). The high support for the position of O. gracilis and the fact that there is no evidence of reticulation (result not shown) are congruent with an autopolyploid origin for this taxon. This could, however, reflect processes of diploidization and genome downsizing in an older allotetraploid, as shown in Nicotiana allopolyploids (Example 1); our result could also be affected by low taxon sampling within Orobanche (i.e., only one parent present). This might explain the difference between our placement and the one based on plastomes.

39Following the tree of, our tree, with section Ervilia, now proposed to be genus Ervilia (Vicia hirsuta, V. sylvatica, and V. ervilia) as the outgroup, shows a closer relationship between V. ervilia and V. sylvatica. Vicia tetrasperma is the sole representative of section Ervum (proposed genus Ervum), which is sister to both Pisum and Lathyrus (Fig. 3e). The type section of Vicia is not included in this analysis. Pisum and Lathyrus are sister taxa, and L. sativus is sister to L. vernus and L. latifolius, a grouping that is incongruent with the weakly supported results of the DNA sequence data, in which L. sativus and L. latifolius are more closely related (Fig. 3e; ). The MP analysis generally has lower support than in other examples, but this broadly mirrors levels of support found in results based on nuclear nrITS and plastid markers. The ML analysis recovers some of the main groupings; however, many nodes collapse, and overall the ML tree is unresolved (Online Appendix 2). An alternative explanation for the lower levels of support found in this example is that perhaps in this case we are looking at a group of more distantly related taxa (comprising several relatively distant related genera), which may be approaching the limits of phylogenetic utility for repetitive DNA. The other examples may well be more closely related, but this is an area in which more investigation is needed to clarify the issues. Compared to the DNA sequence tree, the repeat tree has shorter internal branches and longer external branches, indicating that clustering information may be limited. However deep coalescence and/or extensive reticulation combined with many unsampled taxa may be the underlying problem in this complicated group of legumes. This may explain the poor resolution in as well as difficulties in our analyses.

40Based on resampling the diploid Nicotiana data, several aspects of method performance were evaluated. At low GPs, 0.005–0.040%, trees lack resolution, and groupings are often inconsistent with those inferred from sequence data (Fig. 4a). Additionally, variance is greater at lower levels of GP, making the method less reliable. Above a GP of 0.1%, clustering and phylogenetic inference appear to be consistent and robust in comparison to trees derived from DNA sequence data. With these GPs, the number of clusters necessary to resolve the tree with high support is approximately 150 (Fig. 4b). With lower numbers of clusters (e.g., 5–45), the tree is either topologically inconsistent with trees inferred from sequence data or simply unresolved.

41The matrices were then explored using partitions of 150 clusters to test how the phylogenetic signal varied across a large range of cluster abundances, with CL1 being the most abundant (Fig. 4c). There is no or little variance in phylogenetic signal across partitions when a suitably high GP has been used (e.g., 2%; Fig. 4c). However, when lower GPs are used, the signal degrades rapidly and randomly, particularly with low-abundance elements below cluster number 2000 (Fig. 4c). At a GP of 0.07% the trees produced become highly stochastic, either unresolved or showing inconsistent relationships with low support. At a lower GP the signal degrades more quickly as the data essentially become more quickly “coded” into presence/absence characters—the phylogenetic signal in the actual abundance of repetitive elements is entirely lost.

42Range analysis showed remarkable similarity between trees produced from three independent samplings of 0.32% (as shown in the GP analysis above). The tree produced from the range of the mean cluster abundance ±1SE of the mean had a topology identical to the three samples as expected, and there are only minor branch length differences between all trees. Thus, there seems to be no significant advantage to undertaking such an analysis, although it may be advisable to do so given the ease of effort and proposed advantage of avoiding spurious groupings where species are not actually statistically different (i.e., they have overlapping normal distributions of cluster abundance).

43Repeat types differed in their relative informativeness between datasets, but not in a consistent manner. Overall DNA transposons appeared to have less consistent phylogenetic signal than other types of repeats (Fig. 5). The relative informativeness of retrotransposons seems to be more taxon-specific, e.g., in some datasets Ty1/Copia were more informative than Ty3/Gypsy, whereas in others the opposite is true. The importance of including unclassified repeats is also highlighted by the informativeness of this category in all datasets examined.

44Using one insect and five angiosperm examples that vary in genome size by ∼400-fold (from in Drosophila to 75.7 Gb in Fritillaria), we have shown that using relative abundance of repetitive elements as continuous characters successfully resolves species relationships in a manner similar to that obtained by using DNA sequences from plastid and nuclear ribosomal regions. Using low coverage sequencing of genomic DNA and ≥ 150 cluster abundances, we show that the repetitive DNA-based phylogeny reconstruction method is consistent in resolving expected relationships (i.e., those produced with other, more standard, methods, and data). The method can be seen as an additional source of phylogenetic information from the repetitive, noncoding portion of the genome, which will be a useful comparison to results based on DNA sequences. Gene trees represent the ancestry of particular sequences, and because allele histories may have differences from the species trees plus have coalescent times that differ from each other and species divergence times they often present conflicting topologies. Thus, there has been a recent focus on multiple sequence datasets, which aim to give a genome-wide assessment of species divergence. Here, our method provides data widespread across the genome, and each repeat abundance (cluster) is a marker; a matrix of such abundances likely represents many independently evolving characters. Furthermore at an appropriate level of GP there is little variance in phylogenetic signal across the dataset, showing remarkably consistency, although this requires further testing with additional datasets that include a larger number of taxa.

45Our ML results on the whole agree with our maximum parsimony results, although for many groups the tree is only partially resolved, e.g., Fritillaria and tribe Fabeae (Online Appendix 2). We are unsure as to why ML does not perform as well as MP here, but the large state space of the character coding might provide insufficient information to accurately inform site likelihoods. We will develop further likelihood and Bayesian models for inferring trees from repeat abundances, but these approaches are limited by our understanding of repeat evolution, which is in its infancy. It is possible that horizontal transfer occurs for several repeats, but we believe the overall impact of this on our results is low, due to the large proportion of the genome covered by these analyses. Additional work is needed to further model the evolution of repeat populations, but we believe these results and others (e.g.,,;012) provide evidence that repeats evolve primarily in accordance with random genetic drift; they therefore contain useful phylogenetic signal.

46For reliable estimates of species relationships, we recommend using for clustering and subsequent phylogenetic analysis based on at least 150 cluster abundances. If a lower GP is used we suggest using 1000 cluster abundances in the analysis for reliable detection of the phylogenetic signal present. With larger (and more repetitive genomes) it is possible to resolve relationships with lower levels of GP as read depth will likely still be sufficient at GP < 0.1%, as we observed here for Fritillaria.

47There is an extensive literature on the evolution of repeats, see e.g., recent reviews by,,, and. All genomes contain tandem repeats and transposable elements, predominantly retrotransposons and derivative repeats—in plants these are typically long-terminal repeat (LTR) Ty1/Copia and Ty3/Gypsy elements. Copy number of these elements is highly variable and can change rapidly, contributing a large effect on genome size and architecture. Genome content of these repeats, as a whole is the result of mechanisms of repeat expansion (e.g., retrotransposition, repeat recombination) and contraction (e.g., recombination-based deletion). In this study we show that the variable abundance of different repetitive elements contains phylogenetic signal, i.e., one reflecting the evolutionary history of these species, which would be expected, given the premise that repeats are an inherent structural feature of the genome and in fact underlie much of the evolution of large, complex eukaryotic genomes.

48Previously, it has been shown that the frequency of short sequence repeats provide genomic signatures that can be used to reconstruct phylogenetic relationships (e.g.,; ). Our approach builds on this insight to show that a genomic signature containing phylogenetic information extends to many larger repeat classes (Fig. 5), based on sequence similarity and graphical clustering. Potentially therefore it is likely to be robust to particular features of individual genomes being analyzed, which may or may not be rich in certain categories of repeat. Researchers using genome skimming approaches to assemble high-copy DNA features (plastomes, rDNA cistron) already have these repeat data available (e.g., ), and it provides another data source of noncoding nuclear DNA from which to infer phylogenetic hypotheses.

49This method is genome-wide, and there is no need to attempt to distinguish paralogs, as essentially each cluster represents a homologous family of repeats (as determined by their highly conserved sequence, which is how they are clustered in the first place), the genomic abundance of which is used as phylogenetic data. In contrast, using low-copy nuclear markers requires paralogs to be distinguished accurately from one another. Low-copy nuclear genes are becoming increasingly popular due to often-higher variability when compared with plastid or mitochondrial markers (but not always, ), particularly for recent radiations and population studies. High variability can be seen as a consequence of having long coalescent times (compared to plastid or mitochondrial markers), but this in fact confounds their use, as when two or more young taxa share ancient alleles only due to the absence of fixation in ancestral populations rather than unique descent. This means that many such markers are needed to provide evidence for which ones are providing spurious or conflicting results. Our repeat method, however, shows particular promise for these cases, as repeats evolve rapidly, and there is neither the added complication of comparing markers with different coalescent times nor effect of longer coalescent times and incomplete lineage sorting.

50This method proved useful for inferring parents of some allopolyploid taxa in which extensive diploidization has occurred (i.e., replacement of repeats typical of one parent with those of the other) and repeat abundance is much more similar to one of the parents (here observed in groups of Nicotiana allopolyploids). In other groups, in which homoploid and polyploid hybridization occurs, further conflict will occur in construction of strictly bifurcating trees, but this can be analyzed more thoroughly using networks or pruning of putative parental or hybrid taxa. Using each cluster abundance (a single character) to make a separate estimate and then building a network based on these input trees, we might be able to extract in the repeat data further evidence of reticulation. Nonetheless, in our analyses repeat abundances have proven useful for determining one of the putative parents of allopolyploid taxa, but perhaps with some refinement in methods we could be expected to demonstrate clear evidence for both parents, especially in recently synthesized hybrids (as in the N. tabacum analyses when we looked at individual trees from the bootstrap replicates). Additionally, this method has not yet been tested in a species or species complex of wide geographical range, which would be a useful further test for how this method performs at the intraspecific level. If homogenization of repeats occurs as a result of gene flow, which holds back formation of sequence variants, we predict that divergence occurs quickly once gene flow ceases, as might occur for isolated populations of a widespread species.

51There is an additional caveat that should be mentioned; in order for the method to work reliably it is best that genome size is first estimated (via flow cytometry). If genome size is not used to standardize repeat abundances, then the resultant trees may reflect genome size differences more than shared evolutionary history, as repeat estimates reflect sampling bias rather than true abundance. It should be noted that the Asclepias example did not include genome size standardization as these data were not available. Additionally, we simulated lack of genome size information in the Fritillaria example, using the same number of reads (40 000) for each sample, despite >2-fold genome size differences (30–75 Gb), but results from this analysis closely mirror those presented in the results section, again with high levels of support (Online Appendix 3). Thus when genome sizes are unknown, the method still has the potential to produce a reliable phylogenetic result. It should therefore be possible to include data from various genomic sources, including available genome sequence data, although we do not know if this would vary relative to sequencing effort/coverage used in the clustering.

52The proposed method may also prove useful for sequencing DNA from herbarium specimens, where genomic DNA is often degraded to a greater extent. Highly degraded DNA will be expected to contain intact copies of high-copy regions (i.e., particularly shorter repeats) even when other regions are largely eliminated. Herbarium specimens provide an invaluable source of plant material, often collected from remote regions and for taxa that may have since become rare or even extinct. Utilizing this resource will be a continued focus of research in plant systematics, and bridging the gap from Sanger sequencing of short markers to NGS of genome-wide markers is one current focus of research. Here, we provide one possible solution to this problem, as short NGS reads of low-coverage gDNA give us an invaluable insight into repetitive DNA proportions, a method that is advantageous because repeats are present in high-copy number and distributed across the genome. Unlocking and mining data from museum collections will understandably be a future focus of systematic studies (e.g., ).

53The clustering pipeline utilized here has been shown to be effective in characterizing repetitive elements across various groups of eukaryotes including, for instance, bats. Phylogenetic trees based on repeat abundances estimated with RepeatExplorer could be produced in other groups of animals, including mammals, and fungi. This makes the methodology particularly useful for those carrying out genome evolutionary studies in various types of organisms.

54We conclude that this methodology is quick and easy to implement from the initial stage of DNA extraction and Illumina sequencing through to clustering of reads and tree building, utilizing the RepeatExplorer pipeline and TNT program (both freely available). We predict that as the cost of NGS continues to decrease in coming years (currently the major cost is in the library preparation), the overall cost of using this method will also decrease.

55This method has proven successful in resolving species relationships as previously hypothesized by analyses of DNA sequence data (plastid and nuclear trees) and morphological circumscription in five diverse groups of angiosperms and one insect. There were only a few instances where the results are incongruent with trees derived from DNA sequences, for which evidence about cause should be sought; there are good reasons why in some cases plastid DNA and rDNA could be misleading, so these discrepancies should not be ignored or be thought of as a fault of using repeats as phylogenetic characters. This method does provide an important extension of molecular systematics methods and should be useful for comparative phylogenomics. At the same time as providing data for robust phylogenetic reconstruction in diploid species, this method provides abundant information for understanding genome evolution in the context of repetitive DNA. Indeed, this has already been done for a number of the datasets/partial datasets used in this study (e.g.,,; ).

56Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.vn0gc.