Timing and ecological priority shaped the diversification of sedges in the Himalayas

Study area and taxa selection

For this study, 55 specimens representing 19 species were collected from Western Himalayas range (lies in Pakistan) between 2011 and 2017 (Table S1A). The collecting areas were selected as mentioned in Flora of Pakistan by Kukkonen (Kukkonen, 2001) for the genus. All specimens collected by the lead author (Uzma) were submitted to the herbaria of The Morton Arboretum (MOR), Lisle, Illinois, USA and Pakistan Museum of Natural History (PMNH), Islamabad, Pakistan. Additionally, 27 specimens of 21 species received from E (Royal Botanic Garden Edinburgh), MO (Missouri Botanical Garden) and MSB (Botanishche Staatssammiung München) herbaria were also included (Table S1B), comprising a total of 40 species from the Himalayas from which the new sequences were obtained for this study. In addition, we retrieved sequences for 944 species distributed worldwide from the Global Carex Group (2016) data matrix (available at Dryad: http://dx.doi.org/10.5061/dryad.k05qb). The final dataset included 56% (105 species) of Himalayan taxa out of the total 189 Himalayan species (Table S2). Four outgroup species (Eriophorum vaginatum, Scirpus polystachyus, Trichophorum alpinum and Trichophorum caespitosum) were included from the tribe Scirpeae (Léveillé-Bourret et al., 2014) for phylogenetic analyses. All data and scripts for the analyses are deposited in GitHub (https://github.com/uzma-researcher/Himalayan-Carex-Diversification; https://github.com/uzma-researcher/Himalayan-Carex-Climatic-niche-Evolution).

Molecular methods

Total genomic DNA was isolated from silica-dried leaves of collected specimens and dried leaves of herbarium specimens using DNeasy Plant Mini Kit (QIAGEN, Valencia, California, USA, catalog # 69106) in the laboratory at The Morton Arboretum, Lisle, USA. Amplifications of two nrDNA regions, the internal transcribed spacer (ITS) and 5′  end of the external transcribed spacer (ETS), were performed using primer pairs ITS-IF and ITS-4R (Urbatsch, Baldwin & Donoghue, 2000; White et al., 1990) and ETS-1F and 18S-R (Starr, Harris & Simpson, 2003) respectively. The chloroplast (cpDNA) region matK was amplified using two primer pairs: matK-2.1F and matK-5R (Kew Royal Botanic Garden, http://www.kew.org/barcoding/protocols.html) in the first step and matKF-61 and matKR-673 (nested primers; Global Carex Group, 2016) in the second step with slightly increased in annealing temperature (from 45. 0 °C to 51.0 °C). The PCR reaction mixture of 25µL contained: 2.5 µL 10X MgCl₂-free Taq buffer, 2.5 µL MgCl₂, 1.25 µL DMSO, 0.25 µL BSA, 0.25 µL of each primer at 20 mM, 0.25 µL Taq DNA polymerase (1.25 units), and 1 µL of genomic DNA as a template in amplifications of each region. These regions were selected based on suitability for wide-scale as well as fine-scale phylogenetics in the genus (Starr, Harris & Simpson, 2003; Starr, Naczi & Chouinard, 2009; Global Carex Group, 2016). Amplification cycles for ITS and ETS regions followed the conditions mentioned by Hipp et al. (2006) with minor adjustments in annealing temperature to get appropriate amplicons. However, two-step amplification of matK region involved PCR conditions: 95.0° for 1:00; 30 cycles of: 95.0° for 0:45, 45.0° (51.0° for internal region) for 0:45, 72.0° for 1:30; 72.0° for 3:00. Amplified regions were cleaned and then sequenced following the conditions as described by Begley-Miller et al. (2014) at Pritzker DNA laboratory, The Field Museum, Chicago, USA.

Phylogenetic analyses and time-calibrated molecular phylogeny

The resulting new sequences for ETS, ITS and matK regions were edited and assembled in Geneious version 9.1.6 (Biomatters, Auckland, New Zealand, available from http://www.geneious.com). Sequences were aligned for each region with Global Carex Group (2016) sequences using MUSCLE (Edgar, 2004) as implemented in Geneious. The matrices for each region were trimmed or N-filled to maintain equal characters in all sequences. The best-fit models of molecular evolution were estimated based on the Akaike information criterion implemented in jModelTest2 v.2.1 (Darriba et al., 2012; Table 1). The three DNA regions (ETS, ITS and matK) were first analyzed separately using maximum likelihood (ML) method as implemented in RAxML HPC-PTHREADS-SSE3 version 8.2.4 (Stamatakis, 2006) with model of evolution GTR+I+G (General Time Reversible model with gamma distributions and invariant sites) and node support through 1000 rapid (fast) bootstraps to assess congruence among these nuclear and chloroplast DNA regions. The three matrices (ETS, ITS and matK matrices) were then concatenated into a single matrix (with gaps or missing data). The combined matrix was analyzed using ML in RAxML with model of evolution GTR+I+G and 1,000 fast bootstraps to assess the phylogenetic relationship of the Himalayan Carex. These phylogenetic analyses were performed using the CIPRES Science Gateway v3.3 platform (Miller, Pfeiffer & Schwartz, 2010).

The Himalayan Carex species were non-monophyletic within each of three major clades (Vignea, core Unispicate, core Carex) in the resulting phylogenetic tree. Therefore, the Shimodaira-Hasegawa (SH) test (Shimodaira & Hasegawa, 1999) was performed as implemented in RAxML (version 8.2.4) to evaluate how strongly monophyly of Himalayan lineages was rejected, using five separate constraints that represent alternative scenarios for partial or complete monophyly of Himalayan lineages (Fig. S1; Table S3). Our null hypothesis (H₀) for each test was polyphyly of Himalayan lineages within the major clade (or taxon set) being tested, while H_a was monophyly of Himalayan lineages within that clade. The ML tree generated under each constraint was compared to the unconstrained ML tree.

Divergence times were estimated on the ML tree (constructed on the concatenated dataset) after excluding outgroups using penalized likelihood (Sanderson, 2002) as implemented in treePL (Smith & O’Meara, 2012), which is designed for large phylogenies (Spalink et al., 2016a). We followed Jiménez-Mejías et al. (2016) and employed Carex colwellensis Chandler which was the oldest reliable known fossil ascribable to Carex, dated back to the Priabonian, late Eocene. The crown node was fixed at 37.8–33.9 Ma according to previous reconstructions by Mıguez et al. (2017), an age compatible with Carex colwellensis (Jiménez-Mejías & Martinetto, 2013; Jiménez-Mejías et al., 2016). Penalized likelihood calibration was performed after optimizing priming and smoothing values using a χ² test and cross validation. We assessed run convergence using the “thorough” option and inspected the tree visually using FigTree v1.4.2. The SH test and treePL analyses were performed using supercomputer facility at Research Center for Modeling and Simulation (RCMS), National University of Sciences and Technology (NUST), Islamabad, Pakistan.

Biogeographical distribution and coding

For biogeographical analyses, we considered division of species into two groups: Himalayan and non-Himalayan, according to their presence or absence in the Himalayas (Table S2). Species distributions were based on the World Checklist of Selected Plant Families (Govaerts et al., 2018+), Flora of Pakistan (Kukkonen, 2001) and Flora of China (Dai et al., 2010). However, many species present predominantly in the Himalayas are also present in adjacent regions, making it difficult to study diversification patterns within the mountains. Therefore, we assessed sensitivity of our analyses to alternative biogeographic coding strategies by coding taxa in two different ways: “narrow sense” biogeographic coding, which treats the Himalayas as following strict Himalayan boundaries; and “broad sense” biogeographic coding, which includes the adjoining mountains (Karakoram, Hindu Kush, Tibet plateau, and Hengduan mountains) as part of a broadly construed Himalayas (Table S2). Secondly, to estimate ancestral ranges at the global scale, we scored taxa according to ten ecozones/biogeographical realms based on published distributions. However, our coding diverges from traditional coding (Udvardy, 1975) in two regards: (1) we split the Palearctic region into Western and Eastern to separate Europe, north Africa and western Asia from eastern and central Asia; and (2) we treat the Himalayas in the broad sense as a tenth region (biogeographic coding in Table S4).

Georeferenced data for each Carex species available at GBIF data portal (http://www.gbif.org/) were retrieved using species synonyms to gather all records under their current taxonomic names (Global Carex Group, 2015), listed in Table S5. We downloaded GBIF data using the rgbif package (Chamberlain, 2017) in R (R Core Team, 2013). The largest number of records (17,841) was retrieved for species Carex lasiocarpa, followed by Carex canescens with 16,825. 41 species had only one record and 324 species had fewer than 50 records. We cleaned data to remove specimens georeferenced outside the reported range of the species (according to Govaerts et al., 2018+) and eliminated duplicates using the dismo (Hijmans et al., 2011), ape (Paradis & Schliep, 2018) and magrittr (Bache & Wickham, 2016) packages in R. Species data were also thinned to exclude records within 1.0 to 1.6 km of each other (following Hipp et al., 2018). The post-thinning dataset comprised 965,556 occurrence records for 850 unique species, each represented by an average of 1,136 specimens.

Geographic-dependent diversification and extinction

The GeoSSE model (Goldberg, Lancaster & Ree, 2011) as implemented in the diversitree R package (FitzJohn, 2010) was used to assess historical differences in the rate of speciation, dispersal and extinction in the Himalayas versus non-Himalayan regions in narrow sense (total 966 species of which 10 are endemic to the Himalayas, 861 are endemic to non-Himalayas, 95 are present in both regions) and broad sense (51 endemic to Himalayas, 838 to non-Himalayas, 77 species present in both regions) biogeographic coding (using Table S2 for geographic states coding to narrow and broad sense coding and resulting dated phylogeny). Diversification rate was estimated in three states: A representing species endemic to the non-Himalayas, B representing species endemic to the Himalayas, and AB representing species present in both regions. In this model, we estimate speciation rates of species in region A (sA) and B (sB), as well as sAB, speciation rates for taxa that give rise to two daughter species, one in each region. Likewise, geographical range expansion from A or B to AB were estimated with rates of dispersal dA and dB respectively and range contraction with rates of extinction xA and xB.

Full and constrained alternative GeoSSE models were tested: the full model, in which rates of speciation, extinction and dispersal differed among two regions; a constrained model in which speciation in taxa distributed across two regions was set to zero (sAB ∼ 0); and a model in which rates of speciation and extinction were constrained to be the same between regions (sA ∼ sB; xA ∼ xB). Models were compared using the Akaike Information Criterion. The best-fit model was also fitted using Bayesian Markov chain Monte Carlo (MCMC) for 100,000 generations and an exponential prior to estimate parameter distributions.

Biogeographical reconstruction and ancestral range estimation

Biogeographic history was investigated as transitions among the ten regions initially coded, excluding Antarctica and Oceania for computational reasons. The remaining regions were the Himalayas, Indo-Malaya, Eastern Palearctic, Western Palearctic, Nearctic, Afrotropic, Neotropic and Australasia. The DEC model was utilized as implemented in Lagrange (Ree et al., 2005; Ree & Smith, 2008) and BioGeoBEARS (Matzke, 2014) using BSM (Biogeography stochastic mapping). The DEC method implements the maximum-likelihood approach and allows vicariance, range expansion (dispersal) and range contraction (extinction) processes with inclusion of different parameters. The BSM analysis is based on the Bayesian MCMC approach and simulates the biogeographical history of the events (anagenetic and cladogenetic events) along the branches of the tree. In the BSM analysis, in total 1,000 stochastic mapping replicates with 50,000 maximum trees per branch were conducted on ML tree. The alternative biogeographic model implemented in BioGeoBEARS was not considered, as the inclusion of the jump parameter is not directly testable relative to non-jump models (Ree & Sanmartín, 2018) and introduces complexities that are not necessary to explain our biogeographic scenarios.

Climatic niche modelling

To characterize climatic envelopes for the species for which specimen data were available, we extracted 19 bioclimatic variables from WORLDCLIM data using the R package raster (WorldClim v1.4; Hijmans et al., 2005). The average values for species occurrence data and 19 bioclimatic variables from a total of 965,556 data records for 850 unique species were estimated. The obtained data for 838 species were further proceeded after removing outliers from WORLDCLIM data. For each bioclimatic variable, ancestral character states were estimated assuming a Brownian motion trait evolution process using the ‘fastAnc’ function in phytools (Revell, 2012). Then non-metric multidimensional scaling (NMDS) ordination was employed using Euclidean distances normalized to unit variance on bioclimatic data for the tips as well as the internal nodes of the tree. The stress from K = 1 to K = 10 was calculated and plotted against dimension to estimate how much additional information was extracted from the bioclimatic data with each additional NMDS axis.

Data were plotted with internal nodes of interest, representing ancestors of Himalayan radiations, colored for identification. The Himalayan taxa of the core Unispicate clade are represented by two major subclades, with Kobresia (40 species) sister to a clade dominated by Uncinia (30 species). Placement of the Kobresia ancestors were compared with ancestors of all other Himalayan radiations to evaluate whether the Kobresia ancestors may have represented an ecologically specialized (preadapted) species to Himalayan climatic conditions. Alternatively, the earlier origin of Kobresia species in the Himalayas would support clade age (priority effect and possibly competitive exclusion) as an explanation for success of Kobresia relative to other more recent Himalayan radiations. Niche analyses were conducted using the phytools and vegan (Oksanen et al., 2017; Revell, 2012) packages in R.

However, the analysis above presupposes a Brownian motion model of niche evolution, which is unrealistic in the face of natural selection (Butler & King, 2004; Hansen, 1997). Underparameterizing the model of niche evolution risks inflating Type II error (incorrectly failing to reject the null hypothesis of no difference among clades) if clades tend to adapt to new adaptive regimes to which species migrate. To address this limitation, we also analyzed our niche data using a multivariate Ornstein–Uhlenbeck (O-U) model implemented in the PhylogeneticEM package of R to test transitions in selective regime in an information theoretic framework (Bastide, Mariadassou & Robin, 2017; Bastide et al., 2018). Analyses were conducted using a scalar O-U process and default search settings for the PhyloEM function.

Phylogenetic relationship of Himalayan species and estimation of Himalayan clade divergence in tribe Cariceae

New sequences obtained for this study are deposited in GenBank (Table S6), and details on each of the DNA matrices are reported in Table 1. The matK tree topology (File S1 for matrix) showed low support for most relationships (File S2). Similarly, the phylogeny obtained from ITS (File S3 for ITS matrix) showed phylogenetic incongruence with previous studies (Global Carex Group, 2015; Waterway, Hoshino & Masaki, 2009; Waterway et al., 2016), with low support for the major clades (File S4). The ETS region (File S5 for matrix) was mostly well-defined for major clades with high to moderate supports (File S6). The three matrices with least missing data for one individual per species were concatenated into a single matrix of 970 sequences (outgroups included) with 2,080 sites (File S7), containing 48.7% missing data coded as gaps. The topology of the tree based on combined DNA datasets was highly congruent with previous studies (Global Carex Group, 2015; Waterway, Hoshino & Masaki, 2009; Waterway et al., 2016) and supported at most major clades. All the four major clades (Schoenoxiphium, Vignea, core Unispicate and core Carex) sister to Siderostictae clade were recovered (File S8). The Himalayan species (105) fall into three major clades of Carex: (1) Vignea, (2) core Unispicate and (3) core Carex. The Himalayan taxa (16 species) in the Vignea clade were largely spread out in the clade, not sister to other Himalayan taxa, though two sister species pairs were observed. Himalayan species in core Unispicate clade (30 species) mainly belong to Kobresia, except for two unrelated species formerly classified under section Leucoglochin, arrayed in two clades: Kobresia clade 1, which exhibited only three species (minor clade), and Kobresia clade 2 (major clade) with 27 Himalayan species (File S8). Himalayan taxa (59 species) in core Carex occurred as small clades or isolated (single) species dispersed throughout the clade. Our topology may be poorly supported in some clades, but the groups and branching we retrieved was fully compatible with the Global Carex Group (2016) tree, where the clades retrieved were more strongly supported.

Date calibrations (Fig. S2) are also congruent with previous reports (Escudero et al., 2012; Waterway, Hoshino & Masaki, 2009; Waterway et al., 2016), which is expected given that they are based on the same sources (performed on ML tree). The earliest Himalayan diversification is the Kobresia radiation in the core Unispicate clade at the early Miocene (around 20.6 Ma), followed by diversifications in core Carex (early Miocene, 14.4 Ma) and Vignea (late Miocene, 7.2 Ma; Table 2) clades. Although all Himalayan taxa diversified in three major clades in the epoch Miocene, the timing of diversification of Himalayan Carex was different in three major clades, which suggested multiple origins and also non-monophyletic lineages. Our time calibrations do differ from recent work (Spalink et al., 2016a; Spalink et al., 2016b) due to the fact that we use a calibration based on a more recent date, compatible with the fossil record (Jiménez-Mejías et al., 2016).

Monophyly of the Himalayan taxa within each major clade was strongly rejected at the 0.01 level based on the Shimodaira-Hasegawa test, except for the Kobresia species in the core Unispicate clade, where monophyly cannot be rejected (Table 3). We focused an additional test on the monophyly of Kobresia, in which Kobresia was constrained to be monophyletic and only the Unispicate taxa included. In this test, monophyly of Kobresia was not significantly rejected (ΔlnL = −1.32). Thus, while Kobresia was found to be polyphyletic in previous work (Global Carex Group, 2016) and in the current study, we considered this result poorly supported by the data, pending additional study.

Ancestral range reconstruction

Initial biogeographic reconstruction under the DEC model was conducted on the “narrow sense” (File S9) biogeographic coding, which failed to recover any clade endemic to the Himalayas, which is at odds with our observation that numerous species are predominantly Himalayan in origin. We thus present here analysis of data coded in the “broad sense” (see methods; File S10), using the stochastic mapping (BSM) method on the DEC model as implemented in BioGeoBEARS. We generated 1000 stochastic maps in every 50,000 trees per branch in BSM analysis. In each stochastic map, total event counts of cladogenetic events were higher which were based on vicariance (5.6% events) and sympatric (94.4% events) processes than that of anagenetic which were based on only dispersal events (Table S7). However, an interesting finding in BSM results was that none of the events showed “range switching dispersal”, while all the observed dispersal was of “range expansion”.

Under the broad sense biogeographic coding through BSM, 128 Himalayan taxa were distributed among three major clades: (1) Vignea (16%), (2) core Unispicate (29%) and (3) core Carex (54.6%) (Fig. 1). The highest number of dispersal events into the Himalayas was obtained from Eastern Palearctic (Mean = 39.42; s.d. = 5.23) followed by Nearctic (Mean = 21.44; s.d. = 4.31) and Western Palearctic (Mean = 14.48; s.d. = 3.63) regions. The lowest number of dispersal events was obtained from three regions into the Himalayas; Neotropic (Mean = 3.60; s.d. = 2.47), Australasia (Mean = 3.64; s.d. = 1.83) and Indo-Malaya (Mean = 5.22; s.d. = 2.16). The in situ diversification within the Himalayas that precedes dispersal from the Himalayas region correspondingly into the Eastern Palearctic with the maximum average dispersal events (Mean = 35.26; s.d. = 4.11), however, the following two regions were Indo-Malaya and Western Palearctic (Means = 12. 88 and 12.48, s.d. = 3.31 and 3.72, respectively) (Table S7).

The oldest colonization into the Himalayas, at the crown of the Kobresia clade (comprising 37 Himalayan species out of 128 species in total based on broad sense), occurred during the early Miocene (ca. 20 Ma). Taxa from the Vignea clade originate from the Nearctic and Western Palearctic and Eastern Palearctic (Fig. 1). All origins are reconstructed as arising from relatively recent colonization of the Himalayas (started in late Miocene, ca. 7.2 Ma and last colonization appeared very recently at ca. 0.09 Ma and ca. 0.35 Ma by Carex physodes and Carex canescens, respectively). In the core Unispicate and core Carex clades, by contrast, the Himalayan species cluster into around 12 small clades of 2 to 33 species (File S10). The Himalayan species from the core Unispicate clade derived mostly from Nearctic ancestors with one radiation deriving from the Eastern Palearctic (Fig. 1, File S10). The colonization out of the Himalayas were observed into the Eastern and Western Palearctic, Indo-Malaya, Neotropic, and Nearctic regions around 10 million years ago. While in core Carex, the Himalayan species are arranged in small sub-clades which arose around the early Miocene, there was more recent (ca. 8.5 Ma) in-situ diversification within the Himalayas that served as a source for other regions (File S10). The ancestral ranges of the Himalayan taxa were predominantly in the Eastern Palearctic, Western Palearctic and Nearctic regions. However, in core Carex, Himalayan taxa arise from these three regions as well as the Afrotropic region (Fig. 1). Further, the core Carex dispersals were inferred from the Himalayas into these four regions and particularly Australasia. However, our sampling was somehow biased towards the Nearctic region, this was a direct consequence of using the Global Carex Group (2016) dataset, which is the largest and most complete Carex dataset to date (>50% of the total diversity).

Diversification dynamics of Himalayan species

Among the three GeoSSE models evaluated with AIC values, the best-fit model (Table 4) for the narrow (ΔAIC = 49.94) and broad (ΔAIC = 76.08) sense biogeographic coding was the constrained model in which speciation in taxa distributed across two regions was set to zero (sAB ∼ 0). Both biogeographic coding (narrow and broad) strongly reject the model constrained to have speciation rates (sA-sB) equal in the Himalayas and non-Himalayas and extinction rates also equivalent (xA-xB). Given our finding that the broad-sense biogeographic reconstruction recovers Himalayan lineages (see biogeography results above), we considered the broad-sense biogeographic coding to be better suited to estimate diversification rates within versus outside of the Himalayas (Table 4). Although the Himalayas showed positive net diversification rate (0.16 events per million years) there was no real difference in diversification rates between the two areas (0.06 events per million years). We utilized Bayesian MCMC on the broad sense dataset to estimate parameter uncertainty (Table 4, Fig. S3). Figure 2 represents posterior probability differences between the lineages of the two regions and diversification rates. It suggests the rate of speciation is indeed lower in the Himalayas than non-Himalayas with difference (0.34 events per million years). A similar finding was obtained with the rate of extinction (0.28 events per million years). The rate of dispersal was higher in Himalayan than non-Himalayan lineages with difference events per million years (0.27 events per million years); while not significant, this result suggests that independent origins of Himalayan taxa from non-Himalayan ancestors may be less common than colonization of non-Himalayas from Himalayan radiations (Fig. S3).

Climatic niche evolution

Plotting stress against number of dimensions in a set of initial NMDS ordinations that include data for both the tip states and ancestral reconstructions (Fig. S4) shows significant decreases in stress up to K = 5. However, as visual inspection of the 5-dimensional ordination showed no qualitative differences from the 2-dimensional ordinations, ordination results here were presented for the K = 2 NMDS analysis. For ordination of both the inferred ancestral states (Fig. 3A) and the tip states (Fig. 3B), eight BIOCLIM variables correlate strongly (—r—>0.7) with MDS axis 1 (BIO1, BIO4, BIO6, BIO9, and BIO11, which relate to temperature; and BIO12, BIO13, and BIO16, which relate to precipitation) and three with MDS axis 2 (BIO14 and BIO17, precipitation of the driest month and quarter respectively; and BIO15, precipitation seasonality, estimated as the coefficient of variation). In the ordination of inferred ancestral states, the two Kobresia clade ancestors were not significantly differentiated from the ancestors of the remaining Himalayan taxa. While the maximum likelihood estimator underlying these reconstructions was biased against detecting differentiation among tips closer to the base of the tree (because the estimator under a Brownian motion process is the weighted mean of all tips in the tree; (Felsenstein, 1985), the ordination of tip states (Fig. 3B) was not biased in this way, and should show differentiation if the Kobresia clades are descended from an ancestor that was uniquely adapted to a different climatic niche. To the contrary, the 95% C.I. for the Kobresia clade climatic niche was almost entirely contained within the 95% C.I. for the remaining Himalayan taxa and exhibited a narrower variance, which was expected given that it reflected a narrower phylogenetic diversity (Fig. 3B). Both groups were also not significantly differentiated from the remainder of Carex.

Analysis using the multivariate scalar O-U model, on the other hand, recovered two major shifts in selective regime: one at the base of the Uncinia clade, which is strongly represented in the Neotropics and New Zealand; and one at the base of a clade comprising primarily of the traditional Indicae and Decorae sections, groups of primaril tropical distribution (Fig. 4) for the given data set employed in this study. These clades are characterized respectively by higher temperatures, higher precipitation, and lower seasonality (Uncinia); and higher temperatures, lower precipitation, lower temperature seasonality, and higher precipitation seasonality (Indicae/Decorae clade). Himalayan taxa overall are widespread in climatic niche and exhibit no consistent trend (Fig. S5).

Himalayan orogeny and the origin of the Kobresia clade

Carex taxa of the Himalayas exhibited high phylogenetic diversity both among clades—three of the five major Carex clades were represented in the Himalayas—and within clades, where the Himalayan Carex are mostly highly polyphyletic. The divergence time for the largest Himalayan Carex clade, the Kobresia clade, falls in the early Miocene (ca. 20.6 Ma). This postdates the initial uplift of the Himalayas by about 30–39 million years (see introduction above), during which time the Indian Tibetan continent collided with Asian plates. It corresponds well, however, with the second stage of Himalayan uplift (25–20 Ma), at which time accretion of the upper layer of the Indian continental crust further raised the Himalayas (Molnar & Stock, 2009,Van Hinsbergen et al., 2011; Van Hinsbergen et al., 2012). This combination of uplift events introduced the summer monsoon precipitation regime and alpine climates at altitudes of 5,000–6,500 m (Xie et al., 2014), which might have contributed to migration of floristic elements into the region. Our interpretation is that the crown diversification of the core Unispicate clade that includes the Kobresia clade was roughly simultaneous with the second Himalayan uplift event, and that the dominance of the Kobresia clade among Himalayan Carex is likely due in part to historical contingency: Kobresia simply arrived first to the high alpine Himalaya, excluding potential competitors rather than exhibiting particularly strong preadaptations.

The Kobresia radiation is singular in its diversity: the remaining Himalayan Carex species arose from 43 migrations into the Himalayas from two other major clades, in which all Himalayan species in Vignea clade lacked Himalayan sisters except two clades of 2 species each; while core Carex Himalayan species also largely lacked Himalayan sisters (25 species), except 10 clades of 2–11 species. Based on age alone, we would expect Kobresia to exhibit higher species diversity than the other Himalayan radiations (Wiens et al., 2009). However, our analyses of multiple-regime Ornstein–Uhlenbeck models suggested that there are not special attributes of the Kobresia clade that make it successful in the Himalayas, but rather that its early arrival to the unique Himalayan climates enabled its success. Further, Kobresia exhibited no innate attributes that allowed it to take over the Himalayas. Therefore, our analyses suggested that possibly the contingency (resulting possibly due to historical abiotic events (Fukami, 2015), e.g., uplift of the region during orogeny of the Himalayas) in Kobresia assembly with the combination of clade age (priority effect) and competitive exclusion (Abrams, 1983; Brown & Wilson, 1956; Grant & Grant, 2006; Monroe & Bokma, 2017; Schluter, 2000) likely explain the success of the Kobresia clade and the floristic dominance of its species in the Himalayas when compared to other Himalayan Carex groups.

Long-distance migration and multiple ancestral ranges

Previous biogeographic studies in Carex (e.g., Escudero et al., 2009; Hoffmann, Gebauer & Rozycki, 2017; Miguez et al., 2017; Spalink et al., 2016b; Villaverde et al., 2015b) have focused on specific clades or broader taxonomic scales. Our study utilizes the broadest phylogenetic sampling of the genus Carex to date (966 species, ca. 50% of the extant diversity; Fig. 1) to address the origins of a highly polyphyletic regional flora, the Carex of the Himalayas. The earliest origin of the Himalayan Carex flora was the Kobresia clade, which arose in the early Miocene from a Nearctic ancestor (Fig. 1), presumably by bird-mediated long-distance dispersal (cf. Villaverde et al., 2017b). Such long-distance migrations are well-documented in the genus under diverse scenarios (Escudero et al., 2009; Jiménez-Mejías, Martín-Bravo & Luceño, 2012; Miguez et al., 2017; Villaverde et al., 2015a; Villaverde et al., 2015b; Villaverde et al., 2017a; Villaverde et al., 2017b). There were few radiations out of the Himalayas in this clade to the adjacent Eastern Palearctic and Indo-Malaya as well as to more distant regions (Western Palearctic, Neotropic, and Nearctic). We detected no radiation into the Himalayas from the most closely adjacent region, Indo-Malaya (Fig. 1). Very likely this is due to the fact that Indo-Malaya is dominated by tropical and subtropical dry and moist broadleaf forests, and Carex is a predominantely temperate group (Escudero et al., 2012; Spalink et al., 2018). All these findings were based on the subset of taxa utilized in this study; inclusion of additional taxa might of course influence these results, but we have no reason to suspect that they are biased.

In the core Carex clade, by far the largest Carex clade globally (Global Carex Group, 2015; Global Carex Group, 2016), Himalayan species cluster into small clades (2–11 species) and 25 individual species that arose from outside the Himalayas. Ancestral origins for these clades and species are remarkably disparate (File S10). Thirty Himalayan species derive from the traditional sections Graciles, Decorae, Indicae, Setigerae, Aulocystis, Thuringiaca, Polystachyae, Clandestinae and Radicales derived primarily from Eastern Palearctic ancestors. Twelve Himalayan species from sections Racemosae, Aulocystis, Vesicariae and Phacocystis derived from the Nearctic, while a smaller number of Himalayan species (4) from the traditional Thuringiacae, Spirostachyae, Ceratocystis, and Hallerianae derived from the Western Palearctic. The Afrotropics seemingly contributed only C. obscuriceps to the Himalayas (Fig. S3), although this inference may be due to undersampling of the large traditional section Vesicariae. It is similarly striking how many radiations are inferred to have given rise to geographically distant lineages. Dispersals from the Himalayas to the Afrotropic and Australasian region in core Unispicate and core Carex clades bear further investigation.