Phylogeny of Hawaiian Melicope (Rutaceae): RAD-seq Resolves Species Relationships and Reveals Ancient Introgression

Abstract

Hawaiian Melicope are one of the major adaptive radiations of the Hawaiian Islands comprising 54 endemic species. The lineage is monophyletic with an estimated crown age predating the rise of the current high islands. Phylogenetic inference based on Sanger sequencing has not been sufficient to resolve species or deeper level relationships. Here, we apply restriction site-associated DNA sequencing (RAD-seq) to the lineage to infer phylogenetic relationships. We employ Quartet Sampling to assess information content and statistical support, and to quantify discordance as well as partitioned ABBA-BABA tests to uncover evidence of introgression. Our new results drastically improved resolution of relationships within Hawaiian Melicope. The lineage is divided into five fully supported main clades, two of which correspond to morphologically circumscribed infrageneric groups. We provide evidence for both ancestral and current hybridization events. We confirm the necessity for a taxonomic revision of the Melicope section Pelea, as well as a re-evaluation of several species complexes by combining genomic and morphological data.

Introduction

Oceanic islands have long been a focal point of evolutionary studies, as they represent a microcosm for examining the process of speciation. This microcosm is shaped by a combination of factors: (1) islands are geographically small and discrete units, sometimes far removed from continental landmasses; (2) colonizations or secondary arrivals are relatively infrequent, and thus, gene flow between the source areas and island systems is restricted; and (3) islands often have dynamic geological histories that give rise to extensively varying landscapes with numerous ecological niches (Emerson, 2002; Price and Wagner, 2018). These factors can often lead to high levels of endemism, which is often the result of adaptive radiation of a limited number of colonizers (Price and Wagner, 2004; Losos and Ricklefs, 2009; Keeley and Funk, 2011). Synthesizing the unique aspects of island evolution and extrapolating results to larger scales may allow us to better uncover general patterns and processes in evolution. Such phenomena include identifying factors affecting successful colonization and adaptive radiation (Carlquist, 1967; Carlquist, 1974; Paetzold et al., 2018), morphological or ecological shifts (e.g., “insular woodiness”; Carlquist, 1974; Lens et al., 2013), the spatiotemporal origins of lineages (Appelhans et al., 2018a), reconstructing colonization events (Harbaugh et al., 2009), and studying co-evolution (Roderick, 1997). These insights may result in further questions regarding taxonomy, species richness, medicinal or technical applications, and conservation (e.g., Francisco-Ortega et al., 2000).

Adaptive radiations on islands are of special interest for connecting changes in morphology and ecology through time (Givnish, 1998) but require well-resolved phylogenies to do so. In the Hawaiian Islands, phylogenetic studies based on morphology and taxonomy have sometimes overestimated the number of colonization events, because high levels of morphological diversity led researchers to overestimate lineage diversity and the number of colonization events (Price and Wagner, 2018). In contrast, molecular phylogenetic studies have revealed that many enigmatic Hawaiian plant radiations colonized the islands only once followed by adaptive radiation: the Hawaiian lobeliads (Campanulaceae; Givnish et al., 2009), Psychotria (Rubiaceae; Nepokroeff et al., 2003), Silene (Caryophyllaceae; Eggens et al., 2007), Touchardia/Urera (Urticaceae; Wu et al., 2013), and Melicope (Harbaugh et al., 2009; Appelhans et al., 2014a). Polyploidization and hybridization events were also discovered to predate colonization and radiation in several island lineages, including the Hawaiian silverswords (Asteraceae; Baldwin and Sanderson, 1998; Barrier et al., 1999) and mints (Lamiaceae; Roy et al., 2015) along with the Pan- Pacific sandalwoods (Santalaceae; Harbaugh, 2008), suggesting evolutionary success in young hybrid or polyploid colonists (Carr, 1998; Paetzold et al., 2018).

Time-scaled phylogenies have revealed that most Hawaiian radiations are ≤5 Myr old, which corresponds to the age of the oldest current main islands, Kauaʻi and Niʻihau. This suggests a bottleneck for dispersal from older (and now largely submerged) leeward islands to the current main islands. However, there are several known exceptions of lineages older than 5 Myr, including Drosophila, damselflies, lobeliads, Zanthoxylum (Rutaceae), as well as Melicope (Price and Clague, 2002; Keeley and Funk, 2011; Appelhans et al., 2018a; Appelhans et al., 2018b). Most phylogenetic studies of Hawaiian flora, however, have relied on few sequenced loci and have thus lacked sufficient power to resolve recent rapid radiations where hybridization, incomplete lineage sorting (ILS), and polyploidy may be common. Newer genomic tools are likely to provide more accurate estimates that may transform our understanding of island radiations.

The genus Melicope comprises about 235 species of shrubs and trees distributed throughout SE Asia and Australasia, extending to the Mascarene Islands and Madagascar in the West and most of the Pacific Archipelagos in the East (Hartley, 2001). There are 54 species of Melicope endemic to the Hawaiian Islands (Appelhans et al., 2017; Wood et al., 2017), 41 of which are single island endemics (Stone et al., 1999). Hawaiian Melicope were initially placed in the genus Pelea together with species from the Marquesas Islands (Stone, 1969; Stone et al., 1999) but later incorporated into Melicope, forming the majority of the section Pelea (Hartley, 2001). Hawaiian Pelea was divided into four sections based mainly on the grade of carpel connation: Apocarpa, Cubicarpa, Megacarpa, and Pelea. Since the incorporation of the genus Pelea into Melicope, these sections have not been formally recognized within the larger infrageneric taxonomy for Melicope as recognized by Hartley (2001) but are still being used informally as species groups (Appelhans et al., 2014a), and we refer to them as Stone’s sectional species groups (Stone’s sections) from here on. The most current and comprehensive taxonomic treatment of Hawaiian Melicope was considered “provisional” by the authors (Stone et al., 1999), as species boundaries are difficult to define in some cases. Examples include three described species complexes, where the incorporated species vary from each other primarily in the degree of fruit pubescence; the Melicopeelliptica complex based mainly in Oʻahu (six species), the Hawaiian-based Melicopevolcanica complex (four species), and the Kauaʻi-based Melicopekavaiensis complex (five species) (Stone et al., 1999).

In contrast to other successful island radiations, the colonization of the Hawaiian Archipelago in Melicope was not preceded by a recent polyploidization event. In general, the genus Melicope shows a uniform chromosome number (Paetzold et al., 2018). To date, phylogenetic relationships in Hawaiian Melicope have been investigated in four molecular studies (Harbaugh et al., 2009; Appelhans et al., 2014a; Appelhans et al., 2014b; Appelhans et al., 2018a), with a combination of up to six nuclear and plastid genomic regions amplified using polymerase chain reaction. Hawaiian Melicope was shown to be derived from a single colonization event (Harbaugh et al., 2009). The origin of the lineage was dated to the Mid or Late Miocene (Appelhans et al., 2018a), predating the age of Kauaʻi and Niʻihau (Price and Clague, 2002). In addition, the Hawaiian endemic genus Platydesma is nested within Melicope as a monophyletic sister group to the Hawaiian species and has since been reduced (Appelhans et al., 2017). Statistically supported incongruences between individual genomic regions were not observed, yet the resolution of relationships within and among the clades was in general medium to poor (Harbaugh et al., 2009; Appelhans et al., 2014a; Appelhans et al., 2014b; Appelhans et al., 2018a). However, two independent Hawaiian origins of the Marquesan Melicope radiation, which encompasses seven species, were inferred (Appelhans et al., 2014a; Appelhans et al., 2014b; Appelhans et al., 2017).

Restriction site-associated DNA sequencing (RAD-seq; Miller et al., 2007; Baird et al., 2008) is among the most frequently used reduced representation methods employed in plant systematics. To date, most phylogenetic RAD-seq studies have focused mostly on populations or closely related species (Ree and Hipp, 2015; Díaz-Arce et al., 2016; Hodel et al., 2017). However, a simulated RAD investigation in Drosophila revealed the method to be potentially applicable in groups aged up to 60 Myr (Rubin et al., 2012). Since then, application to deeper species-level relationships has increased (e.g., Eaton and Ree, 2013; Hipp et al., 2014; Eaton et al., 2017), facilitated by the development of RAD-seq assembly pipelines targeted at phylogenetic research (Eaton, 2014).

Incongruence between datasets has been a long-standing occurrence in molecular phylogenetic inference, traditionally manifesting as incongruences between different gene trees. The advance of next-generation sequencing (NGS) technology has shown that the issue is not solved by merely incorporating more data (Jeffroy et al., 2006). There are three possible categories of confounding information in a phylogenetic study: noise, systematic error, and an underlying biological signal. Noise is an effect of the inherently stochastic nature of sequence evolution and leads to a deterioration of phylogenetic signal over time. As such, noise most heavily impacts very small datasets and deep nodes (Misof et al., 2014). Incongruence may also reflect a true biological signal, for example, the presence of ILS or non-tree-like evolution, i.e. introgression, hybridization, or recombination (Misof et al., 2014; Salichos et al., 2014). Effects of hybridization range from introgression of individual alleles, to organelle capture, to hybrid speciation (Currat et al., 2008; Stegemann et al., 2012; Twyford and Ennos, 2012). Either of these processes will result in discordant gene trees, and several approaches have been proposed to unravel them. Based on the distributions of conflicting phylogenetic patterns in the genome, it is possible to distinguish the more stochastic signal of ILS from the directional and asymmetric signal of hybridization (Durand et al., 2011).

Here, we apply RAD-seq to Hawaiian Melicope, a lineage with a crown age of ca. 10 Myr (Appelhans et al., 2018a). We use RAD-seq to infer species-level relationships in the lineage, in a phylogenetic context of several colonization events of individual islands, multiple possible bottlenecks, and adaptive radiations within a lineage. The taxonomic implications of our phylogenetic results are discussed within the framework of evidence for both ancient and current introgression.

Materials and Methods

Taxon Sampling

Table 1 details the identity and origin of the 101 samples of this study: 6 outgroup and 95 ingroup specimens representing 41 Hawaiian species (81% of the lineage). Two samples represent the two independent colonization events to the Marquesas Islands (28% of Marquesan species). Taxonomic treatment follows species recognized in Wood et al. (2016) plus a recently described species (Wood et al., 2017) and including Platydesma (Appelhans et al., 2017). Additionally, morphologically divergent specimens of Melicope barbigera (KW16722 and KW16718) and Melicope ovata (KW16762, KW17082, and MA663) were included (Table 1, asterisk) to elucidate whether these might represent separate taxa. We also included two specimens, KW17111 and KW15733, which correspond closely, though not entirely, to the description of Melicope wawraeana as delimited by Stone et al. (1999). Even the Oʻahu populations that were considered the core of M.wawraeana are variable, suggesting that it is a potentially artificial taxon (Stone et al., 1999). Since the morphology of the two specimens did not correspond entirely to the Oʻahu populations considered to be M.wawraeana, we included them here as Melicope sp. (Table 1).

RAD Library Preparation

DNA was extracted from silica-dried material using the Qiagen DNeasy Plant Mini Kit^® (Qiagen, Hilden, Germany) as per the manufacturer’s instructions with incubation in lysis buffer elongated to 2h. DNA concentration was measured using the Qubit^® fluorometer and the Qubit^® dsDNA BR Assay Kit (Thermo Fisher Scientific, Darmstadt, Germany) and adjusted to 30 ng/µL. Floragenex Inc. (Portland, Oregon, USA) generated RAD libraries using the restriction enzyme SfbI. With a method following Baird et al. (2008) being employed, including the use of sample-specific barcodes, the samples were sequenced on an Illumina^® GAIIx platform to produce 100-bp single-end reads.

RAD Locus Assembly

Quality of raw reads was checked using FastQC (Andrews, 2010). The program ipyrad v.0.7.21 was used to demultiplex raw reads allowing a mismatch of 1 bp. Raw reads were trimmed using cutadapt v.1.9.1 (Martin, 2011) as implemented in ipyrad by removing adapter sequences, trimming bases with Phred scores <30 and removing reads shorter than 35 bp after trimming. Trimmed reads were assembled de novo using the ipyrad pipeline. The software attempts to evaluate orthology by scoring alignments of reads or sequences, as opposed to assessing purely sequence identity (Eaton, 2014). The alignment score is the user-determined clustering threshold to be met. To reduce the risk of introducing assembly error to our dataset, we performed a modified clustering optimization approach (Paris et al., 2017). We iterated over core clustering parameters and plotted assembly matrices (cluster depth, heterozygosity, number of putatively paralogous loci, number of single-nucleotide polymorphisms (SNPs)) to identify parameters introducing excessive assembly errors (Paetzold et al., unpublished results; Paris et al., 2017). In addition, we optimized the clustering of reads within each individual sample and the clustering of consensus sequences across loci separately, reasoning that the divergence found within each individual genome might be significantly different from the ca. 10 Myr of divergence (Appelhans et al., 2018a) within the lineage as a whole. Thus, the assembly was generated using a clustering threshold of 95 for in-sample clustering and 90 for between-sample clustering. The final filtering of loci was performed for values 10, 32, 50, 67, and 85 as the minimum numbers of samples per locus.

Phylogenetic Inference and Quartet Sampling

Phylogenetic inference was performed on all resulting alignments using maximum likelihood (ML) and Bayesian inference (BI). As individual loci are very short and may comprise a high fraction of missing data, a partitioned analysis is neither computationally feasible nor expected to produce reliable results. Thus, all datasets were analyzed solely concatenated. ML was performed using ExaML v3.0.2 (Kozlov et al., 2015) using the new rapid hill-climbing algorithm, a random number seed, the gamma model of rate heterogeneity, and the median for discrete approximation of rate heterogeneity. For datasets containing minimum numbers of 10, 32, and 50 samples, the memory saving option for gappy alignments was activated (-S). Parsimony starting trees were generated using RAxML v8.2.4 (Stamatakis, 2014). RAxML was also used to generate 100 bootstrap replicate alignments and their corresponding parsimony starting trees. ExaML searches were run on every replicate alignment with the above-mentioned settings.

BI was performed using ExaBayes v 1.5 (Aberer et al., 2014). Four independent runs were carried out with a convergence stopping criterion (split frequencies average <5% in three subsequent generations) and for a minimum of 100,000 generations sampling every 100th generation under the GTR+I+G model. Majority rule consensus trees were drawn on topologies of all four runs combined after the first 25% was discarded as burn-in.

Analysis of large-scale, concatenated datasets can result in erroneous relationships with high bootstrap support because of a failure to model the effects of ILS (Gadagkar et al., 2005; Kubatko and Degnan, 2007; Seo, 2008). These effects can be driven by only a few loci (Shen et al., 2017) and especially pertain to short branches (Kumar et al., 2012). On the other hand, a simulation study has shown that concatenated analysis of datasets containing loci with anomalous gene trees will more likely result in unresolved species tree topologies, rather than highly supported false ones (Huang and Knowles, 2009).

Methods implementing the multispecies coalescent (MSC) model explicitly incorporate gene tree conflict into species tree inference and are thus more robust to ILS than concatenation approaches (Kubatko and Degnan, 2007) but are often intractable for large datasets (Liu et al., 2015). Summary methods of species tree inference under the MSC, for example, ASTRAL (Mirarab et al., 2014) or NJst (Liu and Yu, 2011), are based on the analysis of individual gene trees and have become popular due to their comparative speed and accuracy. However, the limited information content of individual RAD loci often limits their application for gene tree inference, which may negatively impact species tree estimation (Salichos and Rokas, 2013; Mirarab et al., 2016). Alternatively, site-based methods avoid estimation of gene trees, instead using SNP data directly, and so are expected to be well suited to short, low-variability loci (Molloy and Warnow, 2018). We employed the SVDQuartets method, which infers quartet trees from SNPs using phylogenetic invariant patterns under the coalescent model and then infers the species tree by quartet joining of the subtrees using algebraic statistics (Chifman and Kubatko, 2014). We converted the SNP datasets into nexus format using the Ruby script convert_vcf_to_nexus.rb (Matschiner, 2018). The SVDQuartets analysis was computed as implemented in the software PAUP*4.0a (Swofford, 2002; Swofford, 2018). We analyzed 250,000 randomly selected quartets and assessed statistical support using 100 nonparametric bootstrap support replicates. For ambiguous positions in the SNP matrix, we chose the “Distribute” option, as these positions represent heterozygous sites.

To estimate the robustness of resolved relationships, we employed the Quartet Sampling method, which aims to measure branch support in large sparse alignments (Pease et al., 2018). As each internal branch divides all samples within a phylogeny into four non-overlapping subsets, the method randomly samples one taxon per subset to produce a quartet phylogeny. The topology of each quartet is either concordant with the tree topology or discordant. Discord is measured and quantified to produce four metrics—quartet concordance (QC), quartet differential (QD), quartet informativeness (QI), and quartet fidelity (QF)—allowing effective assessment of branch-related (QC, QD, and QI) and taxon-related (QF) discordance in the dataset (Pease et al., 2018). The method is implemented in the python script quartet_sampling.py (https://www.github.com/fephyfofum/quartetsampling). We performed Quartet Sampling on all datasets and the respectively resolved topologies using 500 replicates per branch with a minimum required overlap of 300,000 bp in the min10, min32, min50, and min67 concatenated datasets. The minimum overlap was lowered to 140,000 bp in the min85 concatenated dataset, as otherwise five samples would have been excluded from the analysis.

Test for Introgression

The D-statistics (Durand et al., 2011) is a site-based test for introgression. In a four-taxon topology (((P1, P2), P3), O), a derived allele in the P3 lineage is expected to occur also in either P1 or P2 with equal frequency, giving rise to either an ABBA or BABA discordant site pattern (Durand et al., 2011). A statistically significant imbalance in these site pattern frequencies provides evidence of introgression, while equal frequencies are associated with neutral processes like ILS. Unfortunately, this test is not well suited for deeper evolutionary timescales, where the P3 lineage has diverged into multiple sub-lineages, and it also does not allow inference of direction of introgression. Partitioned D-statistics is a system of multiple four-taxon D-statistics in a symmetric, five-taxon phylogeny with the ingroup taxa forming two pairs (P1, P2) and (P3₁, P3₂) and an outgroup taxon (O) (Figure 1) (Eaton and Ree, 2013). The partitioned D-statistics identifies sites, in which either or both of the P3 lineages share a derived allele with either P1 or P2, but not both (Figure 1) (Eaton and Ree, 2013).

Figure 1

Eaton and Ree, 2013).

We used partitioned D-statistics to infer whether discordant relationships inferred between major clades (see below) are caused by ILS or introgression. We defined entire clades as lineages and tested all combinations obeying the symmetric topology.

Results

Raw Data and Assembly

Illumina Sequencing yielded an average of 10,439,082 reads per sample (342,914–34,663,109). After quality trimming, an average of 10,327,562 reads per sample (271,257–34,542,777) were left. The assembled dataset contained a total of 786,169 clusters prior to filtering by sample coverage. Filtering reduced the number of loci by over 90% (Table 2). The final datasets contained between 7,266 (min85) and 59,041 (min10) loci. The number of variable sites (SNPs) ranged from 529,045 (min10) to 82,760 (min85) (Table 2).

Phylogenetic Inference

All five final datasets were used for phylogenetic inference in concatenated BI, ML, and SVD Quartets analyses. Statistical support for inferred relationships was assessed using posterior probabilities (PPs), nonparametric bootstrap (NBS) (ML-NBS and SVD-NBS), and Quartet Sampling. Analyses of the five datasets resulted in mostly congruent relationships, with few exceptions (see below). NBS and PP values are very high across the trees. QI values are high for all nodes (>0.9), and QF scores are average between 0.83 and 0.88 across datasets. Figure 2 shows the result of phylogenetic inference in the concatenated min32 dataset.

Figure 2

Hawaiian Melicope are divided into five main clades corresponding to those previously resolved by Appelhans et al. (2014b). These five clades are fully supported by all statistical methods. The former genus Platydesma represents the earliest diverging lineage (clade V; Figure 2). Clade IV corresponds to Stone’s section Pelea, characterized by whorled leaves. The remaining Stone sections appear to be non-monophyletic. Species ascribed to Stone’s section Apocarpa are resolved as two independent lineages (Clades II and III). Clade I comprises all species of Stone’s sections Cubicarpa and Megacarpa intermingled (Figure 2). Relationships of clade III were resolved incongruently between datasets and analyses. BI and ML analyses resolved clade III as sister to clade IV, and the resulting monophyletic lineage again in a sister–group relationship to clades I + II with maximum PP and high ML-NBS support in four of the datasets (min10, min32, min50, and min85), yet with some discord detected by Quartet Sampling (Figures 2 and 3, Supplemental Figures 1, 2, and 4). The concatenated min67 dataset resolves clade III as sister to clades I + II, and clade IV as sister to clades I + II + III (Supplemental Figure 3) with medium statistical support. Coalescent-based SVDQuartets analysis of SNP datasets resolved a third alternative topology. Here, clade II is resolved as sister to clade III, and the resulting lineage is sister to clades I + IV. This topology receives medium-to-low SVD-NBS support across all SNP datasets, as well as medium-to-high negative QC values, indicating substantial counter-support for this relationship (Supplemental Figures 5–9). The relationship of clade III is highly discordant over quartet replicates (Supplemental Figure 3). Across all datasets, the discord detected by Quartet Sampling for the ancestral branch is skewed favoring one of the tested alternative quartet topologies (QD; Figure 3, Supplemental Figures 1–9).

Figure 3

The remaining relationships within individual clades are fully resolved, improving resolution to the species and intraspecies levels (Figure 2). The majority of all Hawaiian Melicope are resolved in clade I, and relationships among species show many nodes with notable discord and very short branches (Figure 2). Most of the nodes show low QC and medium-to-low QD values (Figure 3). Three samples show incongruent relationships between datasets. This pertains to the Marquesan Melicope hivaoaensis, which is resolved in clade I as either sister to the remaining species (Supplemental Figures 3–9) within the clade or diverging prior to Melicope lydgatei (Figures 2 and 3, Supplemental Figures 1 and 2) as well as to M. kavaiensis and Melicope sp. KW15773 (Figure 3, Supplemental Figures 1–9). In all datasets, QC values show high discord or even counter-support for the placement of these three specimens. However, while QD and QF values are high for M. hivaoaensis, for both M. kavaiensis and Melicope sp. KW15773, QD values are low and QF scores are below average (0.47–0.6 for M. kavaiensis) (Figure 3, Supplemental Figures 1–9). The remaining relationships in clade I are congruent among all concatenation-based analyses. Site-specific coalescence analysis, however, resolved largely incongruent relationships for taxa in this clade, especially pertaining to the most recent divergences. The inferred relationships receive medium-to-very-low SVD-NBS values and show a high amount of discord in Quartet Sampling (Supplemental Figures 5–9).

Clades III and IV are subdivided into two subclades each. Most species sampled with multiple accessions are resolved as monophyletic with high support and no discord detected in Quartet Sampling. Exceptions are Melicope clusiifolia, Melicope haupuensis, Melicope knudsenii, and Melicope feddei. M. clusiifolia is resolved paraphyletic with respect to Melicope haleakalae, which is nested within clade IVB with high-to-maximum support. Specimens of M. haupuensis are resolved as polyphyletic within clade IIIB. The relationships among the three sampled taxa are not resolved consistently across datasets and poorly supported. Quartet Sampling reveals a high level of discord and below-average QF scores (Figure 3, Supplemental Figures 1–9). M. knudsenii is also resolved as polyphyletic with two Maui specimens (MA629 and H41610) monophyletic in clade IIIA, while the third sample from Kauaʻi (KW17119) is resolved as sister to M. barbigera in clade IIIB (Figure 2). Either relationship is virtually uncontested (Figures 2 and 3, Supplemental Figures 1–9). M. feddei is paraphyletic with respect to one of the Kauaʻi M. wawraeana-like specimens (KW17111). The three individuals form a fully supported, monophyletic unit (Figures 2 and 3, Supplemental Figures 1–9).

None of the three species complexes (M. elliptica, M. kavaiensis, and M. volcanica complexes) are resolved as monophyletic. Species of both the M. kavaiensis and M. volcanica complexes are resolved in clade I (Figure 2) in proximity to each other, but not sister to each other. Species of the M. elliptica complex are resolved in different subclades of clade III (Figure 2). Both M. barbigera and M. ovata were resolved as monophyletic, and the morphologically divergent specimens (Table 1, asterisk) are resolved as sister clades to the samples with the typical morphology of the respective species with high support (Figure 2).

The species from the Marquesas Islands are deeply nested within the Hawaiian clade.Melicope inopinata is resolved in clade III as sister to the rest of subclade IIIA. M. hivaoaensis represents a group of six morphologically similar species that form a highly supported monophyletic clade (Appelhans et al., 2014b; Appelhans et al., 2018a) and is nested within clade I here (Figure 2).

Test for Introgression

The min32 dataset was used for the ABBA-BABA test, since it produced the highest number of fully supported nodes. The tree topology in Figure 2 was chosen to represent the species tree topology, as it was recovered by the majority of analyses. The D-statistics was only used to test the incongruent position of clade III, as for incongruent species within clade I, the sampling of the respective populations is not sufficient to draw reliable conclusions. Samples within clades were pooled, and SNP frequencies were used for D-statistic calculations (Durand et al., 2011). All possible relationships complying with the D-statistic assumptions were tested. A total of 24,673 loci covered at least one-third of all samples per clade and, thus, contributed to the test results. Table 3 summarizes the tested topologies and inferred partitioned D-statistics. When clades III and IV are tested as donors for introgressed loci, values for D₁₂ are small and not significant (Z₁₂ < 2.55), while values for D₁ and D₂ are significant, respectively. For tests with either of clade I or II designated as P3 lineages, D₁₂, and D₁ and D₂, are all significant (Table 3). For all tested configurations, the dataset exhibits more than 3,000 discordant site patterns (Table 3).

Discussion

Phylogeny and Introgression

Analysis of ipyrad assemblies consistently resolved five major clades within Hawaiian Melicope (Figure 2). However, the relationships of clade III were incongruent among the five datasets and analysis methods (Figure 2, Supplemental Figures 1–9). Incongruence between datasets may be caused by one of three factors: noise, ILS, or non-tree-like evolution. As noise is expected to impact small datasets and deep nodes most severely (Misof et al., 2014), it is unlikely a sufficient cause of the incongruence observed here, since our RAD-seq alignments are substantial in size (Table 1) and the remaining deep nodes are not affected.

The QD values of the branch illustrate that one of the discordant topologies is inferred significantly more often (0.0–0.4; Figure 3, Supplemental Figures 1–9), which indicates non-tree-like evolution as the cause for the discord. Thus, we used the partitioned D-statistics to test for signals of ancient introgression between clades I through IV with all clades tested as putative donor (P3) lineages. In all cases, values for D₁ and D₂ were each significant, yet values for D₁₂ were only significant when clades I and II were defined as P3 (Table 3). Positive values of D₁ represent introgression between P2 and P3₁, while negative values indicate introgression between P1 and P3₁, and values for D₂ represent events analogous for P3₂ and P2 (Eaton and Ree, 2013; Pease and Hahn, 2015). The significant values for D₁ and D₂ indicate introgression between the respective ancestors of clades I and IV as well as between respective ancestors of clades II and III. Significant values for D₁₂ represent shared ancestral alleles from the clade I + II progenitor introduced into the respective ancestor of clades III and IV (Figure 2, Table 3). All taxa in clades II and III have apocarpous fruits, while all taxa in clades I and IV have syncarpous fruits (Stone et al., 1999), providing a morphological connection between either of the two pairs, which might be linked to introgressed information. However, we interpret these result cautiously, as D-statistic results are sensitive to confounding signals from multiple introgressive events due to phylogenetic non-independence of tests (Eaton et al., 2015).

The origin of the Hawaiian Melicope lineage predates the rise of the current high islands (Appelhans et al., 2018a). Thus, the inferred introgressive events are associated with a time when the ancestral species were still relegated either to refugial areas on small, low islands or shortly after they colonized the young island of Kauaʻi. The time frame under consideration presents a “bottleneck” scenario, where the ancestral lineages were likely in close spatial proximity. Additionally, increased volcanic activity of the Hawaiian hot spot coincided with the rise of Kaua'i (Price and Clague, 2002). This volcanic activity could have produced lava flows, earthquakes, tsunamis, and other catastrophic events, which may have additionally promoted hybridization (Stuessy et al., 2014). The ancestral hybridization events may even have promoted subsequent adaptive radiation on the islands (Kagawa and Takimoto, 2018). Estimation of divergence times in Hawaiian Melicope will be needed to infer the time frame for hybridization events in ancestral lineages. While there is strong evidence for ancient hybridization events within Hawaiian Melicope, the nature of de novo RAD-seq data currently limits our analytic methods. Further information may be obtained through gene tree-based approaches applied to target capture or whole genome-sequencing data (Meng and Kubatko, 2009) or by examining SNP-based patterns, as they vary spatially along a reference genome (Martin et al., 2013).

Bootstrap and PP support values were generally high across trees inferred from different datasets but generally increased with dataset size. Lenient filtering in RAD-seq data is often practiced, as there is a correlation between the size of a data matrix and resolution and support of relationships (Wagner et al., 2013; Hodel et al., 2017). RAD locus dropout is expected to increase with increasing divergence times, as enzyme cut sites will be lost or gained through mutation (Cariou et al., 2013). Loci with a small amount of missing data are therefore expected to represent the conserved spectrum of genomic sites and, thus, provide a limited capacity of resolution. On the other hand, sparse loci are expected to iincrease resolution of relationships despite also introducing noise, as they are assumed to represent the more rapidly evolving genomic fractions (Cariou et al., 2013; Wagner et al., 2013; Eaton et al., 2015). However, including all loci is not advisable either, as there seems to be a point at which inclusion of increasingly more sparse loci might start to decrease support. At this point, noise, due to missing data introduced by the inclusion of more sparse loci, will overpower the informative value these loci provide. However, the Quartet Sampling method seems an adequate approach to evaluating the reliability of the dataset, as the QC value showed the same trend in all datasets regardless of size and offer the QI score to assess the amount and impact of missing data.

We detected some discord between relationships resolved by concatenation and site-specific coalescence-based methods (Figure 2, Supplemental Figures 1–9). The evaluation of the performance of different species-tree inference methods is a matter of ongoing research, especially with regard to genomic datasets. Concatenation-based ML inference can be statistically inconsistent under some conditions in the MSC, that is, ILS causing gene trees to differ from the true species tree (Kubatko and Degnan, 2007). However, the limits of the concatenated approach are poorly understood (Molloy and Warnow, 2018), and the performance of concatenated Bayesian analysis has yet to be formally assessed. Some simulation studies show that concatenated RAD-seq data are robust to gene tree/species tree discord when inferring relationships among taxa (Rivers et al., 2016). In addition, concatenated approaches potentially offer hidden support as a feature overriding gene tree/species conflict (Gatesy and Springer, 2014; Rivers et al., 2016), although hidden support has not been addressed in plant phylogenomic research yet. Coalescence-based methods are statistically consistent under the MSC. Bayesian co-estimation of gene trees and the species tree under the MSC is currently considered the most effective approach, yet computationally very demanding and thus less applicable to large datasets. Hence, summary and site-specific MSC methods have become popular, and several algorithms implementing the concepts do exist (Liu et al., 2015). However, the assessment of the performance of these methods under empirical and simulated conditions is still a matter of active research. For example, gene tree methods have proven to be statistically inconsistent if the cause of gene tree discord is horizontal gene transfer, instead of ILS (Solís-Lemus et al., 2016; Fernández-Mazuecos et al., 2018). Several recent simulation studies compared the accuracy of multiple summary and site-based coalescent methods, including SVDQuartets, as well as concatenated ML under varying levels of ILS and gene tree estimation error (GTEE). Concatenated ML was at least competitive with MSC methods under most conditions and outperformed SVDQuartets under all tested conditions, including high GTEE (Chou et al., 2015; Mirarab et al., 2016; Molloy and Warnow, 2018). The latter would be expected in RAD-seq datasets and should also be present herein.

With respect to species relationships inferred for Hawaiian Melicope and considering the observed lower accuracy of SVDQuartets compared with concatenation-based approaches under conditions typically characterizing RAD datasets, we suggest that the results from concatenated BI and ML are probably more accurate than those based on SVDQuartets and will be discussed below. However, we do stress that none of the approaches have proven to be statistically consistent under conditions observed herein, that is, ILS, GTEE, and horizontal gene transfer (Figure 2).

Taxonomic Implications

The former small genus Platydesma and Stone’s section Pelea are each monophyletic (Figure 2), while the three remaining sections of Stone, comprising the majority of all Hawaiian Melicope species, are not. Apocarpa is divided into two lineages with the majority of species resolved in Apocarpa 1 (Figure 2). The three species of the Apocarpa 2 clade share a number of morphological traits, though none of them is either exclusive or inclusive. All species of Apocarpa 2 occur in mesic forests only and, with the exception of Melicope stonei, share a sprawling, shrubby habit (Stone et al., 1999; Wood et al., 2017). Finally, in all Apocarpa 2 species, both endocarp and exocarp are glabrous and inflorescences are few-flowered, though both of these traits also appear outside of this group (Stone et al., 1999; Wood et al., 2017). In a previous analysis, apocarpous species were resolved in three different clades (Appelhans et al., 2014b), one of which, consisting of M. elliptica only, could not be sampled in this study. Further research will be necessary to identify morphological character combinations distinguishing these lineages. Stone’s sections Cubicarpa and Megacarpa are paraphyletic with respect to each other (Figure 2) with species of each resolved intermingled throughout the clade. The two groups differ by the degree of carpel connation, with carpels “connate from base up to 2/3 of their length” (Stone et al., 1999) characterizing Megacarpa and carpels “nearly to completely” connate (Stone et al., 1999) characterizing Cubicarpa. Carpel connation clearly represents a continuum and not two discrete units. As there is no pattern to the degree of carpel fusion apparent in clade I, the separation of these two of Stone’s sections seems artificial.

Interspecies relationships within clade I are less well supported than in the remaining clades, and Quartet Sampling reveals measurable discord at nearly every branch in the backbone of this clade. For many of the nodes with low QC values, QD values are high (Figure 3, Supplemental Figures 1–9), which characterizes ILS and corresponds to the shortness of these branches. On the other hand, many branches show low QD values, indicating widespread introgression between these lineages. Unfortunately, sampling herein is not sufficient to test individual relationships.

Of the 24 species represented by multiple accessions, 20 were resolved as monophyletic, while four species were either paraphyletic or polyphyletic. M. clusiifolia is the most widespread and morphologically diverse of all Hawaiian Melicope (Stone et al., 1999), and it is paraphyletic with both of the other species of Stone’s section Pelea, M. haleakalae and Melicope waialealae (clade IV, Figure 2). Several attempts have been made to subdivide M. clusiifolia into varying constellations of subspecies, varieties, and forms (St. John, 1944; Stone, 1969). In the most recent taxonomic treatment, Stone et al. (1999) synonymized all subdivisions of the species, arguing that the variable characters seem to represent a continuum rather than distinguishable, discrete units. However, the authors also issued the recommendation that the overall pattern of variability in M. clusiifolia should be studied in detail (Stone et al., 1999). M. haleakalae is characterized as differing from M. clusiifolia, mainly in its persistent sepals (Stone et al., 1999). Considering that M. haleakalae is nested deeply within M. clusiifolia (Figure 2, clade IV), the two might be regarded as conspecific and included in an overall evaluation of the complex. M. waialealae differs from M. clusiifolia mainly in leaf shape (Stone et al., 1999). However, since the leaf shape of M. clusiifolia is highly variable, M. waialealae might represent one end of a continuum across both taxa rather than one of two distinct states. On the other hand, these three species might represent a case of speciation in progress. In this case, the deep nesting, especially of M. haleakalae, within M. clusiifolia, would represent speciation following a progenitor-derivative scenario (Crawford, 2010). The widespread, morphologically variable M. clusiifolia would meet all criteria of the progenitor (p) species. The persistent petals in M. haleakalae and the leaf shape in M. waialealae would represent a variable, morphological feature in the parent being fixed in the respective derivative (d) species. Identification of a true p–d relationship is difficult and rare. However, several candidate species pairs do exist (Crawford, 2010). The p–d species pair Layia glandulosa (Hook.) Hook. & Arn. and Layia discoidea D. D. Keck (Asteraceae) show not only a shift in morphology between progenitor and derivative species but also geographic isolation due to a shift in habitat (Baldwin, 2005). This could be the same for M. waialealae, which is restricted to bogs, whereas the putative progenitor M. clusiifolia occurs in mesic to wet forests (Stone et al., 1999). Unfortunately, there are no data available regarding breeding system or pollinator communities in these species, creating potential barriers to gene flow. Detailed studies of morphological characters, gene flow, and abiotic habitat factors are necessary to determine whether these taxa are separate p–d species pairs or conspecific, as already indicated in previous studies (Appelhans et al., 2014b).

M. knudsenii, delimited by Stone et al. (1999) as the only species occurring on non-adjacent islands, was resolved as polyphyletic, with three samples resolved as two distinct lineages within clade III. Appelhans et al. (2014b) already showed that this taxon is polyphyletic, consisting of three taxa. One of these was recently described as M. stonei (Wood et al., 2017). Our results confirm the previously resolved pattern with the two specimens of M. knudsenii from Maui resolved as sister to Melicope hawaiensis and the specimen from Kauaʻi as sister to M. barbigera (clade III, Figure 2). We confirm that these specimens clearly represent different species. The Maui species will be resurrected under one of the names used in an earlier treatment by Stone (1969), wherein he adopted a narrower species concept than in the later classification (Stone et al., 1999), leaving M. knudsenii restricted to only populations on Kauaʻi.

The three specimens of M. haupuensis included in this study are resolved as paraphyletic. Moreover, they are the only species resolved with incongruent topologies of the individual samples associated with the different datasets (compare Figure 2, Supplemental Figures 1–9). Quartet Sampling shows strong discord for either of the inferred relationships with medium QD values (Figure 3, Supplemental Figures 1–9), indicating the possibility of introgressed sites. Moreover, QF scores for the three specimens are considerably lower than the average, indicating a rogue behavior (Aberer et al., 2013; Pease et al., 2018) of the three taxa. Additionally, several specimens in the field were observed presenting morphologically intermediate forms between M. haupuensis and M. barbigera (personal observation K.R. Wood). QD values for the latter are also low (Figure 3, Supplemental Figures 1–9). Both the morphological intermediates and the incongruence associated with different datasets indicate potential hybridization between these species. However, conclusively identifying putative hybridization events would require sampling at the population level, including any morphological intermediates.

Multiple samples of M. ovata and M. barbigera were included in our analyses, representing both the typical morphology and a deviating morphotype. For either species, the morphologically deviating samples were resolved as the sister group to the samples with the typical habit. Variant morphotypes of M. ovata displayed a pubescent lower leaf surface, whereas leaves are typically glabrous in this species. M. barbigera usually has few-flowered inflorescences (Stone et al., 1999). In contrast, the variant morphotype has inflorescences with a considerably larger number of flowers. Genomic divergence is comparable with that of other species pairs within the lineage. Both groups might be another case of speciation in progress within Hawaiian Melicope. In both cases, detailed morphological studies will be necessary to investigate if the morphologically divergent populations of the two species should be recognized as separate taxa.

The two M. wawraeana-like specimens are resolved in clade I, but not closely related to each other. One specimen (KW17111) is nested within the two samples of M. feddei with high support (Figures 2 and 3). M. wawraeana is very similar to M. feddei and differs mainly in pedicel length (Stone et al., 1999). The present results suggest that some populations might be conspecific with M. feddei, while others (e.g., from the herein unsampled type location) are not. The relationships of the second M. wawraeana-like specimen (KW15733) are resolved incongruently among datasets, as are the relationships of the sampled specimen of M. kavaiensis. The two samples are resolved either as sister groups (Figure 3, Supplemental Figures 1, 4, 5, and 7–9) or as consecutive sister clades within clade I (Supplemental Figures 2, 3, and 6). There is a substantial amount of discord in the dataset for either of the resolved relationships. QD values are low, indicating the possibility of introgression between these morphologically distinct species. Additionally, QF scores for either of the specimens are low corresponding to the rogue behavior of the samples.

The rogue behavior of the aforementioned samples (M. kavaiensis, Melicope sp. KW15733) might also be related to the incongruent placement of M. hivaoaensis, as the three taxa are inferred as closely related, regardless of the relation to the remainder of clade I. For this specimen, QC and QD values are low; however, QF is high (Figure 3, Supplemental Figures 1–9). M. hivaoaensis represents an adaptive radiation of five species endemic to the Marquesas Islands, whose predecessor colonized from the Hawaiian Islands (Appelhans et al., 2014a; Appelhans et al., 2018a). Successful island colonizations have been associated with recent hybridization or polyploidization events (Paetzold et al., 2018). There was no polyploidization event immediately prior to the colonization of the Hawaiian Islands itself (Paetzold et al., 2018), making a polyploidization event prior to the colonization of the Marquesas Islands unlikely. Chromosome counts for Marquesan species are not available for a conclusive answer. However, results herein indicate the presence of several hybridization events within the lineage. Thus, a hybridization event might have predated the colonization of the Marquesas Islands as well. As the incongruent position of M. hivaoaensis seems to correspond to the rogue behavior of Melicope sp. KW15733 and M. kavaiensis, the latter two might represent the parental lineages of the Marquesan Melicope radiation. A conclusive answer to the question is contingent on a thorough sampling of all concerned lineages as well as a prior revision of the M. wawraeana species concept.

We confirm previous results showing that Hawaiian Melicope colonized the Marquesas Islands twice independently, negating the hypothesis that the remote Hawaiian Islands constitute a dispersal sink (Harbaugh et al., 2009; Appelhans et al., 2014a; Appelhans et al., 2018a). The nesting of Marquesan species in different Hawaiian clades is corroborated by fruit morphology (Hartley, 2001), since M. hivaoaensis and its close relatives from the Marquesas Islands have syncarpous fruits as do the species in clade I, while M. inopinata has apocarpous fruits like the species in clade III.

The present study provides unprecedented insight into the relationships of Hawaiian Melicope. Several previous findings could be corroborated and firmly supported by genome-wide data, including the non-monophyly of most of Stone’s sections, which cannot be held up as delimited (Stone, 1969; Stone et al., 1999). The lineage is in need of a taxonomic revision. Understanding the relationships of Hawaiian Melicope would be enhanced by some formal recognition of the subclades with corresponding morphological features. However, the creation of novel formal subgroups within Melicope section Pelea must also include the extra-Hawaiian members of the section. The former genus Platydesma is the most distinctive group within Melicope sect. Pelea and should receive some level of formal recognition. Apocarpa species need to be split into two groups, one of which would include the Marquesan species M. inopinata. However, conclusive treatment of Apocarpa should be adjourned until an improved understanding of the separation within the M. elliptica complex is attained. Delimitations of species within the Pelea group, M. barbigera, M. ovata, and M. haupuensis, may need revision, but levels of hybridization should also be investigated as part of that process. M. wawraeana requires revision as well as a prerequisite to test the putative hybrid character of the Marquesan radiation. Furthermore, the other six Melicope species endemic to the Marquesas Islands would need to be included in a novel taxonomic recognition of Stone’s former sections Megacarpa and Cubicarpa.

Funding

This project was financially supported by the German Science Foundation (DFG; Grant AP 251/3-1 to MA). The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References

• 1

  AbererA. J.KobertK.StamatakisA. (2014). ExaBayes: massively parallel Bayesian tree inference for the whole-genome era. Mol. Biol. Evol.31, 2553–2556. doi: 10.1093/molbev/msu236

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AbererA. J.KrompassD.StamatakisA. (2013). Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice. Syst. Biol.62, 162–166. doi: 10.1093/sysbio/sys078

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AndrewsS. (2010). FastQC: a quality control tool from high throughput sequence data. https://www.bioinformatics.babraham.ac.uk/projects/fastqc

  • Google Scholar
  • View reference in article

• 4

  AppelhansM. S.ReicheltN.GroppoM.PaetzoldC.WenJ. (2018b). Phylogeny and biogeography of the pantropical genus Zanthoxylum and its closest relatives in the proto-Rutaceae group (Rutaceae). Mol. Phylogenet. Evol.126, 31–44. doi: 10.1016/j.ympev.2018.04.013

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  AppelhansM. S.WenJ.DurettoM.CraynD.WagnerW. L. (2018a). Historical biogeography of Melicope (Rutaceae) and its close relatives with a special emphasis on Pacific dispersals. J. Syst. Evol.56, 576–599. doi: 10.1111/jse.12299

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  AppelhansM. S.WenJ.WagnerW. L. (2014a). A molecular phylogeny of Acronychia, Euodia, Melicope and relatives (Rutaceae) reveals polyphyletic genera and key innovations for species richness. Mol. Phylogenet. Evol.79, 54–68. doi: 10.1016/j.ympev.2014.06.014

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  AppelhansM. S.WenJ.WoodK. R.AllanG. J.ZimmerE. A.WagnerW. L. (2014b). Molecular phylogenetic analysis of Hawaiian Rutaceae (Melicope, Platydesma and Zanthoxylum) and their different colonization patterns: molecular phylogeny of Hawaiian Rutaceae. Bot. J. Lin. Soc.174, 425–448. doi: 10.1111/boj.12123

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  AppelhansM. S.WoodK. R.WagnerW. L. (2017). Reduction of the Hawaiian genus Platydesma into Melicope section Pelea (Rutaceae) and notes on the monophyly of the section. PhytoKeys91, 125–137. doi: 10.3897/phytokeys.91.21363

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BairdN. A.EtterP. D.AtwoodT. S.CurreyM. C.ShiverA. L.LewisZ. A.et al. (2008). Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS One3, e3376. doi: 10.1371/journal.pone.0003376

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BaldwinB. G. (2005). Origin of the serpentine-edemic herb Layia discoidea from the widespread L. glandulosa Compositae. Evolution59, 2473–2479. doi: 10.1111/j.0014-3820.2005.tb00956.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BaldwinB. G.SandersonM. J. (1998). Age and rate of diversification of the Hawaiian silversword alliance (Compositae). Proc. Natl. Acad. Sci. U.S.A.95, 9402–9406. doi: 10.1073/pnas.95.16.9402

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  BarrierM.BaldwinB. G.RobichauxR. H.PuruggananM. D. (1999). Interspecific hybrid ancestry of a plant adaptive radiation: allopolyploidy of the Hawaiian silversword alliance (Asteraceae) inferred from floral homeotic gene duplications. Mol. Biol. Evol.16, 1105–1113. doi: 10.1093/oxfordjournals.molbev.a026200

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  CariouM.DuretL.CharlatS. (2013). Is RAD-seq suitable for phylogenetic inference? An in silico assessment and optimization. Ecol. Evol.3, 846–852. doi: 10.1002/ece3.512

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  CarlquistS. (1967). The biota of long-distance dispersal. V. Plant dispersal to Pacific Islands. Bull. Torrey Bot. Club94, 129–162. doi: 10.2307/2484044

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  CarlquistS. (1974). Island biology. New York, NY: Columbia University Press.

  • Google Scholar
  • View reference in article

• 16

  CarrG. D. (1998). “Chromosome evolution and speciation in Hawaiian flowering plants,” in Evolution and speciation of island plants. Eds. StuessyT. F.OnoM. (Cambridge, UK: Cambridge University Press), 97–119.

  • Google Scholar
  • View reference in article

• 17

  ChifmanJ.KubatkoL. (2014). Quartet inference from SNP data under the coalescent model. Bioinformatics30, 3317–3324. doi: 10.1093/bioinformatics/btu530

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  ChouJ.GuptaA.YaduvanshiS.DavidsonR.NuteM.MirarabS.et al. (2015). A comparative study of SVDquartets and other coalescent-based species tree estimation methods. BMC Genomics16, S2. doi: 10.1186/1471-2164-16-S10-S2

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  CrawfordD. J. (2010). Progenitor-derivative species pairs and plant speciation. TAXON59, 1413–1423. doi: 10.1002/tax.595008

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  CurratM.RuediM.PetitR. J.ExcoffierL. (2008). The hidden side of invasions: massive introgression by local genes. Evolution62, 1908–1920. doi: 10.1111/j.1558-5646.2008.00413.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  Díaz-ArceN.ArrizabalagaH.MuruaH.IrigoienX.Rodríguez-EzpeletaN. (2016). RAD-seq derived genome-wide nuclear markers resolve the phylogeny of tunas. Mol. Phylogenet. Evol.102, 202–207. doi: 10.1016/j.ympev.2016.06.002

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  DurandE. Y.PattersonN.ReichD.SlatkinM. (2011). Testing for ancient admixture between closely related populations. Mol. Biol. Evol.28, 2239–2252. doi: 10.1093/molbev/msr048

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  EatonD. A. R. (2014). PyRAD: assembly of de novo RADseq loci for phylogenetic analyses. Bioinformatics30, 1844–1849. doi: 10.1093/bioinformatics/btu121

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  EatonD. A. R.ReeR. H. (2013). Inferring phylogeny and introgression using RADseq data: an example from flowering plants (Pedicularis: Orobanchaceae). Syst. Biol.62, 689–706. doi: 10.1093/sysbio/syt032

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  EatonD. A. R.HippA. L.González-RodríguezA.Cavender-BaresJ. (2015). Historical introgression among the American live oaks and the comparative nature of tests for introgression. Evolution69, 2587–2601. doi: 10.1111/evo.12758

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  EatonD. A. R.SpriggsE. L.ParkB.DonoghueM. J. (2017). Misconceptions on missing data in RAD-seq phylogenetics with a deep-scale example from flowering plants. Syst. Biol.66, 399–412. doi: 10.1093/sysbio/syw092

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  EggensF.PoppM.NepokroeffM.WagnerW. L.OxelmanB. (2007). The origin and number of introductions of the Hawaiian endemic Silene species (Caryophyllaceae). Am. J. Bot.94, 210–218. doi: 10.3732/ajb.94.2.210

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  EmersonB. C. (2002). Evolution on oceanic islands: molecular phylogenetic approaches to understanding pattern and process. Mol. Ecol.16, 951–966.

  • Google Scholar
  • View reference in article

• 29

  Fernández-MazuecosM.MellersG.VigalondoB.SáezL.VargasP.GloverB. J. (2018). Resolving recent plant radiations: power and robustness of genotyping-by-sequencing. Syst. Biol.67, 250–268. doi: 10.1093/sysbio/syx062

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  Francisco-OrtegaJ.Santos-GuerraA.KimS.-C.CrawfordD. J. (2000). Plant genetic diversity in the Canary Islands: a conservation perspective. Am. J. Bot.87, 909–919. doi: 10.2307/2656988

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  GadagkarS. R.RosenbergM. S.KumarS. (2005). Inferring species phylogenies from multiple genes: concatenated sequence tree versus consensus gene tree. J. Exp. Zool. B Mol. Dev. Evol.304B, 64–74. doi: 10.1002/jez.b.21026

  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  GatesyJ.SpringerM. S. (2014). Phylogenetic analysis at deep timescales: unreliable gene trees, bypassed hidden support, and the coalescence/concatalescence conundrum. Mol. Phylogenet. Evol.80, 231–266. doi: 10.1016/j.ympev.2014.08.013

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  GivnishT. J. (1998). “Adaptive radiation and molecular systematics: issues and approaches,” in Molecular evolution and adaptive radiation. Eds. GivnishT. J.SytsmaK. J. (Cambridge, UK: Cambridge University Press), 1–54.

  • Google Scholar
  • View reference in article

• 34

  GivnishT. J.MillamK. C.MastA. R.PatersonT. B.TheimT. J.HippA. L.et al. (2009). Origin, adaptive radiation and diversification of the Hawaiian lobeliads (Asterales: Campanulaceae). Proc. R. Soc. Lond. [Biol.]276, 407–416. doi: 10.1098/rspb.2008.1204

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  HarbaughD. T. (2008). Polyploid and hybrid origins of Pacific Island sandalwoods (Santalum, Santalaceae) inferred from low-copy nuclear and flow cytometry data. Int. J. Plant. Sci.169, 677–685. doi: 10.1086/533610

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  HarbaughD. T.WagnerW. L.AllanG. J.ZimmerE. A. (2009). The Hawaiian Archipelago is a stepping stone for dispersal in the Pacific: an example from the plant genus Melicope (Rutaceae). J. Biogeogr.36, 230–241. doi: 10.1111/j.1365-2699.2008.02008.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  HartleyT. G. (2001). On the taxonomy and biogeography of Euodia and Melicope (Rutaceae). Allertonia8, 1–328.

  • Google Scholar
  • View reference in article

• 38

  HippA. L.EatonD. A. R.Cavender-BaresJ.FitzekE.NipperR.ManosP. S. (2014). A framework phylogeny of the American oak clade based on sequenced RAD data. PLoS One9, e93975. doi: 10.1371/journal.pone.0093975

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  HodelR. G. J.ChenS.PaytonA. C.McDanielS. F.SoltisP.SoltisD. E. (2017). Adding loci improves phylogeographic resolution in red mangroves despite increased missing data: comparing microsatellites and RAD-Seq and investigating loci filtering. Sci. Rep.7. doi: 10.1038/s41598-017-16810-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  HuangH.KnowlesL. L. (2009). What is the danger of the anomaly zone for empirical phylogenetics? Syst. Biol.58, 527–536. doi: 10.1093/sysbio/syp047

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  JeffroyO.BrinkmannH.DelsucF.PhilippeH. (2006). Phylogenomics: the beginning of incongruence? Trends in Genet.22, 225–231. doi: 10.1016/j.tig.2006.02.003

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  KagawaK.TakimotoG. (2018). Hybridization can promote adaptive radiation by means of transgressive segregation. Ecol. Lett.21, 264–274. doi: 10.1111/ele.12891

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  KeeleyS. C.FunkV. A. (2011). “Origin and evolution of Hawaiian endemics: new patterns revealed by molecular phylogenetic studies,” in The biology of island floras. Eds. BramwellD.Caujape-CastellsJ. (Cambridge, UK: Cambridge University Press), 57–88. doi: 10.1017/CBO9780511844270.005

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  KozlovA. M.AbererA. J.StamatakisA. (2015). ExaML version 3: a tool for phylogenomic analyses on supercomputers. Bioinformatics31, 2577–2579. doi: 10.1093/bioinformatics/btv184

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  KubatkoL. S.DegnanJ. H. (2007). Inconsistency of phylogenetic estimates from concatenated data under coalescence. Syst. Biol.56, 17–24. doi: 10.1080/10635150601146041

  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  KumarS.FilipskiA. J.BattistuzziF. U.Kosakovsky PondS. L.TamuraK. (2012). Statistics and truth in phylogenomics. Mol. Biol. Evol.29, 457–472. doi: 10.1093/molbev/msr202

  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  LensF.DavinN.SmetsE.del ArcoM. (2013). Insular woodiness on the Canary Islands: a remarkable case of convergent evolution. Int. J. Plant. Sci.174, 992–1013. doi: 10.1086/670259

  • CrossRef
  • Google Scholar
  • View reference in article

• 48

  LiuL.YuL. (2011). Estimating species trees from unrooted gene trees. Syst. Biol.60, 661–667. doi: 10.1093/sysbio/syr027

  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  LiuL.WuS.YuL. (2015). Coalescent methods for estimating species trees from phylogenomic data. J. Syst. Evol.53, 380–390. doi: 10.1111/jse.12160

  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  LososJ. B.RicklefsR. E. (2009). Adaptation and diversification on islands. Nature457, 830–836. doi: 10.1038/nature07893

  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  MartinM. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal17, 10–12. doi: 10.14806/ej.17.1.200

  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  MartinS. H.DasmahapatraK. K.NadeauN. J.SalazarC.WaltersJ. R.SimpsonF.et al. (2013). Genome-wide evidence for speciation with gene flow in Heliconius butterflies. Genome Res.23, 1817–1828. doi: 10.1101/gr.159426.113

  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  MatschinerT. (2018). Convert_vcf_to_nexus.rb. GitHub repository, https://github.com/mmatschner/tutorials/species_tree_inference_with_snp_data/.

  • Google Scholar
  • View reference in article

• 54

  MengC.KubatkoL. S. (2009). Detecting hybrid speciation in the presence of incomplete lineage sorting using gene tree incongruence: a model. Theor. Popul. Biol.75, 35–45. doi: 10.1016/j.tpb.2008.10.004

  • CrossRef
  • Google Scholar
  • View reference in article

• 55

  MillerM. R.DunhamJ. P.AmoresA.CreskoW. A.JohnsonE. A. (2007). Rapid and cost-effective polymorphism identification and genotyping using restriction site associated DNA (RAD) markers. Genome Res.17, 240–248. doi: 10.1101/gr.5681207

  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  MirarabS.BayzidM. S.WarnowT. (2016). Evaluating summary methods for multilocus species tree estimation in the presence of incomplete lineage sorting. Syst. Biol.65, 366–380. doi: 10.1093/sysbio/syu063

  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  MirarabS.ReazR.BayzidM. S.ZimmermannT.SwensonM. S.WarnowT. (2014). ASTRAL: genome-scale coalescent-based species tree estimation. Bioinformatics30, i541–i548. doi: 10.1093/bioinformatics/btu462

  • CrossRef
  • Google Scholar
  • View reference in article

• 58

  MisofB.MeusemannK.von ReumontB. M.KückP.ProhaskaS. J.StadlerP. F. (2014). A priori assessment of data quality in molecular phylogenetics. Algorithms Mol. Biol.9. doi: 10.1186/s13015-014-0022-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 59

  MolloyE. K.WarnowT. (2018). To include or not to include: the impact of gene filtering on species tree estimation methods. Syst. Biol.67, 285–303. doi: 10.1093/sysbio/syx077

  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  NepokroeffM.SytsmaK. J.WagnerW. L.ZimmerE. A. (2003). Reconstructing ancestral patterns of colonization and dispersal in the Hawaiian understory tree genus Psychotria (Rubiaceae): a comparison of parsimony and likelihood approaches. Syst. Biol.52, 820–838. doi: 10.1080/10635150390251072

  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  PaetzoldC.KiehnM.WoodK. R.WagnerW. L.AppelhansM. S. (2018). The odd one out or a hidden generalist: Hawaiian Melicope (Rutaceae) do not share traits associated with successful island colonization. J. Syst. Evol.56, 621–636. doi: 10.1111/jse.12454

  • CrossRef
  • Google Scholar

• 62

  ParisJ. R.StevensJ. R.CatchenJ. M. (2017). Lost in parameter space: a road map for stacks. Methods Ecol. Evol.8, 1360–1373. doi: 10.1111/2041-210X.12775

  • CrossRef
  • Google Scholar

• 63

  PeaseJ. B.HahnM. W. (2015). Detection and Polarization of introgression in a five-taxon phylogeny. Syst. Biol.64, 651–662. doi: 10.1093/sysbio/syv023

  • CrossRef
  • Google Scholar

• 64

  PeaseJ. B.BrownJ. W.WalkerJ. F.HinchliffC. E.SmithS. A. (2018). Quartet Sampling distinguishes lack of support from conflicting support in the green plant tree of life. Am. J. Bot.105, 385–403. doi: 10.1002/ajb2.1016

  • CrossRef
  • Google Scholar

• 65

  PriceJ. P.ClagueD. A. (2002). How old is the Hawaiian biota? Geology and phylogeny suggest recent divergence. Proc. R. Soc. Lond. [Biol.]269, 2429–2435. doi: 10.1098/rspb.2002.2175

  • CrossRef
  • Google Scholar

• 66

  PriceJ. P.WagnerW. L. (2004). Speciation in Hawaiian angiosperm lineages: cause, consequence, and mode. Evolution58, 2185–2200. doi: 10.1111/j.0014-3820.2004.tb01597.x

  • CrossRef
  • Google Scholar

• 67

  PriceJ. P.WagnerW. L. (2018). Origins of the Hawaiian flora: phylogenies and biogeography reveal patterns of long-distance dispersal. J. Syst. Evol.56, 600–620. doi: 10.1111/jse.12465

  • CrossRef
  • Google Scholar

• 68

  ReeR. H.HippA. L. (2015). “Inferring phylogenetic history from restriction site associated DNA (RADseq),” in Next-generation sequencing in plant systematics. Eds. HörandlE.AppelhansM. S. (Königstein, Germany: Koeltz Scientific Books), 181–204.

  • Google Scholar

• 69

  RiversD. M.DarwellC. T.AlthoffD. M. (2016). Phylogenetic analysis of RAD-seq data: examining the influence of gene genealogy conflict on analysis of concatenated data. Cladistics32, 672–681. doi: 10.1111/cla.12149

  • CrossRef
  • Google Scholar

• 70

  RoderickG. K. (1997). Herbivorous insects and the Hawaiian silversword alliance: coevolution or cospeciation? Available at: http://scholarspace.manoa.hawaii.edu/handle/10125/3219 (Accessed March 14, 2019).

  • Google Scholar

• 71

  RoyT.ColeL. W.ChangT.-H.LindqvistC. (2015). Untangling reticulate evolutionary relationships among New World and Hawaiian mints (Stachydeae, Lamiaceae). Mol. Phylogenet. Evol.89, 46–62. doi: 10.1016/j.ympev.2015.03.023

  • CrossRef
  • Google Scholar

• 72

  RubinB. E. R.ReeR. H.MoreauC. S. (2012). Inferring phylogenies from RAD sequence data. PLoS One7, e33394. doi: 10.1371/journal.pone.0033394

  • CrossRef
  • Google Scholar

• 73

  SalichosL.RokasA. (2013). Inferring ancient divergences requires genes with strong phylogenetic signals. Nature497, 327–331. doi: 10.1038/nature12130

  • CrossRef
  • Google Scholar

• 74

  SalichosL.StamatakisA.RokasA. (2014). Novel information theory-based measures for quantifying incongruence among phylogenetic trees. Mol. Biol. Evol.31, 1261–1271. doi: 10.1093/molbev/msu061

  • CrossRef
  • Google Scholar

• 75

  SeoT.-K. (2008). Calculating bootstrap probabilities of phylogeny using multilocus sequence data. Mol. Biol. Evol.25, 960–971. doi: 10.1093/molbev/msn043

  • CrossRef
  • Google Scholar

• 76

  ShenX.-X.HittingerC. T.RokasA. (2017). Contentious relationships in phylogenomic studies can be driven by a handful of genes. Nat. Ecol. Evol.1, 0126. doi: 10.1038/s41559-017-0126

  • CrossRef
  • Google Scholar

• 77

  Solís-LemusC.YangM.AnéC. (2016). Inconsistency of species tree methods under gene flow. Syst. Biol.65, 843–851. doi: 10.1093/sysbio/syw030

  • CrossRef
  • Google Scholar

• 78

  St. JohnH. (1944). Diagnoses of Hawaiian species of Pelea (Rutaceae)—Hawaiian plant studies, 13. Lloydia7, 265–274.

  • Google Scholar

• 79

  StamatakisA. (2014). RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics30, 1312–1313. doi: 10.1093/bioinformatics/btu033

  • CrossRef
  • Google Scholar

• 80

  StegemannS.KeutheM.GreinerS.BockR. (2012). Horizontal transfer of chloroplast genomes between plant species. PNAS109, 2434–2438. doi: 10.1073/pnas.1114076109

  • CrossRef
  • Google Scholar

• 81

  StoneB. C. (1969). The genus Pelea A. Gray. (Rutaceae, Evodiinae.) A taxonomic monograph. 3rd ed. Stuttgart, Germany: J. Cramer.

  • Google Scholar

• 82

  StoneB. C.WagnerW. L.HerbstD. R. (1999). “Rutaceae,” in Manual of the flowering plants of Hawaiʻi, Revised Edition. Eds. WagnerW. L.SohmerS. H. (Honolulu, HI, USA: University of Hawaii Press and Bishop Museum Press), 1174–1216.

  • Google Scholar

• 83

  StuessyT. F.TakayamaK.López-SepúlvedaP.CrawfordD. J. (2014). Interpretation of patterns of genetic variation in endemic plant species of oceanic islands: genetic variation in island plants. Bot. J. Lin. Soc.174, 276–288. doi: 10.1111/boj.12088

  • CrossRef
  • Google Scholar

• 84

  SwoffordD. L. (2002). Phylogenetic analysis using parsimony (* and other methods). Version 4.0b10. Sunderland, MA, USA: Sinauer Associates.

  • Google Scholar

• 85

  SwoffordD. L. (2018) PAUP:* Phylogenetic analysis using parsimony (and other methods). version 4a165. Available at: https://people.sc.fsu.edu/∼dswofford/paut_test/.

  • Google Scholar

• 86

  TwyfordA. D.EnnosR. A. (2012). Next-generation hybridization and introgression. Heredity108, 179–189. doi: 10.1038/hdy.2011.68

  • CrossRef
  • Google Scholar

• 87

  WagnerC. E.KellerI.WittwerS.SelzO. M.MwaikoS.GreuterL.et al. (2013). Genome-wide RAD sequence data provide unprecedented resolution of species boundaries and relationships in the Lake Victoria cichlid adaptive radiation. Mol. Ecol.22, 787–798. doi: 10.1111/mec.12023

  • CrossRef
  • Google Scholar

• 88

  WoodK. R.AppelhansM. S.WagnerW. L. (2016). Melicope oppenheimeri, section Pelea (Rutaceae), a new species from West Maui, Hawaiian Islands: with notes on its ecology, conservation, and phylogenetic placement. PhytoKeys69, 51–64. doi: 10.3897/phytokeys.69.8844

  • CrossRef
  • Google Scholar

• 89

  WoodK. R.AppelhansM. S.WagnerW. L. (2017). Melicope stonei, section Pelea (Rutaceae), a new species from Kaua’i, Hawaiian Islands: with notes on its distribution, ecology, conservation status, and phylogenetic placement. PhytoKeys83, 119–132. doi: 10.3897/phytokeys.83.13442

  • CrossRef
  • Google Scholar

• 90

  WuZ.-Y.MonroA. K.MilneR. I.WangH.YiT.-S.LiuJ.et al. (2013). Molecular phylogeny of the nettle family (Urticaceae) inferred from multiple loci of three genomes and extensive generic sampling. Mol. Phylogenet. Evol.69, 814–827. doi: 10.1016/j.ympev.2013.06.022

  • CrossRef
  • Google Scholar