DOI: 10.6078/d1ft5g
title: Directionally biased habitat shifts and biogeographically informative cytonuclear discordance in the Hawaiian silversword alliance (Compositae)

description: 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. 
description: <b><i>Population sampling</i>—</b>All 45
 minimum-rank taxa (in 33 species) of <i>Argyroxiphium</i>,
 <i>Dubautia</i>, and <i>Wilkesia</i> recognized in
 the <i>Flora of the Hawaiian Islands</i> (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’-<i>ndh</i>F
 gene, 5’-<i>trn</i>K intron, <i>rpl</i>16 intron,
 and <i>psb</i>A-<i>trn</i>H intergenic spacer).
 The likely extinct <i>A. virescens</i> was represented by a
 putative natural hybrid (<i>A</i>. “Pu‘u ‘Alaea”) with
 <i>A. sandwicense</i> subsp. <i>macrocephalum</i>
 (Carr and Medeiros, 1998). <i>Dubautia kenwoodii</i>, also
 probably extinct, was sampled from the type specimen. Individuals
 representing extensive morphological diversity in
 <i>Dubautia</i> 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,
 <i>Arnica</i> (Baldwin and Wessa, 2000), were also sampled. A
 representative of tribe Madieae sensu Baldwin (Baldwin and Wessa, 2000)
 from outside the <i>Arnica</i>–Madiinae clade, <i>Hulsea
 algida</i>, served as the outgroup (see Appendix 1 for voucher
 information).  <b><i>DNA
 sequencing</i>—</b>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: <i>trn</i>K-3914F (Johnson and Soltis, 1994) and
 884R (Panero and Crozier, 2003) for the 5’<i>trn</i>K intron
 and 5’ end of <i>mat</i>K, 1587MADIA and 607R (Panero and
 Crozier, 2003) for the 3’end of <i>ndh</i>F and
 <i>ndh</i>F–<i>ycf</i>1 intergenic spacer,
 <i>psb</i>AF and <i>trn</i>HR (McGlaughlin and
 Friar, 2011) for the <i>psb</i>A–<i>trn</i>H
 intergenic spacer, and F71 (Jordan et al., 1996) and R1516 (Kelchner and
 Clark, 1997) for the <i>rpl</i>16 intron. ITS amplification products from 10
 samples across <i>Dubautia</i> (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.,
 <i>BGB 667</i> (<i>D. imbricata</i> subsp.
 <i>imbricata</i>), <i>BGB 1386</i> (<i>D.
 kalalauensis</i>), <i>Carr 1047</i> (<i>D.
 knudsenii</i> subsp. <i>knudsenii</i>), <i>BGB
 675</i> (<i>D. latifolia</i>), <i>Carr
 1044</i> (<i>D. microcephala</i>), <i>BGB
 668</i> (<i>D. pauciflorula</i>), and <i>BGB
 530</i>, <i>1184</i>, <i>1213</i>,
 <i>1267 </i>(both subspecies of <i>D. scabra</i>),
 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).             <b><i>Phylogenetic
 and chromosomal analyses</i>—</b>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).<b>  </b>             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 <i>Dubautia</i> genomes
 1–5 and <i>Argyroxiphium</i> genomes 1 and 2.
 <i>Dubautia</i> genome 2 was treated as structurally identical
 to the <i>Wilkesia</i> 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
 <i>Argyroxiphium</i> 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.
 <b>            <i>Assessment of natural
 hybridization</i>—</b>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 <i>Dubautia</i> (<i>D. imbricata</i>
 subsp. <i>acronaea</i>, <i>D. laxa</i> subsp.
 <i>hirsuta</i>, <i>D. paleata</i>, and <i>D.
 waialealae</i>) 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     .             
 <b><i>Estimation of biogeographic
 history</i>—</b>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).
             <b><i>Estimation of growth-form and
 habitat transitions</i>—</b>Each Hawaiian population sampled
 for the phylogenetic analyses was assigned to one of five growth-form
 states <b>— </b>cushion-plant, liana, rosette-plant, shrub, or
 tree <b>— </b>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
 <b>— </b>dry, wet, or bog <b>—</b> 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
 <i>Dubautia scabra</i>, 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 <b>—
 </b>dry, mesic, wet, or bog <b>— </b>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.