Multiple origins of mountain biodiversity in New Zealand's largest plant radiation

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.


| INTRODUC TI ON
Mountains are global hotspots of biodiversity (Antonelli et al., 2018;Hughes & Atchison, 2015;Hughes & Eastwood, 2006), but how they accumulate high levels of species richness remains poorly understood Rahbek, Borregaard, Antonelli, et al., 2019;. Central to answering this question is the extent to which mountains facilitate diversification versus provide opportunities for colonization from elsewhere (Antonelli, 2015;Igea & Tanentzap, 2019). High in situ cladogenesis can accompany mountain uplift because of adaptation to novel niches created by increasing topographic, edaphic and climatic heterogeneity compared with lowland habitats (Antonelli et al., 2018;Hoorn et al., 2013;Quintero & Jetz, 2018;Schwery et al., 2015). For example, adaptive radiation in novel mountain habitats has been invoked to explain the diversity of life-forms in Andean and Rocky Mountain Lupinus Hughes & Atchison, 2015) and niche shifts of New Zealand rockcresses (Joly et al., 2014). Cladogenesis can also be promoted by vicariance associated with range fragmentation that accompanies mountain uplift (Gittenberger, 1991;Steinbauer et al., 2016;Winkworth et al., 2005). By contrast, colonization of open alpine niches by lowland or preadapted lineages (Crisp et al., 2009) may be equally or more important, as in the Himalayas, Hawaii or Mt. Kinabalu, Borneo (Antonelli, 2009;Merckx et al., 2015;Price et al., 2014;Xing & Ree, 2017). The conditions that favour different relative levels of in situ cladogenesis or colonization in generating biodiversity in mountains remain unclear.
The importance of diversification versus colonization for generating the standing diversity found today in the world's mountains may depend on the characteristics of a given mountain range. These characteristics can include the age, elevation, latitude and ruggedness of a given range, as well as its surrounding habitat and biodiversity, spatial and temporal continuity and glaciation history. If uplift accelerates diversification, colonization may be more important only in mountains with slower or older uplift. For example, while in situ cladogenesis is high within the young (~8 Myr old) and rapidly uplifting Hengduan Mountains, diversity in the much older neighbouring Himalayas (up to 40 Myr old) has resulted largely from repeated dispersal of lineages from diverse continental surroundings (Xing & Ree, 2017; though the uplift timing of the Hengduan Mountains is debated, e.g. Su et al., 2019). However, in younger mountains, uplift may be so recent that there has not been sufficient time for extensive in situ cladogenesis. Biodiversity accumulation patterns can similarly vary with the age of habitat even within a mountain range. More recently emerged habitats, such as the high alpine zone, may have fewer species than older lower elevation bands due to having less time for speciation, though this pattern may be disrupted if colonization of the high alpine zone is common (Quintero & Jetz, 2018). Thus, it is worth examining mechanisms of diversity accumulation at finer elevational and habitat scales and for entire mountain ranges. Finally, colonization may remain dominant if mountains provide a weak barrier to gene flow, for example, due to low topographic relief (Antonelli et al., 2018), or if extinction is high, for example, due to high fragmentation or lack of glacial refugia (McCulloch et al., 2010;Winkworth et al., 2005). Other factors influencing colonization rates include the strength of evolutionary priority effects (Tanentzap et al., 2015;Uzma et al., 2019), immigration pressure (Xing & Ree, 2017) and environmental filtering (Edwards & Donoghue, 2013).
As mountain assemblages are complex with multiple lineagespecific histories, focusing on the biogeographic and evolutionary history of a diverse clade can provide insight into the processes generating biodiversity. Here, we study the largest plant radiation in New Zealand, the hebes or Veronica sect. Hebe (Plantaginaceae). This group comprises 124 extant species, of which 96 occupy mountain habitats (Bayly & Kellow, 2006). The core of the radiation (hereafter, core hebe subclade) includes approximately 80 species formerly classed as genus Hebe, alongside several smaller, morphologically and genetically distinct subclades, including the snow, speedwell, sun and semi-whipcord hebes (Meudt et al., 2015). Veronica sect.
Here, our aim is to test processes and understand patterns of mountain diversity in Veronica. We ask three questions: (1) What was the relative importance of in situ cladogenesis and colonization for generating mountain diversity in Veronica at different stages of mountain uplift? (2) Was colonization of mountain habitats associated with niche divergence relative to species arising from in situ 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.

K E Y W O R D S
biodiversity, cladogenesis, climate niche, colonization, historical biogeography, niche filling cladogenesis? (3) What were the relative roles of adaptation to novel ecological niches and vicariance in promoting mountain cladogenesis? Our approach combines historical biogeographical and bioclimatic niche modelling with tests of niche disparity and range overlap using a new, time-calibrated phylogeny of Veronica sect.
Hebe. Together, our work provides new insight into the accumulation of mountain biodiversity by studying a large clade from relatively young mountains.

| 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% (n = 41 loci), HybPiPer paralog warnings for >10% of taxa (n = 42 loci), alignment length <150 bp and <20 parsimony informative sites or average percent identity <65.5%. Further discussion of paralogs resulting from polyploidy is in Thomas et al. (2021). The filtering process left 239 loci from which the final subset was selected as described below.
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 (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 Plantaginaceaewide tree that was calibrated with Plantago fossils along with geomorphological and palaeobotanical data (Meudt et al., 2015). This calibration estimated the posterior distribution of the crown age of New Zealand Veronica as 5.3-10.3 Ma. Accordingly, we used a normal prior (mean = 7.5 Ma, standard deviation = 1 Ma) for the most recent common ancestor of New Zealand Veronica sect. Hebe, excluding the Australian outgroup. Although this prior simplifies the calibration date of Meudt et al. (2015), for example, by rounding and assuming symmetrical errors, the differences are small so are unlikely to bias the estimation of the posterior distribution of ages. We ran a Markov Monte Carlo chain for 1.8 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 R package 'rwty' v1.0.2 (Warren et al., 2017) showed convergence for the species tree topology with ESS >1000.

| Biogeographical modelling
To compare the relative importance of in situ cladogenesis and colonization, we first estimated ancestral habitats of New Zealand Veronica using BioGeography with Bayesian (and likelihood) Evolutionary Analysis in R Scripts ('BioGeoBEARS') v1.1.2 (Matzke, 2013). We performed two analyses, one with two elevational habitat bands, mountain and lowland, and the other with the mountain band subdivided further into subalpine (below treeline) and alpine (above treeline). Species were assigned to elevation bands based on elevational data from herbarium records, field guides and floras (Bayly & Kellow, 2006;NZPCN, 2021), and expert advice (P. Garnock-Jones, pers. comm.). Generally, these bands corresponded with >400 m a.s.l. for mountain, 400 to 1000-1500 m for subalpine, and >1000-1500 m for alpine, with the boundaries between bands varying with latitude (Wardle, 1991). Based on elevations extracted from georeferenced records, the mean elevation (± standard error) was 1214 (±30) m for mountain species, 181 (±30) m for lowland species and 592 (±43) m for species present in both habitats. We excluded from downstream analyses 7 species endemic to outlying islands due to a lack of records and isolation from upland environments. Excluding these species may change some absolute rate calculations in the tree, but should not change relative differences and our main conclusions. BioGeoBEARS was run with six standard models: DEC (dispersal-extinction-cladogenesis), DIVALIKE (dispersalvicariance, likelihood version), BAYAREALIKE (BayArea, likelihood version) and each with jump dispersal (Matzke, 2013). We set a time constraint for the availability of mountain habitat based on uplift of the Southern Alps. Mountain and subalpine habitat were available from 4 Ma, and alpine habitat was available from 1.9 Ma (Heenan & McGlone, 2013). Model support was compared with the Akaike information criterion corrected for small sample sizes (AICc).
Using the best supported BioGeoBEARS model (see Appendix S1 in Supporting Information, Table S1.1-2), we reconstructed in situ cladogenesis and colonization at different stages of mountain uplift given the estimated ancestral habitats.
We generated 100 biogeographical stochastic mappings (BSMs) that estimated colonization and cladogenesis events across habitats in our phylogeny while incorporating model uncertainty (Dupin et al., 2017). Cumulative number of colonization and in situ cladogenesis events for each habitat band were calculated by binning the events from the BSMs into 0.5-Myr intervals.
We calculated rolling per-capita rates for colonization and in situ cladogenesis after Igea and Tanentzap (2019) and Xing and Ree (2017) by dividing the median number of events in a time bin for a given habitat band and event type by the number of lineages present in that habitat in the previous time bin. We accounted for phylogenetic and node age uncertainty by running this analysis for 100 randomly sampled trees from the StarbeaSt2 posterior distribution.
To test explicitly the association between mountain uplift and diversification rates in New Zealand Veronica, we fitted four models that estimate speciation and extinction rates given palaeo-elevation via the R package 'RPANDA' v1.9  and compared these to four simple models dependent on time alone using AICc (Appendix S2; Table S1.3). Palaeo-elevation in the Southern Alps was estimated by interpolating linearly between radiometric dating of rocks over time by Tanentzap et al. (2015).

| Physiological niche modelling
To quantify niche occupancy and the importance of niche adaptation, we first quantified 24 physiological climate traits. We estimated the traits with a process-based species distribution model (SDM) derived from the Thornley Transport Resistance (TTR) model of plant growth (Higgins et al., 2012). The TTR SDM estimates carbon and nutrient uptake and allocation given environmental data associated with species occurrences. The resulting niche dimensions for each species are defined by upper and lower environmental limits on physiological processes (see Appendix S2, Table S2.4 for details). The parameters defining these niche dimensions represent niche traits in downstream analyses. With the R package 'ttr.sdm' v0.5.1, we fitted the SDM for each species using geo-referenced herbarium records supplemented by records from the Global Biodiversity Information Facility (see Supplementary Methods in Appendix S2). A total of 22,652 occurrences were used, with 11-2391 records for each species (mean: 231.1, median: 118; Table S1.6). Environmental data used to fit the models included temperatures from WorldClim v1.4 (Hijmans et al., 2005), soil nitrogen (Shangguan et al., 2014) and soil moisture and solar radiation (Trabucco & Zomer, 2010). We generated random pseudoabsence points equal in number to occurrences after Larcombe et al. (2018). Pseudoabsence points were distributed proportionally across environmental zones classified from the environmental variables with the CLARA algorithm in the R package 'cluster' v2.0.7.
To reconstruct niche occupancy, adaptation and vicariance throughout time, we estimated ancestral niche traits from the TTR parameters. First, to produce a computationally tractable dataset for ancestral niche modelling, we reduced the multivariate niche space to nine parameters using phylogenetic principal component analysis (PCA) in the R package 'phytools' v0.7.47 (Revell, 2012; Appendix S2, Supplementary Methods and Tables S2.4, S2.5).
We could then estimate the ancestral niche space across the tree by fitting multivariate evolutionary models using the R package 'mvMORPH' v1.1.4 (Clavel et al., 2015). We compared the fit of five macroevolutionary models using AICc. Brownian motion (BM, neutral evolution), Ornstein-Uhlenbeck (OU, trait variance constrained around an optimum) and Early Burst (accelerated rates early followed by deceleration) were tested with a single rate parameter or optimum per trait. Branches were assigned to habitats based on the best supported BioGeoBEARS model using BSMs converted to stochastic character maps (Bollback, 2006). Two more complex models allowing multiple BM rates (BMM) and OU optima (OUM) for mountain, lowland and mountain generalist branches were run for 50 BSMs. Ancestral traits could then be estimated at each node with the best supported model (Table S1.7; Appendix S3, Figure S3.1).
Finally, we used the ancestral trait reconstructions to estimate the niche differences among species through time. First, we constructed a multidimensional niche space by calculating a pairwise distance matrix among the estimated ancestral traits for the internal nodes and tips of the tree using the R package 'Claddis' v0.6.3 (Lloyd, 2016). The matrix was subjected to classical multidimensional scaling ordination, retaining the full dimensions of the pairwise matrix so that distance calculations in different subsets of the tree over time would be made within the same trait space (after Guillerme & Cooper, 2018). We then identified the most recent ancestral nodes of all the branches present along a 'time slice' for each node in the phylogeny using the R package 'dispRity' v1.6 (Guillerme, 2018). We only included ancestral nodes in our calculations that were in the same habitat band as the focal node according to the best supported historical biogeographic model. At each time slice, we calculated the multidimensional Euclidean distance between the centroid of that assemblage's niche space and (i) the focal node and (ii) all nodes in the assemblage (overall trait disparity).

| Niche occupancy of colonizing species
To test whether colonizing mountain species have distinct niches from in situ mountain species, we compared niche distances from resident mountain assemblages for focal nodes resulting from colonization versus in situ cladogenesis. We assigned each focal node a status of colonization or in situ cladogenesis based on whether it occupied a different habitat than its most recent ancestral node.
We calculated the difference between median niche distance for the two groups and tested its statistical significance by randomly shuffling the internal nodes 100 times to generate a null model. We derived a p-value by calculating the percentile of the observed difference relative to the distribution of differences calculated from the null model. We repeated this analysis for 50 BSMs. To account for phylogenetic uncertainty, we repeated the analysis for 20 trees from the posterior distribution with 10 BSMs each.

| Adaptation versus vicariance
We first tested if adaptive divergence was important in clade diversification. We estimated disparity through time (DTT) as a proxy for ecological divergence among mountain species, expecting to see an increase in disparity with mountain uplift. DTT was calculated with two complementary methods. First, we derived overall trait disparity from ancestral niche traits estimated at time slices across the phylogeny as described above. We did this by calculating the mean of the trait distances of all nodes in a time slice from their centroid and plotting these values over time. Second, we used a method that does not rely on ancestral state estimation, average subclade DTT.
This approach instead calculates the average disparity of descendant subclades of each node in the tree, standardized by the disparity of all species in the tree (Harmon et al., 2003). Partitioning disparity either among or within subclades can help distinguish between the relative contributions of adaptive or allopatric speciation over time, as can the morphological disparity index (MDI) that measures area between the DTT curve and a null model (Harmon et al., 2003).
Specifically, higher disparity among subclades than within subclades, with DTT values closer to 0 and negative MDI, suggests vicariance has likely been more important for speciation because closely related species separated by complex topography have retained similar niches . Higher disparity within subclades, with values closer to 1 and positive MDI, suggest that niche adaptation has accompanied speciation (Harmon et al., 2003;Liu et al., 2016).
We calculated average subclade disparity through time with the R package 'geiger' v2.0.7 (Pennell et al., 2014) and compared the value to 100 simulations under Brownian motion (Harmon et al., 2003) for the MCC tree and 100 posterior trees.
To compare the importance of ecological adaptation more directly to vicariance in driving cladogenesis, we estimated niche and range overlap between species. Species niche overlap was estimated by first projecting environmental suitability across New Zealand for each species from the fitted TTR SDM. Using Schoener's D, a standard proportional niche overlap index ranging between 0 and 1 (Warren et al., 2008), we calculated the mean overlap between descendants of each node in the phylogeny with the niche.overlap function from the R package 'phyloclim' v0.9.5. The node-based overlap calculations considered phylogenetic relationships by weighting individual pairs by their phylogenetic distance (Fitzpatrick & Turelli, 2006). This produces phylogenetically informed 'node averages' of niche overlap across time. To focus on mountain habitat, we pruned the phylogeny to only species in the mountains or in both mountains and lowlands (86 species) before calculations. To calculate spatial range overlap, we first delineated the observed range for each species by merging 50-km buffers around occurrence points and reducing the outer perimeter of the merged area by 40 km to reduce range overestimation (after Li et al., 2018). Range polygons were rasterized at 5-km resolution, and node-average overlap was calculated with Schoener's D.
To test whether allopatric (vicariant) or sympatric speciation dominated historically, we correlated node-averaged range overlap against node ages (age-range correlation, Fitzpatrick & Turelli, 2006). Dominant allopatric speciation would be supported by a small intercept and statistically significant positive slope (increasing towards the past), reflecting range shifts that have caused higher overlap among older species than younger species (Fitzpatrick & Turelli, 2006). A large, positive intercept and negative trend would suggest dominant sympatric speciation. We also applied this technique to niche overlap to test the effect of speciation mode on the level of adaptive divergence. Compared to age-range correlation, the reverse relationships might be expected. However, a high rate of range drift or niche evolution can mask the temporal signal estimated by these methods (Phillimore et al., 2008). We therefore also compared median node-averaged overlap to a null model expectation for niches and ranges. Further interpretations of this comparison are detailed in Appendix S2, Supplementary Methods.

| Patterns of diversity
We found evidence that both colonization and in situ cladogenesis

| Processes generating diversity
Overall, estimated biogeographic events for Veronica suggest in situ cladogenesis has generated more species in the mountains than colonization, but that its rate slowed over time. After the mountains began to uplift 4 Ma, the median number of mountain cladogenesis events estimated by the two-habitat biogeographic model rose steadily and overtook lowland cladogenesis by ca. 2.5 Ma (Figure 2a).
Colonization events remained low until the last 0.5 Ma, when they increased sharply but still only produced a median of 25 total events F I G U R E 1 One major colonization of mountain habitat by New Zealand Veronica. Reconstructed biogeographic history of New Zealand Veronica in mountain and lowland habitats based on the best supported biogeographic model (dispersal-extinction-cladogenesis; Table S1.1). Tips represent species. Pie charts show probabilities of ancestral habitats for each node. Thicker branches have posterior probability support ≥0.9. Dashed line marks when mountain habitat first became available. Black dots indicate presence on North Island (N) and/or South Island (S). Inset map: elevation contour line at 400 m indicating the division between mountain and lowland regions as defined in this study (projection: New Zealand Transverse Mercator).  (Figure 2a). As a rolling percapita rate, median cladogenesis was similarly higher than colonization in the mountains until the last half-million years (Figure 2b). The three-habitat model was generally similar to the two-habitat model in the subalpine zone ( Figure S3.4b). In the alpine zone, colonization events and rates strongly outpaced cladogenesis ( Figure S3.4a), reflecting colonization from the subalpine. The temporal decline in cladogenesis in both models was also supported by a negative correlation between uplift and speciation rate inferred from a palaeoelevation-dependent model of diversification ( Figure S3.5).
Colonization tended to introduce species with more divergent physiological traits into mountain niches in Veronica. The median niche distance from existing ancestral mountain nodes was, on average, 21.1% (95% confidence interval: 11.4-30.5%) greater for mountain taxa arising from colonization than from in situ cladogenesis under the best supported model of niche evolution with the MCC tree ( Figure 3). Median niche distance was statistically different between in situ and colonizing taxa compared to a null model for 88% Stages of species accumulation 1 2 3

FI G U R E 3 Niches of new mountain species of New Zealand
Veronica arising from colonization are more distant from existing assemblages than species diversifying in situ. Histogram shows median niche distance from ancestral taxa for nodes that appear in mountains either by colonization or in situ cladogenesis. Distances were calculated for each of 50 BSMs from the best supported biogeographic model and were statistically significant for 44 BSMs (p < 0.05).

Non-adapƟve vicariance
ii) High range High niche

Non-adapƟve vicariance + 2°sympatry
iii) Low range Low niche Non-adapƟve vicariance + adapƟve divergence iv) High range Low niche AdapƟve sympatric speciaƟon to 2.5 Ma) resulted from lowland ancestors colonizing newly available subalpine habitat and undergoing high rates of in situ cladogenesis. Novel habitat expansion has been suggested as an important trigger for radiation of lineages that successfully colonize vacant niches Hughes & Atchison, 2015;Moore & Donoghue, 2007;Schwery et al., 2015). The second stage in our study (ca. 2.5 to 0.5 Ma) was characterized by in situ cladogenesis dominating over colonization but slowing over time. The saturation of niche trait disparity suggests that the decrease in cladogenesis rate may be due to niche filling as competition increased between close relatives Gavrilets & Losos, 2009;Price et al., 2014;Rabosky, 2013;Roquet et al., 2013). An alternate explanation for this pattern could be undetected speciation and/or extinction of earlier colonizers. Finally, colonization became important again as cladogenesis continued to slow (ca. 0.5 Ma to present).
During this stage, the emergence of high alpine habitat increased rates of alpine colonization but not in situ cladogenesis, perhaps due to the extinction risk of increasing specialization, competition from generalists and/or lack of time (Heenan & McGlone, 2013). The importance of cladogenesis in the first and second stages of our proposed diversity model ultimately requires novel habitat for species radiation. As the Southern Alps formed, the resulting rocky, exposed habitat expanded providing new opportunity for in situ cladogenesis in Veronica (Heenan & McGlone, 2013). Montane habitat below treeline has proven similarly important for diversification elsewhere, including in Lupinus and Ericaceae Schwery et al., 2015). Allopatric speciation is likely a common precursor to ecological adaptation to novel habitats, even in radiations with sympatric species (Aguilée et al., 2011;Rundell & Price, 2009), for example, due to secondary sympatry. Our results in Veronica support such a process of vicariance followed by range shifting contributing to high cladogenesis given the combination of higher-than-expected niche and range overlap. In the Southern Alps, repeated glaciation could have caused higher elevation populations to descend independently into isolated valleys (refugia) in diverse habitats unable to support forests. These populations could then have speciated allopatrically and rejoined in secondary sympatry at higher elevations during glacial retreat (Winkworth et al., 2005). High average subclade disparity through time also points to high climate niche lability in ancestral Veronica species (Harmon et al., 2003), such as for adaptation to a fluctuating climate (Hua & Wiens, 2013;Ogburn & Edwards, 2015). This lability may help explain how Veronica persisted in the mountains despite climate instability and extinction risk for specialized species (Hua & Wiens, 2013). Lability is also evident in the recolonization of the lowland by a mountain ancestor in the North Island, followed by diversification. Several members of this North Island subclade are highly endemic, allopatric coastal species, suggesting that vicariance may have also accompanied lowland speciation.
Deceleration of in situ cladogenesis during the second and third stages of our radiation shows that cladogenesis is facilitated but not always prolonged by geological uplift. Orogeny has been associated with increasing diversification in groups such as Andean orchids (Pérez-Escobar et al., 2017) and Lupinus  and larger correlative studies Xing & Ree, 2017).
By contrast, the slowdown in Veronica occurred despite increased uplift and the emergence of the alpine zone ( Figure S3), more resembling the classic pattern of an adaptive radiation where niche filling constrains cladogenesis (Price et al., 2014;Skeels & Cardillo, 2019).
Other mountain radiations also display this slowdown despite their initial response to uplift, suggesting that lineages with accelerating diversification may be in earlier stages of radiation characterized by high turnover of species that soon go extinct (Hughes & Atchison, 2015;Linder, 2008;Roquet et al., 2013). Speciation slowdowns can be attributed to factors other than niche filling, such as filling of geographic space (Moen & Morlon, 2014); however, the niche filling hypothesis is congruent with the saturation of climate niche disparity we observed in Veronica (Aristide & Morlon, 2019).
The high alpine zone, rather than providing new opportunities for in situ cladogenesis, chiefly promoted range expansion of generalist subalpine species (Heenan & McGlone, 2013). The main exception to this pattern was the alpine specialist snow hebe subclade, which comprises several cushion species, possibly reflecting the importance of this key innovation for inhabiting the harshest alpine environments (Boucher et al., 2016;Roquet et al., 2013). The dominance of generalist species further suggests that established species with wide thermal niches were highly competitive even as new habitat expanded (Tanentzap et al., 2015). Extinction may thus be higher for range-limited specialists in the frequently disturbed, fragmented alpine, but this is unlikely to explain the overall slowdown, as it began before the alpine emerged. Another explanation is that some cladogenesis is too recent to detect due to incomplete reproductive isolation, leading to deflated rates towards the present (Egan & Crandall, 2008;Etienne & Rosindell, 2012). However, species in smaller radiations have been dated to within the last million years in the Southern Alps (Heenan & McGlone, 2013;Joly et al., 2014).
During the third stage of our model for generating mountain diversity, colonization may have become important partly because it introduced species with more divergent traits into niches that were progressively filled. Although niche filling typically reduces invasion potential (Stegen et al., 2012), divergent colonizing species were likely able to exploit marginal habitat made available by a lag in reexpansion by in situ species from glacial refugia (Dullinger et al., 2012;Guo et al., 2015), or by upward shifts of in situ generalists into the alpine zone. Alternatively, in situ species may have been replaced by colonizing species within filled niches, as predicted by taxon cycle dynamics (i.e. lineages wax and wane in dominance even as overall species richness remains constant; Reijenga et al., 2021). Another possibility is that earlier successful colonizations were obscured by later extinction or migration during glacial cycles, which current biogeographical models cannot capture. Future modelling may simulate more complex dispersal processes (Hackel & Sanmartín, 2021), and clarify what drives mountain colonization.

| Limitations and prospects for future research
Our conclusions are subject to several assumptions and potential limitations. First, the phylogenetic relationships and ages underlying our analyses are based on a simplified model of evolutionary rates, which is generally unavoidable using Bayesian methods for a dataset of this size. This reflects a trade-off between more sophisticated estimation methods that better incorporate uncertainty (e.g. StarBEAST2) and the number of assumptions modelled. Second, this model does not explicitly consider possible complexities introduced by polyploidy or introgression. Polyploidy has arisen repeatedly in Veronica sect. Hebe (Albach et al., 2008). Although potential paralogs were removed and phylogenetic uncertainty was considered in downstream analyses, some low-support relationships may have introduced noise into the analysis. This would have not affected the trend observed in biogeographical analysis, but may have reduced the power of ancestral reconstructions. Finally, the discrete time periods assigned to habitat categories are an obvious simplification, as mountain uplift varied both spatially and temporally. Nevertheless, we found that exposure to increased relief after a certain time was still important at a broad scale. Advances in phylogenetic and biogeographical methods that increase computational efficiency and incorporate more realistic levels of complexity will be important for further refining these hypotheses Bromham et al., 2018).
Other factors besides adaptation along macroclimatic axes or vicariance may have contributed to Veronica cladogenesis. For example, biotic interactions such as competition with other taxa excluded in our study may have filled niches in Veronica habitat, constraining colonization and cladogenesis or increasing extinction (Tanentzap et al., 2015). Future studies combining multi-lineage biogeography and trait evolution with biotic interactions would thus be valuable Weber et al., 2017). Ecological conditions specific to New Zealand, such as the absence of grazing mammals and fire regimes through evolutionary time (Antonelli et al., 2011), could also have promoted diversification. Veronica may have been consumed by now-extinct mega-herbivores (Wood et al., 2020), which could have driven differences in extinction or defence mechanisms.
Pollination and seed dispersal interactions, by contrast, are unlikely to have been important, as Veronica pollinators tend to be generalists and seeds are wind or gravity dispersed (Bayly & Kellow, 2006).
Cladogenesis can also be promoted by whole genome duplication (Meudt et al., 2015). Polyploidy has been associated with greater island diversity generally , and in Veronica, potentially by promoting range expansion into new biomes (Liddell et al., 2021). However, extrinsic factors such as habitat area tend to be more important than intrinsic traits in explaining diversification rates (Vamosi & Vamosi, 2010). Local abiotic niche axes such as soil type, microclimate and microtopography (e.g. aspect and slope) at scales much smaller than that of widely available environmental data are also known to influence alpine plant distributions and extinction dynamics (Kulonen et al., 2018). While process-based niche modelling such as used here captures important limits on plant growth, its accuracy for alpine niches could be improved with very-high resolution climate downscaling (Carlson et al., 2013). Finally, advances in integrating genomic data with landscape change are a promising avenue for further teasing apart cause and effect in mountain diversification (Dolby et al., 2022).
Through past climate oscillations, both in situ cladogenesis and colonization maintained diversity in the mountains. In the future, as anthropogenic climate change advances, both processes will likely continue to operate on a clade level, but their relative roles may be harder to predict. Extinctions or range contractions may create open niches, reducing the effect of niche filling (Ricklefs, 2010), but whether in situ cladogenesis or colonization respond will depend on how these processes generate species with adaptive traits in pace with climate change. Our results suggest that alpine specialists may not keep up with colonizers responding to climate change, especially generalists already present in montane habitats but not previously seen in the alpine zone (Alexander et al., 2015). Consequently, overall species diversity, functional diversity, species interactions and evolutionary potential may all decline (Thomas et al., 2004;Trisos et al., 2020;Urban et al., 2012). Future work incorporating novel competitive interactions alongside niche traits and biogeographical processes (Jackson et al., 2009;Svenning et al., 2014;Wisz et al., 2013) will be crucial for predicting the fate of alpine biodiversity in a changing world. War Memorial Museum Herbarium (AK).

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest. to preparation of sequence data and alignments) and https://github. com/annet homas/ veron ica-biogeo (for scripts related to biogeographical and species distribution modelling analysis).