Utility of a large, multigene plastid data set in inferring higher-order relationships in ferns and relatives (monilophytes)^†

Taxonomic and genomic sampling

The final matrix considered here includes 64 representatives of the major land-plant lineages (Appendix S1; see Supplemental Data at http://www.amjbot.org/cgi/content/full/ajb.0900305/DC1; and see Appendix 1 for a list of newly generated sequences and associated GenBank accession numbers). Of these, four represent the major lineages of lycophytes (Selaginella uncinata and Huperzia lucidula sequences were obtained from GenBank; accession numbers AB197035 and NC_006861, respectively). We included 16 gymnosperms and 6 angiosperms to represent the deepest splits in seed-plant phylogeny (48; 12). We also included four bryophyte outgroups (GenBank accession numbers for Anthoceros formosae, Marchantia polymorpha, and Physcomitrella patens are NC_004543, NC_001319, and NC_005087, respectively; see Appendix 1 for Sphagnum sp.). The remaining 34 taxa sample the major clades of monilophytes according to 42, 1) and 58) and were chosen, as far as possible, to sample the crown-clade root nodes (= most recent common ancestors, MRCAs) of relevant clades. They include eight eusporangiate monilophytes and 26 leptosporangiate ferns (Adiantum capillus-veneris, Angiopteris evecta, and Psilotum nudum sequences were obtained from GenBank, accession numbers NC_004766, NC_008829, and NC_003386, respectively; see Appendix 1 for the remainder). We sampled all seven orders and 21 of the 33 families of leptosporangiate ferns recognized by 62; unsampled families are mainly small ones in Cyatheales and Polypodiales (largely comprising one to five genera; 62).

We surveyed 17 genes and associated noncoding regions considered by 15 and 47, 36). However, we excluded an intergenic spacer between rps7 and ndhB from consideration. This region is present in all seed plants surveyed for our study but is apparently not present as a contiguous region in most leptosporangiate ferns, owing to a large inversion that involves a substantial portion of the inverted repeat (49; 71).

DNA extraction, amplification, and sequencing

We extracted genomic DNA using the protocol of 8), as modified in 47, from fresh, silica-dried, and herbarium specimens. Amplification and sequencing follows 15 and 47, 36). We designed a set of 16 new fern-specific primers to facilitate amplification and sequencing (Appendix 2). With a few minor exceptions, all products were completely sequenced in both forward and reverse directions. The new data were added to a previously generated alignment (48) using criteria outlined in 16 and 48. Regions in some taxa that we could not amplify or sequence were coded as missing data in the final matrix (see Appendix 1).

The aligned regions considered for analysis included all of the coding regions (Table 1) in addition to unambiguously aligned noncoding regions from two introns (rpl2 and ndhB) and seven intergenic spacer regions (i.e., 3′-rps12-rps7, and three each in the psbE-psbF-psbL-psbJ and psbB-psbT-psbN-psbH clusters). We dealt with hard-to-align regions, often limited to individual taxa, by staggering them in the alignment (e.g., 53), coding the gap cells as missing data. Offset unique regions are effectively ignored in parsimony analysis and should have only minor effects in maximum-likelihood analysis (e.g., on the estimation of base frequencies). This approach contributes significantly to the large length of the full concatenated alignment (36139 bp, based on ∼10 kb per taxon in monilophytes; Table 1). For the full alignment, 2573 bp are variable and uninformative across vascular plants and 7479 bp are parsimony informative.

Phylogenetic analyses

We performed a heuristic maximum likelihood (ML) search using PhyML 2.4.4 (20) with a BIONJ starting tree, NNI (nearest-neighbor-interchange) branch swapping, and model parameters estimated from the data in each case. We chose a DNA substitution model for ML analysis using the hierarchical likelihood ratio test (hLRT) and the Akaike information criterion (AIC) in Modeltest 3.7 (41). The optimal model chosen in each case was GTR + Γ + I (general-time-reversible [GTR] rate matrix with proportion of invariable sites [I] considered, and among-site rate variation accounted for using the gamma [Γ] distribution as described by the shape parameter alpha (α)]. We also performed a heuristic maximum-parsimony (MP) search using PAUP* 4.0b10 (67), with all characters and character-state changes equally weighted, and using TBR (tree-bisection-reconnection) branch swapping, with 100 random-addition replicates, and otherwise using default settings. We performed nonparametric bootstrap analysis (11) using the same search criteria, with 100 bootstrap replicates (and a single random-addition replicate per parsimony bootstrap replicate). We use “weak,” “moderate,” and “strong” in reference to clades that have bootstrap support values <70%, 70–89%, and ≥90%, respectively (e.g., 13; 47, 36).

Inference of nucleotide rate classes and exploration of the effect of long branches

We partitioned the data into nine rate classes using HyPhy (40; and see 1; 48). We used the best MP tree topology to assign each of the aligned sites to its most likely individual rate category (using the GTR model and eight discrete rate classes). We then excluded the two fastest rate classes (RC7 and RC8) and performed a maximum-likelihood search and bootstrap analysis on the remaining characters (RC0–RC6), using the search criteria outlined above. The sequence for Selaginella is highly divergent compared with the rest of the vascular plants (see Results; also evident from a visual examination of the alignment). There is evidence of elevated plastid genome evolution for the family, reported for rbcL by 33; 1), which may be a consequence of the extensive RNA editing observed in Selaginella (69). We therefore recalculated rate classes with Selaginella removed from consideration and repeated the ML search (again for RC0–RC6). Overlaps in the nucleotides estimated to be in each rate class, with Selaginella either included or deleted, were determined by including or excluding the relevant rate-class character sets in PAUP*.

Phylogeny of ferns and relatives

We have reported on results for seed plants elsewhere, and so we focus here on phylogenetic relationships involving monilophytes and lycophytes. The majority of the vascular-plant clades inferred here were strongly supported by both ML and MP analysis. These include all seven monilophyte families sampled for two or more taxa and 18 of the 27 multifamily clades recovered in the best ML tree (which are labeled as clades a–aa; Fig. 1, Table 2). The latter clades comprise lycophytes, monilophytes, Psilotopsida (= Ophioglossaceae-Psilotaceae), all five multifamily orders in 62; i.e., Cyatheales, Gleicheniales, Polypodiales, Salviniales, Schizaeales), leptosporangiate ferns, leptosporangiate ferns excluding Osmundaceae, core leptosporangiate ferns (Salviniales+Cyatheales+Polypodiales), eupolypods, eupolypods I, and several others of interest (clades labeled as: a, j, m, r, and v in Fig. 1 and Table 2). The eupolypods II clade (= clade z) was strongly supported by ML and moderately well supported by MP, and two clades were moderately strongly supported by both ML and MP (clades l and p). Thus, 78% of clades (= 21 of 27 taxon partitions) seen in the ML tree that involved two or more families of ferns or lycophytes simultaneously had moderate to strong support from ML and MP bootstrap analysis.

Two clades were moderately strongly supported by ML analysis but had <50% support from MP (clade g, with Marattiaceae sister to leptosporangiate ferns, and clade c = euphyllophytes). The latter clade, the euphyllophytes, was also contradicted by MP (clade ad = lycophytes+seed plants, which had 94% support from parsimony). Two clades had weak support from ML but were moderately well supported by MP (clade aa, and clade e = monilophytes excluding Equisetaceae). The remaining clades recovered for the shortest MP or ML trees conflicted between phylogenetic inference methods but were also poorly supported (i.e., clades u and x, seen in the ML tree [Fig. 1], and clades ab, ac, and ae, observed in the single MP tree but not in the best ML tree [Table 2]). Only one substantial conflict was observed in total (clade c vs. ad, mentioned above).

Rate-class analyses

When we filtered the most rapidly evolving rate classes (RC78), with Selaginella included or excluded during site-rate assignment, and reran the ML analyses with the remaining more conservative characters (RC0–RC6), 15 clades had approximately the same ML support values before and after removal of rapidly evolving characters, considering both filtering methods. The ML bootstrap values changed (fell) by 10% or more bootstrap support for both filtering methods for five clades total in the best ML tree (i.e., clades k, l, m, o, aa; see Table 2). We considered >10% change in bootstrap support as noteworthy, following 57: fig. 2). Considering the two filtering cases here (fast sites removed, with Selaginella included or excluded during the rate-class assignments), distinctly different patterns of change in ML bootstrap support were observed in six cases, for comparisons of bootstrap support values before and after filtering of fast sites (Table 2). In four cases (clades b, c, d, and h), the ML bootstrap support fell by 10% or more compared with the unfiltered analysis when Selaginella was included during rate-class assignment but was essentially static for the same clades when Selaginella was excluded. All four clades are early branches of monilophyte or vascular-plant phylogeny (i.e., lycophytes as a whole, euphyllophytes, monilophytes as a whole, and leptosporangiate ferns as a whole). Support for clade a (Isoëtaceae+Selaginellaceae) fell by more than 10% when Selaginella was included in rate-class inference, with no corresponding value available for when Selaginella was excluded (by definition).

A possible sister-group relationship between Dennstaedtiaceae and Pteridaceae (clade x) was poorly supported in the unfiltered ML analysis (<50% support), but this increased to moderate support when Selaginella was excluded from rate assignments (71% ML bootstrap support); a similar pattern was observed for the first branch in monilophyte phylogeny: Equisetaceae was inferred to be the sister to all other monilophytes with weak support in the unfiltered ML analysis (clade e, 58% support) but moderate support with fast sites removed and Selaginella excluded from rate assignments (71% support; Table 2).

We compared how the ML assignments of individual sites to rate classes differed when Selaginella was included or deleted. We observed substantial overlap in rate-class assignments with or without this problematic taxon (Table 3), and most differences in how rate classes were assigned involved neighboring rate classes (either one rate class higher or lower). Thus, rate-class assignments were in general similar whether Selaginella was present or deleted. However, with either filtering approach (Selaginella included vs. deleted), the alternative case tended to place more characters in the next-highest class, rather than the next-lowest class, particularly for the lower rate classes (considering displacements off the diagonal in Table 3; note the relative size and direction of change in rate assignments moving out from the diagonal). On the whole, the subset of characters that were assigned differently by the two methods tended to be placed in faster rate classes when Selaginella was included in rate assignments. This effect was most pronounced for the zero-rate class. When Selaginella was included, 349 characters were assigned to RC2–RC8 that were assigned to RC0 with Selaginella excluded. By contrast, no characters assigned to RC0 with Selaginella excluded were assigned to higher rate classes when this taxon was included (Table 3).

Congruence and discordance with other phylogenetic studies of ferns and relatives

Phylogenetic congruence, the ability to corroborate phylogenetic relationships using different gene and taxon sets, is a key pillar of systematic biology (e.g., 39; 24; 13; 42; 56) because it builds confidence that the recovered topologies reflect evolutionary history rather than arbitrary groupings due to stochastic error. Conversely, discordance among studies may also allow us to flag and investigate possible instances of systematic error, such as long-branch attraction (9). Our analysis of vascular-plant relationships using a large multigene plastid data set is, in general, highly congruent with earlier multigene studies, corroborating clades in common across studies (Fig. 1, Table 2). In addition, our ML bootstrap support values for major multifamily clades either are within a few percentage points of other studies or are equal to them or better (Table 2). Even at the currently relatively limited taxon density, we found moderate to strong ML bootstrap support for ∼85% (23 of 27) multifamily clades in the best ML tree (Fig. 1), including some of the most problematic nodes. This supports the general utility of our approach for the current and future investigations of monilophyte phylogeny. However, several major branches of monilophyte phylogeny are poorly supported or in conflict in all published studies. Here, we focus on three important areas of monilophyte phylogeny that are currently poorly understood: the branching order of the earliest splits in monilophytes, leptosporangiate ferns, and polypod fern phylogeny, respectively.

Effect of long branches and rapidly evolving sites on inference of deep fern phylogeny

Maximum likelihood is expected to be less sensitive to long-branch attraction than parsimony (e.g., 4; 27, 7; 64; 68). Strong disagreement between ML and MP analysis may therefore be a flag for long-branch attraction in phylogenetic inference. However, considering the multifamily clades of monilophytes and lycophytes, our parsimony and likelihood estimates were generally similar, which we view as a reassuring result. Rapidly evolving sites should be accommodated in phylogenetic inference using maximum likelihood, so long as the substitution model is a reasonable fit. When we filtered out high-rate sites (RC78), this generally had little effect on phylogenetic inference beyond reducing some ML bootstrap support values (Table 2), even though the model MP tree used for rate classification likely had an incorrect relationship (i.e., lycophytes + seed plants). These reductions in support values may be consistent with expectations of increased sampling error in smaller data sets. However, we observed an exceptionally long branch for Selaginella (Fig. 1), also noted by 33, 1) for rbcL. When we removed this taxon prior to rate classification, and repeated ML analysis on the filtered data set, four of these clades no longer experienced a reduction in support (and support for two clades surpassed what was seen for the unfiltered ML analysis; clades e and x). These clades include some of the earliest splits in euphyllophyte phylogeny. Examination of how nucleotides were assigned to different rate classes with or without this taxon indicates that some sites behave quite differently when Selaginella is included in the rate assignments (Table 3). This may be at the root of how this taxon may reduce some ML support values in ferns and relatives. However, until this effect is better characterized, it may be useful to exclude Selaginella as an outgroup in future studies of monilophyte phylogeny, at least for plastid-derived sequence data.

One strong conflict between ML and MP involves a well-supported (but likely incorrect) sister-group relationship inferred here between lycophytes and seed plants using MP (with corresponding poor MP support for euphyllophytes; Table 2). Our ML analysis did not recover this relationship, but instead moderately supported euphyllophyte monophyly (81% support). 56 also recovered the euphyllophyte clade in their morphological study, with strong MP bootstrap support (Table 2). However, in other molecular studies a similar pattern of good ML versus poor MP bootstrap support for euphyllophytes has been found, for example by 42, who used a comparable sampling of bryophyte outgroups (Table 1). Subsequent molecular studies by Pryer and colleagues either did not report parsimony bootstrap values for euphyllophytes (70) or did not include bryophytes in the analysis (44), precluding full assessment of euphyllophyte monophyly (Table 2). The best MP and ML support for euphyllophytes among these multigene studies was found by 46, 2), who used a much denser sampling of other land plants in their analyses. The possible negative effect of long bryophyte branches on the inference of early monilophyte phylogeny has not been systematically studied, but it may be ameliorated by substantially increasing the number of taxa sampled. This should be a key sampling consideration for future broad phylogenetic studies of euphyllophyte phylogeny, and we are currently working to substantially increase taxon density for this gene set in bryophytes (Y. Chang and S. W. Graham, unpublished manuscript).