Explosive radiation in high Andean Hypericum—rates of diversification among New World lineages

The páramos, high-elevation Andean grasslands ranging from ca. 2800 m to the snow line, harbor one of the fastest evolving biomes worldwide since their appearance in the northern Andes 3–5 million years (Ma) ago. Hypericum (St. John's wort), with over 65% of its Neotropical species, has a center of diversity in these high Mountain ecosystems. Using nuclear rDNA internal transcribed spacer (ITS) sequences of a broad sample of New World Hypericum species we investigate phylogenetic patterns, estimate divergence times, and provide the first insights into diversification rates within the genus in the Neotropics. Two lineages appear to have independently dispersed into South America around 3.5 Ma ago, one of which has radiated in the páramos (Brathys). We find strong support for the polyphyly of section Trigynobrathys, several species of which group within Brathys, while others are found in temperate lowland South America (Trigynobrathys s.str.). All páramo species of Hypericum group in one clade. Within these páramo Hypericum species enormous phenotypic evolution has taken place (life forms from arborescent to prostrate shrubs) evidently in a short time frame. We hypothesize multiple mechanisms to be responsible for the low differentiation in the ITS region contrary to the high morphological diversity found in Hypericum in the páramos. Amongst these may be ongoing hybridization and incomplete lineage sorting, as well as the putative adaptive radiation, which can explain the contrast between phenotypic diversity and the close phylogenetic relationships.


INTRODUCTION
High altitude mountain regions in the tropics constitute an open grassland vegetation type that is characterized by large rosette and small cushion plants, bunch grasses, and evergreen sclerophyllous shrubs (Luteyn, 1999). These high-elevation grasslands occur above the tree-line and below the upper limits of plant life "along the crests of the highest mountain ranges or on isolated mountaintops, like islands in a sea of forest" (Luteyn, 1999). The páramos of South America are discontinuously distributed along the Andean Cordilleras between ca. 2800-5000 m a.s.l. in Venezuela, Colombia, Ecuador, and northern Peru, with outliers in the Cordillera de Talamanca in Costa Rica and adjacent Panama.
A hyperdiverse high-elevation ecosystem has evolved since the northern Andes were formed in the Pliocene above the modern tree-line that marks the lower limit of páramo vegetation (ca. 3-5 Ma ago;Gregory-Wodzicki, 2000;Hoorn et al., 2010), which makes it one of the youngest of the Neotropical ecosystems (Graham, 2009). Today, the páramos alone comprise around 3500 vascular plant species (Luteyn, 1999;Sklenár et al., 2010;Graham, 2011a), half of which are of temperate origin at the generic level (Smith and Cleef, 1988;Sklenár et al., 2010;Weigend et al., 2010;Antonelli and Sanmartín, 2011a) as in other high elevation tropical floras (Gehrke and Linder, 2009). In contrast, the páramo flora is the richest overall tropical mountain flora and has the largest number of genera and endemic elements (Smith and Cleef, 1988;Myers et al., 2000).
The flowering plant genus Hypericum (St. John's wort, Hypericaceae) is a prominent and often abundant component of the flora with arborescent shrubs in the sub-páramo (H. laricifolium, H. irazuense), dwarf shrubs in pastures or meadows from elfin forest to higher zones (H. andinum, H. mexicanum, H. juniperinum), or on rocky and disturbed places (H. cardonae, H. humboldtianum), and prostrate plants in damp areas in the grass páramo (H. selaginella, H. prostratum). Around 65 of the ∼100 species described in Hypericum in South America are páramo endemics. Overall, Hypericum is of temperate origin and has its main center of species richness in the Old World (Nürk and Blattner, 2010;Meseguer et al., 2013;Nürk et al., 2013). Molecular phylogenetic studies using sequences of the nuclear rDNA internal transcribed spacers (ITS) (Meseguer et al., 2013;Nürk et al., 2013) and chloroplast data (Meseguer et al., 2013) revealed the three large Hypericum sections Myriandra, Brathys and Trigynobrathys sensu Robson (1977Robson ( , 1987Robson ( , 1990Robson ( , 2012 as a monophyletic group, which includes ca. 90% of the Hypericum species native to the New World 1 , and verified the genus Triadenum to be included within Hypericum as sister to the Myriandra+Brathys s.l. clade of Nürk et al. (2013); i.e., clade B in Meseguer et al. [2013;see Ruhfel et al. (2011) for taxonomic implications].
Hypericum sect. Myriandra species are distributed mainly in the Nearctic with some located in Honduras, Bermuda and the Caribbean. The majority of New World Hypericum is classified in sects. Brathys and Trigynobrathys (Table 1). Although predominately a high-elevation group in the Neotropics, a small number of herbaceous Hypericum species has adapted to lower elevations, distributed below 3000 to less than 1000 m in lowland regions of temperate South America (Robson, 1977 onwards).
While shrubs constitute the dominant life form in North America, a multitude of life forms evolved in the Neotropicsup to 6 m high sclerophylous arborescent to small dwarf shrubs, prostrate shrubs, and perennial to annual herbs. Most of this phenotypic diversity can be found in the páramos of the Andes (Figure 1). This diversity is mirrored by the great variety of ecological conditions occupied (for an overview see Crockett et al., 2010) and cytology; chromosome numbers range from (2n=) 8, 12, 16, 18, 22, 24, and 32 Robson), but see Ruhfel et al. (2011). b Only one species is endemic to Central America (H. limosum from Cuba). (Robson, 1987(Robson, , 1990Moraes et al., 2009) (Robson, 1987(Robson, , 1990Moraes et al., 2009). This prevents comparative phylogenetic studies on the extent and distribution of polyploidy across the Neotropical species. The recency of the páramo biome, the extent of ecological and phenotypic diversity (Figure 1), as well as high species richness (summarized in Robson, 2012) suggest that Hypericum could represent an adaptive radiation (Simpson, 1953) in the páramos of South America. Following Glor (2010), adaptive radiation is a response to natural selection and ecological opportunity involving diversification of species and associated adaptations. To define and diagnose adaptive radiation three operational criteria are considered in the phylogenetic context: (a) multiplication of species and common descent, (b) extraordinary diversification, and (c) adaptation via natural selection (phenotype-environment correlation and trait utility; Schluter, 2000;Sudhaus, 2004;Glor, 2010).
Insights into the evolutionary history of Hypericum in the Neotropics first demand phylogenetic hypotheses in an explicit framework to test questions about diversification rates and niche shifts, biogeography, key innovations, or polyploidization events (Harvey and Pagel, 1991;Emshwiller and Doyle, 1998;Pagel, 1999;Emshwiller and Doyle, 2002;Stephens and Wiens, 2003;Wiens and Donoghue, 2004;Donoghue, 2007, 2009;Crisp et al., 2011;Vamosi and Vamosi, 2011). The consideration of divergence times of lineages in the phylogenetic context is essential, as dispersal patterns and species richness are related to time (Ricklefs and Latham, 1992;Stephens and Wiens, 2003;Wiens and Donoghue, 2004), and area availability (Vamosi and Vamosi, 2010).
In a recent study, Meseguer et al. (2013) estimated divergence times and conducted biogeographic analyses based on chloroplast sequence variation (trnS-trnG, trnL-trnF, psbA-trnH) for the genus Hypericum. On the base of a rather reduced sampling of South American Hypericum (5 species, 7 accessions) they revealed a single Neotropical clade that had a mean crown group age of 3.9 Ma (Meseguer et al., 2013), sister to a clade containing North American and Asian species belonging to sect. Trigynobrathys (together called the "Brathys-group"). Furthermore, they suggest merging sect. Trigynobrathys into sect. Brathys sensu Robson (1977Robson ( , 1987Robson ( , 1990Robson ( , 2012, as one species from sect. Trigynobrathys groups within a clade containing mainly species belonging to sect. Brathys, a result also revealed in phylogenetic studies analyzing morphological (Nürk and Blattner, 2010) and rDNA ITS data (Nürk et al., 2013).
In this study, we conduct new phylogenetic analyses and estimate divergence times in a Bayesian framework, using nuclear ITS sequences of a representative sampling of New World Hypericum. We employed a wide sampling across the New World species of Hypericum and a dense sampling within the Neotropics to expand prior phylogenetic hypotheses (Meseguer et al., 2013;Nürk et al., 2013) on Hypericum in South America. Despite potentially problematic issues resulting from the multicopy nature of ITS (Baldwin, 1992;Baldwin et al., 1995) for phylogenetic inference (e.g., possible paralogs, and/or pseudogenes; Álvarez and Wendel, 2003;Bailey et al., 2003;Nieto Feliner and Rosselló, 2007) we used this marker system as it offers expanded species/accession sampling (Blattner, 1999;Nürk et al., 2013) while being aware of the potential challenges (see Discussion).
Based on a dense sampling within Neotropical Hypericum, we investigate phylogenetic relationships of the South American and especially the high-elevation Andean species. In particular, we aim to address the following questions: (1) Are Neotropical Hypericum species monophyletic, i.e., is there only one clade in the Neotropics, as suggested by Meseguer et al. (2013), or is there more than one? (2) Is the polyphyletic nature of sect.
Trigynobrathys as reported in previous studies (Meseguer et al., 2013;Nürk et al., 2013) an artifact of species sampling? (3) Are ITS based divergence time estimations congruent with those previously reported from cpDNA (Meseguer et al., 2013)? (4) How many lineages colonized the páramos and what are the ages of these lineages? And finally, (5) are there differences in diversification rates between clades of New World Hypericum?
The first three questions examine previously published hypotheses and compare ITS age estimates to those inferred by analysis of chloroplast sequence variation. The latter two questions are a first attempt to diagnose adaptive radiation in Andean high-elevation Hypericum, by aiming at two of the three operational criteria that define adaptive radiation as proposed by Glor (2010), (a) multiplication of species and common descent, and (b) extraordinary diversification.

Our approach involved broad sampling within New World
Hypericum to ascertain hypotheses about backbone relationships (major lineages), as well as of closer related species from the Andean highlands. Samples were obtained from herbarium collections (ANDES, BM, GAT, HEID) and freshly collected silica-gel dried material. Forty-five ITS sequences selected from GenBank were included into the final data set, additionally to 135 sequences (representing 56 species) newly generated for this study (see Appendix, voucher), which have been submitted to the European Nucleotide Archive (http://www.ebi.ac.uk/ena/data/ view/HG004646-HG004780, Accession No. HG004646-780). In total, the final data set contained 180 accessions (93 species), including 135 native to South America (56 species); a three-to ten-fold increase in species sampling in the Neotropics compared to previous studies (Meseguer et al., 2013;Nürk et al., 2013).
Genomic DNA was extracted with the Invisorb® Spin Plant Mini Kit (Stratec Molecular GmbH, Berlin, Germany) following the manufacturer's protocol. Amplification of the ITS region (including ITS-1, 5.8S rDNA, and ITS-2) followed the procedure detailed in Nürk et al. (2013), with separate amplification (and sequencing) of ITS-1 and ITS-2 for poorly preserved herbarium exsiccatae (Blattner, 1999). Cleaned amplification products were sent for sequencing to Eurofins MWG Operon (Ebersberg, Germany). Forward and reverse sequences from each template were manually edited and combined into a single consensus sequences with Geneious v5.4 (Biomatters, available from www. geneious.com). Sequences were checked for patterns in the chromatograms, which suggest multiple non-identical ITS copies (paralogs or pseudogenes), and multiple reads per site were coded as ambiguities (present in the spacer regions only, in 23 newly generated sequences with one to four ambiguous sites per sequence and in eight sequences downloaded from GenBank with one to seven ambiguous sites per sequence).

PHYLOGENETIC INFERENCE
Sequences were aligned using the L-INS-I algorithm implemented in the software Multiple Alignment using Fast Fourier Transform (MAFFT) v6.9 (Katoh et al., 2002;Katoh, 2005) and manually adjusted using PhyDE v0.9 (Available online: http:// www.phyde.de). MrModeltest v2.3 (Nylander, 2004) was used to select the appropriate model of sequence evolution and the SYM + model (Yang, 1993(Yang, , 1994Zharkikh, 1994) was chosen according to the Akaike Information Criterion (Akaike, 1974;Posada and Buckley, 2004). Phylogenetic analyses were performed under Maximum likelihood (ML; Felsenstein, 1981) and Bayesian Inference (BI; Mau et al., 1999). Hypericum faurei R.Keller [=Triadenum japonicum (Blume) Makino] belonging to the Triadenum clade was used as outgroup, following Nürk et al. (2013) who showed this clade to be sister to the rest. For ML analysis the RAxML GUI v1.1 (Stamatakis, 2006;Silvestro and Michalak, 2011) was used with the GTRCAT model and clade support was evaluated with 10000 rapid bootstrap pseudoreplicates (Stamatakis et al., 2008). For BI optimization MrBayes v3.2.1 (Ronquist and Huelsenbeck, 2003) was started with 4 independent runs, each with 4 chains for 10 million generations with the appropriate substitution model (SYM + ), setting temperature to 0.01, sampling every 1000 generations, and using the ML tree as a starting tree, but introducing random perturbations into the starting tree to initiate parameter calculation from different priors to enable detection of possible convergence problems (using the command "mcmcp nperts = 5"). Following the results of Meseguer et al. (2013), which showed a more realistic branch length estimation when introducing a lambda parameter correction in Bayesian phylogenetics on ITS sequence data in Hypericum, we used a 'corrected' exponential prior on branch length of 1/λ = 0.1 ["prset brlenspr = Unconstrained:Exp(100)"]. Convergence of the parameters was monitored using Tracer v1.5 . After discarding 25% of the sampled trees as burnin, posterior probabilities were calculated on the BI stationary sample.

DIVERGENCE TIME ESTIMATIONS
We choose the Bayesian tree to test for rate constancy among lineages. The likelihood scores associated with branch length were calculated on this tree in PAUP * (Swofford, 2002) under the optimal model of sequence evolution and associated parameters with and without a strict molecular clock enforced. We followed the approach of Huelsenbeck and Rannala (1997) to assess significance. A global molecular clock was rejected (p < 0.05) for the rDNA sequence data. Therefore, divergence times were estimated under a relaxed molecular clock employing the uncorrelated lognormal (UCLN) model (Drummond et al., 2006) that assumes branch specific substitution rates to be drawn from a single lognormal distribution estimated from the data. Implementation of the UCLN model in BEAST v1.7.2  together with the use of Markov chain Monte Carlo (MCMC) sampling methods estimates both topology and substitution rates and calculates absolute divergence times and confidence intervals when calibrated with external data (fossils, or secondary calibration points like estimated ages revealed in other studies). The following calibration points were considered, relying on the age estimates reported in Meseguer et al. (2013). (1) The age of the root node estimated to (23.82−) 29.30 (−35, 23) Ma was constrained with a normal distribution that had a mean of 29, and a standard deviation of 5. (2) The age of the Myriandra+Brathys s.l. crown node estimated to (16.99−) 21.92 (−27, 33) Ma was constrained with a normal distribution that had a mean of 22, and a standard deviation of 4. (3) The age of the Myriandra crown node estimated to (9.41−) 13.59 (−18.48) Ma was constrained with a normal distribution that had a mean of 13.5, and a standard deviation of 3.
Analyses were performed in two independent runs in BEAST to test for convergence in divergence times, each consisted of 100 million generations, and sampling a tree every 10000 generation. Each run started from the tree obtained by ML search, after performing a semi-parametric method based on penalized likelihood (Sanderson, 2002) in R (R Development Core Team, 2013) with the "chronopl" command as implemented in the package APE (Paradis et al., 2004). The GTR model of nucleotide substitution was applied with the model of site heterogeneity. The birth and death model of speciation considering incomplete sampling (Stadler, 2009) was set as tree prior. Accessions grouping in the BI tree within the Myriandra and Myriandra+Brathys s.l clades were constrained monophyletic. Convergence of the parameters was monitored using Tracer v1.5  and the resulting trees of the two runs were combined in LogCombiner  with a burnin of 50%. Means and confidence intervals were calculated on the remaining 10002 trees in TreeAnnotator  to obtain a final consensus tree [maximum clade credibility tree that has 95% of the highest posterior density (HPD)] for visualization in FigTree v1.3.1 (Rambaut, 2006(Rambaut, -2013.

DIVERSIFICATION RATES
As approximation of net diversification rates (r) we used the simple macro-evolutionary constant rate, pure-birth (Yule model) taxonomic likelihood estimate of Magallón and Sanderson (2001) calculated as r = (ln N 1 − ln N 0 )/t, where N 1 = extant species (standing taxonomic diversity), N 0 = initial species diversity, here taken as 1, and t = inferred clade age (time in Ma). We calculated r based on crown ages using the Bayesian mean and 95% HPD age estimates over the entire phylogeny, and cladespecific r for Triadenum (6 species), Myriandra (29 species), Trigynobrathys s.str. (≤52 species), Brathys (≥97 species), and the Páramo clade (≤67 species) and assigned species numbers according to species richness of the sections given in Robson (2012). The polyphyly of sect. Trigynobrathys, however, complicates the assignment of species richness to the Trigynobrathys s.str. and Brathys clades, as less than 30% of the species assigned to sect. Trigynobrathys are sampled in this study. Thus, we used the numbers given above, as we cannot approximate, which species belonging to sect. Trigynobrathys sensu Robson (1977 onwards) group within the Brathys and Trigynobrathys s.str. clade, respectively. For diversification rate estimation of the Páramo clade we assigned species richness in a conservative way, using only the number of species reported to be native to páramo habitats (Robson, 1987(Robson, , 1990(Robson, , 2012, i.e., used an underestimated number of species belonging to this clade. To test the hypothesis of extraordinary diversification in the páramos, we applied the likelihood-based approach given in Magallón and Sanderson (2001). Specifically, we ask which of the major clades in our phylogeny are unexpectedly speciesrich (or poor), given their age and the estimated net diversification rate for the entire phylogeny (i.e., the background rate for Triadenum+Myriandra+Brathys s.l.). We calculated 95% confidence intervals for the background net diversification rate (r) based on the crown group age using the formula k upper (t) = 1 + logβ 0.025(1 + α)/r(1 − α − β + αβ) + α + 2β − 1 for the upper boundary value, and k lower (t) = 1 + logβ 0.975(1 + α)/r(1 − α − β + αβ) + α + 2β − 1 for the lower, where β = (e rt − 1/e rt − ε) and α = εβ, and assuming (1) no extinction (relative extinction rate ε = 0), and (2) a reasonably high relative extinction rate (ε = 0.9). The 95% confidence interval of the expected number of species is the range of values between k upper and k lower in a semi-log plot for crown group ages of (log) species diversity vs. age at time t after the origin of a clade, under a given r and ε. Those clades that fall outside the confidence intervals are then regarded as being exceptionally species-rich or poor (Magallón and Sanderson, 2001).

PHYLOGENY
The aligned data matrix comprised 755 characters, of which 366 were variable. The 50% majority-rule consensus tree of the Bayesian (BI) stationary sample (n = 30004 trees) (Figure 2) is highly congruent with the ML consensus tree (not shown) in the sense that there were no conflicts between strongly supported clades (>97% Bayesian posterior probability [BPP], >75% bootstrap support [BS]). Moreover, correction of the exponential prior (lambda) of branch length in the Bayesian analysis resulted in average branch length congruent with those inferred under ML. The trees and data sets produced in this study are available from TreeBase (http://www.treebase.org) study number 14179.
Within the Brathys clade, relationships are well-supported at the most basal dichotomy, but are almost without support between the remaining taxa. However, all species native to the páramos group within one clade (Figure 2), although without strong support (0.81 BPP, 69 BS). The species belonging to sect. Trigynobrathys sensu Robson (1977

AGE ESTIMATES AND DIVERSIFICATION RATES
The ultrametric time-calibrated maximum clade credibility tree (chronogram) obtained by the Bayesian relaxed clock analyses is shown in Figure 3A. The crown age estimates for the major clades are summarized in

Triadenum
Hypericum FIGURE 2 | Phylogeny of New World Hypericum inferred from rDNA ITS sequence Bayesian analysis. In the tree sketch top left, the phylogenetic position of the species in this study within the genus Hypericum is marked by clade names in boldface. Bayesian posterior probabilities and ML bootstrap support is given above the branches (BPP|BS). Rooting is that of Nürk et al. (2013). Accession names consist of a section number (Robson, 2012), the species, and a collection/GenBank identifier. Clade names not considered in Figure 3B are given in gray.  and expected (95% confidence interval, lines), assuming no extinction (solid lines, ε = 0) and high rates of extinction (dashed lines, ε = 0.9), based on crown ages. Note the species richness in the Páramo clade of Hypericum is higher than expected, under both high rates of extinction and no extinction.

Frontiers in Genetics
of Brathys s.l. (this node was not constrained for calibration) is highly congruent with the ages reported by Meseguer et al. (2013). The same is true for the Trigynobrathys s.str. and the Brathys clade when compared to the ages revealed in Meseguer et al. (2013). According to our divergence time estimation, the subclade within Trigynobrathys s.str. containing the Neotropical "lowland" species diversified 4.1 Ma ago (2.2-6.5 95% HPD). The Páramo clade of Hypericum is 3.8 Ma old (2.3-5.6 95% HPD), with slightly younger estimates for the core Páramo clade (3.3 Ma, 1.9-4.8 95% HPD).
Based on the ages obtained by Bayesian divergence time analysis, estimation of net diversification rates (r) under a Yule model and considering taxonomic richness reveals a higher rate for the Páramo clade of Hypericum [r = (0.75−)1.10(−1.86)] compared to all other clades (summarized in Table 2). The test of excessive species richness (Figure 3B) shows extraordinary species diversity www.frontiersin.org September 2013 | Volume 4 | Article 175 | 7 in the Páramo clade when assuming a high relative extinction rate (ε = 0.9) for the entire phylogeny. The same is true when extinction is assumed to be zero (ε = 0.0). In the latter, however, all clades within the Brathys s.l. clade, that is, all clades containing Neotropical species, are revealed to be more species-rich than expected when compared to the background net diversification rate of the entire phylogeny.

DISCUSSION
Our analysis confirms the monophyly of Myriandra and Brathys s.l and corroborates the relationships among key lineages recovered in previous studies (Meseguer et al., 2013;Nürk et al., 2013). Four major clades, Triadenum, Myriandra, Trigynobrathys s.str., and Brathys comprise around 184 Hypericum species, which constitute ca. 90% of the species diversity of Hypericum in the New World (Robson, 2012). Nevertheless, it partially rejects the traditional infrageneric classification (Robson, 1977 onwards), revealing sect. Trigynobrathys (Robson, 1990) to be nonmonophyletic 2 , a result reported also in other studies (Nürk and Blattner, 2010;Meseguer et al., 2013;Nürk et al., 2013). The extensive sampling in this study highlights the number of species, which group in a polyphyletic position. What emerges from our analyses is striking evidence that Neotropical species, which group within two distantly related clades, underwent two independent biogeographic histories, as might be expected from their phylogenetic position. Moreover, both lineages exhibit different ecologies and life forms: one containing herbaceous species native to the low-and uplands of South America (Trigynobrathys s.str., subclade Neotropic "lowland"), the other mostly shrubby high Andean species (the Páramo clade). Thus, we assume that dispersal into South America occurred twice in two independent lineages within Hypericum. This is in contrast to Meseguer et al. (2013), who 2 Hence, we refer with Brathys s.l. to sects. Brathys + Trigynobrathys sensu Robson (1977 onwards), and with Trigynobrathys s.str. to the East Asian-North American species (H. japonicum, H. gramineum, etc.) and the species native to lowland and upland areas of South America (H. brasiliense, H. connatum, etc.). Consequently, Brathys refers to sect. Brathys including sect. Trigynobrathys pro parte. revealed a single dispersal to South America (which is due to restricted sampling of New World species in their analysis). Further conclusions demand formal biogeographic analyses on the base of a sound phylogeny with comprehensive species sampling, to reveal the spatiotemporal evolution of species and source areas of dispersals. Similarly, the inferred monophyly of the páramo species of Hypericum, and thus, the hypothesis of secondary evolution of "lowland" species (H. silenoides; placed in derived positions within the Páramo affinis clade) needs to be tested by ancestral area reconstructions incorporating altitudinal variation. Until then, monophyly of páramo Hypericum is the parsimonious hypothesis and the conditional diversification rates (see below) are conservative estimates.
Bayesian relaxed clock analyses of rDNA ITS revealed age estimates, which are in congruence to those inferred from analyses of cpDNA sequence divergence (Meseguer et al., 2013). To the best of our knowledge, this is the first report of divergence times in Hypericum, which have been estimated from nuclear data. According to these estimates ( Figure 3A; Table 2), the stem lineages of the major clades, Myriandra, Trigynobrathys s.str., and Brathys, diverged during the Miocene (5.3-23.0 Ma ago), with an early split of Triadenum from Myriandra+Brathys s.l. in the Oligocene (23.0-33.9 Ma ago), and Myriandra from Brathys s.l. at the Oligocene/Miocene boundary (ca. 23 Ma ago). In contrast, the main diversification within the major clades is revealed to have taken place from the Pliocene (5.3 Ma) onwards (Figure 3). Within the Páramo clade of Hypericum, the inferred time of the onset of species divergence (ca. 3.3-3.8 Ma ago) coincides with the final uplift of the Andes, and thus, with the (early) emergence of the páramos (Gregory-Wodzicki, 2000;Hooghiemstra et al., 2006;Graham, 2009Graham, , 2011b. Moreover, the fossil record for Hypericum in Andean high valleys (pollen fossils; van der Hammen et al., 1973;Wijninga and Kuhry, 1990;Wijninga, 1996) reported for the Late Pliocene (2.5-3.6 Ma ago) is in good agreement with the revealed age estimates. Nevertheless, the used secondary calibration approach (calibrating the nodes of the tree with the estimated ages of another study) is problematic in the sense that the uncertainty produced in the original study is accumulated in our analysis.
Although topology and age estimates inferred from ITS sequence variation is highly congruent to those reported from chloroplast sequence analyses (Meseguer et al., 2013), the multi-copy nature of the 18S-ITS1-5.8S-ITS2-26S nuclear ribosomal cistron potentially confounds species tree reconstruction due to the possible presence of paralogs (derived from gene duplication) or pseudogenes (non-functional copies; Álvarez and Wendel, 2003;Bailey et al., 2003). Intra-individual rDNA polymorphism in Hypericum has been documented by Nürk et al. (2013), suggesting incomplete concerted evolution. If not artifactual, evidence for polymorphic ITS types has also been found in this study (indicated by double peaks in the chromatograms). On the other hand, the polymorphisms observed within an individual do not exceed that expected for a heterozygous individual (e.g., Muir et al., 2001). The direct sequencing approach used here, however, does not permit full investigation of ITS copy variation. While comparison with cpDNA phylogenetics (Meseguer et al., 2013) indicates no cases of ITS paralogy across major clades, the position of several species in multiple subclades of the Páramo clade (e.g., H cardonae, H. juniperinum, H. prostratum; Figure 2), could be evidence that divergent ITS lineages are present within individuals (deep paralogy sensu Bailey et al., 2003; duplication and divergence prior to speciation). Multiple ITS types within species could also result from hybridization and allopolyploidization Doyle, 1998, 2002;Blattner, 2004;Soltis et al., 2008;Kiefer and Koch, 2012). The sympatric occurrence of species (e.g., H. mexicanum with H. juniperinum, or H. strictum with H. tetrastichum, H. prostratum, and H. selaginella) provides the background in which hybridization is possible. Hence, gene trees and reticulate evolution could confound species relationships within the Páramo clade of Hypericum. Also, the possible presence of non-functional ITS copies (pseudogenes) could result in long branch attraction (Felsenstein, 1978), or repulsion (Siddall and Whiting, 1999). If present, pseudogenes affect divergence time estimations by accumulating substitutions that cause increased branch length estimates, which, in turn, result in earlier divergence time estimates. Comparison to future studies analyzing chloroplast sequence divergence or single/low copy genes and employing a comprehensive species sampling will reveal the influence of possibly included pseudogenes and the accuracy of the age estimates provided in the present study. Hence, ages reported here can be used as a conservative estimate, when implemented in diversification rate estimations.
We calculated net diversification rates (r), in which taxonomic richness of clades is considered ( Table 2). The taxonomic likelihood approximation (Magallón and Sanderson, 2001) reveals a two-fold higher net diversification rate for the Páramo clade (mean: 1.1 speciation events per million year), compared to the Trigynobrathys s.str. clade (0.5 sp/Ma), i.e., the second group containing Neotropical species. The inferred rate of net diversification of páramo Hypericum is in the range reported from other high Andean plant groups, e.g., Halenia with 1.0 sp/Ma von Hagen and Kadereit, 2003), Gentianella with 1.7 sp/Ma (von Hagen and Kadereit, 2001), or Lupinus with exceptional 1.9-3.7 sp/Ma (Hughes and Eastwood, 2006). When minimum and maximum rates of Páramo Hypericum are considered (0.75-1.86 sp/Ma), it is similar to the average net diversification rate of the páramos estimated as 1.36 sp/Ma (Madriñán et al., unpublished data).
That is, the Neotropical radiation in Hypericum seems to be related to the emergence of high-altitude habitats, as unexpected high species richness is mainly detected in the Páramo clade, and is less pronounced in the second Neotropical clade containing "lowland" species. We hypothesize that the availability of the "new" páramo habitats have had a causal impact on diversification of Neotropical Hypericum, allowing species to adapt into different niches within these newly emerging ecosystem and to rapidly diversify. A further possibility is that Pleistocene climatic fluctuations, repeated fragmentations of the "sky-island"like páramo habitats promoted allopatric speciation leading to increased diversification (Rauscher, 2002;Hughes and Eastwood, 2006;Moore and Donoghue, 2007). Additionally, the diversification analysis is inference of acceleration from method of moments (Magallón and Sanderson, 2001), as the Yule model assumes constant rates over time. Punctual extinction or decreasing extinction rates might produce a similar pattern (Antonelli and Sanmartín, 2011b;Crisp and Cook, 2011;Stadler, 2013) that might not be associated to an increase in speciation in Hypericum in the páramos. These hypotheses need to be tested and compared to the influence of morphological, physiological, and spatiotemporal patterns to reveal causal cohesive motives underlying the observed extraordinary species richness (de Aguiar et al., 2009;Moore and Donoghue, 2009;Crisp et al., 2011;Stadler, 2011).
Further investigations on the impact of traits and/or events on diversification of páramo Hypericum demands a supported phylogeny for this group, which is not provided in the ITS data. Likewise, the amount and influence of hybridization and/or incomplete lineage sorting needs to be investigated by incorporation of both further nuclear and chloroplast markers (Rieseberg and Soltis, 1991;Jakob and Blattner, 2006;Carine et al., 2007). Comparison of nuclear low copy genes and cpDNA marker will offer to investigate the amount of reticulation, and incomplete lineage sorting/introgression (Peters et al., 2007;Nosil et al., 2009;Brassac et al., 2012) involved in diversification of Hypericum in the Andes.
To summarize, we conclude that only one lineage in Hypericum dispersed and diversified in the páramos. The age estimate for the Páramo clade of Hypericum correlates with the early emergence of high-elevation Andean grasslands. Based on these age estimates, extraordinary diversification is inferred for páramo Hypericum, as species numbers within these high-Andean grasslands are excessively rich. Great phenotypic diversity has evolved in a short time frame (3.3-3.8 Ma) in páramo habitats, although low genotypic differentiation is observed in the nuclear rDNA. Thus, keeping in mind the limitations discussed above, we propose that adaptive radiation-ecological and phenotypic diversity driven by intraspecific selection (on regulatory divergence rather than protein structure; Schluter, 2000;Losos, 2010;Mariac et al., 2010)-under strong ecological pressure caused the morphological diversity. That is, the rapid radiation in páramo Hypericum has likely been promoted by the uplift of the Andes, via adaptive radiation and/or (allopatric) speciation induced by the orogeny and topography of the northern Andean Cordilleras.