Phylogeny of the plant genus Pachypodium (Apocynaceae)

Taxon sampling

We generated new ITS and trnL-F sequences from 56 Pachypodium samples representing all 27 minimum-rank taxa (species and subspecies) in the most recent revision of the genus (Lüthy, 2004; Tables 1 and 2). An additional ITS sequence was generated for Funtumia africana—a close relative of Pachypodium (Livshultz et al., 2007)—for use in rooting the ITS tree. Pachypodium and Funtumia tissues for DNA analysis were taken from greenhouse or garden plants (Appendix 1). Tissues were obtained by D. Burge, Walter Röösli, Nicholas Plummer, and Anurag Agrawal. Specimens were identified by W. Röösli, N. Plummer, or D. Burge according to the taxonomic revision of Lüthy (2004) and subsequent descriptions of new taxa (Lüthy, 2005; Lüthy & Lavranos, 2005). Plants were selected based on geographic distribution, with a larger amount of sampling for widespread taxa. Between one and eight populations of each taxon were used (Tables 1 and 2). Additional non-Pachypodium trnL-F sequences, for rooting trees, were obtained from GenBank (F. africana [EF456206], Holarrhena curtisii [EF456122], Kibatalia macrophylla [EF456119], Malouetia bequaertiana [EF456243], and Mascarenhasia lisianthiflora [EF456174]). These taxa were selected on the basis of their close relationship with Pachypodium (Livshultz et al., 2007).

Molecular methods

Total genomic DNA was extracted from silica-dried leaves or seeds using the DNeasy Plant Mini Kit (Qiagen, Germantown, MD) according to the manufacturer’s instructions. For seeds, up to three excised embryos from a single parent plant were pooled prior to DNA extraction (Burge & Barker, 2010). DNA was extracted from seeds when silica-dried material for the same plant was not available, or proved recalcitrant to extraction of high quality DNA. All polymerase chain reactions were performed using Qiagen Taq DNA Polymerase. Amplifications were performed using an initial incubation at 94°C for 10 min and 30 cycles of three-step PCR (1 min at 94°C, 30 s at 45°C, and 2 min at 72°C), followed by final extension at 72°C for 7 min. PCR was performed on a Perkin Elmer GeneAmp thermocycler. The primers ITS4 (White et al., 1990) and ITSA (Blattner & Kadereit, 1999) were used to amplify the ITS1-5.8S-ITS2 region of the nuclear ribosomal DNA. Primers ‘c’ and ‘f’ (Taberlet et al., 1991) or a combination of these with internal primers ‘d’ and ‘e’ were used to amplify the trnL-F chloroplast region. For some plants, sequencing of ITS was problematic as a result of variation in length among ITS copies present in individual plants. Consequently, cloning of the ITS region was required for some plants. Cloning was carried out using the pGEM-T Easy Vector kit (ProMega, Madison, WI) according to the manufacturer’s instructions. NIA inserts were amplified directly from up to four positive colonies using the PCR protocol described above. For all PCR reactions, excess primer and dNTPs were removed using exonuclease I (New England Biolabs, Ipswich, MA [NEB]; 0.2 units/µl PCR product) and antarctic phosphatase (NEB; 1.0 unit/ µl PCR product) incubated for 15 min at 37°C followed by 15 min at 80°C. For sequencing we used Big Dye chemistry (Applied Biosystems, Foster City, CA) according to the manufacturer’s instructions. Sequences were determined bidirectionally on an Applied Biosystems 3100 Genetic Analyzer at the Duke University Institute for Genome Science and Policy Sequencing Core Facility.

Sequence editing and alignment

All sequences were assembled and edited in Sequencher 4.1 (Gene Codes Corporation). In the case of the five plants for which ITS was cloned, we assessed sequence variation using an alignment of cloned sequences (hereafter ‘isolates’). Two plants yielded pools of identical isolates (P011 and P021, Appendix 1), one yielded four different types of isolate (P053), and two were represented by a single successfully cloned isolate (P046 and P048). For the plant with more than one isolate type (P053), we included all four isolates in the phylogenetic analyses of ITS; for the plants with identical isolates, we selected a single isolate to represent each plant. New ITS and trnL-F sequences for Pachypodium were deposited in GenBank (Appendix 1).

The 60 new ITS and 55 new trnL-F sequences, along with additional outgroup sequences from GenBank, were used to create separate alignments for the two regions (Table 3; Alignments S2 and S3). Sequences were aligned in MUSCLE (Edgar, 2004) under default settings. For ITS, several indel- and repeat-rich regions (54 bp total) were excluded due to alignment ambiguity. A portion of trnL-F not available for some taxa (the 3′ trnL intron) was recoded as missing data. Indels were not recoded for analysis.

Following individual alignment of ITS and trnL-F, we endeavored to create a combined alignment. Preliminary analyses showed that for the single Pachypodium sample represented by more than one cloned ITS sequence (P053; Table 1), the four sequences formed a monophyletic group. Thus, a single sequence from this group was selected at random. For the final combined alignment (Alignment S2), the entire trnL-F region was coded as missing data for the two Pachypodium samples from which trnL-F was not obtained (P. bispinosum A049 and P. brevicaule subsp. leucoxanthum, P066). To test for conflict between the nuclear (ITS) and chloroplast (trnL-F) portions of the alignment, we used the incongruence length difference test (Farris et al., 1995), implemented in PAUP* v 4.0 (Swofford, 2002) as the partition homogeneity test. The test used 1000 random repetitions of the parsimony analysis described below (see Phylogenetic analyses). Results showed significant disagreement between ITS and trnL-F (P = 0.047; 953/1000 trees). To account for this conflict, we ran all of our phylogenetic analyses on the separate trnL-F and ITS alignments, noting any well-supported conflicts between the results, and compared these to results from the combined alignment (see Discussion).

Phylogenetic analyses

Trees were reconstructed using Bayesian, maximum likelihood (ML), and maximum parsimony (MP) techniques. Bayesian analyses were carried out based on the best fit model of evolution from jModelTest 2, under default parameters (Posada & Crandall, 1998; Guindon & Gascuel, 2003; Darriba et al., 2012; ITS: GTR +I +G; trnL-F: GTR +I). Bayesian sampling was performed in MrBayes v 3.2.1 (Ronquist & Huelsenbeck, 2003), using the models of sequence evolution identified by jModelTest 2; all other parameters of MrBayes were left at default values; for the combined tree, no rate or model constraints were imposed between the two partitions. Analyses were carried out as follows: (1) three separate runs of 1 × 10⁷ MCMC generations, sampling every 1000 generations, (2) examination of run output for convergence (standard deviation of split frequencies nearing 0.001) (3) removal of the first 1000 samples (10%) as burnin after visual inspection of likelihood score plots, (4) comparison of consensus trees for each run, and (5) combination of post-burnin samples from all three runs to compute a 50% majority-rule consensus tree (conducted in PAUP* v 4.0 (Swofford, 2002)). A partitioned model of sequence evolution was used for the analysis of the combined data.

Maximum likelihood analyses were carried out in GARLI v 2.0 (Zwickl, 2006). For each alignment, two search replicates were performed in a single execution. Models of evolution were the same as those described for Bayesian analyses, with a partitioned model applied to the combined alignment. Other parameters were kept at default. Statistical support was inferred with 100 replicates of bootstrap reweighting (Felsenstein, 1985), implemented as in the tree search.

Maximum parsimony analysis was conducted in PAUP* v 4.0 (Swofford, 2002). An initial heuristic search of 100 random taxon addition replicates was conducted with tree-bisection-reconnection branch swapping (TBR) and MULPARS in effect, retaining only ten trees from each replicate. A strict consensus of these trees was then used as a constraint tree in a second heuristic search using the similar parameters as above, but with 1000 random sequence addition replicates, and retaining 100 trees per addition replicate. We used this method due to the excessive number of trees generated by unconstrained searches. This strategy checks for shorter trees than those found by the initial search, demonstrating that the final consensus tree reflects all of the most parsimonious trees (Catalán, Kellogg & Olmstead, 1997). We also ran searches on the three alignments using an unconstrained search with the nearest neighbor interchange (NNI) swapping algorithm, which produced trees of exactly the same length as the constrained searches. In the interest of brevity, we present results only for the constrained searches. We estimated Bootstrap support (Felsenstein, 1985) for our parsimony trees using 100 pseudoreplicates and the same search setting as described above, including use of a constraint tree. We treated gaps as missing data for all phylogenetic analyses.

Topology testing

We used Templeton’s nonparametric test (1983), as implemented in PAUP* v 4.0 (Swofford, 2002), to evaluate several key phylogenetic relationships. Templeton’s test compares pairs of topologies, measuring relative statistical support for the trees within a sequence dataset (alignment). For these tests, we compared the best tree from the original parsimony tree search (see above) to the best tree from a search using a constraint (e.g., African Pachypodium constrained as monophyletic). For more on these tests, see below (Results).

Alignments

The ITS region had an aligned length of 658 bp (Alignment S1). Of the 156 (included) variable positions within the ingroup, 110 were parsimony informative (Table 3). The trnL-F region had an aligned length of 961 bp (Alignment S2). Of the 33 variable positions within the ingroup, 18 were parsimony-informative (Table 3). The combined alignment contained 61 terminals, with an aligned length of 1619 bp (Alignment S3). Of the 184 (included) variable positions in the ingroup, 114 were parsimony informative.

Phylogenetic trees

The Bayesian 50% majority-rule consensus tree for ITS contained 13 internal nodes with a posterior probability (PP) of 1.0 (Treefile S4A; Fig. 2). By contrast, the trnL-F-based Bayesian tree contained only five internal nodes with a PP of 1.0 (Treefile S5A; Fig. 2). The combined ITS and trnL-F tree contained 17 internal nodes with a PP of 1.0 (Treefile S6A; Fig. 3).

Maximum parsimony searches based on ITS data alone resulted in 4851 trees of 324 steps (Table 4; Treefile S4B); a total of 12 internal nodes had bootstrap (BS) support greater than or equal to 95% (Treefile S4C). Searches using trnL-F data alone resulted in 8 trees of 71 steps (Table 4; Treefile S5B); only one internal node had BS support greater than or equal to 95% (Treefile S5C). Searches on the combined ITS and trnL-F data resulted in 4582 trees of 394 steps (Table 4; Treefile S6B); a total of 8 internal nodes had BS support greater than or equal to 95% (Treefile S6C; Fig. 3). In all cases, use of a constraint tree failed to find any trees of equal or shorter length that contradicted the respective consensus trees.

Maximum likelihood (ML) analyses support similar relationships to those indicated by maximum parsimony and Bayesian analyses. The best ML tree for ITS alone contained 14 internal nodes with BS support greater than or equal to 95% (Treefiles S4D and S4E). By contrast, the trnL-F-based ML tree contained only one internal node with BS greater than or equal to 95% (Treefiles S5D and S5E). The best ML tree based on ITS combined with trnL-F contained 10 internal nodes with BS support greater than or equal to 95% (Treefiles S6D and S6E; Fig. 3).

Pachypodium is recovered as monophyletic in the trnL-F tree (Fig. 2A), but lack of broad outgroup sampling for ITS prevents assessment of Pachypodium monophyly based on nuclear DNA; support for Pachypodium monophyly in the combined tree is driven by trnL-F. Six of the 11 minimum-rank Pachypodium taxa (species and subspecies) represented by more than one sampled plant (Table 1) are monophyletic in the combined tree, four with strong support (PP 1.0; MP bootstrap ≥ 95%; P. baronii, P. decaryi, P. rosulatum subsp. rosulatum, and P. windsorii; Fig. 3). The following multi-taxon clades are also recovered with high levels of support in the combined tree (PP = 1.0; MP BS ≥ 95%): (1) the Malagasy P. decaryi, P. rutenbergianum, and P. sofiense, (2) the African P. lealii and P. saundersii, (3) the African P. namaquanum, P. succulentum, and P. bispinosum, (4) an 11-taxon group corresponding to section Gymnopus (Table 2), and (5) a smaller group nested within Gymnopus comprising P. brevicaule subsp. brevicaule, P. densiflorum, P. eburneum, P. inopinatum, and P. rosulatum subsp. bicolor.

Topology test results

Based on the results from our initial tree searches (Figs. 2 and 3), we were interested to know whether the data could reject (1) monophyly of African Pachypodium, (2) monophyly of Malagasy Pachypodium, and (3) reciprocal monophyly of African and Malagasy Pachypodium. These tests were done by comparing the most parsimonious tree from the original heuristic tree search to the most parsimonious tree from a search in which one of the above groups was used as a constraint. We carried out these analyses for ITS and for the combined data. Because the trnL-F region was not sampled for one of the African species (P. bispinosum), we were not able to evaluate these hypotheses on the basis of chloroplast DNA alone. For ITS, the shortest tree compatible with the first constraint (monophyletic African Pachypodium) was four steps longer (328 steps) than the unconstrained tree (324 steps), which was judged not to be significant based on a Templeton test (P = 0.25). A similar result was obtained for the combined data (396 steps in the constrained tree versus 394 steps in the unconstrained tree; P = 0.64). For the second constraint (monophyletic Malagasy Pachypodium), the shortest ITS tree compatible with the constraint was only one step longer than the unconstrained tree, which was also not significant based on a Templeton test (P = 0.71); again, the combined data were in agreement (both trees 394 steps; P = 1.0). Finally, for the third constraint (reciprocal monophyly of African and Malagasy Pachypodium), the shortest ITS tree compatible with the constraint was five steps longer than the unconstrained tree, which was not a significant difference (P = 0.1); the combined data support this result 398 steps in the constrained tree versus 394 steps in the unconstrained tree; P = 0.29.

Conflict

Our study identified significant conflict between the nuclear and chloroplast datasets, based on the incongruence length difference test (see Materials and Methods). However, we elected to combine the datasets for further analysis. Our choice to unite the conflicting datasets is a conditional combination approach (Bull et al., 1993; Huelsenbeck, Bull & Cunningham, 1996), based on the lack of conflict between well-supported internal nodes (also called “hard conflict”) in the trnL-F and ITS trees (Fig. 2). Our combined approach should be treated as tentative, despite the lack of clearly conflicting internal nodes in ITS versus trnL-F trees.

Phylogenetic relationships

Our trnL-F trees suggest that Pachypodium is monophyletic, based on sampling of closely related genera. However, because of a lack of appropriate outgroups for the nuclear region (ITS), we were unable to evaluate the hypothesis of Pachypodium monophyly on the basis of both genomes. Nevertheless, the monophyly of Pachypodium is generally uncontroversial, and is supported by other molecular phylogenetic research (Livshultz et al., 2007), as well as a suite of morphological characters, including alternate phyllotaxy (most Apocynaceae have opposite leaf arrangement), a horseshoe-shaped retinacle (the connection between the anther and the style head), loss of colleters associated with the calyx, and stem succulence (Sennblad, Endress & Bremer, 1998).

Overall, our data do not provide sufficient phylogenetic resolution to draw conclusions concerning the monophyly or non-monophyly of African and Malagasy Pachypodium. Despite the recovery of several well-supported lineages in both African and Malagasy Pachypodium, the basal branching relationships among these lineages is not well resolved by ITS, trnL-F, or the combined data (Figs. 2 and 3). However, it should be noted that trnL-F provides some evidence for the cohesiveness of African Pachypodium (Fig. 2B); lack of ITS data for reliably vouchered P. bispinosum makes it impossible to test this hypothesis using trnL-F, although sequence data for samples of P. bispinosum of unknown wild origin (horticultural strains) do group with other African species in trnL-F trees (D. Burge and A. Agrawal, unpublished data). In general, there are four mutually exclusive hypotheses on the relationship between African and Malagasy Pachypodium, each of which may represent a valid interpretation of our results: (1) reciprocally monophyletic African and Malagasy Pachypodium, (2) monophyletic Malagasy Pachypodium derived from within a basal grade of African Pachypodium, rendering African Pachypodium paraphyletic, (3) monophyletic African Pachypodium arising from a basal grade of Malagasy Pachypodium, with Malagasy Pachypodium paraphyletic, and (4) neither African nor Malagasy Pachypodium monophyletic. Topology tests could not reject any of these hypotheses.

A recent estimate of 37–64 Ma for the divergence of stem Apocynaceae from its closest relatives among the Gentianales (Bell, Soltis & Soltis, 2010) implies that the crown age of Pachypodium is probably more recent than the ∼80 Ma timing for the isolation of Madagascar from Africa (Yoder & Nowak, 2006). In fact, a recent review of Madagascar biogeography suggests that most of Madagascar’s biotic connections are best explained by long-distance dispersal during the Cenozoic, rather than ancient Gondwanan vicariance (Yoder & Nowak, 2006). Thus, if Pachypodium did not originate in Madagascar, it must have arrived on the island via long-distance dispersal. However, the lack of phylogenetic resolution among major African and Malagasy lineages of Pachypodium prevents preventing reliable reconstruction of geographic range evolution, including dispersal-vicariance scenarios between Africa and Madagascar.

Additional molecular phylogenetic work will be required to obtain better support for basal-branching relationships in Pachypodium, particularly the relationship between African and Malagasy species. This work will likely require the sequencing of additional loci, from both the chloroplast and nuclear genome. Resolution of relationships among species from Lüthy’s (2004) section Gymnopus will also require additional work. In Gymnopus, a number of widespread species (e.g., P. densiflorum and P. brevicaule) are non-monophyletic. The lack of phylogenetic cohesiveness among populations in such species is consistent with both hybridization following initial divergence, as well as incomplete lineage sorting (retention of ancestral polymorphisms; Pamilo & Nei, 1988; Maddison & Knowles, 2006), a phenomenon that often occurs during rapid diversification). For future studies on section Gymnopus, rapidly evolving genetic markers such as low-copy nuclear genes may help to discern species-trees from gene-trees, while population genetic markers such as AFLPs and microsatellites might also help to decipher complex relationships, especially in regions of geographic overlap among species.

Testing classification

Our exhaustive sampling of Pachypodium species and subspecies (Table 1) has provided the opportunity to test existing morphology-based hypotheses on infrageneric relationships. Our results support the most recent infrageneric classification of Pachypodium proposed by Lüthy (2004; Table 2). Lüthy’s (2004) shrubby, predominantly yellow-flowered section Gymnopus is clearly monophyletic (Fig. 3, PP 1.0; MP & ML BS 100%), as is the shrubby, red-flowered section Porphyropodium (Fig. 3, PP 0.98; MP and ML 98%). Our results also indicate a very close relationship between Porphyropodium and Gymnopus (Fig. 3, PP 0.96; MP and ML BS >86%), a relationship not emphasized by past classifications. Finally, the third section recognized in Lüthy’s (2004) classification, the mostly arborescent, white-flowered Leucopodium, is marginally supported in the combined phylogenetic tree (Fig. 3, PP 0.94; ML BS 71%). Overall, our results also support the tradition of using corolla color as a basis for circumscription of taxa within Pachypodium (Fig. 3; Poisson, 1924; Pichon, 1949; Perrier de la Bâthie, 1934; Lüthy, 2004). Nonetheless, we agree with Lüthy (2004) that an ideal infrageneric classification should use multiple morphological characteristics to define groups.

Below the section level, previous classifications of Pachypodium are not well supported by our molecular phylogenetic results. One clear exception is Lüthy’s (2004) series Contorta (Table 2), which was defined on the basis of seed morphology to include the arborescent P. rutenbergianum and P. sofiense, as well as the limestone-endemic P. decaryi. Our results show that this group is strongly monophyletic (Fig. 3, PP 1.0; MP and ML BS 100%), confirming the detailed work of Lüthy (2004). However, this contrasts with most previous opinions. Pichon (1949), for example, allied P. decaryi with another limestone endemic, P. ambongense.

Within section Gymnopus, Lüthy’s (2004) series Densiflora (Table 2) roughly corresponds to a clade that we recover nested inside Gymnopus (Fig. 3, Clade A, PP 1.0; MP BS 76%). However, Clade A includes P. rosulatum subsp. bicolor and P. brevicaule subsp. brevicaule, both considered members of series Ramosa by Lüthy (2004). Our results indicate that the floral characters used by Lüthy (2004) and others to define groups within Gymnopus (Table 2) are homoplasious.

Most past classifications of Pachypodium have dealt in very sparse detail, if at all, with the distinctive and morphologically heterogeneous African members of the genus. As discussed above (see Phylogenetic relationships), our results suggest that African Pachypodium comprises two distinctive lineages, one containing the morphologically similar P. lealii and P. saundersii (Rapanarivo et al., 1999), and a second containing the bizarre monopodial tree P. namaquanum and the tuberous shrubs P. bispinosum and P. succulentum. The close relationship between P. lealii and P. saundersii (Fig. 3, PP 1.0; MP BS 95%) has been noted for some time, as indicated by a reduction to synonymy under P. saundersii that was undertaken by Rowley (1973). The close relationship of P. namaquanum to P. bispinosum and P. succulentum was less expected (Fig. 3, PP 1.0; MP and ML BS 100%). Vorster & Vorster (1973) did propose a close relationship between P. namaquanum and P. bispinosum based on corolla shape. However, these authors also proposed that the asymmetrical flowers of P. succulentum linked this species to P. lealii and P. saundersii more than to P. bispinosum. Our results clearly show that P. bispinosum and P. succulentum are one another’s closest relatives, sister to P. namaquanum.

Conservation

Conservation planning for threatened flora and fauna must take into consideration the evolutionary potential of populations and taxa (Forest et al., 2007). Ignoring evolutionary potential will lead to losses of diversity that compromise the ability of these groups to adapt and survive in the long-term. In the case of Pachypodium, phylogenetic results presented here show that several species and groups of species are strongly divergent from other Pachypodium (e.g., P. decaryi and most African Pachypodium; Fig. 3). These groups represent important islands of phylogenetic diversity within Pachypodium, the loss of which would drastically reduce the overall diversity of the genus. Many members of the Gymnopus section of Pachypodium, by contrast, are very shallowly divergent based on our results (Fig. 3). The members of Gymnopus are adapted to a great variety of habitats, and therefore may contain much ecological diversity in terms of local adaptation (Lüthy, 2004). However, in comparison to highly divergent taxa such as P. decaryi, each individual Gymnopus taxon represents a very small proportion of the total phylogenetic diversity of Pachypodium. In light of the always-limited resources available for conservation, an effort should be made to prioritize the protection of phylogenetically divergent lineages of Pachypodium as well as the overall genetic diversity of the genus. We recommend stronger conservation measures—including greater restrictions on the trade of wild-collected plants—for very narrowly distributed species having Bayesian PP of 1.0 in the combined ITS and trnL-F tree (Fig. 3). This includes the Malagasy P. baronii, P. windsorii, and P. decaryi. The highly divergent African species are not included in this list due to their relatively wide geographic distributions.