Aim: How mountains accumulate species diversity remains poorly understood, particularly the relative role of in situ cladogenesis compared with colonization from lower elevations. Here, we estimated the contributions of in situ cladogenesis and colonization in generating biodiversity of a large mountain plant radiation and determined the importance of niche adaptation and divergence in these processes. We expected cladogenesis would accompany novel habitats formed by mountain uplift but colonization would become more important with time as dispersal opportunities accrue.

Location: New Zealand, Southern Alps

Taxon: Veronica sect. Hebe (Plantaginaceae)

Methods: We estimated the most complete time-calibrated phylogeny to date for Veronica sect. Hebe to quantify rates of in situ cladogenesis and colonization of mountain habitat based on historical biogeographical models. We used environmental niche modeling to quantify species’ climate niches and estimate niche disparity and divergence over time.

Results: In situ cladogenesis generated more species in the mountains than colonization from lowlands. Whereas cladogenesis slowed over time, colonization increased, especially in the alpine zone. Both adaptive ecological speciation along climate niche axes and non-adaptive, vicariant speciation contributed to cladogenesis. However, climate niche disparity through time became saturated, suggesting competition for niche space was important. Colonization brought more divergent species into mountain niches.

Main Conclusions: We suggest mountain diversity accumulates through three main stages: high cladogenesis after initial colonization, decreasing cladogenesis with increasing competition, and increasing colonization after niches saturate, likely promoted by niche divergence. Combining lineage and mountain uplift trajectories, these stages provide a conceptual model to understand how diversity accumulates elsewhere. Assuming these deep-time findings apply to anthropogenic conditions, alpine specialists could struggle to outcompete colonizers facilitated by climate change, especially from generalist clades. Considering novel competitive interactions alongside niche traits and biogeographical processes will be crucial for predicting the fate of alpine biodiversity in a changing world.

Phylogeny

We sampled 115 of 124 New Zealand species of Veronica by adding 38 newly sequenced species to the 77 in Thomas et al. (2021). We followed Thomas et al. (2021) for sample collection, DNA extraction, and sequencing. Briefly, we used target capture enrichment on recent herbarium and fresh collected samples using the Angiosperms353 bait kit that targets 353 low-copy nuclear genes (Johnson et al., 2018; McDonnell et al., 2021). Samples were sequenced on an Illumina NextSeq 500 with 2×150-bp paired-end reads. We processed sequences as in Thomas et al. (2021), retrieving target loci and flanking intron regions with HybPiper v1.2 (Johnson et al., 2016). We trimmed and filtered alignments using the “full” filtering scheme described in Thomas et al. (2021), which consisted of removing markers with taxon occupancy <50%, alignment length less than 150 bp and fewer than 20 parsimony informative sites, or average percent identity <65.5%.

To minimize gene tree conflict observed during phylogenetic reconstruction with maximum likelihood methods by Thomas et al. (2021), we used a Bayesian multispecies coalescent (MSC) model to estimate a time-calibrated phylogeny with StarBEAST2 (Ogilvie et al., 2017). We first used Sortadate (Smith et al., 2018) to reduce the Angiosperms353 dataset to the 50 most informative supercontigs (loci plus their flanking introns), as the full dataset was computationally intractable. Genes were selected by averaging the ranks for gene tree length (as a heuristic for information content) and congruence with a maximum likelihood species tree estimated by IQTREE v1.6 (Nguyen et al., 2015; see Thomas et al., 2021). We partitioned the MSC model by locus and used a HKY site model with empirical base frequencies, strict clock with an estimated rate, and the Calibrated Yule speciation model (Heled & Drummond, 2012), all with default priors. Because there are no accepted Veronica fossils in New Zealand, we relied on secondary calibration from a Plantaginaceae-wide tree with Plantago fossils (Meudt et al., 2015). This calibration estimated the crown age of New Zealand Veronica as 5–10 Ma. Accordingly, we used a normal prior (mean = 7.5 Ma and standard deviation = 1 Ma) for the most recent common ancestor of the New Zealand Veronica sect. Hebe, excluding the Australian outgroup. We ran a Markov Monte Carlo chain for 1.6 billion generations in the CIPRES Science Gateway computing facility (Miller et al., 2015) and confirmed all relevant parameters had effective sample sizes (ESS) > 200 in Tracer v1.7.1 (Rambaut et al., 2018) after a burn-in of 1 billion generations. The ‘RWTY’ R package (Warren et al., 2017) showed convergence for the species tree topology with ESS > 1000. We sampled 1000 trees from the posterior distribution and calculated the Maximum Clade Credibility tree in TreeAnnotator v2.6.3 (Bouckaert et al., 2019).

Alignments: text files; can be viewed with Geneious or Integrative Genomics Viewer

Tree files: text files; can be manipulated with R package ape or viewed with FigTree; posterior tree distribution can be visualized with DensiTree

Beauti XML: xml file; can be loaded with Beauti or run with BEAST2 software