Detecting the phylogenetic signal of glacial refugia in a bryodiversity hotspot outside the tropics

Glacial refugia have likely been important in shaping diversity gradients outside the tropics. Many taxa that have high extratropical diversity in the present day, such as mosses, may have persisted in glacial refugia. However, the biogeographical histories of most species within refugia remain unclear.


| INTRODUC TI ON
The diversity of life is unevenly distributed across the globe.
Perhaps the most well-recognized biodiversity pattern is the latitudinal gradient in species richness, which describes the tendency for richness to be highest in equatorial regions and taper towards the poles (see Hawkins et al., 2003;Hewitt, 2000). However, while uncommon, some taxa show countergradients, with highest richness at high latitudes (Kindlmann et al., 2007); for example lagomorphs, arboreal ants, molluscs and aphids are all more diverse outside of the tropics (Morales-Castilla et al., 2020). Understanding these 'exceptions to the rule' can provide unique insights into the mechanisms shaping biodiversity gradients. Vascular plants demonstrate a typical latitudinal gradient in species richness, with highest diversity in the warm and moist tropics (Brown, 2014). The biogeography of non-vascular plants has been less well explored; while there is some evidence that liverworts and hornworts also exhibit a general latitudinal richness gradient (Wang et al., 2017), mosses may exhibit inverse or no latitudinal diversity gradient (Mateo et al., 2016;Shaw et al., 2005). Additionally, several hotspots of moss species richness are found outside of the tropics, including in British Columbia, the northern Andes, Japan, Madagascar, the East African Highlands, the Himalayan region, central Europe and Scandinavia (Geffert et al., 2013). Here, we use a regional phylogeny of extant mosses from the bryophyte diversity hotspot and putative glacial refugia of Haida Gwaii (formerly known as the Queen Charlotte Islands), an island archipelago off the north-west coast of British Columbia, Canada, to explore the historical mechanisms that have shaped this exceptionally rich northern latitude bryophyte flora.
The moss diversity of Haida Gwaii is remarkable, with over 380 species within an area of 10,180 km 2 (Golumbia & Bartier, 2004), including the endemic Carey's small limestone moss (Seligeria careyana; Vitt & Schofield, 1976), which has only ever been observed from three localities. While the majority of moss species (195) (Schofield, 1989). It is possible that rare single long-distance dispersal events contribute to the unique diversity of disjunct species in glacial refugia (Heinrichs et al., 2009;Provan & Bennett, 2008).
However, some bryophytes with disjunct distributions appear to lack mechanisms for effective long-range dispersal, suggesting vicariance might also be important (Shaw, 2001). The occurrence of similarly disjunct non-bryophyte plants indicates a possible common process: Eurasian origin, dispersal to North America and subsequent western North America disjunction (Xiang et al., 1998). The presence of these disjunct and endemic species on Haida Gwaii and their absence elsewhere in north-western North America, which resembles a more boreal flora compared with the more temperate flora on Haida Gwaii (Alaback, 1996), suggest a unique spatial and temporal history of the Haida Gwaii ecosystem that has sustained a rich bryoflora to the present (Heusser, 1989).

One explanation for the unique bryophyte diversity of Haida
Gwaii is linked to suggestions that the archipelago was a glacial refugium during the Late Wisconsin Glaciation (Mathewes & Clague, 2017). Graham island, the largest island in the Haida Gwaii archipelago, was glaciated at least twice, once >52,000 years ago, then again between 27,500 and 16,000 years ago (Demboski et al., 1999;Mathewes et al., 2015). However, Haida Gwaii may have been partially ice-free as early as 16,000 years ago, when the lowland phase of the Late Wisconsin Glaciation was at its maximum, as indicated from the succession of Cyperaceae captured in pond sediments from 16,830 years BP (Lacourse et al., 2005). Possible locations of refugia include the Queen Charlotte Ranges, the west coast of Graham and Moresby islands, and the shelf beneath Hecate Straight, which was at various times exposed with freshwater lakes (Shafer et al., 2010).
Present-day plant communities located between 900 and 1100 m in elevation potentially represent remnants of the refugial flora (Heusser, 1989), as the climate at these elevations may resemble the climatic conditions of historical refugia.
To date, the biogeographical histories of most species within glacial refugia remain unclear. For a few species, population studies and phylogeographic methods have allowed us to identify continental refugia and patterns of post-glacial range expansion (Allen et al., 2015;Provan & Bennett, 2008), for example in oaks, common beech, black alder and silver fir in Europe (Hewitt, 1999). Commonalities across relictual species can also be found, for example, as reflected in their present-day genetic diversity and geographic distributions-species that survived within several large refugia have greater diversity than those that were restricted to fewer small refugia (Roberts & Hamann, 2015), which more often exhibit disjunct present-day distributions.
Glaciation would have likely had a disproportionate impact on clades with temperate origins and high extratropical diversity. Refugia at higher latitudes may have been important in preserving relictual lineages, more so than for clades following more commonly observed latitudinal biodiversity gradients with greater species richness at lower latitudes.
Here, we reconstruct the regional phylogeny of the Haida Gwaii bryoflora-focusing on mosses-using rbcL and trnL-F DNA barcode markers extracted from field samples and herbarium specimens. We use this phylogeny to explore the macroevolutionary and macroecological imprint of glacial history on the moss flora of Haida Gwaii by examining the relationship between evolutionary distinctiveness, dispersal traits (i.e. spore size and sexuality) and range continuity in their present-day distributions. We test whether species with more disjunct distributions are also phylogenetically isolated, and whether habitats associated with glacial refugia, such as those at higher elevations on Haida Gwaii, support more phylogenetically dispersed bryophyte communities. We suggest that clades exposed to glacial expansion may have experienced a wave of extinctions, making surviving species, such as those in refugia, more phylogenetically distinct. If glaciation also led to a significant contraction of species' distributions, with remnant populations restricted to glacial refugia geographically distant from centres of their evolutionary origin, these species would tend to have few close relatives with which they co-occur and more disjunct present-day distributions.

| Plant material and molecular sequencing
We used the two molecular barcode markers, rbcL and trnL-F, to place the mosses of Haida Gwaii on the robust, multi-gene backbone phylogenetic reconstruction of major moss lineages proposed by Liu et al. (2019). Using the species list (with corrections to conform to current nomenclature conventions) provided in Golumbia and Bartier (2004)-a comprehensive baseline inventory report for the bryophytes of Haida Gwaii-we downloaded the available rcbL and trnL-F sequences from the GenBank database (https://www.ncbi.nlm.nih.gov/genba nk/) for each species. Species that did not have either existing rbcL or trnL-F sequences were targeted for field collection and sequencing.
During July and August 2019, we collected voucher specimens representing 109 species from Haida Gwaii using standard field collection methods, as outlined in Schofield (1992). The 21 collection sites (Figure 1), consisting of two alpine locations, four hyper-oceanic locations, eleven low-to mid-elevation forest locations (four along creeks) and four low-elevation open wetlands, were selected to encompass as many as possible of the different habitat types described by Golumbia and Bartier (2004). All microhabitat and substrate types within each site were targeted for collection. Location description, latitude and longitude, microhabitat and other environmental characteristics were recorded for each voucher. Each specimen was photographed in situ, and accessioned at the University of British Columbia Herbarium (UBC). Taxonomic identification and nomenclature of the specimens conform to volumes 27 and 28 of the Flora of North America (Flora of North America Bryophyte Editorial Committee, 2007Committee, & 2014. In total, 99 field samples and an additional 86 species previously vouchered at the UBC Herbarium were sampled and sequenced for both molecular barcode regions. Samples were ground at 60 Hz for 60 s using TissueLyser 24 (Shanghai Jinxin Industrial Development Co. Ltd.), and total genomic DNA for each sample was extracted using the HP Plant DNA Kit (Omega Bio-Tek) following the manufacturer's protocol. A nested PCR system was optimized based on the recommended protocol from Hentschel et al. (2006) and, using the primers in Feldberg et al. (2016): rbcL1F and rbcL1390R for the first round, and rbcL210F and rb-cL1200R for the second round of PCR for rbcL; and trnL-F and trnL-R for the one round of PCR for trnL-F. The amplification products were F I G U R E 1 Location and topographical map of Haida Gwaii, BC. Red symbols indicate sampling locations sequenced at the Beijing Genomic Institute (BGI). The generated sequences for each species were then validated against the NCBI database using the BLAST algorithm. For some species, it was only possible to sequence one of the two barcode genes.

| Phylogenetic reconstruction
We reconstructed a dated molecular phylogeny of Haida Gwaii mosses (n = 300 species) using a backbone constraint (described below), and the combined dataset of the newly generated rbcL and trnL-F sequences and sequences downloaded from GenBank. Sequences were aligned using Geneious Prime (v. 2020.0.3; Biomatters Ltd.; Kearse et al., 2012). For new barcode sequences, ambiguous ends were trimmed, and a consensus sequence was generated using forward and reverse reads. Then, rbcL and trnL-F sequences were aligned separately using MUSCLE v3.8.13 (Edgar, 2004) and optimized using Gblocks Server (Castresana, 2002) to eliminate poorly aligned regions. The aligned genes were then concatenated into a single matrix.
We reconstructed the phylogeny using RAxML v8.2.12 (Stamatakis, 2014), with the GTR + GAMMA + I model for both rbcL and trnL-F genes, as selected by jModelTest v2.1.10 (Darriba et al., 2012). Phylogenetic information within the two barcode genes is limited, and while useful for differentiating among species, they may be less informative at resolving deeper phylogenetic relationships.
We therefore enforced topological constraints at the family level, assigning genera to families using the topology reflecting the ordinal relationships of mosses proposed by Liu et al. (2019). For the genus Sphagnum, we further enforced subgenera-level constraints using the topology in Shaw et al. (2019). Bootstrap analyses were implemented by GTR + CAT approximation for 100 replicates. We rooted the resulting topology with Takakia lepidozioides as the outgroup.
Branch lengths were made proportional to time using the penalized likelihood method, as implemented in r8s v1.8.1 (Sanderson, 2003).

| Species life history traits, evolutionary distinctiveness and distribution data
We used the habitat and distributional information from Golumbia and Bartier (2004), as well as life history strategies from the bryological literature (Söderström & During, 2005), to generate a matrix of species distribution and habitat preferences, species sexuality, asexual reproduction and spore size (minimum, maximum and mean) (see Table   S1). We also calculated evolutionary distinctiveness (ED), a specieslevel estimate of phylogenetic isolation on a time-calibrated phylogeny (Jetz et al., 2014), using the evol.distinct function in the package picante (Kembel et al., 2010) and the 'equal splits' method (Redding & Moores, 2006) in R version 1.1.463 (R Development Core Team, 2018). When calculated at the regional level, this metric provides an indication of whether a species has close relatives on the Haida Gwaii archipelago (high ED species have few close relatives). We recorded spore size as a continuous variable, and species sexuality (0 = monoicous, 1 = dioicous) and asexual reproduction (0 = no, 1 = yes) as discrete variables using data from volumes 27 and 28 of the Flora of North America (Flora of North America Bryophyte Editorial Committee, 2007Committee, & 2014 augmented with information from the California moss eflora (2020).
We classified species geographic distributions using information from appendix A of Golumbia and Bartier (2004), which presents a continuum of 10 categories (indicated in parentheses, below). First, we broadly classified species from widespread to local on a scale of 1 to 6: 1 (widespread in Northern Hemisphere), 2 (Amphi-Pacific), 3 (western European-western North American), 4 (western Europeanwestern North American showing a widely interrupted pattern in the Northern Hemisphere and Eastern North American disjuncts), 5 (tropical-subtropical disjuncts and species for which Haida Gwaii contains the only known Canadian population) and 6 (western North America and Pacific coast). Second, we generated a discrete classification scheme distinguishing species with continuous distributions (grouping distribution codes 1 and 2), disjunct distributions (grouping distribution codes 3, 4 and 5) and localized distributions (distribution code 6). Third, we used a binary classification to distinguish between species with continuous (distribution codes 1, 2 and 6) vs. disjunct distributions (distribution codes 3, 4 and 5). Where possible, species missing distribution data from Golumbia and Bartier (2004) were classified based on the available distribution information in volumes 27 and 28 of the Flora of North America (Flora of North America Bryophyte Editorial Committee, 2007Committee, & 2014. Taxa on our reconstructed phylogeny without distribution data (n = 13) were excluded from subsequent analyses. Finally, we coded the presence of species within each of the 16 habitat and vegetation types on Haida Gwaii (hereon referred to as habitat types), as described in Golumbia and Bartier (2004). The assignment of habitat associations was derived from the habitat codes in appendix A of Golumbia and Bartier (2004), which were determined through extensive habitat surveys conducted on Haida Gwaii, and broadly describes the species habitat preferences summarized in Schofield and Hong (2002).

| Analysis of species distribution and evolutionary history
We explored whether the composition of mosses across the 16 habitat types listed in Golumbia and Bartier (2004) was characterized by species with disjunct distributions and high ED, as this might suggest a common cause for both species' attributes linked to habitat preference. We first summarized the phylogenetic diversity (PD), the sum of lengths of all branches on the phylogeny that span the members within the habitat type (Faith, 1992), and species richness (SR) for each habitat type. Next, we calculated the standard effect size of phylogenetic diversity (ses.PD) for bryophyte assemblages associated with each of the habitats, using the 'tip-swap' algorithm and 1000 randomizations. Ses.PD quantifies whether observed PD differs from a null expectation derived from sampling species at random from the species pool. Last, we calculated the mean evolutionary distinctiveness (mED) of the different habitat types by taking the mean of the species ED as estimated on the reconstructed regional phylogeny.
We then used three different approaches (linear models, generalized linear models and analysis of variance) to test the relationship between species distribution, life history traits and evolutionary distinctiveness. The phylogenetic non-independence of species violates assumptions of classical statistical methods (Felsenstein, 1985;Harvey et al., 1995); therefore, we fit our models using phylogenetic generalized least squares (PGLS), which is equivalent to GLS with an error structure specified by the phylogeny, phylogenetic ANOVA (phylANOVA) and phylogenetic generalized linear models (phyloglm).
First, we used a phylogenetic generalized least squares regression to model trait correlations in the R package caper (Orme et al., 2012).
For purposes of model fitting, zero-length branches were assigned a nominal unit length. The phylogenetic variance-covariance matrix was included to model phylogenetic non-independence in the data, fitting the maximum-likelihood value of λ. Species distribution, coded on a continuous scale from 1 to 6 (see above), was the response, and separate models were fit with species sexuality, spore size, asexual reproduction and evolutionary distinctiveness as independent predictors, and then with species sexuality, spore size and evolutionary distinctiveness together as grouped predictors. We checked model fits using standard diagnostics using the plot() function in Caper.
Second, to further explore the relationship between ED and distribution, we compared differences in ED between the three distribution classes: continuous, disjunct and regional (see above) with a phylogenetic ANOVA using the phylANOVA function from the R package phytools (Revell, 2012). Because the ANOVA could not be fitted to an incompletely resolved phylogenetic tree, we removed polytomies by randomly dropping all except one taxon per terminal clade using the function thin_terminal_polytomies (Davies et al., 2012), and fit the model on this reduced dataset. p-Values for the post hoc tests were adjusted using Holm's correction for multiple comparisons. We then repeated this thinning procedure 1000 times to generate a distribution of test results.
Third, we compared ED between species with continuous vs. disjunct distributions using phylogenetic logistic regression fit on the full tree using the function phyloglm from the R package phylolm (Ho & Anné, 2014), and with 100 independent bootstrap replications.
Here, species distribution was modeled as a binary response with ED as the explanatory variable.

| Phylogenetic reconstruction of Haida Gwaii moss species
The regional phylogeny of Haida Gwaii mosses (Figure 2) was constructed using 264 trnL-F sequences (145 newly generated sequences and 119 from GenBank; Table S2) and 324 rbcL sequences (127 newly generated sequences and 195 from GenBank; Table S2), and includes 300 out of the approximately 380 species estimated present on the islands. Table S3 includes a complete moss species list for Haida Gwaii with up-to-date taxonomic nomenclature, and indicates the species we were not able to obtain sequences with an asterisk. While we made efforts to collect samples in the field to maximize species coverage, some rare species lack genetic data and thus could not be included in our analyses. However, ED metrics are surprisingly robust to missing taxa (Weedop et al., 2019), and we do not believe that missing taxa would bias the observed relationship between evolutionary distinctiveness and range disjunction.
Due to the limited phylogenetic information within the two genes, major clades were topologically constrained to agree with the backbone phylogeny proposed by Liu et al. (2019), while the sequence data from the barcode genes were used to resolve more fine-scale relationships.
This approach is well established (Joly et al., 2014) and has been shown to help facilitate robust phylogenetic reconstructions even in highly diverse systems, such as tropical forests (see Kress et al., 2010). Following time calibration, some short branches were collapsed to polytomies, resulting in a few unresolved clades within the Hypnales, including in Campyliaceae, Calliergonaceae and Leskeaceae. However, these collapsed nodes are most likely associated with relatively rapid divergence events (where there has been limited time for mutations to accumulate) and thus should not have a large influence on species ED estimates.  (Table 1). Epiphytes, seaside outcrops, disturbed soil and highelevation habitats have the highest proportion of species with disjunct distributions, while forest floor, epiphytes, aquatic and high-elevation habitats have the highest proportion of species with localized western North American distributions ( Table 1).

| Phylogenetic attributes and species range structure across habitats
Phylogenetic structure and species diversity of bryophytes also varied across the different habitat types. Broad evidence for phylogenetic clustering (as indexed by ses.PD) indicates that most habitats capture less phylogenetic diversity (PD) than predicted from their species richness (Table 1). Epiphytes have the lowest ses.PD, indicating that epiphytic taxa are the most phylogenetically clustered. The mean evolutionary distinctiveness of species (mED) is also lower for epiphytes than for other habitat types ( Table 1). In contrast, high-elevation habitats capture somewhat more PD than expected from their species richness (positive ses.PD), and mean species evolutionary distinctiveness is higher than in other habitats (mED = 113.68). Overall, there was a general positive correlation between ses.PD and mED (Pearson's r = .66), while the relationship between ses.PD and total PD or SR was somewhat weaker and negative (Pearson's r = −.11 and −.38, respectively). Thus, while habitats with a less phylogenetically constrained moss flora tended to be less species-rich, the species within them were more evolutionarily distinct.
However, we did not observe any strong relationship between ses.PD and the proportion of species with disjunct distributions within habitats.

| Relationship between species distributions and evolutionary distinctiveness
In our phylogenetic generalized least squares regressions, we found no significant relationship between either dispersal traits (i.e. spore size and sexuality) or evolutionary distinctiveness and species distribution when modeled as a continuous variable assuming six distribution classes. Although there appeared to be some weak evidence for a relationship between sexuality and distribution (slope = 0.66, p = .04) in the multiple regression model, the overall model was not significant (p > .05, model R 2 = .02; Table 2), and this relationship was not supported in the single-trait model (all p > .05, model R 2 = .01; Table 2).
Reclassifying distributions into three discrete classes, localized, disjunct and continuous, revealed strong and significant differences in the evolutionary distinctiveness of species with different distributions (p < .05 from the phylogenetic ANOVA; Table S4). The corrected post hoc pairwise comparisons using the Holm adjustment indicates that species with disjunct distributions are more evolutionarily distinct than species with continuous or localized distributions, but that there was no significant difference between species with continuous and regional distributions ( Figure 3).
We further examined the difference in evolutionary distinctiveness between species with disjunct distribution and those with continuous or regional distributions (Figure 4) Figure S1).

| DISCUSS ION
One of the most universal global diversity patterns is the latitudinal gradient in species richness. However, unlike most vascular plant taxa, mosses exhibit no general latitudinal diversity gradient, and temperate regions are equally as diverse as the tropics (Geffert et al., 2013;Möls et al., 2013). The temperate Northern Hemisphere hotspots of moss diversity, such as in British Columbia, which also contain a disproportionate number of rare species, are particularly noteworthy. Explanations for these unusual diversity gradients frequently assume that extratropical species richness reflects phylogenetic niche conservatism (Stephens & Wiens, 2003;Wiens & Donoghue, 2004) in clades with temperate or high latitude origins.
Mosses consist of a diversified group of plants and likely had a temperate origin (Shaw et al., 2005), with different lineages evolving various strategies to persist across diverse environmental conditions, for example, related to water availability (Chen et al., 2015). TA B L E 1 Comparison between phylogenetic diversity (PD), species richness (SR), standard effect size of phylogenetic diversity (ses.PD), mean evolutionary distinctiveness (mED), and the proportion of species with continuous, disjunct and localized distributions to the total number of species across the 16 habitat types Note: Siliceous rocks and water body margins have the highest overall species richness, while siliceous rocks and calcareous rocks have the greatest frequency of species with continuous distributions, and calcareous rocks and epiphytes have the greatest frequency of both disjunct and localized species. High-elevation sites are the only habitat that have higher PD than predicted by their species richness (positive ses.PD), and the largest mean evolutionary distinctiveness (mED). Epiphytes have the lowest ses.PD, and second lowest mED, and swamp forest has the lowest mED.

TA B L E 2
Phylogenetic generalized least squares regressions assuming a continuous classification of range size on a scale of 1-6 (where 1 = Northern Hemisphere, 2 = Amphi-Pacific, 3 = western European-western North American, 4 = Northern Hemisphere disjuncts, 5 = tropical-subtropical disjuncts, and 6 = localized in western North America) as response variables, and with spore size, evolutionary distinctiveness and sexuality as predictors  1958), and other basal clades, such as Sphagnum, are more phylogenetically diverse in the temperate Southern Hemisphere (Shaw et al., 2005). Extratropical niche conservatism and ancestral adaptation to temperate environmental conditions may therefore help explain the high latitude species richness of mosses. Shaw et al. (2005) suggest the lack of a strong latitudinal gradient in molecular diversity highlights the 'migratory prowess of mosses', and mosses may be more constrained in their distributions by fit to the environment (filtered into temperate climates matching to their evolutionary origins) than dispersal limitation per se. Nonetheless, glacial history likely provides an additional part of the explanation for present-day hotspots of temperate moss diversity.
During glacial cycles, the distributions of many high latitude and temperate species contracted towards the equator, reducing species richness and phylogenetic diversity at higher northern latitudes.
Following glacial retreat, northward recolonization may have been slow for some taxa (Svenning & Skov, 2007), and undoubtedly restructured many species assemblages. Glacial refugia provided longterm stable habitats (Sundaram et al., 2019) that allowed some areas to retain high diversity away from the tropics, and to more rapidly recolonize higher latitudes following glacial retreat (Keppel et al., 2012). Several Northern Hemisphere hotspots of present-day moss diversity might represent such historical refugia. Indeed, historical climatic conditions throughout the different glacial-interglacial cycles may have imposed consistent constraints on the functional diversity of plants, reducing regional bryophyte diversity, while more stable habitats, such as glacial refugia, were able to support and accumulate greater evolutionary and functional diversity through time (Dynesius & Jansson, 2000;Ordonez & Svenning, 2017).

| Phylogenetic structure of glacial refugia
Phylogenetic niche conservatism predicts that newly emerging environments, such as those revealed from receding ice sheets, will select for species with traits preadapted to these novel habitats (Donoghue, 2008). Such environmental filtering of species by postglacial habitats will therefore tend to result in phylogenetic clustering (Baeten et al., 2015;Ma et al., 2016). However, glacial refugia would not have been subject to similar habitat filtering, allowing them to retain a more phylogenetically dispersed species assemblage. Nonglaciated habitats, such as nunataks, were likely tundra, resembling present-day high-elevation habitats, which is consistent with our observation of greater phylogenetic dispersion in mosses from highelevation sites. The phylogenetic overdispersion in high-elevation habitats contrasts starkly with the phylogenetic clustering observed across other habitat types, consistent with habitat filtering.
Phylogenetic overdispersion has often been associated with the process of competitive exclusion at small spatial scales, that is within ecological communities, and phylogenetic under-dispersion is F I G U R E 3 Boxplot of mean evolutionary distinctiveness for Haida Gwaii mosses with (i) continuous distribution throughout N. Hemisphere, (ii) N. Hemisphere disjuncts and (iii) localized in W. N. America distributions. Colour and letters above the boxes represent statistical significance (p < .05) from the phylogenetic ANOVA. Species with disjunct distribution have a significantly higher average evolutionary distinctiveness compared to both species with continuous and localized distributions (see Table S2) F I G U R E 4 Visual comparison of evolutionary distinctiveness and distribution mapped on the reconstructed Haida Gwaii moss phylogeny. Evolutionary distinctiveness, based on the dated phylogeny in millions of years, is shown on the left, mapped as a continuous variable under assumptions of Brownian motion using the contMap function from the phytools R library (Revell, 2012). Blue indicates lower evolutionary distinctiveness, and red indicates higher evolutionary distinctiveness. The figure on the right shows species distributions classified as a binary trait (continuous vs. disjunct) and mapped onto the phylogeny using stochastic character mapping using the make.simmap and densityMap functions from the phytools R library (Revell, 2012). Species at the tip of blue branches have continuous distributions, while species at the tips of red branches have disjunct distributions associated with environmental filtering of the regional species pool (Webb et al., 2002). However, at regional scales, spatial patterns of phylogenetic dispersion are more likely to reflect historical biogeography (Davies & Buckley, 2011). We suggest that present-day highelevation habitats may provide a modern refuge to glacial relicts, disrupting the expected patterns of phylogenetic clustering commonly predicted in strongly filtered communities at high elevations, and contribute to phylogenetic dispersion in glacial refugia (Brooks & Bandoni, 1988;Fryxell, 1962;Shooner et al., 2018). Other species characteristics of the mosses on Haida Gwaii and within these high-elevation habitats, such as the high proportion of evolutionary distinct and disjunct species, further support the glacial refugia hypothesis.

| Glacial history and range disjunction
The because we a priori hypothesized that they might relate to dispersal ability, showed little relationship with species geographic extents.
While propagule size has previously been associated with dispersal ability (Söderström & During, 2005), correlations between spore size and distribution are weak and non-significant in our data. It is possible that, in mosses, neither small nor large spores have a dispersal advantage. Small spore species can produce higher quantities of spores, but large spore species establish more easily within a locality and have better survival (During, 1997). While other traits might explain some of the variation in geographic extent, trait data for mosses are limited, and we suggest that these disjunct distributions more likely reflect the signal of glacial history, a consequence of the process of range expansion and contraction between glacial and interglacial periods (Zhao et al., 2019).
If high-elevation sites best match historical refugia, we might then also predict that they should have a higher proportion of disjunct species. Given differences in the relative abundance of different habitat types across the landscape, such as the wide-

| Glacial refugia support more evolutionarily distinct lineages
We find that the mean evolutionary distinctiveness of the constituent species in high-elevation habitats is much greater than found in any other habitat type, suggesting that species within high-elevation habitats tend to have few close relatives. Although it is also possible that high-elevation habitats encompass a larger number of vegetation types than other habitats, which could drive greater phylogenetic dispersion (there may be large differences in phylogenetic membership between vegetation types), we do not suspect this is the primary factor driving the observed phylogenetic patterns because such spatial structuring cannot explain why species found in high-elevation habitats are more evolutionarily distinct, and evolutionary distinctiveness was quantified at the level of the entire Haida Gwaii archipelago. Here, we suggest that more evolutionarily distinct species with few close relatives are found in refugia (i.e. highelevation habitats) either because their close relatives may have been lost during previous cycles of glacial expansion, and failed to recolonize following glacial retreat, or because they represent isolated populations of previously widespread species that diversified elsewhere ( Figure 5 illustrates one possible scenario).

| Evolutionary distinctiveness as a predictor of range disjunction
We also show that evolutionary distinctiveness is a strong predictor of range disjunction, supporting our hypotheses that glacial history shapes both the phylogenetic structure and range distribution of species. Although our results were significant across models and alternative classifications of range distributions, many additional factors likely influence species' evolutionary distinctiveness, most obviously the underlying diversification processes that have shaped moss phylogeny. Including such additional information would increase model explanatory power and could allow us to more clearly demonstrate the link with range disjunction. Unfortunately, a complete global phylogeny at the species level for mosses would be required, and we are still some way from achieving this goal.

| Bryodiversity hotspots vs. global plant diversity hotspots
While the signature of phylogenetic overdispersion in the mosses of high-elevation habitats on Haida Gwaii may reflect the signature F I G U R E 5 Schematic illustrating one possible shift in regional phylogenetic relationship through glacial cycles. (a) shows the distribution of species ⍺, ⍬, Ω and µ in the pre-glacial landscape. During glaciation, lowland areas become covered by ice (b). Temperatures during this period are also much cooler, and high-elevation areas and north-facing slopes may also be ice-covered or too cold to support plant life. As a result, species µ is lost from the regional pool as its climate niche contracts. Fortunate survivors, species ⍬ and Ω, persisted in warmer mid-elevation refugia habitats that resemble nunataks, while species ⍺ persists in incised valley refugia. In the post-glacial flora (c), species ⍬ and Ω expand their ranges upwards in elevation as the temperature warms; meanwhile, warm-adapted species capable of rapid dispersal and establishment (species ⍺) rapidly recolonize lower elevation areas. In high-elevation habitats, the resulting post-glacial flora (species ⍬ & Ω) is more phylogenetically dispersed and has higher mean evolutionary distinctiveness (mED), relative to the pre-glaciation flora (species Ω & µ). While it is possible that species µ might have persisted in a refugia elsewhere, the disjunct distribution of high-elevation (refugial) habitats will likely restrict its dispersal and ability to recolonize, at least in the short term. Furthermore, the regional mED of high-elevation sites may further increase if new species that are descendants from distant clades colonize the habitat of glacial history in this region, the phylogenetic structure of refugial communities likely differs among taxa and biomes. For example, phylogenetic clustering may be observed in tropical refugia dominated by angiosperms (Costion et al., 2015), which have a broadly tropical evolutionary origin, with many angiosperm lineages diversifying when tropical environments were globally more extensive (Crane & Lidgard, 1989). In contrast, mosses likely diversified extensively outside the tropics (Shaw et al., 2005). Thus, while tropical refugia for flowering plants might be more likely to sample from a regional pool of more recently diversified lineages, temperate refugia for mosses may be more likely to sample older lineages from a more phylogenetically diverse species pool.
In the temperate moss refugia, the persistence of relictual lineages may contribute greater phylogenetic diversity to regional assemblages than the migration of invading lineages from other regions post-glaciation (Harrison & Noss, 2017). Relict species might possess unique traits that allowed them to persist through glacialinterglacial cycles (Tang et al., 2018), or they may simply represent fortunate survivors. The narrow ecological niche preferences among species that persisted within glacial refugia (Frahm, 2007) might have further reinforced the range disjunction of relict species due to the isolated distribution of these habitat types in the present day (Lv et al., 2018). It would be interesting to evaluate whether there is a generalizable relationship between range disjunction and ED across refugia for other temperate clades. We hope our study will inspire further exploration of such relationships.

| Conservation value of refugial habitats
The biodiversity of Haida Gwaii comprises three intertwined features: ecosystem diversity, species diversity and genetic diversity.
Maintaining genetic diversity within species is essential for sustaining healthy populations (Hughes et al., 2008). Similarly, maintaining the phylogenetic diversity of species within ecosystems contributes to sustaining ecosystem processes and integrity (Cadotte et al., 2008). Phylogenetic diversity is also a measure of evolutionary heritage with intrinsic value (Faith, 2015;Mooers et al., 2005), and evolutionary distinctiveness is important when prioritizing habitats for conservation (Costion et al., 2015). On Haida Gwaii, we have shown that high-elevation sites likely represent historical glacial refugia and support relictual species with high evolutionary distinctiveness.
These habitats therefore have high conservation value even though they may not be as species-rich as other habitat types. In contrast, epiphytic habitats have high species richness, but the phylogenetic clustering of species within them captures relatively little evolutionary history. The conservation of refugial habitat, especially highelevation habitats, is increasingly important with current warming trends placing additional stress on cold-adapted species, including many Northern Hemisphere bryophytes (Wu et al., 2018) .
By examining the imprint of range contraction and expansion during glacial cycles from a community-wide species perspective (see also Mastrogianni et al., 2019), our study illustrates how the present-day phylogenetic structure of species composition can reveal the signal of glacial refugia, and help explain why some taxa have higher species richness outside of the tropics. We present an alternative but complementary approach to traditional population genetics and palaeobiological studies, which we hope can contribute to providing a richer picture of the recent biogeographic history of glacial refugia and Northern Hemisphere bryophyte diversity.

ACK N OWLED G EM ENTS
The Gougherty for help with the location map. Finally, special thanks to Steve Joya for his assistance with identification and verification of the sequenced specimens.

CO N FLI C T O F I NTE R E S T
The authors declare that there is no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All newly generated barcode sequence data are available on the