Evolution of Philodendron (Araceae) species in Neotropical biomes

Taxon and gene sampling

We have sequenced new data for 110 extant species of Philodendron and 16 species of Homalomena of the following molecular markers: the nuclear 18S and external transcribed spacer (ETS); and the chloroplast trnL intron, trnL-trnF intergenic spacer, the trnK intron and maturase K (matK) genes. Additionally, 13 outgroup species were analyzed, comprising the genera Cercestis, Culcasia, Colocasia, Dieffenbachia, Heteropsis, Montrichardia, Nephthytis, Furtadoa and Urospatha. Outgroup choice was based on the close evolutionary affinity of these genera to Philodendron, as suggested by previous studies. The complete list of species included in this study, the voucher and GenBank accession numbers were listed in Tables 1 and 2 of the Supplemental Information 1.

Ancestral biome reconstruction is dependent on the estimated phylogeny and the current geographic distribution of sampled species terminals. Thus, taxon sampling may impact the inference of ancestral species distribution along biomes. As indicated in Table S1 we have sampled all P. subg. Meconostigma species; 82 P. subg. Philodendron species and 7 P. subg. Pteromischum species. Our sampling strategy is representative of the current Philodendron diversity. Although ∼75% of the sampled species are P. subg. Philodendron in our analysis, ∼82% of Philodendron species consist of P. subg. Philodendron (Boyce & Croat, 2013; Coelho et al., 2016).

DNA isolation, amplification and sequencing

Genomic DNA was isolated with QIAGEN DNeasy Blood & Tissue kit from silica-dried or fresh leaves. Primers used for amplification and sequencing were listed in Table S3. Sequencing reactions were performed in the Applied Biosystems 3730xl automatic sequencer and edited with the Geneious 5.5.3 software.

Alignment and phylogenetic analysis

Molecular markers were individually aligned in MAFFT 7 (Katoh & Standley, 2013) and then manually adjusted in SeaView 4 (Gouy, Guindon & Gascuel, 2010). We estimated individual gene trees (Fig. 1, SM) for each molecular marker in MrBayes 3.2.2 (Huelsenbeck & Ronquist, 2001; Ronquist & Huelsenbeck, 2003) using the GTR + G substitution model. The Markov chain Monte Carlo (MCMC) algorithm was ran twice for 10,000,000 generations, using four chains. Chains were sampled every 100th cycle and a burn-in of 20% was applied. A supertree was estimated from the tree topologies of each molecular marker using the PhySIC_IST algorithm, available at the ATGC-Montpellier online server (http://www.atgc-montpellier.fr/physic_ist/). Only clades with posterior probability ≥ 85% were considered to estimate the supertree. We have used this approach to avoid the impact of missing data in phylogeny estimation (Scornavacca et al., 2008). As PhySIC_IST calculates non-plenary supertrees, it removes taxa with significant topological conflict and/or with small taxon sampling (Scornavacca et al., 2008). The final supertree was thus composed of 89 terminals, as 50 terminals were discarded due to conflicting resolutions.

In order to assess the stability of the (Philodendron + American Homalomena) clade, we have calculated the log-likelihoods of alternative topological arrangements in PhyML 3.0 (Guindon et al., 2009) using the species sampling of the supertree. We have tested the following topologies: (I) (American Homalomena (P. subg. Philodendron +P. subg. Meconostigma); (II) (P. subg. Meconostigma (P. subg. Philodendron + American Homalomena) and (III) (P. subg. Philodendron (P. subg. Meconostigma + American Homalomena). The significance of the difference in log-likelihoods between topologies was tested with the approximately unbiased (AU) and the Shimodaira-Hasegawa (SH) tests implemented in CONSEL 1.2.0 (Shimodaira & Hasegawa, 2001).

Divergence time inference

Dating Philodendron evolutionary history is difficult mainly because of the scarcity of the fossil record (Loss-Oliveira et al., 2014). For instance, Dilcher & Daghlian (1977), based on fossilized leaves, described a putative P. subg. Meconostigma fossil from the Eocene of Tennessee (56.0–33.9 Ma). However, Mayo (1991) identified the referred fossil as a Peltranda. Thus, we have decided not to use this fossil as calibration information. Alternatively, in order to estimate divergence times, we have assigned a prior on the rate of nucleotide substitution. We were then prompted to infer the evolutionary rates of plastid coding regions of monocots using a large sample of publicly available chloroplast genomes. Nuclear genes were excluded from dating analysis because of the absence of prior information on evolutionary rates.

To estimate monocots substitution rate, we used chloroplast genomes from 154 Liliopsida species retrieved from the GenBank (Table S4). All orthologous coding regions were concatenated into a single supermatrix. Maximum likelihood phylogentic reconstruction was implemented in RaxML 7.0.3 (Stamatakis, 2006) under GTR model. Molecular dating of monocots (Liliopsida) was conducted under a Bayesian framework, using fossil information obtained from Iles et al. (2015) (Table S5). Because the number of terminals used was large, rate estimation was conducted with the MCMCTree program of PAML 4.8 package (Yang, 2007) using the approximate likelihood calculation (Dos Reis & Yang, 2011) and the uncorrelated model of evolution of rates. In MCMCTree, posterior distributions were obtained via MCMC; chains were sampled every 500th cycle until 50,000 trees were collected. We performed two independent replicates to check for convergence of the estimates. Calibration information for Liliopsida was entered as minimum and maximum bounds of uniform priors. The estimated mean substitution rate was inferred at 3.26 × 10⁻⁹ substitutions/site/year (s/s/y). This value is significantly higher than the previous estimate of Palmer (1991), which reported an average substitution rate of 0.7 × 10⁻⁹ s/s/y for angiosperm platids. As the credibility interval of our estimate was large, we adopted a Gaussian prior for evolutionary rates with a 95% highest probability density (HPD) interval that included maximum and minimum values of our estimate and that of Palmer’s.

Dating analysis of Philodendron and Homalomena species was performed in BEAST using a relaxed molecular clock with evolutionary rates modeled by an uncorrelated lognormal distribution; the GTR + G□model of sequence was applied. The MCMC algorithm was ran for 50,000,000 generations and sampled every 1,000th cycle, with a burn-in of 20%.

Biome shifts

To unveil how Philodendron species colonized the Amazon forest, Atlantic forest, Cerrado and Caatinga, we conducted a Bayesian Binary MCMC (BBM) (Yu, Harris & He, 2012; Ronquist & Huelsenbeck, 2003) implemented in Reconstruct Ancestral State in Phylogenies 2.1b (RASP) software (Yu, Harris & He, 2012). The input tree topology was the supertree estimated in PhySIC_IST. BBM chains were ran for 10,000,000 generations and were sampled every 1,000th cycle. State frequencies were estimated under the F81 model with gamma rate variation. Information on the occurrence of each Philodendron species along Neotropical biomes was obtained from Coelho et al. (2016) and from the (Team) CATE Araceae (http://araceae.e-monocot.org).

Phylogenetic relationship between Philodendron and Homalomena

In this study, Asian Homalomena was recovered as sister to the (Philodendron + American Homalomena) clade, and Furtadoa mixta clustered with the Asian Homalomena clade. The evolutionary affinities of American Homalomena, P. subg. Meconostigma and P. subg. Philodendron were not strongly supported. However, the topological arrangement in which Philodendron is a monophyletic genus was statistically significant by the AU and SH tests, suggesting the monophyly of Philodendron.

Previous studies have reported conflicting results concerning the monophyly of Philodendron and the phylogenetic status of American Homalomena (Fig. 5). For instance, Barabé et al. (2002), based on the trnL intron and the trnL-trnF intergenic spacer, proposed P. subg. Philodendron as a paraphyletic group and was unable to solve the (P. subg. Meconostigma + Asian + American Homalomena) polytomy (Fig. 5A). Gauthier, Barabé & Bruneau (2008) recovered the American Homalomena as sister to Philodendron and the Asian Homalomena as sister to the (American Homalomena + Philodendron) clade, although their Bayesian analysis inferred a paraphyletic Philodendron, with P. subg. Pteromischum sister to the American Homalomena (Figs. 5B and 5C, respectively). Alternatively, Cusimano et al. (2011) recovered a monophyletic Philodendron, with Homalomena as sister lineage of Furtadoa (Fig. 5D). Recently, Yeng et al. (2013) estimated the Homalomena phylogeny based on the nuclear ITS marker and also sampled Philodendron species. In the ML and Bayesian trees reported in their study, P. subg. Pteromischum was closely related to the American Homalomena, whereas P. subg. Meconostigma and P. subg. Philodendron were recovered as sister taxa (Fig. 5E).

Discrepancies between previous works and our analysis may be due to different choice of phylogenetic methods, markers and taxon sampling. Gauthier, Barabé & Bruneau (2008) was the only study intended to investigate specifically the systematics of Philodendron genus. When compared to their analysis, our study included a larger sampling of taxa and molecular markers with the aim of estimating the phylogeny of Philodendron and Homalomena species; it is also the first analysis that used a supertree approach to this end.

Our phylogeny characteristically presents short branch lengths within the Philodendron clade. The high frequency of polytomies indicates the genetic similarity among terminals, which is further corroborated by the ease in obtaining artificial hybrids between different species. This corroborates a scenario of low genetic differentiation and low reproductive isolation (Carlsen, 2011).

Philodendron diversification may also consist of several recent rapid radiation events. Phylogenetic reconstruction under this scenario is challenging, because of a significant amount of substitutions is needed to accumulate within short periods of time (Maddison & Knowles, 2006). However, morphological variation of Philodendron is remarkable, which seems contradictory considering the previously discussed features. However, it has been extensively discussed that morphological variation is not a suitable proxy for genetic variation (e.g., Prud’Homme et al., 2011; Houle, Govindaraju & Omholt, 2010). Many environmental and epigenetic factors may can increase phenotypic variation even in the absence of DNA sequence variation (Prud’Homme et al., 2011). Evidently, we cannot rule out the possibility that DNA regions that present significant genetic differences between species were not sampled in this work.

Diversification of Philodendron and Homalomena

Although the chronology of Philodendron divergence was not extensively focused by previous studies, Nauheimer, Metzler & Renner (2012) analyzed the global history of the entire Araceae family based on a supermatrix composed of five chloroplast markers and several well-established calibration points. Their analysis included a single Philodendron species and estimated age of the Philodendron/Asian Homalomena divergence at approximately 20.0 Ma (ranging from 31.0–9.0 Ma). This study, however, also included a single species of Asian Homalomena.

The wide range of the posterior distribution credibility intervals of Nauheimer, Metzler & Renner (2012) hampers the proposition of putative biogeographic scenarios for the evolution of Philodendron, American and Asian Homalomena. Differences between their timescale and the divergence times proposed in this study might therefore be due to methodological differences caused by their reduced taxonomic sampling. Nevertheless, both our estimate of the age of the Philodendron divergence from Asian Homalomena and that of Nauheimer, Metzler & Renner (2012) suggests that this event took place when South America was essentially an isolated continent.

The isolation of the South American continent persisted from approximately 130.0 Ma (Smith & Klicka, 2010) to 3.5 Ma (Vilela et al., 2014), with the rise of the Panamanian land bridge. Therefore, from the Early to Middle Miocene there was no land connection with North America, Asia or Africa (Oliveira, Molina & Marroig, 2010). If dispersal, rather then vicariance, is the most plausible hypothesis to explain Philodendron and American Homalomena colonization of the Neotropics, hypotheses on the possible routes of colonization should be investigated. Based on the continental arrangement during the Miocene, we propose that the dispersal of Philodendron and American Homalomena ancestor could have followed four possible routes (Fig. 6): (1) from Asia to North America through the Bering Strait; (2) from Africa to the Neotropics by crossing the Atlantic ocean; (3) from Asia to Neotropics by crossing Pacific ocean; and (4) from Asia to Neotropics , also by crossing the Atlantic ocean.

The Araceae fossil record is currently assigned to Florida, Russia, Germany, United Kingdom, Venezuela, Yemen, Colombia and Canada (Shufeldt, 1917; Berry, 1936; Bogner, Hoffman & Aulenback, 2005; Chandler, 1964; Dorofeev, 1963; As-Saruri, Whybrow & Collinson, 1999; Wilde & Frankenhauser, 1998; Wing et al., 2009; Stockey, Rothwell & Johnson, 2007). However, as none of the fossil specimens was described as closely related to Philodendron or Homalomena, the Araceae fossil record fails to corroborate any dispersal hypothesis in particular.

Considering route 1, although the Bering Strait have connected Asia to the North America during most of the Cenozoic period (Butzin et al., 2011), there is no evidence of extant Philodendron and Homalomena in North America or North Asia. Route 2 involves long-distance oceanic dispersal through ca. 2,000 km—the minimum distance between Africa and the Neotropics (Oliveira, Molina & Marroig, 2010)—through Atlantic paleocurrents, which were probably stronger than Pacific currents. This hypothesis is congruent with the clustering of Philodendron and American Homalomena into a single clade, assuming Africa as the center of diversification of Asian and American Homalomena, as well as Philodendron. However, we should conisder that the last recent common ancestor of Philodendron and Homalomena was distributed in Africa. On the other hand, this hypothesis is corroborated by the distribution of the extant Philodendron and Homalomena species. Givnish and colleagues (2004) also suggested two long-distance dispersal events through the Atlantic, but in the opposite direction. Their analysis indicated that Bromeliaceae and Rapateaceae arose in the Guayana Shield of northern South America and reached tropical west Africa via long-distance dispersal at ca. 6–8 Ma.

When considering long-distance dispersal events, it is crucial to evaluate their viability as related with the plant’s ability to produce dispersal structures that would tolerate aquatic and saline conditions for long periods of time (Lo, Norman & Sun, 2014). Although such features have not been evualuated for Philodendron and Homalomena, some Homalomena species inhabits swamp forests and open swamps. Thus, features that would favor their survival in waterlogged environments could also influence their maintenance in seawater.

Although route 3 is geographically unlikely due to the 8,000 km distance between Asia and the Neotropics through the Pacific Ocean (Oliveira, Molina & Marroig, 2010), it cannot be completely discarded because it is corroborated by the extant distribution of Homalomena and Philodendron. Finally, route 4 suggests the dispersal through the Atlantic ocean from Asia to the Neotropics. This is also an improbable hypothesis because the African continent would act as a barrier between Asia and the Neotropics, requiring the dispersal through both the Indian and the Atlantic oceans.

The extant distribution of Philodendron and Homalomena species and the scarcity of fossil information challenge the proposition of a scenario for the origin of Philodendron and American Homalomena in the Neotropics. However, the biological and geographical information provided to date indicates a long-distance oceanic dispersal through the Atlantic, as suggested by route 2, as the most plausible hypothesis to explain Philodendron and American Homalomena colonization of the Neotropics.

Early diversification of Philodendron species

According to our analysis, the last common ancestor of Philodendron and the American Homalomena was distributed in the Amazon forest about 8.6 Ma (11.1–6.8 Ma) during the Middle/Late Miocene. Interestingly, this time estimate is very close to the age of the divergence between the (Philodendron/American Homalomena) clade from the Asian Homalomena (Fig. 4, node A). The Middle and Late Miocene were characterized by wetland expansion into western Central Amazonia, which fragmented the rainforest and formed extensive wetlands (Jaramillo et al., 2010). According to our analysis, Philodendron earliest divergence events took place in this scenario. The Amazon forest, from the Late Miocene to the beginning of Pliocene, was composed of a diverse and well-structured forest. The Amazon river landscape was well established; this probably allowed the extensive development of the Amazonian terra firme forest (Jaramillo et al., 2010). This scenario is compatible with the biology of extant species of Philodendron because a well-structured forest would allow the development of epiphyte and hemiepiphyte species, such as Philodendron.

Philodendron diversification along Neotropical biomes

Our results suggest that Philodendron species occurred exclusively at the Amazon forest for ca. 5.0–6.0 Ma. During the Pliocene, as result of the glacial cycles, climate cooling and drying permitted the expansion of the open savanna areas, mostly represented by the ‘dry diagonal’, which is constituted by the Caatinga, Cerrado and Chaco biomes. This consisted of a crucial event, because it resulted in the isolation of the Atlantic forest in the east coast of South America (DaSilva & Pinto-da-Rocha, 2013), which is synchronous to the inferred age of the early diversification of Philodendron in this biome. This also corroborates the hypothesis that the Atlantic forest taxa present a closer biogeographic relationship with the Amazon forest, as proposed by Amorim & Pires (1996) and Eberhard & Bermingham (2005). After the separation between Atlantic and Amazon forests during the Pliocene, species dispersal might have been common through the forest patches (DaSilva & Pinto-da-Rocha, 2013).

Roig-Juñent & Coscarón (2001) and Porzecanski & Cracraft (2005) suggested that the Atlantic rainforest also presents similarities in organismal composition with the Cerrado biome. This association would have been a result of dispersal events through gallery forests. The history of the formation of Cerrado biome is still uncertain (Zanella, 2013; Werneck, 2011), but our analysis indicated that the ancestors of Philodendron clades from the Cerrado were distributed in the Atlantic forest. Therefore, we also corroborate the hypothesis of lineage dispersal from the Atlantic forest to the Cerrado biome. These events would have occurred after the colonization the Atlantic forest by Philodendron species.

Final considerations on Philodendron evolution

Given the significant morphological diversity of Philodendron, its widespread distribution in the Neotropics and the age of the Araceae family (∼140.0 Ma, Nauheimer, Metzler & Renner, 2012), it would be expected that the origin of this genus was older. In sharp contrast, we have estimated phylogenies with very short branch lengths and very recent divergence times. A similar scenario was reported by Carlsen & Croat (2013) for Anthurium, which is the most diverse Araceae genus, and also by Nagalingum and colleagues (2011) for cycads. Therefore, the inferred tempo and mode of evolution of Philodendron species were reported in several plant families.