Repeated range expansion and niche shift in a volcanic hotspot archipelago: radiation of Hawaiian Euphorbia (Euphorbiaceae)

Aim The taxon cycle hypothesis describes the cyclic movement of taxa during range expansion and contraction, accompanied by an evolutionary shift from open and often coastal vegetation to closed, and often inland forest vegetation in island systems. The Hawaiian Archipelago is an ideal system to test this hypothesis given the linear fashion of island formation and a relatively well-understood geological history. Location Hawaiian Islands. Methods We sampled 153 individuals in 15 of the 16 native species of Hawaiian Euphorbia section Anisophyllum on six major Hawaiian Islands, plus 11 New World close relatives, to elucidate the biogeographic movement of this lineage along the Hawaiian island chain. We used a concatenated chloroplast DNA data set of more than eight kilobases in aligned length and applied maximum likelihood and Bayesian inference for phylogenetic reconstruction. Connectivity among islands and habitat types was estimated using BayesTraits. Age and phylogeographic patterns were co-estimated using BEAST. In addition, we used nuclear ribosomal ITS and the low-copy genes LEAFY and G3pdhC to investigate the reticulate relationships within this radiation. Results We estimate that Hawaiian Euphorbia first arrived on Kauai or Niihau ca. 5 million years ago and subsequently diverged into 16 species on all major Hawaiian Islands. During this process Euphorbia dispersed from older to younger islands in a stepping-stone fashion through open, dispersal-prone habitats. Taxa that occupy closed vegetation on Kauai and Oahu evolved in situ from open vegetation taxa of the same island. Consequently, widespread species tend to occupy habitats with open vegetation, whereas single island endemic species predominantly occur in habitats with closed canopy and are only found on the two oldest islands of Kauai and Oahu. Main conclusions The spatial and temporal patterns of dispersal and range shifts in Hawaiian Euphorbia support an intra-volcanic-archipelago version of the taxon cycle hypothesis.

Connectivity among islands and habitat types was estimated using BayesTraits. Age and phylogeographic patterns were co-estimated using BEAST. In addition, we used nuclear ribosomal ITS and the low-copy genes LEAFY and G3pdhC to investigate the reticulate relationships within this radiation.
The taxon cycle hypothesis describes the cyclic pattern of range expansion, contraction, and habitat shifts that occur in island archipelagos (Wilson, 1959(Wilson, , 1961Ricklefs & Cox, 1978;Ricklefs & Bermingham, 2002). According to this hypothesis, taxa that colonize oceanic islands have inherently high dispersal abilities and are typically adapted to disturbed habitats. After the initial colonization, some populations expand to habitats that have closed forest canopies and are less susceptible to disturbance. During this transition, selection favors survival and competitive ability under closed vegetation rather than the dispersal ability that facilitated the initial colonization. Species that remain in the lowland, drier and more disturbed habitats retain their dispersal ability and continue to colonize new islands, whereas species moved to forests with closed vegetation are increasingly specialized and isolated, with a tendency to become narrow endemics.
The taxon cycle hypothesis was originally proposed to explain patterns observed among the ants of Melanesia (Wilson, 1959(Wilson, , 1961. It has been further developed in birds in the Lesser Antilles (Ricklefs & Cox, 1972;Ricklefs & Cox, 1978;Ricklefs & Bermingham, 1999), and supported by genetic evidence in birds of the Indo-Pacific (Jønsson et al., 2014) and in ants across the Old World (Economo et al., 2015). Volcanic hot spot island archipelagos, such as the Hawaiian Islands, present an ideal system to test the taxon cycle hypothesis. The recent, relatively well-understood geological history makes the dynamics of taxon movement more tractable than in older, more complicated continental island systems.  (Dorsey et al., 2013) and is outside the scope of this study. The 16 Hawaiian Euphorbia species occur in all major island habitats, from coastal strand to dry forests, wet forests and bogs, and range in habit from subshrubs and shrubs to trees ten meters tall ( Fig.   1). Four of these species have two or more recognized varieties. Ten species are endemic to a single major island, while the remainder are known from two or more major islands (Table 1).
Previous phylogenetic studies suggest that Hawaiian Euphorbia originated through allopolyploidy, with their closest relatives being small herbs occurring in dry, warm, and cinerascens Engelm ( Fig. 1f; Yang & Berry, 2011). Given the overlapping distribution of the putative sister species in North America, the allopolyploidy event likely happened before longdistance dispersal to Hawaii.
In this study, we use multiple chloroplast and nuclear markers to reconstruct the history of radiation in Hawaiian Euphorbia. We adopt a modified classification of the taxon cycle stages (Wilson, 1959;Ricklefs & Cox, 1978), with stage I consisting of species that are undifferentiated and that occupy multiple islands ("widespread"), stage II consisting of widespread species with considerable intraspecific differentiation and recognized varieties, and stage III consisting of single island endemics. If Hawaiian Euphorbia fits the taxon cycle hypothesis, we predict that: (1) Hawaiian Euphorbia colonized younger islands successively from older islands.
(2) Habitat connectivity is high among open vegetation and low among closed vegetation.
(  (Price & Elliott-Fisk, 2004). Of the 153 DNA accessions, 125 were obtained from the Hawaiian Plant DNA Library (Morden et al., 1996;Randell & Morden, 1999), complemented by 18 additional samples newly collected from the field or cultivated sources (Appendix S1). Forty-three DNA accessions in the DNA Library were collected by M.J. Sporck-Koehler and L. Sack, accompanied by some of Hawaii's most experienced field botanists (see Acknowledgements) as part of an ecophysiological study of the Hawaiian Euphorbia (Sporck, 2011); permits for many of the species were limited to less than ten leaves per plant, with vouchers not permitted for State and Federally listed endangered or very rare and extremely vulnerable taxa. We described in Appendix S1 alternative voucher specimens representing the same population and/or additional locality information for source populations. Genomic DNA extraction, plus PCR amplification and sequencing of both ITS and cpDNA followed the protocols in Yang & Berry (2011). A total of seven chloroplast (cpDNA) noncoding regions were sequenced: rpl14-rpl36 spacer, psbB-psbH spacer, atpI-atpH spacer, psbD-trnT spacer, trnH-psbA spacer, rpl16 intron, and the trnL-F region. For the ITS region, sequences with continuous superimposed peaks were excluded. Two of these excluded PCR products, E.
celastroides var. kaenana 5840 and E. kuwaleana 5700 were cloned following the protocol of Yang & Berry (2011) to evaluate allelic variation. The second intron of the nuclear low-copy gene LEAFY and intron of glyceraldehyde 3-phosphate dehydrogenase subunit C (G3pdhC) were PCR amplified and cloned following the protocol in Yang & Berry (2011), except that at least 24 clones from each PCR product were sequenced. Copy-specific primer pairs were designed for both markers and at least eight clones were sequenced from each copy-specific PCR reaction.
See Supplementary Methods in Appendix S2 for primers used, and PCR and cloning procedures using copy-specific primers. complemented before concatenating all seven cpDNA regions into the first character set of the cpDNA matrix. Indels were scored following the simple gap-coding criterion (Simmons & Ochoterena, 2000) in SeqState v1.4.1 (Müller, 2006) and were treated as the second character set of the cpDNA matrix.
Trees were sampled every 1,000 generations. Parameters were unlinked between the two partitions except tree topologies. The binary indels were subject to "rates=gamma". In addition, a branch length prior "brlenspr=unconstrained:exponential(100.0)" was applied to the nucleotide partition to prevent unrealistically long branches (Marshall, 2010). Diagnostic parameters were visually examined in the program Tracer v1.5 (Rambaut & Drummond, 2007) to verify stationary status. Trees sampled from the first 1 million generations were discarded as burn-in, and the remaining 18,002 trees were used to compute the majority rule consensus (MCC) tree.
Maximum likelihood (ML) analysis was carried out using RAxML v7.4.2 (Stamatakis, 2006), partitioning nucleotides vs. indels. The nucleotide substitution model was set to GTR+γ, 500 rapid bootstrap replicates were performed, followed by a thorough search for the best tree.
M o l e c u l a r d a t i n g The Hawaiian island chain was formed by the Pacific plate moving northwestward over a fixed hot spot (Carson & Clague, 1995). We assumed that a new island was colonized soon after it emerged (Fleischer et al., 1998), and that given the extremely small colonizing population, deep divergence from ancestral polymorphisms in the colonizing population is highly unlikely. We cross-validated these two assumptions with a preliminary analysis estimating the stem ages of Maui Nui and Hawaiian clades by constraining the stem age of the oldest Oahu clade on the cpDNA data set with the time of full development of the Waianae Mountains (a normal prior with mean 3.86 my, standard deviation 0.089 my; Sherrod et al., 2007;Lerner et al., 2011).
A final analysis was carried out by applying the following age constraints to the cpDNA data set: (1) 3.86 ± 0.089 my for the stem age of the Oahu clade; and (2) 2.14 ± 0.117 my for the stem age of Maui Nui clades (the age of Penguin Bank, which formed the past land connection between Oahu and Maui Nui; Carson & Clague, 1995;Price & Elliott-Fisk, 2004;Sherrod et al., 2007;Lerner et al., 2011). The analysis was performed in BEAST v1.7.4 (Drummond et al., 2012), using the concatenated cpDNA data set without coding indels. The substitution model HKY+I+γ was applied as selected by jModeltest v0.1.1 (Posada, 2008), with an uncorrelated lognormal relaxed clock and a pure-birth Yule model. Four independent runs of 60 million generations were carried out, sampling every 10,000 generations starting from a random starting tree. Convergence diagnose parameters were visualized in Tracer, and trees sampled from the first six million generations were discarded as burn-in. A MCC tree was calculated in Discrete phylogeographic analysis (Lemey et al., 2009) was used to reconstruct the pattern of dispersal along the island chain from the cpDNA data set. Phylogeographic analysis was carried out in BEAST using two independent continuous-time Markov chains (CTMCs) by manually editing the xml file generated by BEAUti from the previous molecular dating analysis following Lemey et al. (2009). Most recent common ancestor (MRCA) of all Hawaiian accessions was set to Kauai according to molecular dating results. Convergence diagnostic parameters were visualized in Tracer, and the first six million generations were discarded as burn-in.
We categorized coastal strand, scrub, and dry forest as "open vegetation". These vegetation types are either fully exposed or have relatively open forest canopy coverage, and are generally relatively low in elevation, although the upper elevation limit of lowland dry forest varies from 150 to 1,500 meters depending on the island and the direction of the slope, and E. olowaluana occurs in montane dry forests from 700 to 2,800 m in elevation on Hawaii (Gagné & Cuddihy, 1990;Koutnik & Huft, 1990 Three independent final rjMCMC runs were carried out each with 50 million generations, sampling every 1,000 generations, with rate deviation set to 0.5 to optimize the acceptance rate, and a uniform hyperprior between 0 and 5 to seed the exponential prior. Results were visualized in Tracer, and rates sampled from the first 10 million generations were discarded as burn-in. We obtained sequences of all seven chloroplast non-coding regions from each of the 164 DNA accessions included in this study (Table S2. Both the low-copy nuclear genes LEAFY and G3pdhC showed increased copy numbers among Hawaiian taxa compared to outgroups, but the resolution within each copy was low ( Fig.   S2.2b). Four copies of LEAFY were recovered, of them only one copy was detected from the known outgroups. Similarly, six copies of G3pdhC were detected in Hawaiian Euphorbia, among which three had a clear association with known outgroup taxa.
We used cpDNA only for dating and phylogeographic analyses to track dispersal via seed dispersal or vegetative fragments. Using island age for molecular dating can potentially be biased by delayed arrival long after island formation, multiple dispersal events, local extinction, and ancestral polymorphism. Another consideration is that at the time Kauai formed ca. 5 mya, the adjacent island of Niihau was of similar size and prominence (Price & Clague, 2002). To crossvalidate our assumptions and their potential caveats, a preliminary analysis was carried out constraining the stem age of Oahu clade 1, the most diverse and well supported Oahu clade, by the date at which the Waianae Mountains of Oahu formed ( Fig. 2; 3.86 ± 0.089 my; Sherrod et al., 2007;Lerner et al., 2011 Sherrod et al., 2007;Lerner et al., 2011), whereas Maui Nui clade 3 is a much smaller and younger clade including only a single coastal taxon, and probably represented a more recent dispersal event. As for the Hawaii clade, both its stem age (1.9 my; 1.1-2.9 my) and crown age (1.3 my; 0.7-2.1 my) were much older than the age of the island of Hawaii (≈0.59 my; Sherrod et al., 2007;Lerner et al., 2011). Both taxa in the Hawaii clade, E. multiformis var. microphylla and E. olowaluana, also occur on Maui Nui (Koutnik, 1987). It is likely that the "Hawaii clade" diverged on Maui Nui before dispersing to Hawaii.  The number of overall species per island was highest in Oahu, and decreased towards both older and younger islands, showing a humped trend (Fig. 5a). Stage I species were most numerous on Maui Nui; stage II species were approximately equally numerous on all four islands; and stage III species were only found on Kauai and Oahu, the two oldest islands, and absent from the two younger island groups. Alternatively, breaking down infraspecific taxa by endemic to one island vs. on two or more islands (Fig. 5c) showed a similar trend, with the number of endemic taxa highest on Kauai, and the number of widespread taxa peaking on Maui Nui.
The species-habitat plot (Fig. 5b)  Our analyses supported radiation from a single colonization through allopolyploidy. Both of the nuclear low-copy genes examined in this study, LEAFY and G3pdhC had increased copy numbers in Hawaiian Euphorbia compared to outgroups. Two of the four copies detected in LEAFY, and three of the six copies detected in G3pdhC were not associated with known outgroups previously identified using ITS and chloroplast markers, consistent with previous results from cloning another nuclear low-copy gene, EMB2765 (Yang & Berry, 2011).
Additional support for the allopolyploid origin of Hawaiian Euphorbia comes from chromosome numbers: while the outgroup species E. cinerascens and E. stictospora have chromosome counts of 2n=12 and 2n=32 respectively (Urbatsch et al., 1975) are rounded in shape and fold upwards and is specialized in coastal beach habitats ( Fig. 1c; Koutnik, 1987). We included multiple Kauai, Oahu, and Maui Nui accessions of E. degeneri, and they are placed by cpDNA in multiple, separate clades within Oahu clade 1, and by ITS in a polytomy consisting mostly of open Oahu and younger island accessions. Given its distinctiveness in both morphology and habitat, its polyphyly in the cpDNA tree, and the lack of resolution in the ITS tree, what we recognize morphologically and ecologically as E. degeneri is likely a result of reticulate evolution and/or rapid diversification lacking species monophyly.
A second highly polyphyletic species, E. celastroides (Fig. 1d-e), is more variable in both morphology and habitats and has eight recognized varieties (  Fig. S2.2a).
Despite being highly variable, E. celastroides is morphologically recognizable, with entire, distichous leaves that are oblong to obovate in shape (Fig. 1e). It can be distinguished from the vegetatively similar E. multiformis, also a widespread stage II species, by its erect fruits and appressed cyathial glands, as opposed to recurved fruits and protruding glands in E. multiformis (Koutnick 1985).
In addition to presumed older reticulations, we also found evidence for more recent hybridization events. Euphorbia multiformis var. microphylla 5622 and 5624 were collected from the Pohakuloa Training Area of Hawaii, and they share an almost identical cpDNA haplotype with E. olowaluana accessions that came from the same area (Fig. 2). In the ITS phylogeny, however, both E. multiformis var. microphylla 5622 and 5624 are not part of a monophyletic E. olowaluana (Fig. S2.2a). Pittosporum (Pittosporaceae; Bacon et al., 2011), and Bidens (Asteraceae; Knope et al., 2012).
Together these examples caution against using single representative samples per species, or relying on cpDNA and/or ITS as the sole source for studying recent, rapid radiations. Oahu, Maui Nui, and finally Hawaii, following the "progression rule" from older to younger islands (Hennig, 1966;Funk & Wagner, 1995), with at least one dispersal event in the reverse direction through a stage I species.
Our BayesTraits analysis with outgroup ranges coded as "unknown" reconstructed the most likely origin of Hawaiian Euphorbia as Kauai (52%), followed by Oahu (29%; Table 2).

Molecular dating using ages of Oahu and Maui Nui further supported a Kauai or Niihau origin of
Hawaiian Euphorbia, given these two islands were of similar sizes 5 mya (Price & Clague, 2002). Our estimation of the stem age of Hawaiian Euphorbia at ca. 5 my is also consistent with previous molecular dating based on fossil evidence, which estimated the split between E. hirta and E. humifusa to be ca. 9 mya (Horn et al., 2014), a split much deeper than the stem of Hawaiian Euphorbia (Yang & Berry, 2011).
Following the initial establishment in Kauai, dispersal from Kauai to Oahu occurred at least once (Fig. 3), with coastal Oahu serving as a hub and stepping stone for diversification in Hawaiian Artemisia (Hobbs & Baldwin, 2013) and in flightless alpine moths in Hawaii and Maui (Medeiros & Gillespie, 2011). In Hawaiian violets, however, an nrITS phylogeny recovered a dry clade and a wet clade, each having species from multiple islands (Havran et al., 2009). Since that study relied solely on the nrITS marker in a group of polyploid species, its results may have been biased (Marcussen et al., 2012).  5a-b). In contrast, all single island endemic species (stage III) occurred on Kauai and Oahu, the two oldest islands, and most of them occur in closed habitats. No stage III species occur on Maui or Hawaii, despite their larger sizes and higher elevations than the older islands. Therefore it appears that when a young island emerges, it is first colonized by stage I and stage II species.
Stage III species only arise later through evolution of adaptation to forest understories.
Certain infraspecific taxa in Hawaiian Euphorbia are geographically and morphologically distinctive enough that it is sometimes unclear whether separate species should be recognized (see discussion in Koutnik 1985;Koutnik 1987;Koutnik & Huft, 1990). This kind of taxonomic uncertainty may reflect the process of increased isolation during the transition from stage II to stage III of a taxon cycle. We further broke down the distribution of infraspecific taxa by island and habit respectively (Fig. 5c-d). The trend is very similar to species vs. island and habitat relationships: single-island endemic taxa are most numerous on Kauai and Oahu, and absent from Hawaii, whereas the number of widespread taxa peaks on Maui Nui. Among habitat types, endemic taxa are most numerous in mesic forests, whereas widespread taxa are most numerous in dry forests.
Both dispersal ability and habitat specialization in Hawaiian Euphorbia appear to be associated with seed characters. Hawaiian Euphorbia most likely arrived from North America via tiny seeds that adhered to birds through a mucilaginous seed coat (Carlquist, 1966(Carlquist, , 1980aPrice & Wagner, 2004). A survey of mucilaginous seed coats across Euphorbia sect.
Anisophyllum (Jordan & Hayden, 1992) showed that it is present in most mainland species as well as in E. celastroides, a stage II species and one of the most widespread members of Hawaiian Euphorbia. It is absent, however, in all four stage III Hawaiian species surveyed: E.
clusiifolia, E. halemanui, E. remyi, and E. rockii. Interestingly, E. degeneri, a stage I species occurring on coastal strand of all major Hawaiian islands, lacks a mucilaginous seed coat.
Instead, it is able to float in sea water (Carlquist, 1980b), which likely explains its coastal distribution and offers an alternative dispersal mechanism besides sticking to birds. In contrast, neither E. celastroides (stage II) nor E. clusiifolia (stage III) appear to have floating seeds (Carlquist, 1966). In addition to the difference in dispersal ability between species of different taxon cycle stages, stage III species such as E. rockii and E. clusiifolia have seeds about twice as large in diameter as stage I and stage II species. Such larger, non-sticky, non-buoyant seeds may enhance seedling survival in forest understory while simultaneously reducing their dispersal ability.
By being restricted to localized habitats and having much reduced dispersal ability, stage III species are likely more vulnerable to extinction, as predicted by the taxon cycle hypothesis.
Of the six species and four additional varieties of Hawaiian Euphorbia that are federally listed as endangered, five species plus two additional varieties belong to stage III and predominantly occur in wet or mesic forest (E. deppeana, E. eleanoriae, E. halemanui, E. herbstii, E. remyi var. kauaiensis, E. remyi var. remyi, and E. rockii), with E. kuwaleana being the only lowland scrub species listed; two remaining varieties are localized populations of stage II species (E. celastroides var. kaenana and E. skottsbergii var. skottsbergii) that occur in coastal habitats.
Although the cyclic pattern of dispersal and habitat shift seen in Hawaiian Euphorbia supports aspects of the taxon cycle hypothesis, there are also fundamental differences between the taxon cycles described by Wilson (1959Wilson ( , 1961 and our intra-volcanic-archipelago version of the taxon cycle. First, Wilson described cyclic patterns at a larger scale, with repeated movement from the continent to an archipelago. In contrast, the pattern seen in Hawaiian Euphorbia occurred over smaller distances and times. Because of the highly isolated nature of the Hawaiian Archipelago, colonization from a continent is infrequent. Instead of "colonization waves" from the continent, colonization occurred in waves from older to younger islands, starting soon after a young island formed.
Secondly, in the original taxon cycle concept (Wilson, 1959(Wilson, , 1961, single island endemics formed from range contraction and local extinction from previously widespread species. In the Hawaiian Islands, however, ecological speciation occurred rapidly within islands such that populations occupying closed habitats speciated while the widespread species dispersed to the next new island, producing a high level of polyphyly among coastal taxa. Many stage II species and their infraspecific taxa are likely at various stages of moving to specialized habitats and becoming increasingly isolated. The parallel formation of new species in closed vegetation from lowland open vegetation on both Kauai and Oahu provides a promising case study for underlying mechanisms and processes driving ecological speciation. Finally, the time series of island formation and erosion in volcanic island systems adds another dimension to the dynamics of taxon cycle (Whittaker et al., 2008). This is evident from the hump-shaped curve of overall species versus island age relationship typical in volcanic island systems (Fig. 5), here showing a peak for Oahu. Stage I species lead the curve and peak on Maui Nui; stage II species had a more even distribution among major island groups; and stage III species occurred only on the two oldest islands. The total species number of Hawaiian lobelias, in contrast, peaked on Maui (Givnish et al., 2009), the second youngest major island. That pattern was likely due to an earlier arrival on a former tall island further up the island chain and/or more rapid divergence, making the Hawaiian lobelioids "saturated" with ca. 130 All DNA sequence data were deposited in GenBank (Appendix S1, to be deposited). All alignment files in NEXUS format are in Appendix S3.  Distribution of the 16 Hawaiian Euphorbia species on six major Hawaiian Islands. Habitat types are sorted from wetter habitats generally at higher elevations to lower elevation and dryer ones, whereas ages of islands are ordered left to right from older to younger (Koutnik, 1987;Koutnik & Huft, 1990;Lorence & Wagner, 1996;Morden & Gregoritza, 2005). Taxa with an "*" are federally listed as endangered.
T a b l e 2 Median transition rates between ranges inferred from BayesTraits. Rates that are nonzero in >95% of the post-burn-in values are bold. Shaded cells involve either the source or destination range as being closed vegetation.
Hawaiian Euphorbia (a-e) and their closely related North American species (   Maximum clade credibility (MCC) tree recovered from BEAST phylogeographic analysis. Node labels are mean ages, and node bars are 95% highest posterior density (HPD) intervals. Outgroups are not shown. Following each taxon name is the DNA accession number, island initials for the individual, and vegetation type for the taxon. Superimposed blue lines on the Hawaii clade indicates the most likely scenario inferred from the clade age and historical distribution of E. olowaluana on Maui. Map in the upper left corner shows the inferred dominant pattern of dispersal among islands. Transition rates among all six ranges where Hawaiian Euphorbia occurs. Thickness of arrow is proportional to the median value of the transition rate. For clarity, rates lower than 0.03 are now shown.