Directionally biased habitat shifts and biogeographically informative cytonuclear discordance in the Hawaiian silversword alliance (Compositae)

PREMISE OF THE STUDY: Expanded phylogenetic analyses of the endemic Hawaiian silversword alliance (Argyroxiphium, Dubautia, Wilkesia; Compositae) were undertaken to obtain enhanced resolution of habitat shifts, interisland dispersal, evolutionarily relevant hybridization, and clade diversity in this prominent example of adaptive radiation. METHODS: Samples spanning the geographic and ecological distributions of all recognized taxa were included in phylogenetic and biogeographic analyses of nuclear ribosomal DNA (nrDNA) and cpDNA sequences. Bayesian model testing approaches were used to model habitat evolution and the evolution of nuclear chromosomal arrangements while accounting for phylogenetic uncertainty. KEY RESULTS: Cytonuclear discordance detected previously appears to reflect chloroplast capture, at least in part, with nrDNA trees being largely congruent with nuclear chromosomal structural data and fine-scale taxonomy. Comparison of biogeographic histories estimated from the posterior distributions of nrDNA and cpDNA trees, including inferred chloroplast-capture events, provides additional resolution of dispersal history, including a back-dispersal to Maui Nui from Hawai'i. A newly resolved major nrDNA clade of endemic Kaua'i taxa that mostly were described as new-to-science since the 1980s strengthens the earlier hypothesis that diversification on Kaua'i has not waned since the island began to decline in area through subsidence and erosion. A bias in habitat shifts was estimated, with transitions from dry-to-mesic or -wet and from wet-to-mesic or -bog habitats dominating diversification of the silversword alliance from a dry-exapted tarweed ancestor. CONCLUSIONS: The habitat-transition biases estimated here may indicate more limited pathways of ecological evolution than proposed previously for an adaptive radiation involving such major ecological shifts.Funding provided by: National Tropical Botanical GardenCrossref Funder Registry ID: http://dx.doi.org/10.13039/100011856Award Number: Funding provided by: National Science FoundationCrossref Funder Registry ID: http://dx.doi.org/10.13039/100000001Award Number:Population sampling—All 45 minimum-rank taxa (in 33 species) of Argyroxiphium, Dubautia, and Wilkesia recognized in the Flora of the Hawaiian Islands (Wagner et al., 2005) were sampled from all islands and most volcanoes of known occurrence for sequences of nrDNA (ITS and ETS) and cpDNA (3'-ndhF gene, 5'-trnK intron, rpl16 intron, and psbA-trnH intergenic spacer). The likely extinct A. virescens was represented by a putative natural hybrid (A. "Pu'u 'Alaea") with A. sandwicense subsp. macrocephalum (Carr and Medeiros, 1998). Dubautia kenwoodii, also probably extinct, was sampled from the type specimen. Individuals representing extensive morphological diversity in Dubautia from the summit bog of Wai'ale'ale, Kaua'i were sampled for ITS and ETS sequences to verify their identities, including parentage of putative hybrids (see "Assessment of natural hybridization" below). All DNAs examined earlier for cpDNA restriction sites and ITS sequences by Baldwin et al. (1990) and Baldwin and Robichaux (1995) were included among the samples. North American tarweeds of subtribe Madiinae, including the "Madia" lineage sensu Baldwin (2003) to which the silversword alliance belongs, and the sister group to Madiinae, Arnica (Baldwin and Wessa, 2000), were also sampled. A representative of tribe Madieae sensu Baldwin (Baldwin and Wessa, 2000) from outside the Arnica–Madiinae clade, Hulsea algida, served as the outgroup (see Appendix 1 for voucher information). DNA sequencing—Standard DNA extraction and polymerase chain reaction (PCR) methods were used to amplify nrDNA and cpDNA regions, as detailed in Appendix S1 (see Supplemental Data with this article). The 3' end of the nrDNA ETS, upstream of the 18S gene, was amplified and sequenced using primers ETS-Hel-1 and 18S-ETS (Baldwin and Markos, 1998). The four cpDNA regions were amplified and sequenced using the following primer pairs: trnK-3914F (Johnson and Soltis, 1994) and 884R (Panero and Crozier, 2003) for the 5'trnK intron and 5' end of matK, 1587MADIA and 607R (Panero and Crozier, 2003) for the 3'end of ndhF and ndhF–ycf1 intergenic spacer, psbAF and trnHR (McGlaughlin and Friar, 2011) for the psbA–trnH intergenic spacer, and F71 (Jordan et al., 1996) and R1516 (Kelchner and Clark, 1997) for the rpl16 intron. ITS amplification products from 10 samples across Dubautia (chosen based on nucleotide polymorphism at phylogenetically informative sites or prior phylogenetic evidence for a hybrid history of the taxon) were cloned using the zero blunt TOPO cloning kit (Thermo Fisher Scientific, Waltham, Massachusetts, USA), following the manufacturer's protocol. For these samples, i.e., BGB 667 (D. imbricata subsp. imbricata), BGB 1386 (D. kalalauensis), Carr 1047 (D. knudsenii subsp. knudsenii), BGB 675 (D. latifolia), Carr 1044 (D. microcephala), BGB 668 (D. pauciflorula), and BGB 530, 1184, 1213, 1267 (both subspecies of D. scabra), products were pooled prior to cloning from three replicate PCRs per sample, with only 20 rounds of cycling per reaction to limit possible PCR drift or selection (see Wagner et al., 1994), in an attempt to detect any variation across nrDNA repeats that might reflect a history of hybridization. Sanger sequencing of both DNA strands for all PCR products was performed at the UC Berkeley DNA Sequencing Facility (Barker Hall). For sequencing of the ITS region, ITS5 (White et al., 1990) was used instead of ITS-I (Urbatsch et al., 2000). Opposing DNA strands for each sample were assembled and reconciled using Geneious 6.1.8 (see Data Availability for GenBank accession numbers). Phylogenetic and chromosomal analyses—DNA sequences of each gene region were aligned separately using MAFFT version 7.017, as implemented in Geneious 6.1.8 with default parameters, and refined manually using the similarity criterion of Simmons (2004). Resulting sequence matrices were concatenated into an nrDNA (ETS+ITS) matrix and a four-gene cpDNA matrix after preliminary phylogenetic analyses of separate gene regions (not shown) indicated lack of conflict between ETS and ITS trees or among trees of each of the four cpDNA regions based on comparisons of clade support from bootstrap values and posterior probabilities. Prior to phylogenetic analysis of each concatenated matrix, samples assignable to the same minimum-rank taxon of Carr and colleagues (see Wagner et al., 2005) were merged into a common operational taxonomic unit (OTU) if they were identical at all sites or potentially so with resolution of sequence polymorphisms (hereafter called identical). The cpDNA matrix represented 111 OTUs, with a combined length of 3410 bp plus 28 insertions/deletions (indels) treated as additional characters using simple indel coding (Simmons and Ochoterena, 2000); the nrDNA matrix included 98 OTUs, with a combined length of 1307 bp. For phylogenetic analyses, the combined spacers (ETS, ITS-1, and ITS-2) and the 5.8S gene were treated as two distinct partitions in the nrDNA matrix and the four plastome regions were treated as a single partition in the cpDNA matrix based on results using Partition Finder version 1.1.1 (Lanfear et al., 2012). The best-fitting nucleotide substitution model was chosen by ModelFinder (Kalyaanamoorthy et al., 2017), using the Bayesian information criterion. The GTR substitution model (Tavare, 1987; Rodriguez et al., 1990) with rate variation across sites modeled as a discretized gamma distribution (Yang, 1994) was applied to the ETS, ITS-1, and ITS-2. For the 5.8S gene we used a TrN substitution model (Tamura and Nei, 1993) with a proportion of invariable sites. For the cpDNA matrix we applied the GTR substitution model uniformly to all nucleotide sites with rate variation across sites modeled as a discretized gamma distribution. For the cpDNA indels we used the F81 substitution model (Felsenstein, 1981). Maximum likelihood phylogenetic analyses were conducted using IQ-Tree 1.6.12 for MacOSX (Nguyen et al., 2015). Branch support was obtained with ultrafast (UF) nonparametric bootstrapping (1000 replicates), and a Shimodaira-Hasegawa (SH)-like approximate likelihood ratio test (aLRT; 1000 replicates) (Guindon et al., 2010). For Bayesian inference of phylogeny, Markov chain Monte Carlo (MCMC; Metropolis et al., 1953; Hastings, 1970) analyses were run in RevBayes (Höhna et al., 2016) for 40,000 iterations, with each iteration consisting of 944 independent Metropolis-Hastings steps. Trees were sampled every 10 iterations, and the first 25% of all samples were discarded as burnin. The maximum a posteriori (MAP) tree was then calculated from the 3000 sampled trees. Secondary calibrations for divergence time estimation were obtained from Landis et al.'s (2018) re-analysis of a fossil-calibrated cpDNA dataset for Compositae (Barreda et al., 2015), with inclusion of sequences from additional taxa of subtribe Madiinae and close relatives (see Data Availability for information on accessing sequence matrices and trees). The history of nuclear chromosomal arrangements (Carr and Kyhos, 1986; Carr, 2003) was estimated on both the cpDNA and nrDNA trees, as done earlier for a subset of the samples and sequences studied here (Baldwin, 1997, 2003), to evaluate clade incongruence between cpDNA and nrDNA trees based on the relatively low likelihood of homoplasy in chromosomal rearrangements compared to nucleotide substitutions in general (e.g., O'Grady et al., 2001). The evolution of nuclear chromosomal arrangements was modeled using a seven-state model, where each state represented a different chromosomal arrangement. The seven chromosomal arrangements were characterized by Carr and Kyhos (1986) and Carr (2003a) as Dubautia genomes 1–5 and Argyroxiphium genomes 1 and 2. Dubautia genome 2 was treated as structurally identical to the Wilkesia genome of Carr and Kyhos (1986) and Carr (2003a), consistent with available cytogenetic data. The rate of chromosomal translocations, t, was estimated using an exponential prior with the mean rate of one translocation over the length of the phylogeny. The rate of transitioning between each of the chromosomal arrangements was set to t/n, where n is the number of rearrangements differentiating the two genomes according to Carr and Kyhos (1986) and Carr (2003a). The number of chromosomal rearrangements differentiating Argyroxiphium genomes 1 and 2 from some of the other genomes is uncertain (Carr and Kyhos, 1986). For those transitions, we estimated each rate using a mixture distribution in which the rate could be t, t/2, t/3, or t/4, with equal prior probability. Ancestral chromosomal arrangements and parameters were estimated using MCMC analyses in RevBayes (Höhna et al., 2016). Phylogenetic uncertainty was accounted for by sampling across 3000 trees of the posterior distributions estimated as described above. The MCMC was run for 5000 iterations, with each iteration consisting of eight separate Metropolis-Hastings steps. The first 10% of samples were discarded as burnin. Assessment of natural hybridization—Previous documentation of a wide diversity of natural hybrids in the silversword alliance (Carr and Kyhos, 1981; Carr, 2003b) warranted instances of sampling to detect unconfirmed hybridization, especially in the summit bog of Wai'ale'ale, Kaua'i, where four taxa of Dubautia (D. imbricata subsp. acronaea, D. laxa subsp. hirsuta, D. paleata, and D. waialealae) are sympatric and evidently hybridize (Carr, 1985, 1999). Direct ETS and ITS sequences of putative hybrids were examined for additivity of mutations diagnostic for each of the four sympatric taxa, i.e., for double-peaks in electropherograms corresponding to nucleotide states that differentiate taxa. Pollen stainability of flowering hybrids was also assessed by light microscopy from a sample of 200+ grains per individual as an estimate of fertility to consider intrinsic potential for gene flow and possible chromosomal differentiation between parental taxa of unknown genomic arrangement. Pollen from late-stage floral buds of dried specimens was stained overnight with cotton blue in lactophenol prior to microscopic examination . Estimation of biogeographic history—Biogeographic history was estimated for the nrDNA and cpDNA trees at the scale of island or, for the prehistoric island Maui Nui (including modern-day Lana'i, Maui, and Moloka'i), island group using the Dispersal-Extinction-Cladogenesis (DEC) model (Ree and Smith, 2008), as implemented in RevBayes (Höhna et al., 2016). Probabilities for biogeographic states were obtained from sampling across 3000 trees of the posterior distribution for each of the two sets of BI trees, from nrDNA and cpDNA sequences. All non-Hawaiian taxa were coded as North American. OTUs representing samples from multiple geographic areas were coded so that the tip state probabilities were equal among that OTU's observed geographic areas. Ancestral ranges were estimated during an MCMC analysis that was run for 5000 iterations, with each iteration consisting of 15 separate Metropolis-Hastings steps. The first 10% of samples were discarded as burnin. In addition, isolation-by-distance (Wright, 1943) analyses were performed for the nrDNA and cpDNA datasets using Mantel (1967) tests to compare log transformed phylogenetic distance in molecular branch lengths and log transformed geographic distance among all samples. For each of the two datasets, 1000 Mantel tests were performed, with sampling across the posterior distribution of trees to account for phylogenetic uncertainty. Mantel tests were performed using the ade4 package (Dray and Dufour, 2007) in the R programming language (R Core Team, 2017). Estimation of growth-form and habitat transitions—Each Hawaiian population sampled for the phylogenetic analyses was assigned to one of five growth-form states — cushion-plant, liana, rosette-plant, shrub, or tree — following the habit designations of Blonder et al. (2016), which align closely with those of Robichaux et al. (1990). Each Hawaiian population also was assigned to one of three habitat states — dry, wet, or bog — based on the ecological designations in Table 1 of Robichaux et al. (1990) and Blonder et al. (2016), with "cinder and lava" treated as "dry" except for Dubautia scabra, which likely accesses abundant water at depth (Robichaux, 1984) and was scored only as "wet." Other taxa or populations not scored previously for habitat under the above classification were assigned designations consistent with that system, based on plant associations, local climate, and published descriptions. As an alternative approach to ecological categorization, each Hawaiian population was assigned to one of four habitat states — dry, mesic, wet, or bog — with the bog designation unchanged and the other three designations based on moisture regimes of Gagné and Cuddihy (1990) and mean annual rainfall estimates for collection localities (Giambelluca et al., 2013; http://rainfall.geography.hawaii.edu/). For taxa from Lana'i high montane sites, where fog precipitation is extensive relative to rainfall, habitats were scored as "wet," as reflected by field records and moisture zonation estimates for the Hawaiian Islands (Price et al., 2012; https://hub.arcgis.com/datasets/HiStateGIS::moisture-zones). For North American taxa, habitat designations for both sets of habitat categorizations were based on average annual precipitation estimates from PRISM data (period 1981–2010; http://www.prism.oregonstate.edu/normals/). Growth-form and habitat evolution under both systems of habitat classification were each estimated for the nrDNA trees in RevBayes (Höhna et al., 2016). Reversible‐jump MCMC (Green, 1995) was used to explore a set of continuous‐time Markov models of character evolution and to infer ancestral states while marginalizing over the posterior distribution of trees. The growth-form analyses were restricted to the silversword alliance because of lack of overlap in habit between North American and Hawaiian members of the "Madia" lineage (Baldwin 2003). The set of models of character evolution considered included all models in which each transition rate was either estimated using an exponential prior with the mean rate of one transition over the length of the phylogeny or fixed to zero. The reversible‐jump MCMC sampled from the models in proportion to their posterior probability. This approach enabled model‐fit comparisons through Bayes factors (Kass and Raftery, 1995) and provided the opportunity to account for both phylogenetic and model uncertainty by making model‐averaged ancestral state and parameter estimates (Madigan and Raftery, 1994; Kass and Raftery, 1995; Huelsenbeck et al., 2004; Freyman and Höhna, 2018). The MCMC was run for 5000 iterations, with each iteration consisting of 360 separate Metropolis-Hastings steps. The first 10% of samples were discarded as burnin