Plastid genomes reveal support for deep phylogenetic relationships and extensive rate variation among palms and other commelinid monocots

Taxon sampling, plant material, and deposition of vouchers

Leaf material from various botanical gardens (see Table 1) was silica gel-dried and stored at −20°C. Representatives from all five subfamilies had at least one species from each of the 15 non-Arecoid tribes and three of the 14 Arecoid tribes, following the subfamilial and tribal classifications of Asmussen et al. (2006) and Dransfield et al. (2005), respectively. Six species with previously sequenced plastomes or plastid gene sets were chosen to represent the subfamily Arecoideae, based on Baker et al. (2009, 2011): Areca vestiaria and Veitchia spiralis (tribe Areceae); Cocos nucifera, Elaeis guineensis and E. oleifera (tribe Cocoseae); and Chamaedorea seifrizii (tribe Chamaedoreeae). The focus of sampling was not within the Arecoideae, as relationships in this subfamily have been treated by Baker et al. (2011), and are the subject of a recent plastome-based investigation (Comer et al., 2015). Here, arecoid taxa were represented by ‘placeholders’, in order to focus sampling more intensively on relationships within and among the remaining subfamilies. All four genera of Dasypogonaceae are represented; Baxteria was newly sequenced for this study, while sequences from Dasypogon, Calectasia and Kingia were generated previously (Givnish et al., 2010; Barrett et al., 2013). Nearly complete sets of coding genes from previous plastome-based studies were included in monocot-wide, protein-coding matrices (see later; Givnish et al., 2010; Steele et al., 2012; Barrett et al., 2013, 2014a,b; GenBank numbers listed in respective papers).

DNA isolation, sequencing and read trimming

Silica-dried leaf material was ground in liquid nitrogen, and total DNAs were isolated using a modified CTAB protocol (Doyle & Doyle, 1987). The DNA concentration was quantified using a NanoDrop spectrophotometer (Thermo Scientific, Carlsbad, CA, USA). Isolations with concentrations > 30 ng μl⁻¹ were chosen for Illumina sequencing. The total gDNAs were shipped to Cold Spring Harbor Laboratory (Cold Spring Harbor, Woodbury, NY, USA) or to the University of Missouri-Columbia for library preparation and Illumina sequencing. The samples were sequenced over a 4-yr period (2009–2013), including single end reads of 71 or 96 bp, and paired end reads of 100 bp (Table 1). Sequences were trimmed and adaptors removed with Trimmomatic v.0.33 (Bolger et al., 2014). Trimming was completed from the 5′ and 3′ ends of each read, removing bases below PHRED 30 within a sliding window of four bases, keeping only reads/read-pairs of 51 bases or longer.

Plastome assembly and annotation

A multistep approach was taken to build plastomes, including reference guided assembly, de novo assembly and remapping of reads, to ensure high-quality assemblies. Reference-guided assemblies were conducted using Yasra (Ratan, 2009; http://www.bx.psu.edu/miller_lab), based initially on plastomes of Phoenix dactylifera (Arecaceae, Genbank GU811709) and Typha latifolia (Typhaceae; GenBank GU195652). The resulting plastomes were used as subsequent references. De novo assemblies were constructed in Velvet (Zerbino & Birney, 2008), using kmer lengths of 51–81 bp. Contigs from each assembly method were merged in Sequencher v.5.2 (GeneCodes, Ann Arbor, MI, USA). Reads were then remapped to references for each taxon to check for misassemblies or rearrangements in Geneious v.6 (http://www.geneious.com; Kearse et al., 2012) under stringent parameters (3% mismatch to reference; maximum gap size 1000 bp), and reads matching the draft reference were assembled de novo, also in Geneious under suggested parameters. Inverted repeat boundaries were determined in Sequencher, and verified by remapping reads in Geneious.

The UNIX function ‘grep’ was used to close contig gaps by searching among read files. Finished plastid genomes were annotated via BLASTX in Dogma (Wyman et al., 2004). Smaller exons (< 30 bp) were manually annotated. Lastly, primers were developed with Primer3 (Untergrasser et al., 2012) to close low-coverage gaps between contigs (for a few single end datasets, Supporting Information Fig. S1). Amplifications and Sanger sequencing were conducted as in Davis et al. (2004).

Plastid gene/genome alignment and data matrix construction

Phylogenetic analyses

Of the ‘coding’ matrices, mono132 and comm58 were analysed with unpartitioned, codon-partitioned and gene-partitioned GTR + Γ models, using 25 rate categories. PartitionFinder v.1.1.1 (Lanfear et al., 2012, 2014) was used for the mono132 matrix, to find an improved partition modelling strategy among 75 genes and three codon positions, using GTR models of RAxML, linked branch lengths, the Bayesian information criterion (BIC), and the ‘hcluster’ algorithm due to the large size of the dataset. Whole plastome alignments were analysed in RAxML (Stamatakis, 2006), under an unpartitioned GTR + Γ model. Analyses were conducted with 3–5 replicate maximum likelihood (ML) searches, and branch support was assessed with 1000 standard bootstrap pseudoreplicates. Partition models were compared using the corrected Akaike information criterion (AICc, Burnham & Anderson, 2002), which identifies the best fit model while penalizing overparameterization. The number of free parameters for each model (K) = P(2N-3) + (P × S), where N is the number of taxa, S is the number of GTR model parameters and P is the number of partitions. The AICc was calculated as −2log_eL + 2K(n/(n−K−1)), where log_eL is the likelihood, K is the number of free parameters, and n is the number of characters.

Additional ML analyses were conducted separately based on codon positions 1 + 2 and position 3 (unpartitioned GTR + Γ, as described earlier), respectively, of the mono132 matrix to assess differences in tree topology and lineage specific distances resulting from these two partitions. Furthermore, an analysis was run using ‘RY’ coding (with the unpartitioned ‘BINGAMMA’ binary model but otherwise as above), whereby pyrimidines and purines are coded as such, to explore the potential effect of systematic biases (e.g. GC content or transitional biases among clades).

Parsimony searches were conducted in TNT (Goloboff et al., 2008) with 100 random addition sequences, holding up to 10 000 trees, with tree bisection–reconnection branch swapping. Support was estimated via 1000 Jackknife pseudoreplicates, with 37% deletion probability using the same search parameters as described earlier. Branch support in all analyses is considered strong/robust if > 90, and weak if < 80.

Plastid substitution rate analyses, inference of rate changes

To assess variation in substitution rates among commelinid monocot orders, patristic distances from the common ancestral node of each individual order to each tip were calculated for the ‘mono132’ dataset, using the ‘ape’ and ‘geiger’ packages in R (Paradis et al., 2004; Harmon et al., 2008) and plotted in Past v.3.1 (Hammer et al., 2001). Statistical tests for differences in branch lengths among clades were not attempted due to issues of nonindependence of the data, (i.e. shared ancestral branches). Nearly all commelinid monocot families were represented; noncommelinid orders were not included in comparisons due to some limited sampling (e.g. for Alismatales, Dioscoreales, Pandanales, Liliales).

To test the null hypothesis that monocots evolve via a ‘Global Clock,’ analyses were conducted in the ‘baseml’ module of Paml v.4.8 (Yang, 1997, 2007). Specifically, different ‘branch models’ were tested, allowing rates to vary in prespecified regions of the tree corresponding to monocot orders, as opposed to a ‘background’ rate. Model M0 specified a global clock for all monocots; model M1 allowed Arecales to evolve via a local clock; M2 allowed Arecales and Poales to have independent local clocks; M3 allowed local clocks for all commelinid orders. Nested models were tested via likelihood ratio tests, where Λ = 2(log_eL₁−log_eL₂), which is chi-square distributed, and also using the corrected Akaike Information Criterion (Burnham & Anderson, 2002), where AICc = −2log_eL + 2K(n/(n–K−1)), in which n is the number of sites in the alignment, and K is the number of free parameters.

To investigate more explicitly the variation in substitution rates across the commelinid monocots in a phylogenetic context, a Bayesian approach was taken using the random local clock model in Beast v.1.81 (Drummond and Rambaut, 2007; Drummond & Suchard, 2010; Drummond et al., 2012). The RLC model has the advantage of requiring no a priori designation of the location of rate shifts, by assigning a prior on the number of rate changes throughout the tree, ranging from a global clock to a model in which all branches have independent rates. The objective was not to use fossil-calibration points for divergence time estimation, but to explore the probable number and locations of major shifts in relative substitution rate. Two matrices were analysed: the ‘comm58’ matrix with all coding loci, and a reduced matrix (comm58-3gene) that was much more computationally tractable than the former. Analyses were run using an unpartitioned GTR + Γ model (four rate categories, Yule speciation prior, Poisson rate change prior, and clock rate using CTMC rate reference (Ferreira & Suchard, 2008) with initial value of 1.0).

For the comm58matrix, three sets of 30 million generations from random seeds were run, and for the comm58-3gene matrix, seven sets of 100 million generations were run, with burn-in values from 20 to 85 million generations, respectively. Major ‘deep’ nodes were constrained based on relationships and support from previous analyses (Givnish et al., 2010; Barrett et al., 2013; this work) including: all commelinids, Zingiberales (Z), Commelinales (C), Zingiberales + Commelinales (ZC), Poales (P), ZC + P, Arecales (A), Dasypogonaceae (D), and A + D. Stationarity within and convergence between runs was verified using Tracer v.1.6 (Rambaut & Drummond, 2009), with effective sample sizes > 200 for parameters. Runs were combined with LogCombiner, trees modified in TreeAnnotator, and then visualized in FigTree (Rambaut, 2009), which allowed rate changes to be indicated among branches.

Mesquite v.3.03 (Maddison & Maddison, 2015) was used to assess the influence of third codon positions on topology and lineage specific substitution rates. The maximum likelihood tree based on codon positions 1 + 2 was compared with that from position 3 in a patristic distance correlation between the two trees. Analyses were run on a 12-core workstation with 64 Gigabytes of RAM (Linux).

Read coverage, plastome structure in the palms and alignments

The mean coverage depth of plastomes ranged from 19.5× (Mauritia) to 1918.1× (Salacca); these two accessions also had the lowest and highest total numbers of reads (1294 602 and 115 114 590, respectively; Table 1). Plastome structure is conserved across palms and Dasypogonaceae, with two exceptions. Eugeissona tristis (Calamoideae) has lost ndhF, based on 55.1× mean coverage of the plastome (excluding ndhF) vs 2.5× mean coverage of ndhF, and on contiguity of paired reads across the deleted region. Tahina spectabilis (Coryphoideae: Chuniophoeniceae) has lost nearly all of Inverted Repeat B (IR_B), except for a 159 bp fragment of ndhB (Fig. 1). Remapping reads to the annotated plastome of Tahina did not suggest the typical ‘doubling’ of read depth over the region corresponding to the IR observed in other plastomes (with IR_A removed), suggesting major excisions of the IR.

In addition, the plastome of Tahina has a 1944 bp inversion in the region corresponding to the Large Single Copy in other palms, with breakpoints between the 5′ exon of rps16/trnQ-UUG and trnS-GCU/5′ exon of trnG-UUG. Boundaries were verified by aligning the plastomes of Tahina and the closely related Chuniophoenix nana, and by mapping paired reads back to Tahina and visually examining the regions corresponding to the breakpoints of the inversion to check that each was spanned by several pairs of reads.

Phylogenetic placement of the palms among the commelinid monocots

Branch support based on the mono132 dataset is generally strong (Fig. 2), with 100 ML Bootstrap and 100 maximum parsimony (MP) Jackknife values (Fig. S6) for the monophyly of each order of monocots (except Pandanales, with one representative plastome). As well as for nearly all of the ‘deep’ branches relating orders of the monocots. Within the commelinids, there is moderate to strong support for (Arecales, Dasypogonaceae) (AD clade), and strong support for (Zingiberales, Commelinales, Poales; ZCP clade). Support is > 80 and 90 under all GTR partition models for AD and ZCP, respectively, but there is variation in support for these two clades (Fig. 2). Model partitioning affects support for these two nodes, particularly for AD; ML bootstrap support ranges from 91 (unpartitioned model) to 81 (codon partitioned) for this node, indicating that selection of model partition is significant. Comparison of the four ML model partitions used for the mono132 matrix indicates the model with the highest number of free parameters (gene partition) has the lowest AICc score, followed by the model identified by PartitionFinder (Table 2). Analysis of the RY-coded, mono132 matrix yielded identical major relationships among commelinid orders, with similar support values compared with those from standard GTR analyses based on nucleotide sequences (83 for AD and 86 for ZCP; Fig. S7c). The few weak to moderately supported nodes near the base of Poales are mostly associated with the placement of Typhaceae (Typha, Sparganium), Bromeliaceae (Brocchinia, Neoregelia, Puya), Mayacaceae (Mayaca) and Eriocaulaceae (Syngonanthus); under an unpartitioned model the latter two families form a clade, while they do not form a clade under other partitioning schemes. In all cases support is weak (< 60). Results in this region differed slightly for MP analyses, also with weak support (Fig. S6).

Phylogenetic relationships among the palms based on protein coding genes and complete plastomes

There is 100% support for the monophyly of each palm subfamily (excluding the monotypic Nypoideae) based on coding loci, as well as for the ‘deep’ branching pattern among subfamilies (Fig. 3). In particular, there is strong support for the position of Nypa. Support is generally strong for all relationships in Calamoideae, but the placement of Eugeissona is problematic. Depending on the analysis, Eugeissona is either weakly supported as sister to the rest of the calamoids (codon, gene partitioned models, and PartitionFinder), or (Eremospatha, Mauritia), then Eugeissona are successively sisters to the rest of the calamoids (tribe Calameae; unpartitioned model). Subfamily Coryphoideae comprises two strongly supported principal clades, each with four tribes. In one clade, Chuniophoeniceae (including Tahina) is sister to (Caryoteae (Corypheae, Borasseae)), all with 100% support in all analyses and for all clades. In the other clade, ((Sabaleae, Cryosophileae), (Phoeniceae, Trachycarpeae)), all have at least 99% support.

Overall, topology based on whole aligned plastomes (including spacers, introns, tRNAs and rRNAs; Fig. 4) is identical to that based only on protein coding genes, except for the placement of (Eremospatha, Mauritia) as sister to Calamoideae vs Eugeissona. In general, relationships and average support did not change drastically based on the different percentages of the full plastome alignment (moving from the full dataset to that containing only the most conserved sites as determined by Gblocks; Table 3). Average bootstrap support in the ML analysis was highest for the ‘Relaxed 1’ dataset (93% of total 159 074 aligned positions), with the lowest SD among mean values; on the other hand, ‘Strict 1’ (77% of full dataset) had the lowest mean support with highest SD (Table 3; Fig. 4). Interestingly, moving from more to fewer characters, overall support seemed to increase, then drop, followed by another increase, but overall gains/losses in support were small, ranging only 2.79%.

Notable differences in a few places, based on different datasets, include the position of the clade (Eremospatha, Mauritia) as sister to all other calamoids in all datasets except ‘Strict 2,’ which contains the fewest but most conserved characters; there was low support for relationships at the base of subfamily Calamoideae for all datasets. Within Coryphoideae, Colpothrinax is sister to Pritchardia (Trachycarpeae), collectively sister to a clade of (Washingtonia, ((Brahea, Chamaerops), (Licuala, (Acoelorrhaphe, Serenoa)))) for all datasets but ‘Strict 2’ with 100% support. For the ‘Strict2’ dataset, Pritchardia is sister to the latter clade, with Colpothrinax sister to all of these collectively. Within Trachycarpeae, Licuala is sister to (Acoelorrhaphe, Serenoa) in all but the ‘Full’ dataset; support for this relationship generally increases as data are removed.

Molecular evolutionary rate variation among clades of the monocots

Figures 2(a) and 5 show the striking differences in branch lengths estimated under the GTR + Γ model in RAxML across orders of the commelinid monocots. Among the commelinids, Arecales have relatively short branch lengths, while Poales have the longest, but with extensive variation. Within Poales, earlier diverging families have shorter branches (namely, Bromeliaceae and Typhaceae; Fig. 2a), with a trend of increasing branch length moving from the root of the order. Nucleotide composition also varies among plastomes, with higher mean GC content in Commelinales and Poales, but with extensive variation among plastomes within these two orders (Fig. 5a); Arecales and Dasypogonaceae have similar, relatively low GC profiles, but ML analysis of the RY-encoded mono132 matrix does not reveal major differences in topology or support compared with that of the standard nucleotide matrix (Fig. S7c). Patristic distances between ML trees based on codon positions 1 + 2 and position 3 had a correlation of 1.0, and relationships among major clades of the commelinids were identical between trees based on these two partitions (Fig. S7a,b).

Branch models in baseml/PAML (comm58; Table 4) indicate a significant departure from a ‘global clock’ (all rates equal among clades). A model allowing a local clock for Arecales (M1) has a significantly better fit than the global clock (M0); as does a model in which Arecales and Poales both have separate local clocks (M2). The best fitting model is one in which all commelinid orders evolve via separate local clocks (M3), as indicated by AICc comparison and Bonferroni-corrected likelihood ratio tests (Table 3).

Bayesian analyses under the random local clock (RLC) model in Beast of the comm58 dataset seemed to reach stationarity based on ESS values > 200 for parameters of interest, but did not converge between runs for several parameters (including the number of rate changes, the parameter of most interest), and also required months of computation time. When data were reduced to atpB, matK and rbcL (‘comm58-3gene’ dataset), four of seven runs converged on a median value of 22 substitution rate changes with 95% Highest Posterior Density (HPD) of 20–24 (Figs 6, 7). Rate accelerations occur in Poales, Commelinales and Zingiberales, and a few relatively minor rate changes occur in Dasypogonaceae and Arecales, with some rate decelerations in Arecales (Fig. 6). Comparison of posterior densities of the four converged runs vs runs without data (priors only) suggests that the 95% HPD of 20–24 rate changes is significantly different from the 95% HPD of 0–2 for prior only (Fig. 7).

Structure of the plastome in Arecales

The loss of ndhF from the plastome in Eugeissona tristis is surprising, but has been documented in other clades (Wicke et al., 2011; Barrett et al., 2014a). The ndh complex is encoded by nuclear and plastid genes, and functions in electron cycling, especially in environments of variable or high light intensity, ameliorating the effects of oxidative stress (Peng et al., 2011). Losses or pseudogenes have been observed in both photosynthetic and heterotrophic plants (Braukmann et al., 2009; Blazier et al., 2011; Wicke et al., 2013; Barrett et al., 2014a), in which selection on photosynthetic function may have been altered. Reads mapped to the ndhF region of the closely related Calamus, suggesting a copy of ndhF in either the nuclear or mitochondrial genome. Coverage of 2.5× for ndhF (based on Eugeissona reads mapped to Calamus) is more similar to that of the mitochondrial genome than to that expected for the nuclear genome (< 1×), suggesting a functional copy in the mitochondrial genome.

Additional sequencing or long-range PCR may allow a more definitive answer as to the location of ndhF. However, a mitochondrial-encoded ndhF would not likely be adaptive (assuming the ndh complex is functional in Eugeissona); the ndh complex would thus require protein products from all three genomes. If ndhF is encoded in the mitochondrial genome and still functional, it would possibly serve an alternative function or may be pleiotropic, functioning in both the mitochondrion and chloroplast. Fragmented mitochondrial copies of ndh genes have been observed in the orchid Erycina pusilla but these are not likely expressed (Lin et al., 2015).

The lack of an IR in Tahina spectabilis (except a 159 bp, duplicate fragment of ndhB) is intriguing, since the ‘canonical’ structure of nonpoalean and nonheterotrophic monocots is generally conserved. Loss of the IR has been observed, for example, in papilionoid legumes (Palmer et al., 1987; Lavin et al., 1990; Liston, 1995) and some gymnosperms (Strauss et al., 1988; Wu et al., 2011; Wu & Chaw, 2014), which is hypothesized to have led to a destabilization of the gene order in the plastome, evidenced by rearrangements. The simplest mechanistic explanation for IR loss in Tahina is that two massive boundary shifts of the single copy/IR junctions occurred, diminishing the extent of the IR. According to the IR stabilization hypothesis, a complete or nearly complete loss of the IR would be followed by one or more gene order rearrangements. Tahina has a 1994 bp inversion, which is in line with that expectation, assuming the inversion occurred after IR reduction. The fact that changes in plastome structure were detected among 38 representative taxa suggests that deviations from canonical structure may be more prevalent than previously hypothesized, given the high degree of palm plastome conservation. It also suggests that structural changes do not occur deep in the palm plastid phylogeny as synapomorphies but unite more exclusive clades or may be present idiosyncratically in species or populations.

Phylogenetic placement of the palms among the commelinid monocots

The placement of the palms as sister to Dasypogonaceae, with moderate to strong support, is in agreement with recent plastid-based studies (e.g. ML analysis of Givnish et al., 2010; Barrett et al., 2013; but see Davis et al., 2013), as well as the placement of this clade as sister to a strongly supported clade of (Poales, (Zingiberales, Commelinales)). The increased taxon sampling in the current study based on nearly complete plastid gene sets has resulted in slightly stronger support for AD and ZCP (e.g. comparing support values for Fig. 2b and those in the ML analyses of Barrett et al., 2013), some of the most recalcitrant nodes in the monocot Tree of Life. Dense taxon sampling can ameliorate problems caused by the ‘node density artifact’ (Fitch & Bruschi, 1987; Venditti et al., 2006), specifically, that a lack of dense taxon sampling may result in extremely short internodes and a resulting lack of support. Since relationships are identical in this and earlier plastome-based studies among the major clades of the commelinids, there was therefore not a strong node-density artifact in earlier studies with sparser sampling – or at least the effect was not strong enough to affect topology. However, the observed slight increase in support in the current study may be due to the effects of denser taxon sampling for the two relatively short branches leading to AD and ZCP (Fig. 2).

Model based partitioning schemes affect support values for AD and ZCP; for the former, support values vary substantially, by 10. For both clades, codon partitioned models give lower support than gene or unpartitioned models, while those based on PartitionFinder give similar values to codon based models (Fig. 2b). Thus, methods that improve fit to the data based on gene × codon partitioning (e.g. PartitionFinder) are crucial to conclusions of support, especially for complex, heterogeneous datasets such as mono132.

Phylogenetic relationships among the palms based on protein coding genes and whole plastomes

There is complete resolution and strong support for relationships among the five palm subfamilies, reflecting relationships in previous multigene analyses (Figs 3, 4; Asmussen et al., 2006; Baker et al., 2009). Of particular interest is support for the placement of Nypa fruticans, which comprises the subfamily Nypoideae. Whereas earlier molecular studies were inconclusive on the positions of Nypa and Calamoideae relative to other palms (Uhl et al., 1995; Lewis & Doyle, 2001), later studies based on a few plastid loci but more complete taxon sampling suggested that Calamoideae are sister to Nypa + remaining Arecales. Here, the position of Nypa is supported as sister to (Coryphoideae, (Ceroxyloideae, Arecoideae)), and the calamoids sister to this entire clade, with 100% support in all analyses (Figs 3, 4). Asmussen et al. (2006) recovered strong support for this topology based on four plastid loci and nearly complete generic level sampling, while Baker et al. (2009) found strong support using all available data and generic level sampling in supertree analyses (e.g. > 90% of input trees supported Nypa in this position).

Asmussen et al. (2006) and Baker et al. (2000a,b, 2009; supertree analysis, but see supermatrix), recovered Eugeissona as sister to the rest of the subfamily Calamoideae. Here, Eugeissona is either sister to the remaining calamoids, or a clade of (Eremospatha, Mauritia) and then Eugeissona are successively sister to the remaining calamoids (Figs 3, 4; depending on the analysis). Findings of a strongly supported sister relationship between Eremospatha and Mauritia have implications for the monophyly of the tribe Lepidocaryeae, which was paraphyletic in the supermatrix analysis of Baker et al. (2009; but see supertree analysis). Sampling of additional calamoid genera, including Raphia, would likely yield stronger support in this part of the tree.

Support is strong for relationships among the tribes of Coryphoideae (Figs 3, 4). The current findings (Figs 3, 4) closely resemble the relationships of Baker et al. (2009) based on a supertree analysis using irreversible weights matrix representation with parsimony. Baker et al. (2009) recovered varying coryphoid relationships among several supertree and supermatrix analyses. Both their supertree and the trees recovered here contained two principal clades within Coryphoideae, each represented by four tribes; the only difference is the relative placement of Caryoteae and Chuniophoeniceae within the ‘syncaropus’ clade, which also contains Borasseae and Corypheae (Figs 3, 4; Dransfield et al., 2008; Baker et al., 2009). However, relationships among coryphoid tribes generally did not have strong support in the supertree analysis of Baker et al. (2009), whereas here they received strong support. Asmussen et al. (2006) recovered a weakly supported clade of (Trachycarpeae, Phoeniceae) as sister to (Chuniophoeniceae, (Caryoteae, (Corypheae, Borasseae))). Plastid genomes appear to resolve tribal relationships among the coryphoids with strong support.

The diverse and widespread tribe Trachycarpeae has been problematic from a phylogenetic perspective. Based on a multilocus approach with dense taxon sampling at the generic/species level, Bacon et al. (2012) were able to resolve many relationships among the genera of this tribe; our findings are similar, but the relative placement of Washingtonia, Pritchardia and Colpothrinax differs. Here (Figs 3, 4), the clade of (Pritchardia, Colpothrinax) is sister to (Washingtonia, remaining Trachycarpeae) with strong support, while Bacon et al. (2012) recovered (Washingtonia, Pritchardia) as sister to (Colpothrinax, remaining Trachycarpeae), but with weak support (< 75 in ML Bootstrap analysis). Support for a sister relationship between (Pritchardia, Colpothrinax) is moderate based on coding loci (Fig. 3), but strongly supported based on whole aligned plastomes (Fig. 4).

Except for the placement of Eugeissona among the calamoids, palm relationships are identical between coding and whole plastid matrices. This suggests that at least for slowly evolving groups, alignment of whole plastomes does not introduce enough homoplasy (i.e. among highly variable spacers and introns) as to drastically change relationships, and in many places increases support relative to more reduced matrices. The effect of data removal (relative to the ‘full’ plastome matrix) generally is a slight decrease in support (Fig. 4), with the exception of the placement of Licuala as sister to (Serenoa, Acoelorrhaphe).

Variation in plastid substitution rate among clades of the monocots

It has been hypothesized that monocots display extreme substitution rate heterogeneity, based on earlier plastid and nuclear studies (Bousquet et al., 1992; Gaut et al., 1992, 1996; Morton & Clegg, 1995), but these studies mainly included representatives of the palm and grass families and were based on few loci. Here we included sampling across all monocot orders, with representative plastomes of nearly all families of the commelinids. This study represents the largest genome-scale phylogeny of commelinids (e.g. compared with the 44 commelinids in Givnish et al., 2010; and 27 in Barrett et al., 2013) and reevaluation of a classic case of rate variation in the genomic era. Our findings (Figs 5-7) provide a genome-wide perspective on extreme plastid rate variation, support previous hypotheses, and further indicate that earlier findings are not an artifact of locus choice (e.g. rbcL, ADH), but that extreme substitutional rate variation is pervasive based on nearly complete coding complements of the plastome. Moreover, patristic distances are highly correlated among ML trees (coefficient = 1.0) based on codon positions 1 + 2 and position 3, suggesting that substitutional rate variation among clades is not driven exclusively by rate changes at codon position 3, but by changes across all codon positions. Based on Bayesian RLC analysis, we find strong evidence for c. 22 rate shifts among the commelinids (Fig. 7); additional sampling of taxa and loci may increase this estimate (but see earlier regarding issues with the convergence of the RLC model with all 75 genes). Accelerations mostly occur in the ZCP clade, and are prevalent among Commelinales, Zingiberales and Poales, while decelerations occur in Arecales.

These findings have phylogenetic implications, as extreme variation in substitution rates may cause issues in phylogenetic analysis (discussed in Rodriguez-Ezpeleta et al., 2007; Philippe & Roure, 2011), which is exacerbated by sparse taxon sampling. Here, we have sampled representatives of nearly all families of the commelinids, and observed highly similar large scale relationships (A, D), (ZC, P) relative to previous studies (Givnish et al., 2010; Barrett et al., 2013), with strong to moderate support, suggesting that the sparse sampling and rate heterogeneity in earlier studies did not severely obscure plastid relationships (Givnish et al., 2010; Barrett et al., 2013).

What are the potential correlates/drivers of such rate variation among commelinid clades? Some researchers have identified correlates across plants such as species diversity, generation time, growth form or environmental factors (Barraclough & Savolainen, 2001; Davies et al., 2004; Smith & Donoghue, 2008; reviewed in Burleigh, 2012; Lanfear et al., 2013). Lanfear et al. (2013) accounted for both phylogeny and a number of the aforementioned factors, and found that a strong relationship remained between plant height and substitution rate across angiosperms: taller angiosperms having slower substitution rates. This relationship was observed for both nuclear and plastid DNA. The authors explain their findings in terms of the ‘rate of mitosis’ (ROM) hypothesis. In plants, different tissues derived from the apical meristem may develop into reproductive tissues (accruing heritable, somatic mutations), unlike in animals, in which the germline is sequestered early in development. Furthermore, some plants may experience a deceleration of the rate of mitosis in the apical meristem as they near their maximum height (Lanfear et al., 2013), but whether this occurs in palms is unclear due to the need for replacement of ephemeral but essential organs such as leaves. Palms generally are the tallest of the monocots, and palm shoot tissues often derive from a single apical meristem (although many palms are sympodial, with clustered or branching stems). The ROM hypothesis is a good candidate for future investigation as a potential explanatory factor for diminished plastid substitution rates in the palms and rate accelerations in some lineages of Poales, Commelinales and Zingiberales (though some bamboos are > 10 m).

Different factors may affect substitution rates in different ways across the commelinids or the importance of multiple factors may also vary among clades. For example, accelerated substitution rates among clades within Poales (e.g. Poaceae) may be a combined product of high ROM (and short stature) and other life history traits such as C₄ photosynthesis, potentially leading to high rates of speciation and widespread ecological dominance, but these remain as hypotheses to be tested. The potential correlates/drivers of rate variation among monocot orders warrant more explicit investigation, incorporating trait, environmental, nuclear and mitochondrial data, along with increased taxon sampling, and the use of multivariate, phylogenetic comparative methods (Lanfear et al., 2010).

Conclusion

Plastid gene sets generated by Illumina NGS, and representing nearly all families of the commelinid monocots, provide robust branch support and evidence for extreme heterogeneity of plastid substitution rates. Analysis of whole plastomes in the palms further resolves and provides strong support for relationships among the five subfamilies, and for relationships among tribes of the diverse subfamily Coryphoideae. Model based analyses suggest that commelinid orders deviate significantly from a global clock in terms of plastid substitution rates, and also that c. 22 rate shifts have occurred across the commelinids. These findings have implications for the systematics of the commelinid orders, which represent immense species diversity and ecological/economic importance, and lay the groundwork for future comparative tests of the drivers and correlates of substitution rates.