Overcoming Deep Roots, Fast Rates, and Short Internodes to Resolve the Ancient Rapid Radiation of Eupolypod II Ferns

. —Backbone relationships within the large eupolypod II clade, which includes nearly a third of extant fern species, have resisted elucidation by both molecular and morphological data. Earlier studies suggest that much of the phylogenetic intractability of this group is due to three factors: (i) a long root that reduces apparent levels of support in the ingroup; (ii) long ingroup branches subtended by a series of very short backbone internodes (the “ancient rapid radiation” model); and (iii) significantly heterogeneous lineage-specific rates of substitution. To resolve the eupolypod II phylogeny, with a partic- ular emphasis on the backbone internodes, we assembled a data set of five plastid loci ( atpA , atpB , matK , rbcL , and trnG-R ) from a sample of 81 accessions selected to capture the deepest divergences in the clade. We then evaluated our phylogenetic hypothesis against potential confounding factors, including those induced by rooting, ancient rapid radiation, rate heterogeneity, and the Bayesian star-tree paradox artifact. While the strong support we inferred for the backbone relationships proved robust to these potential problems, their investigation revealed unexpected model-mediated impacts of outgroup composition, divergent effects of methods for countering the star-tree paradox artifact, and gave no support to concerns about the applicability of the unrooted model to data sets with heterogeneous lineage-specific rates of substitution. This study is among few to investigate these factors with empirical data, and the first to compare the performance of the two primary methods for overcoming the Bayesian star-tree paradox artifact. Among the significant phylogenetic results is the near-complete support along the eupolypod II backbone, the demonstrated paraphyly of Woodsiaceae as currently circum- scribed, and the well-supported placement of the enigmatic genera Homalosorus , Diplaziopsis , and Woodsia . [Moderate data; outgroup rooting; Phycas ; phylogeny evaluation; rate heterogeneity; reduced consensus; star-tree paradox; Woodsiaceae.]

A classic problem in phylogenetics is the reconstruction of "ancient rapid radiations," broadly defined as evolutionary histories where long branches are intercalated among a series of short backbone internodes (see Fig. 1; Whitfield and Lockhart 2007;Jian et al. 2008). Accurately resolving such topologies is a welldocumented challenge for phylogenetic inference (Gaut and Lewis 1995;Huelsenbeck 1995;Jackman et al. 1999;Anderson and Swofford 2004;Wang et al. 2009) and is also of considerable practical importance-this ancient rapid radiation model is a prominent feature of many phylogenetic problems (Whitfield and Lockhart 2007). Furthermore, the ancient rapid radiation pattern rarely exists unaccompanied; rather, it tends to coincide with other well-recognized analytical challenges. First, the phylogenetic root is often long with respect to ingroup branches (Fig. 1;Bergsten 2005;Schuettpelz and Hoot 2006). Because signal deteriorates along phylogenetic branches (in a likelihood framework), long branches are less likely than short ones to strongly affix to any single point in the topology (Wheeler 1990;Swofford et al. 1996;Huelsenbeck et al. 2002). Furthermore, although the monophyly of the ingroup and of all ingroup relationships may be fully supported, uncertainty in the placement of the root may nonetheless reduce apparent support for relationships among ingroup clades when one uses consensus-based measures to assess support (Wilkinson 1996;Roberts et al. 2009). Second, lineage-specific heterogeneity in rates of substitution is common, making "fast" taxa particularly difficult to place (Fig. 1;Felsenstein 1978;Hillis and Bull 1993;Soltis et al. 1999;Takezaki and Gojobori 1999; but see Ho and Jermiin 2004;Nickrent et al. 2004;Drummond et al. 2006). Finally, the presence of both very short and very long branches-regardless of their topological arrangement-poses additional challenges. While long-branch attraction has been well characterized (Felsenstein 1978;Anderson and Swofford 2004;Bergsten 2005), other branch-length related inconsistencies are just beginning to attract attention Yang and Rannala 2005;Yang 2008;Marshall 2009;Roberts et al. 2009;Brown et al. 2010).
One option for tackling problems associated with reconstructing ancient rapid radiations is to amass character-rich (often genome-scale) data sets (e.g., Pereira and Baker 2006;Hallstrom et al. 2007;Jian et al. 2008;Wang et al. 2009;Regier et al. 2010). However, the specific challenges inherent to this sort of phylogenetic problem are not necessarily amenable to resolution by greatly expanded character data (Philippe et al. 2011). Rather, increasing character data can yield increasingly strong support for erroneous relationships, especially 490 2012 ROTHFELS ET AL.-ANCIENT RAPID RADIATION OF EUPOLYPOD II FERNS 491 FIGURE 1. Challenges inherent in resolving the eupolypod II phylogeny. Eupolypods II phylogram modified from , in (a) unrooted and (b) rooted form. (i) Outgroup taxa are on long branches. (ii) Backbone internodes are very short, suggesting an "ancient rapid radiation." (iii) The ingroup is marked by significant heterogeneity in rates of evolution, with the members of Aspleniaceae on much longer branches than other eupolypod II taxa.
in cases of branch-length variation such as is inherent in the long-root, short internode, and rate heterogeneity features common under the ancient rapid radiation model (Gaut and Lewis 1995;Soltis et al. 2004;Bergsten 2005;Philippe et al. 2005;Steel and Matsen 2007;Whitfield and Lockhart 2007;Rannala and Yang 2008;Susko 2008;Yang 2008). Here, we focus on resolving an ancient rapid radiation, that of the eupolypods II clade, using moderate amounts of character data but a strongly expanded taxon sample (for a recent application of this "moderate data/many taxa" approach, see Parfrey et al. 2010). This fern clade has resisted elucidation by both morphological and molecular data (Ching 1964a;Sledge 1973;Smith 1995;Sano et al. 2000a;Smith et al. 2006;Wei et al. 2010;Kuo et al. 2011), and previous molecular studies indicate that it exhibits all of the analytical challenges outlined above (see Fig. 1).
The Eupolypods II, together with its sister group, Eupolypods I, comprise the large eupolypod clade, which encompasses two-thirds of living fern species (Fig. 2;Pryer et al. 2004;Schneider et al. 2004b;Smith et al. 2006). The ancestors of Eupolypods I and II diverged from each other in the Early Cretaceous (Pryer et al. 2004;Schneider et al. 2004b;Schuettpelz and Pryer 2009). The eupolypod II clade started to diversify shortly thereafter (its crown group is approximately 100 million years old; Schuettpelz and Pryer 2009) and has subsequently grown into a lineage-rich clade comprising nearly 30% of extant fern diversity. Eupolypods II includes some of the most familiar groups of ferns (the lady ferns, ostrich fern, sensitive fern, marsh ferns, and spleenworts), as well as some of the most speciesrich genera: Thelypteris s. lat. (∼950 species); Asplenium (∼700 species); Diplazium (∼350 species); Athyrium (∼220 species); and Blechnum s. lat. (∼150 species).
The eupolypod II clade is cosmopolitan in distribution, with the subgroups primarily temperate to tropical, and the larger subclades each well represented in both areas. However, many of the phylogenetically enigmatic taxa in this clade are limited to the Himalayas or Southeast Asia, and critical members of several genera are rare and/or infrequently collected. This pattern of rarity and endemism, in conjunction with the richness and geographical breadth of the clade as a whole, is undoubtedly a contributing factor to the incomplete sampling of these ferns in previous phylogenetic studies.
Not surprisingly, given the clade's size and age, eupolypod II taxa are morphologically disparate and seemingly incohesive. However, early workers did tend to recognize the close affinities among many of the taxa in this clade, although frequently with members of Eupolypods I interdigitated among them (Holttum 1947;Sledge 1973;Mickel 1974;Tryon and Tryon 1982). The cohesiveness of the Eupolypods II started to become apparent with the earliest applications of molecular phylogenetic techniques to ferns Hasebe et al. 1995;Sano et al. 2000a) and has been strongly supported in recent broad studies (Schneider et al. 2004b;Kuo et al. 2011). None of these studies, however, found support for the backbone relationships within Eupolypods II, and only Kuo et al. (2011) attempted to sample its major lineages. It remains one of the few areas of the fern tree-of-life where the backbone relationships remain elusive (Smith et al. 2006;. Our approach to resolving the eupolypod II phylogeny couples a considerably expanded taxon sample with moderate character sampling. Our objectives include identifying well-supported major (approximately "family-level") clades and determining the backbone relationships among these clades. Given the anticipated 492 SYSTEMATIC BIOLOGY VOL. 61 FIGURE 2. Broad phylogeny of ferns. Approximate number of species per clade is given in parentheses. Modified from Smith et al. (2006). phylogenetic challenges and potential for artifacts in our data, we explicitly evaluate our phylogenetic hypothesis against these analytical pitfalls, placing strong emphasis on the use of the reduced consensus technique (Wilkinson 1996) to isolate the effects of signal weakness from those of signal conflict (e.g., Wiens 2003;Cobbett et al. 2007). Our study aims for a comprehensive and well-supported phylogeny of this important group of ferns, and for novel inferences about the behavior of our choice of methods, gleaned from their performance on this data set.

MATERIALS AND METHODS
Taxon Sampling We selected an ingroup of 67 species, intended to maximize our capture of the deep divergences (Zwickl and Hillis 2002) within Eupolypods II. Decisions for inclusion were based on data from previous molecular (Gastony and Ungerer 1997;Murakami et al. 1999;Sano et al. 2000aSano et al. , 2000bSmith and Cranfill 2002;Tzeng 2002;Cranfill and Kato 2003;Wang et al. 2003;Schneider et al. 2004a;Kuo et al. 2011) and morphological studies (Brown 1964;Kato 1975aKato , 1975bKato , 1977Kato , 1979Kato , 1984Kato and Darnaedi 1988;Sano et al. 2000b;Wang et al. 2004;Wang 2008). While 67 species is sparse coverage of the approximately 2600 species in the clade, our utilization of past results (both molecular and morphological/taxonomic) in selecting our taxon sample allows us a high degree of confidence that we have captured a great majority of the deepest branches, if not all of them. Most unsampled taxa are known to be deeply nested in crown clades, especially in the large genera Asplenium, Athyrium, Blechnum, Diplazium, and Thelypteris (sensu lato).
Wherever possible, we included generic and familial types, to facilitate future taxonomic revisions. Based on data from  and Liu et al. (2007), our broad outgroup sample included 10 representatives from the sister group to the Eupoly-pods II (Eupolypods I, see Fig. 2). To better evaluate the effect of uncertainty in outgroup placement on the ingroup topology and to better understand the divergence between Eupolypods I and II, we also included two representatives from each of the two potentially successive sister groups to the Eupolypods (Notholaena and Cryptogramma from Pteridaceae; Dennstaedtia and Pteridium from Dennstaedtiaceae; see Fig. 2 and Schuettpelz and Pryer 2007). Our total sample has 81 terminal taxa (Appendix 1).

Amplification and Sequencing
DNA was extracted from silica-dried or herbarium material, using either (i) a modified Carlson-Yoon protocol (<0.01 g dried plant material, silica beads, 750 μL Carlson buffer, and 20 μL mercaptoethanol added to a 2-mL tube and ground for 45 s using a Mini-Beadbeater (BioSpec Products), followed by incubation at 65˚C for 45 min; Yoon et al. 1991) or (ii) the protocol of Pryer et al. (2004) or (iii) the protocol of Kuo et al. (2011). For material extracted under the Carlson-Yoon protocol, the extracted DNA was purified by Illustra GFX PCR DNA and Gel Band Purification Kit (GE Healthcare).
Five plastid loci were selected for analysis: atpA, atpB, matK, rbcL, and the trnG-trnR intergenic spacer (henceforth "trnG-R"). All loci, except for matK, were amplified according to either the "standard" or "difficult" reaction protocols (below) depending on the source of the material (standard for most extractions; difficult for those from herbarium specimens greater than 10 years old), using the primers listed in Table 1. The standard amplification reaction used standard Taq polymerase with the following cycle: a 3 min initial denaturation at 95˚C, followed by 35 cycles of 30 s denaturation at 95˚C, 1 min annealing at 54˚C, and 2 min elongation at 72˚C, followed by a final elongation of 10 min at 72˚C. The difficult amplification reaction, using Phusion High Fidelity DNA Polymerase (Finnzymes), was 1 min initial denaturation at 98˚C, followed by 35 cycles of 10 s denaturation at 98˚C, 30 s annealing at 58˚C, and 1 min elongation at 72˚C, followed by a final elongation of 8 min at 72˚C. Amplification of all matK sequences followed the protocol of Kuo et al. (2011).
PCR products from Carlson-Yoon extractions were purified using MultiScreen Plates in a vacuum manifold (Millipore) and sequenced by Macrogen Inc. (South Korea). For material extracted under the protocol of Pryer et al. (2004), each PCR product was cleaned using 0.5 μL of exonuclease I and 1 μL of Shrimp Alkaline Phosphatase (USB, Cleveland, OH); reaction tubes were incubated at 37˚C for 15 min and then heated to 80˚C for 15 min to inactivate the enzymes, prior to sequencing on a Applied Biosystems 3730 xl at the Duke IGSP Sequencing Facility (Duke University). Material extracted under the protocol of Kuo et al. (2011) was sequenced at Genomics (Taipei, Taiwan). We completed our sampling with an additional 100 previously published sequences from GenBank (Appendix 1).

Alignment and Tree Search
Sequences were manually aligned in Mesquite 2.72 (Maddison and Maddison 2009). Indels (limited to matK, trnG-R, and the ends of the atpA alignment) were assessed by eye, and ambiguously aligned areas were excluded prior to phylogenetic analysis. Any gaps associated with unambiguous indel regions were treated as missing data. In one rapidly evolving region of the trnG-R alignment, we were unable to confidently align the Pteridaceae sequences to those of the other taxa. In order to retain this otherwise unambiguous region, we excised those portions of the Pteridaceae sequences, replacing them with question marks.
To evaluate congruence among our loci, we performed maximum likelihood (ML) tree searches on 1000 bootstrap data sets for each locus individually, under a GTR+I+G model using the default settings in Garli v1.0.695 (Zwickl 2006, see SI Table 1, available from http://www.sysbio.oxfordjournals.org). The majority rule bootstrap consensus trees from each locus were manually compared and examined for strongly supported conflicts (Mason-Gamer and Kellogg 1996), after which we concatenated the full data with abioscript (Larsson 2010), producing a single annotated fivelocus data set, with excluded regions removed. This alignment is largely complete (361 of the possible 405 sequences are present, for an average of 4.5 loci per terminal taxon) and contains 13.3% missing data and 6595 characters, of which 3641 are variable (Table 2). Our alignment is available on TreeBase (accession number S11464); the full length of all newly generated sequences (including any portions excluded prior to analysis) are deposited in GenBank (see Appendix 1).
To obtain a point estimate of the phylogeny, we performed 10 tree bisection-and-reconnection heuristic searches of the concatenated (unpartitioned) data, each from a different random-addition-sequence starting tree, under ML using a GTR+I+G model as implemented in PAUP* 4.0b10 for Unix (Swofford 2002). The values for the exchangeability parameters, base frequencies, gamma shape parameter, and proportion of invariant sites were fixed at their ML values as optimized under a Garli 0.951 (Zwickl 2006) tree search, using default genetic algorithm and termination settings.  Notes: Missing data include both uncertain bases (e.g., ?, R, Y) and gaps (-). Support values are listed as MLBS or Bayesian PPs.
We assessed support using ML bootstrapping and Bayesian inference. For the ML bootstrapping, we performed 5000 replicates in each of PAUP* 4.0b10 for Unix (Swofford 2002), Garli 1.0.695 (MPI parallel version; Zwickl 2006), and RAxML v7.2.6 (Stamatakis 2006). The PAUP* runs were performed with the parameters optimized as above, reconnection limit set to eight ("reconlim = 8"), and with only a single random-additionsequence per bootstrap replicate. In Garli, we ran 5000 bootstrap replicates on the concatenated data, under a GTR+I+G model using the default genetic algorithm and termination settings. In RAxML, we ran 5000 bootstrap replicates on the data partitioned by locus, with each locus assigned a GTR+G model. For Bayesian inference, we ran four runs of four chains each (one cold; three heated), for 15 million generations, under a partitioned GTR+I+G model in the parallel version of MrBayes v3.1.1 (Ronquist and Huelsenbeck 2003;Altekar et al. 2004). Each of the five loci was assigned its own partition, with substitution parameters unlinked among partitions, and branch lengths linked (with a proportionality parameter to account for rate heterogeneity among partitions); the posterior was sampled every 1000 generations. Visual inspection in AWTY (Wilgenbusch et al. 2004;Nylander et al. 2008) revealed that the runs converged within the first 500,000 generations. To be conservative, we excluded the first 2 million generations of each run as burn-in prior to summarizing the posterior. The posterior thus comprised 52,000 samples (13,000 post-burn-in samples from each of the four runs).

Phylogeny Evaluation
As stated above, earlier studies (e.g.,  indicate that the eupolypod II phylogeny is likely to include several key challenges for phylogenetic inference, specifically a series of long branches among very short backbone internodes (an ancient rapid radiation), lineage-specific rate heterogeneity, and a distantly related outgroup. Given these concerns, we sought to explicitly evaluate our topology and support values against these potential artifacts, with particular emphasis on the support values along the backbone of the tree.
Our approach to phylogeny evaluation involved permutations of both the models and the implementation of those models (i.e., programs). The models were deliberately selected according to their varying degrees of susceptibility for each of the risk factors in question, in an attempt to isolate potential model-based biases. The choice to additionally vary the programs used was in part due to constraints of implementation-no single program offered all the models we wished to compare. This approach has the added benefit of demonstrating a further level of robustness: if our phylogenetic results are insensitive to both the differing models and the myriad of incompletely quantified differences among programs, we can be all the more confident in our conclusions. Additionally, varying both the models and their implementation more closely matches the options available to empirical phylogeneticists seeking to resolve ancient rapid radiations. This approach suffers a clear liability, however, in that the effects of model differences and implementation differences are conflated. In the event of differing results, we may not be able to isolate the effects of one from the other; therefore, the added value to empirical phylogeneticists comes at the cost of reduced utility of our results to program developers and theorists.
The specific evaluations performed are described more thoroughly in the Results section. Computationintensive analyses were run on either the Duke Shared Cluster Resource (https://wiki.duke.edu/display/SC SC/DSCR) or the Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX). When appropriate, multiple tree files were summarized onto a target phylogram with SumTrees 2.0.2 (Sukumaran and Holder 2010) for subsequent inspection or manipulation in FigTree 1.3.1 .

Data and Topology Point Estimate
Tree-wide mean ML bootstrap support (MLBS) values (summed bootstrap values on the ML tree divided by the number of internodes in that tree) for the individual loci ranged from 74% (atpB) to 84% (matK). The concatenated data have a mean MLBS of 92% and strongly support (i.e., have 70% MLBS and 0.95 posterior probability [PP]) 90% of the partitions (Table 2). Across data sets, ML and Bayesian inference consistently inferred strong support for a comparable number of bipartitions (Table 2), a result that offers further empirical corroboration for the approximate equivalence of 70% MLBS and 0.95 PP (Hillis and Bull 1993;Alfaro et al. 2003). pruned from trees, MLBS with Aspleniaceae removed from analysis, and posterior support from BEAST. c) Controlling for the Bayesian star-tree paradox artifact using the polytomy prior in Phycas. The four PPs listed for each internode are from: MrBayes 3.1.1 (susceptible to the star-tree artifact), Phycas with C = 1, Phycas with C = e, and Phycas with C = 10. d) Controlling for the Bayesian star-tree paradox artifact using the Yang branch-length prior. The four PPs listed for each internode are from: MrBayes 3.1.1 (susceptible to the star-tree artifact), MrBayes 3.1.2 with branch-length prior mu1/mu0 = 0.01, MrBayes 3.1.2 with branch-length prior mu1/mu0 = 0.001, and MrBayes 3.1.2 with branch-length prior mu1/mu0 = 0.0001.
There are two well-supported conflicts among the individual-locus ML trees. The first involves a tip relationship (matK unites Deparia pterorachis with D. unifurcata, with 72% MLBS, whereas rbcL places D. pterorachis as sister to the rest of the genus, with 75% support) that is peripheral to the focus of this study. The second is deeper in the tree: matK has 80% MLBS for a clade uniting Thelypteridaceae with the athyrioids, Blechnaceae, and Onocleaceae, to the exclusion of Woodsia and allies, whereas both atpA and atpB place Woodsia and its allies within that clade (92% in atpA; 71% in atpB). Given that we confirmed this conflict to not be attributable to laboratory or identification errors, and because the loci involved are linked and the taxa are long-diverged, we do not ascribe this conflict to differences in evolutionary history and proceeded to concatenate all the data for subsequent analyses.
Each of our ten ML best-tree searches of the concatenated data in PAUP* (from different random-additionsequence starting trees) inferred the same tree (Fig. 3a), suggesting that tree space is unimodal for our data set. Most partitions in this topology point estimate were strongly supported by both ML bootstrapping and Bayesian PPs (Table 2); the different ML programs (PAUP*, Garli, and RAxML) inferred very similar MLBS levels (data not shown).

Bayesian Star-Tree Paradox Artifact
For certain branches, we observed very high Bayesian PPs from the MrBayes analysis, but much lower levels of support from the ML bootstrapping (Fig. 4a); these support discrepancies are disproportionately represented among short branches (Fig. 4b). This pattern is consistent with artificially high Bayesian support due to the star-tree paradox artifact-most implementations of Bayesian phylogenetic inference do not consider polytomies among the option set for the MCMC sampler and thus can return high PPs for branches that are unsupported by the data ; Yang and VOL. 61 To ensure that this "star-tree paradox" artifact was not influencing our assessment of support, we compared the results from our original MrBayes 3.1.1 analysis (Ronquist and Huelsenbeck 2003; MrBayes 3.1.1 is potentially vulnerable to the star-tree paradox artifact) with those of a non-Bayesian analysis (ML bootstrapping in PAUP* from our initial assessment of support), as well as with two implementations of Bayesian inference that use different approaches to reduce their vulnerability to the star-tree paradox artifact.
First, we analyzed our data with Phycas 1.1.2-r (Lewis et al. 2010). Phycas uses reversible-jump MCMC to allow the sampling of incompletely resolved topologies, controlled via the incorporation of a polytomy prior, "C." A value of C = 1 means that unresolved and fully resolved topologies are weighted equally; under a value of C = 10 a trichotomy is 10 times more likely, a priori, than either of its fully resolved resolutions . We performed three runs of our full concatenated data, under a GTR+I+G model, using a branchlength hyperprior (default values), for 200,000 cycles (sampling from the posterior every 10 cycles; note: Phycas makes proposals to each free parameter in each cycle, and thus Phycas cycles are not comparable to Mr-Bayes generations, which include a proposal to only a single parameter), with C = 1, C = e (2.718), and C = 10, respectively. Inspection of the AWTY-type plots (see Nylander et al. 2008) and sojourn plots (see  revealed that the runs converged before 40,000 cycles; to be conservative, we excluded the first 50,000 cycles (5000 samples) as burn-in. These trees are available as SI Figures 1-3.
We then reanalyzed our data with a modified version of MrBayes 3.1.2 that incorporates exponential priors on internal and external branch lengths (Yang 2007(Yang , 2008. These "Yang branch-length priors" allow the concentration of the prior mass on topologies with very short internal branches, for an intended effect similar to that of the polytomy prior (above); only if the data strongly support a branch can the short internal branch prior be overcome. We performed three runs of four chains each under the settings used in the initial MrBayes analysis but with the addition of the branch-length priors. The mean of the external branch-length prior (mu1) was set to 0.1, and the mean of the internal branch-length prior (mu0) to 0.00001, 0.0001, and 0.001, successively. As in the initial analyses, 2 million generations from each run were discarded as burn-in (trees available as SI Figs. 4-6).
The results of these analyses show that four of the backbone internodes (A, B, C, and D; Fig. 3a) were largely insensitive to either the polytomy prior (Fig. 3b) or the Yang branch-length prior (Fig. 3c); their PPs stayed at or within three percentage points of 1.0 for all seven analyses (the original MrBayes 3.1.1 analysis, three Phycas runs with increasingly strong polytomy prior values, and three MrBayes 3.1.2 runs with increasingly strong Yang branch-length prior values). Interestingly, the two approaches (polytomy prior vs. branch-length prior) had very different effects on the other backbone internodes, despite the approaches being designed to overcome the same shortcoming in Bayesian phylogenetic inference. For example, the only backbone internode that was not well supported by the original ML and Bayesian analyses (internode F; Fig. 3a) exhibited increased support under weak versions of either the polytomy prior or the branch-length prior. However, the norms of reaction for the two priors were opposed: as the polytomy prior increased in strength (C = 1, e, and 10), the posterior support for internode F decreased (PP = 0.93, 0.89, 0.83; Fig. 3b), whereas as the branch-length prior strength increased (mu0/mu1 = 0.01, 0.001, 0.0001), the PP on internode F also increased (PP = 0.97, 1.0, 1.0; Fig. 3c).
The remaining three internodes (E, G, and H; Fig. 3a) were well supported by the original ML and Bayesian analyses but showed some sensitivity to either the polytomy or the branch-length prior, again in opposing ways. Internode E was largely insensitive to the branchlength prior (Fig. 3c) but was strongly weakened by the polytomy prior (Fig. 3b), whereas internodes G and H were largely unaffected by the polytomy prior but were unsupported under strong values of the branch-length prior (Fig. 3b,c).

Lineage-Specific Rate Heterogeneity
To investigate whether the rapid rate of evolution for the Aspleniaceae (Figs. 1 and 3) was biasing tree reconstruction, we attempted to isolate the effects of this rate heterogeneity in three ways. First, we pruned the Aspleniaceae from 1000 full-data Garli ML bootstrap trees prior to building the consensus tree and evaluating support. This "reduced consensus" approach (Wilkinson 1996;Burleigh et al. 2009) removes any effects due solely to uncertainty in the placement of these long-branch taxa. If the remaining relationships are well supported, then overall support values will appear low in the standard consensus but will be restored under the reduced consensus. Second, we reran the Garli ML bootstrap analysis on a data set where the Aspleniaceae had been removed prior to analysis, to eliminate any effect that these taxa might have on the optimization of model parameters, and to allow the model to better fit the remaining data (the "reduced data" approach). Third, we ran the full data set in BEAST 1.5.4 , incorporating a relaxed-clock model that explicitly models lineage-specific rate variation  and thus should be less sensitive to any artifacts induced by the strongly heterogeneous rates in our data. We ran three independent runs on the full concatenated data set, each for 20 million generations (sampling the posterior every 1000 generations), with the following settings: birth-death tree prior; lognormal uncorrelated relaxed clock; and GTR+I+G substitution model. Priors were left at their default values, with the exception of those for six time-to-most-recent-common-ancestor (TMRCA) age parameters, which were each given normal distributions with a mean equal to the inferred age estimated for that clade by Schuettpelz and Pryer (2009) and a standard deviation equal to 10% of that mean. The relevant taxon sets, and their TMRCA prior means, are: tree root (165.6 MA), Dennstaedtiaceae (119.3 MA), Eupolypods (116.7 MA), Pteridaceae (110.8 MA), Eupolypods II (103.1 MA), and Eupolypods I (98.9 MA). None of the taxon sets was constrained to be monophyletic. The use of secondary constraints such as these is clearly inferior to the use of fossil data for divergence time dating (Magallón 2004), but as no such data are available, and our interest is more in the relative than absolute divergence times, this approach seemed best. Visual inspection in Tracer  demonstrated that the runs converged before 1 million generations; to be conservative, we excluded the first 3 million generations of each run prior to analyzing the pooled posterior of 51,000 samples (17,000 from each run; SI Fig. 7). For this sample, the effective sample size for each parameter was above 300.
None of these attempts to mitigate potential effects of the increased rates of molecular evolution associated with Aspleniaceae strongly affected support values along the backbone. Support values from the full taxonsample consensus data (Fig. 3d, first values) differed from those from the reduced consensus (Fig. 3d, second values) by at most one percentage point. Removing Aspleniaceae from the data set prior to the bootstrap tree searches had a larger effect (up to a five percentage point change in support; Fig. 3d, third values), but in no case resulted in an internode moving from well supported ( 70% MLBS) to poorly supported or vice versa. The support values from the BEAST analysis (Fig. 3d, fourth values) were concordant with those of the ML runs, especially in that the internode uniting Rhachidosorus with the diplazioids, Hemidictyum, and Aspleniaceae (internode F; Fig. 3d) was the only one without strong support (it had a PP of 0.90).

Rooting Uncertainty
To evaluate any effects that uncertainty in root-branch placement might have on apparent levels of support within the ingroup, we compared ingroup backbone MLBS values from the analysis of our complete data (full outgroup) with those from each of six different variations in outgroup composition: (i) ingroup only; (ii) ingroup + Dryopteris; (iii) ingroup + Dryopteris and Didymochlaena; (iv) ingroup + Dryopteris and Notholaena; (v) ingroup + Dryopteris and Notholaena and Pteridium; and (vi) ingroup + our full eupolypod I sample (Fig. 5). This outgroup sampling regime was selected to successively bisect the longest outgroup branches, with a particular emphasis on breaking the proximate root branch (the branch connecting the ingroup to the first outgroup node).
We evaluated support for each of the six outgroup sampling regimes using both the reduced consensus approach (full data included in the analysis, but outgroup pruned down to the desired sample prior to forming the consensus tree; Fig. 5, first values; Wilkinson 1996;Burleigh et al. 2009) and a reduced data approach (outgroup reduced to the desired sample prior to the analyses; Fig. 5, second values). The former approach controls for uncertainty in outgroup placement alone (i.e., it offers a metric of the signal strength), whereas the latter approach additionally accounts for model fit. All analyses were based on 1000 ML bootstrap replicates of the concatenated data in Garli 1.0.695 (MPI parallel version;Zwickl 2006), under a GTR+I+G substitution model, using the default genetic algorithm and termination settings. The results of these rooting comparisons demonstrate that our initial concerns-that the outgroup would wander and thus reduce support measures within the ingroup-were largely unfounded. The reduced consensus support values were minimally different from those with the full outgroup (Fig. 5, first values). When a "consensus interference" effect did appear (first values > 0 in Fig. 5b,c,f), it was correlated with the maximum root length rather than with proximate root length, i.e., it is the long Notholaena branch that wanders rather than the outgroup as a whole.
In stark contrast, outgroup composition had a strong effect on backbone support if the outgroup was changed prior to the tree searching steps. When we reduced our outgroup sample and reran the ML bootstrapping (the reduced data approach), backbone internode support values changed from their full-outgroup values by up to 32 percentage points (Fig. 5, second values). The largest of these changes are reductions in support for branches proximate to the root (internodes E, F, and G; Fig. 3a) and are due to uncertainty in the ancestral state of the smaller outgroup sample (as demonstrated by the reduced consensus values from each of the smaller data sets; data not shown).

Eupolypod II Phylogeny
Our results demonstrate that the vast majority of internodes in the ML tree are strongly supported by both ML bootstrapping and Bayesian PPs ( Fig. 6 and Table 2), and these support values proved robust to our phylogeny evaluations. In particular, the ML tree has 10 highly supported major (approximately family unit) ingroup clades: Cystopteris/Gymnocarpium; Rhachidosorus; Diplaziopsis/Homalosorus; Hemidictyum; Aspleniaceae; Thelypteridaceae; Woodsia and allies; Onocleaceae; Blechnaceae; and the athyrioids (Fig. 6b). Of these, Cystopteris/Gymnocarpium is sister to the remaining eupolypod II taxa, followed by the loosely supported assemblage of Rhachidosorus with Diplaziopsis/Homalosorus + Hemidictyum + Aspleniaceae. Blechnaceae is sister to Onocleaceae, and they together are successively sister, in a pectinate pattern, to the athyrioids, then to Woodsia and allies, and finally to the Thelypteridaceae (Fig. 6). This broad phylogeny is in general agreement with earlier molecular phylogenetic studies that included members of the Eupolypods II (Gastony and Ungerer 1997;Murakami et al. 1999;Sano et al. 2000aSano et al. , 2000bSmith and Cranfill 2002;Tzeng 2002;Cranfill and Kato 2003;Wang et al. 2003;Schneider et al. 2004a;Wei et al. 2010; see particularly Kuo et al. 2011). However, the backbone of the phylogeny is strongly supported for the first time; the only backbone internode lacking such support is the one attaching the Rhachidosorus branch to the rest of the tree (Fig. 6, internode F). Additionally, we are finally able to confidently place the enigmatic genera Cheilanthopsis, Diplaziopsis, Homalosorus, Protowoodsia, and Woodsia.

Bayesian Star-Tree Paradox Artifact
Although one would anticipate that internodes across a topology would differ in their sensitivity to the startree paradox approaches (not all short branches are inferred equal), it is unclear what is driving the different responses in our data-neither the original Bayesian posteriors nor the MLBS levels correlate with the behavior of a given internode under the additional priors (Fig. 3). This study is the first to examine the performance of these star-tree paradox methods on empirical data; their nonparallel effects were perhaps the most surprising result of this portion of the analyses. However, while they were developed for the same function, the methods differ strongly in their approach and should not be expected to result in similar behavior. The branch-length prior, in effect, flattens the posterior for topologies. As the mu0/mu1 ratio decreases, the relative influence of any data supporting an internal branch is reduced, and the external branches come closer to being randomly arranged. However, each topology sampled from the posterior must still be fully resolved and thus any reduction of support for a particular topology must be accompanied by increased support for some other one. In this sense, the branch-length prior is less a measure of intrinsic support for a given internode than it is a measure of whether that node is better supported than all alternative resolutions. Under circumstances of low support for an entire set of relationships, the branchlength prior favors the best of a bad lot. The polytomy prior, on the other hand, allows the direct comparison between a given resolution and a polytomy. Strong values of the branch-length prior lead to many trees, each with a low posterior, whereas strong values of the polytomy prior lead to a star tree, with a high posterior.
In interpreting the performance of these two methods on our data, it is important to stress that we did not attempt to tightly isolate the effects of the Yang branch-length prior and those of the polytomy prior. Rather, each was bundled with other elements of its host program (MrBayes 3.1.2 and Phycas 1.1.2-r, respectively), which differ from each other in both their models, and their implementation of those models. In particular, important model differences include data partitioning in MrBayes (the Phycas runs were on unpartitioned data) and the incorporation of a branch-length hyperprior in Phycas (there is no such hyperprior in MrBayes); important implementation differences include the limitation of Phycas to Larget-Simon moves (Larget and Simon 1999), whereas MrBayes utilizes a broader suite of topology proposals.
Regardless of the different performance of the two methods, the backbone support levels in our data were generally robust to the star-tree paradox artifact approaches (Fig. 3b,c), suggesting that the high Bayesian support values for these internodes are valid. Even under extreme values of the polytomy prior, for example (C = 10, or trichotomies 10 times more likely, a priori, than their fully resolved alternatives), the posterior consensus tree still resolved each of the eight critical backbone internodes, and only one fell below 0.95 PP (internode E; Fig. 3b). The differences between Bayesian PP and MLBS values in our data (Fig. 4), therefore, reflect something other than the failure of the original Bayesian analyses to include polytomies in the option set; these differences may simply be due to Bayesian inference being more sensitive to small amounts of data than is bootstrapping and thus more likely to support short internodes (Alfaro et al. 2003).

Lineage-Specific Rate Heterogeneity
The absence of any effect of lineage-specific rate heterogeneity on our topology estimation or support levels is particularly interesting in light of recent questions ) about the general applicability of the unrooted model (aka "no clock" model; Yang and Rannala 2006;Wertheim et al. 2010) in phylogenetic inference. Given the dramatic lineage-specific rate heterogeneity that is present in our data set, one might expect the unrooted and relaxed-clock models to fit very differently, and given that the fast lineages in our data are intercalated among short internodes, our topology would be expected to be sensitive to such model differences. However, no effects are seen; our data, at least, do not support concerns about the application of the unrooted model in phylogenetic inference, a result that provides empirical support to the simulation results of Wertheim et al. (2010).

Rooting Uncertainty
The effects of differing outgroup compositions on support levels for branches phylogenetically distant from the root were unexpected and may reflect a combination of both stochastic variation in ML bootstrapping and of factors of model optimization on the different data sets. Neither of these explanations is heartening. The latter-the "model-mediated" effect-requires that changes in outgroup composition have strong and somewhat idiosyncratic effects on support levels on parts of the tree phylogenetically distant from the root itself (for a similar case, see Roberts et al. 2009). These effects are not due to the outgroup itself changing position (that possibility is eliminated by comparison with the reduced consensus values) and must instead be mediated through nontopological factors. The former explanation-stochastic variation in bootstrap support values-would suggest that 1000 pseudoreplicates are insufficient to get accurate support estimates for these data. Regardless of the precise mechanism by which the outgroup affects support values, these results emphasize the wisdom of including a broad outgroup sample, particularly when the outgroup is distantly related to the taxa of interest (Swofford et al. 1996;Graham and Iles 2009).

Eupolypod II Phylogeny: Major Clades
The affinities of Cystopteris s. lat. (including Acystopteris, e.g., Blasdell 1963) and Gymnocarpium have been the object of considerable taxonomic disagreement. Both genera, individually or in tandem, have been thought to be allied with the Dryopteridaceae (in Eupolypods I) or the Athyriaceae; in either position they were inevitably highlighted as being anomalous (see Sledge 1973). Early molecular studies supported Cystopteris and Gymnocarpium as sister species and demonstrated their lack of close affinity to either Dryopteridaceae or Athyrium but were unable to pinpoint their phylogenetic position Hasebe et al. 1995). Recent studies Kuo et al. 2011) were the first to support a sister group relationship between a Cystopteris/Gymnocarpium clade and the rest of Eupolypods II, a result corroborated and strengthened by our data (Fig. 6, internode E), the first to include multiple accessions of Acystopteris and Cystopteris s. str.
Historically, arguments about Rhachidosorus focused on its validity as a genus, distinct from either Athyrium or Diplazium (Ching 1964a;Kato 1975b). Early molecular phylogenies Tzeng 2002;Wang et al. 2003) provided the first evidence that Rhachidosorus might not be closely related to either, a result further emphasized by the three-gene results of Kuo et al. (2011). In our study, the two included Rhachidosorus species form a tight clade phylogenetically distant from any other taxon; their closest relatives appear to be Diplaziopsis, Homalosorus, Hemidictyum, and Aspleniaceae. While our data do not strongly support a precise position for Rhachidosorus (Fig. 6, internode F), of note is the 100% MLBS and 1.0 PP for internode D (Fig. 6), which separates Rhachidosorus from the athyrioids. Thus, our data very strongly reject a close relationship between Rhachidosorus and its presumed allies, the athyrioids, an unanticipated conclusion based on morphological data (Kato 1975b). Indeed, our data suggest that the two groups last shared a common ancestor nearly 100 MA (SI Fig. 7).
As with Rhachidosorus, Homalosorus and Diplaziopsis were long thought to be allied with the athyrioids, where they are typically treated as members of Diplazium (Ching 1964b;Kato 1975aKato , 1977Kato and Darnaedi 1988;Wang et al. 2004). The first indication that this placement might be inaccurate came from the study of Sano et al. (2000a), which strongly supported Homalosorus (a monotypic genus) as sister to their lone Diplaziopsis accession and placed the two genera distant from Diplazium, a result corroborated by Wei et al. (2010) and Kuo et al. (2011). Our study includes two Diplaziopsis species, which are strongly supported as sister to each other, and together are sister to Homalosorus. These two genera form a clade that is strongly supported, for the first time, as sister to Hemidictyum + Aspleniaceae (Fig. 6, internode G).
Woodsia has been underrepresented in molecular phylogenetic studies to date; no study has included more than one species, and none has been able to strongly infer the position of that species, either. Here, we establish that Woodsia s. lat. is likely to be monophyletic (Fig. 6, seven species included in our analysis) and we demonstrate that two of the three segregate genera (Cheilanthopsis and Protowoodsia) recognized by Shmakov (2003) are nested within Woodsia s. str.; only Hymenocystis is as-yet unsampled. Additionally, our study finds strong support for the position of Woodsia s. lat. to be far from the other Woodsiaceae genera (Athyrium, Acystopteris, Cornopteris, Cystopteris, Deparia, Diplazium, Gymnocarpium, Rhachidosorus, Diplaziopsis, Homalosorus, Hemidictyum), as circumscribed in the most recent family level fern classification (Smith et al. 2006); compare Figure 6a with Figure 6b.
The "athyrioids" have been a source of great disagreement in fern systematics (e.g., Ching 1940;Alston 1956;Ching 1964a;Sledge 1973;Tryon and Tryon 1982). Molecular data confirmed their distant relationship to the dryopteroid ferns (Dryopteridaceae, in Eupolypods I), but uncertainty regarding their delimitation and affinities has persisted until very recently. Sano et al. 2000a were the first to extensively sample the athyrioids, and they provided the initial evidence that the group, as then understood, was strongly heterogeneous. Our data corroborate the results of earlier studies Wang et al. 2003; in revealing three major clades within the athyrioids s. str.: one containing Athyrium and close allies ("athyriids"); one containing Diplazium s. lat. ("diplaziids"); and one containing Deparia s. lat. ("depariids"). Our novel finding is the well supported, early diverging position of Athyrium skinneri with respect to the other athyriids included in our sample. This species belongs to a small group of predominantly Mexican taxa, none of which had been included in previous phylogenetic studies. Its position as sister to the rest of the included athyriids (including Cornopteris and Pseudocystopteris) emphasizes the paraphyly of Athyrium as currently circumscribed and has important implications for our understanding of the evolution of both the athyriids and the diplaziids. Our study provides additional novel 502 SYSTEMATIC BIOLOGY VOL. 61 support for the placement of the athyrioids as phylogenetically distant from Rhachidosorus, Cystopteris, Gymnocarpium, Woodsia, Diplaziopsis, Hemidictyum, and Homalosorus, a topology that is in conflict with the recent classifications of the group (Wang et al. 2004;Smith et al. 2006); both Athyriaceae sensu Wang et al. (2004) and Woodsiaceae sensu Smith et al. (2006) are shown here to be strongly paraphyletic (Fig. 6).

Eupolypod II Phylogeny: Morphological Stasis and Disparity
A striking pattern in our phylogeny is its incongruence with previous morphology-based hypotheses of relationship, particularly with respect to the position of the genera of Woodsiaceae sensu Smith et al. (2006): Acystopteris, Cystopteris, Diplaziopsis, Gymnocarpium, Hemidictyum, Homalosorus, Rhachidosorus, Woodsia and allies, as well as the athyrioids (Fig. 6). Some of these groups have been historically difficult to place and thus their isolation from Woodsia or the athyrioids (the bulk of Woodsiaceae sensu Smith et al. (2006) is in the athyrioids) is not particularly surprising. Smith et al. (2006) themselves noted that their Woodsiaceae might prove to be not monophyletic. The placement of three genera, however, was utterly unanticipated by morphological data: Diplaziopsis, Homalosorus, and Rhachidosorus. These taxa have not only been considered closely related to the athyrioids, they have been nearly universally considered to be members of the large genera Diplazium (first two) or Athyrium (Rhachidosorus). Their phylogenetic position, deeply isolated from their presumed relatives, underscores the complex patterns of morphological evolution in Eupolypods II; further morphological investigations are necessary to determine whether the apparent similarities between these three genera and the athyrioids are due to convergence or symplesiomorphy.
This trend of shared morphological syndromes across very deep splits in the tree by some members of the "Woodsiaceae" is in contrast to the interdigitation, among those same taxa, of a series of distinct morphologically unique groups, including the Aspleniaceae, Blechnaceae, Onocleaceae, and Thelypteridaceae. The coarse picture of eupolypod II morphological evolution, then, is marked by two seemingly opposing patterns. On the one hand are the autapomorphy-rich clades, whose individual phylogenetic coherence is strong, but whose deep relationships were obscure based on morphological data. And, on the other, the morphologically consistent yet phylogenetically incoherent members of the "Woodsiaceae".
Although not the focus of this study, our phylogeny contains rich information on relationships closer to the tips of the tree, within the approximately family unit clades. For example, within the athyrioids and Blechnaceae, morphological evolution is complex, and nonmonophyletic generic concepts are common. Generic delimitation within these families is in need of much further study. In addition, a cursory comparison between the Onocleaceae and their sister group, the Blechnaceae, is revealing. Both clades have approximately the same crown ages (SI Fig. 7) yet exhibit strikingly different patterns of diversification. The Onocleaceae branch is marked by few well-spaced divergences leading to the five extant species. Conversely, the Blechnaceae branch features multiple, very short internodes; this family includes approximately 200 extant species.

Phylogeny Evaluation
Despite the presence in our data set of each of the anticipated challenges to robust phylogenetic inference (long outgroup branch; strong lineage-specific rate heterogeneity; ancient rapid radiation model; Figs. 1 and 3a), we were able to infer a phylogeny with strong backbone support (Fig. 6), and our various evaluations gave no indication that the support for the internodes in our ML tree is due to artifacts. However, different approaches to controlling for the Bayesian star-tree paradox artifact, and different outgroup sampling regimes all influenced support levels; only lineage-specific rate heterogeneity had a negligible effect.
These effects give further weight to arguments for rigorously evaluating phylogenies against potential artifacts. While specific vulnerabilities may be data set dependent, the core elements of our analysis regime are broadly applicable, including the inspection of preliminary phylogenetic hypotheses for potential confounding factors, the investigation of those factors through scrutinizing the performance of multiple models and multiple implementations of those models, and the utilization of the reduced consensus approach to isolate topological effects of signal weakness from those of signal conflict. Although this study is focused on the post-data set steps, preanalysis components (taxon sampling, character sampling, character evaluation) are also vital. In particular, in our case, the use of a broad taxon sample with moderate character data proved effective.