Testing which axes of species differentiation underlie covariance of phylogeographic similarity among montane sedge species

Co‐distributed species may exhibit similar phylogeographic patterns due to shared environmental factors or discordant patterns attributed to the influence of species‐specific traits. Although either concordant or discordant patterns could occur due to chance, stark differences in key traits (e.g., dispersal ability) may readily explain differences between species. Multiple species’ attributes may affect genetic patterns, and it is difficult to isolate the contribution of each. Here we compare the relative importance of two attributes, range size, and niche breadth, in shaping the spatial structure of genetic variation in four sedge species (genus Carex) from the Rocky Mountains. Within two pairs of co‐distributed species, one species exhibits narrow niche breadth, while the other species has broad niche breadth. Furthermore, one pair of co‐distributed species has a large geographical distribution, while the other has a small distribution. The four species represent a natural experiment to tease apart how these attributes (i.e., range size and niche breadth) affect phylogeographic patterns. Investigations of genetic variation and structure revealed that range size, but not niche breadth, is related to spatial genetic covariation across species of montane sedges. Our study highlights how isolating key attributes across multiple species can inform their impact on processes driving intraspecific differentiation.

A primary objective of phylogeography is to identify the processes that shaped the spatial structure of species' genetic variation (Avise 2000). By focusing on multiple taxa within a comparative phylogeographic framework, for example, species that share many life history similarities but differ in specific attributes, it is possible to isolate the individualistic effects of these attributes on the geographic distribution of genetic structure. In some cases, phylogeographic patterns among species may be concordant because environmental influences, such as the cycling of glacial and interglacial periods during the Pleistocene, overwhelm any species-specific differences (Bemmels et al. 2016;Liggins et al. 2016;Gil-López et al. 2017). In other circumstances, species-specific differences are so influential that they overwhelm signatures of shared history and cause species to react to pressures in profoundly different ways (i.e., phylogeographic discordance; Myers et al. 2017, reviewed in Sork et al. 2016). As such, comparative phylogeographic methods can facilitate hypothesis generation in relation to how species' attributes affect genetic structure across members of a community, and/or how traits interact with one another to shape phylogeographic patterns (Papadopoulou and Knowles 2016;Hodel et al. 2018).
Two attributes of species that shape population genetic structure are range size and niche breadth; both theory and empirical evidence support their effect on genetic structure (Gaston 1994), which we define as the degree of differentiation between populations within a species, as measured by metrics such as F ST . The magnitude and direction of the relationship between niche breadth, which we consider the range of environmental variation that determines the fundamental abiotic requirements of a species (sensu Hutchinson 1957), and a species' genetic structure is not universal. Hereafter, we use the term niche breadth to refer to the range of species-level environmental factors, as opposed to population-level or individual-level descriptions of niche, that characterize habitat specialization. Species may have broad niche breadths due to local adaptation to diverse environments in their ranges (e.g., to thermal environments), which may lead to high genetic structure (Campbell-Staton et al. 2016). Alternatively, broad niche breadths may be associated with little genetic structure (e.g., Svanbäck and Schluter 2012) in species that display high phenotypic plasticity (Finch et al. 2019, Parker et al. 2003. High genetic structure can also arise in species with narrow niche breadths, such as when habitat specialists persist across a landscape where environmental heterogeneity is high, and hence suitable habitats are isolated (e.g., Afonso Silva et al. 2017). Narrow niche breadth species may also exhibit low genetic structure among populations, especially if the suitable habitat is rare and contiguous. Finally, niche breadth can have little impact on genetic structure; for example, broad-niche habitat generalists and narrow-niche habitat specialists can both exhibit a similar lack of genetic structure (Reece et al. 2011).
Similar to niche breadth, the relationship between range size and genetic structure is difficult to generalize (Duminil et al. 2007;Pelletier and Carstens 2018). Good dispersers may have large ranges and high connectivity among populations, which can contribute to low genetic structure (Brown et al. 1996). However, in some cases, dispersal is not related to range size, and when decoupled, species with larger ranges can have moderate to high levels of genetic structure (Lester et al. 2007). Another possibility is that gene flow from the center of a species' range reduces genetic structure and limits adaptation in peripheral populations, preventing the range size from increasing (Kirkpatrick and Barton 1997). Rather than being linked to dispersal per se, range size may impact genetic structure when it is positively correlated with opportunities for isolation. In such circumstances, the geographic structuring of genetic variation may be driven by events such as glacial cycling or the formation of landscape barriers (Letsch et al. 2016).
The combined effect of a species' range size and niche breadth on the observed genetic structure is complex. Because both factors may affect genetic variation, it is difficult to assess their relative importance in the absence of an empirical frame-work designed to disentangle one from the other (Slatyer et al. 2013). Niche breadth and range size are frequently positively correlated with one another (Brown 1984;Slatyer et al. 2013) and in many cases, geographic range and niche breadth are used as proxies for one another (Carscadden et al. 2020). However, many studies do not integrate appropriate controls to isolate their individual effects. In cases where range size and niche breadth are not positively correlated with one another (Williams et al. 2006;Kambach et al. 2019), life history characteristics may play an important role (Luna and Moreno 2010). Even when niche breadth and range size appear to covary, there are clearly other important factors involved (e.g., environmental heterogeneity and life stage (Finch et al. 2019) or species' morphological traits (Beck and Kitching 2007)).
In this study, the phylogeographic patterns of four high elevation sedge species (Carex spp., Cyperaceae) from the western United States are inferred to test the influence of niche breadth versus range size on genetic structure. The four sedge species investigated (C. bella L.H. Bailey, C. chalciolepis Holm, C. epapillosa Mack., and C. pelocarpa F.J. Herm.) likely share a Pleistocene origin ). In addition, traits that commonly impact genetic patterns (e.g., dispersal ability/mechanism and pollination syndrome) are similar in the four species (A. Reznicek, pers. comm., Ball and Reznicek 2002). Due to differences in their geographic distributions and niche breadths, these species form a natural experiment with two axes of variation (Fig. 1). Within the small-range pair of taxa, one species (C. bella) typically occupies wetter habitat (i.e., mesic meadows and subalpine forest), whereas the other species (C. chalciolepis) prefers drier montane habitats (i.e., rocky alpine slopes; Murray 1969, Massatti andKnowles 2014). Similarly, within the large-range pair of taxa, one species (C. epapillosa) shows an affinity to meadows in lower elevation habitats (similar to C. bella), and one species (C. pelocarpa) prefers higher elevation habitats (similar to C. chalciolepis) (Murray 1969). The species were selected for their observed differences in habitat preference, which were confirmed to correspond to differences in niche breadth via analyses in the present study. We implement a comparative phylogeographic approach to investigate the roles of niche breadth and range size in influencing genetic patterns in co-distributed species. The goal of the study is not to assess whether the pairs of co-distributed species are phylogeographically concordant (i.e., co-distributed species display spatially consistent geographic breaks in genetic structure), or to discover whether isolation by distance (IBD) or isolation by environment (IBE) describe phylogeographic patterns. Rather, we use the above metrics to identify similarities within species pairs united by either niche breadth or range size. The similarities in patterns between species are then used to make inferences about the relative importance of two orthogonal axes of species differences-niche breadth versus range size ( Fig. 1).

SAMPLE COLLECTION
Samples were field-collected from 10 sites that represented the majority of the range for each species (see maps in Figs. 2 and 3; Table S1). In our experimental design, the small range species were sampled from 10 identical locations and the large range species were also sampled from 10 identical locations. We used four species so that we could generate ample genomic data to sufficiently characterize genetic structure between 10 sampling locations per species; it would be ideal to include additional species, but the sampling effort in the field and lab required to obtain sufficient genetic data to confidently infer phylogeographic patterns for each species prevented us from adding more species to this study. Across all species, an average of 7.8 individuals were collected per sampling location. Individuals were sampled from the same sites for the two small-range species (C. bella and C. chalciolepis) across Colorado, New Mexico, and Utah (a range of approximately 250,000 km 2 ) and for the two large-range species (C. epapillosa and C. pelocarpa) across Colorado, Idaho, Montana, Nevada, Oregon, Utah, and Wyoming (a range of approximately 1,000,000 km 2 ). Due to the close proximity of different habitat types in montane regions, we were able to sample species with different habitat preferences from the same sampling sites. The average distance among sampled plants at a locality was 300m, and the minimum distance between samples was 35m to minimize the chance of sampling siblings or close relatives. For each sampled individual, leaf tissue was collected and placed directly into silica gel desiccant. Leaf tissue was kept in the dark and stored at room temperature until DNA extraction. DNA was extracted with DNeasy Plant Mini Kits (Qiagen, Hilden, Germany) using the manufacturer's standard protocol.

PROCESSING
We used a restriction site associated DNA sequencing (RAD-seq) approach to generate thousands of anonymous loci using six single-end Illumina HiSeq 2500 sequencing lanes following the protocol of Peterson et al. 2012. Briefly, two restriction enzymes, EcoRI and MseI, were used to digest genomic DNA, and Illumina adaptor sequences and unique 10 bp barcodes were ligated to the restriction sites. Pooled ligation products in each library were PCR amplified for 12 cycles and 400-500 bp fragments size selected using a Pippin Prep (Sage Science). All libraries were sequenced at The Centre for Applied Genomics (Hospital for Sick Children, Toronto, Canada). We used ipyrad 0.7.20 (Eaton 2014) to demultiplex and filter reads from each library, allowing one mismatch per barcode and using the most stringent filtering to identify and exclude contamination from adapters and/or primers. As some individuals were represented in more than one library, every library was initially demultiplexed and filtered separately (ipyrad steps one and two). Reads corresponding to the same individual were then combined and all reads were trimmed to a length of 41 bp using the fast-x toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). We used a 90% clustering threshold for de novo assembly; all other parameter settings that affect steps three through seven were set as default (see SI Methods for additional details). Across all species, between two and 10 individuals were sequenced per sampling location, for a total of 89 C. bella individuals, 97 C. chalciolepis individuals, 71 C. epapillosa individuals, and 55 C. pelocarpa individuals (SI Table S1, S2).

OCCUPIED BY EACH SPECIES
We used environmental variables that quantified soil characteristics, temperature, and precipitation as input for a PCA of the environmental space occupied by each species. This included 19 'bioclim' layers (Hijmans et al. 2005) and nine soil variables (bulk density, coarse fragment percentage by volume, percentage silt, percentage sand, percentage clay, pH in water, pH in KCL, cation exchange capacity, organic carbon content) (Hengl et al. 2017). We used the R packages 'dismo' (Hijmans et al. 2017)  and 'raster' (Hijmans 2016) to extract the environmental values at a spatial resolution of 50m 2 from each of the 28 variables using the GPS coordinates for hundreds of georeferenced specimen occurrences per species from vetted herbarium specimens (SI Table S3). The first three principal components, which explain a total of 63.0% of the variation, were used to define the environmental space occupied by each species. Discrepancies in how niche breadth changed over time among species could represent an uncontrolled axis of variation. To investigate how niche breath changed through time, we built niche models using MAXENT v3.3.3 (Phillips et al. 2006) for each species and hindcast them backwards in time to the mid-Holocene. We used bioclimatic layers BIO1, BIO2, BIO3, BIO7, BIO15, BIO17, BIO18 (Hijmans et al. 2005) with 30s resolution after removing other highly correlated layers (cutoff r > 0.7). Niche models for the present were constructed using 10,000 background points and 100 cross-validated replicates. Niche models were hindcast to the mid-Holocene using the CCSM4 climate projections (SI Fig. S1). Normalized Levin's (1968) B, a measure of niche breadth, was calculated for each species for the present and the mid-Holocene using the 'raster.breadth' function in the R package ENMTools (Warren et al. 2019). The measure of B implemented in ENMTools is interpreted as a quantification of the smoothness of geographic distribution of the habitat suitability scores. We do not consider B values to be precise estimates of the fundamental niche of each species; rather we use ENM analyses to infer the relative change of niche breadths of co-distributed species through time to act as a control against an unintended additional axis of variation (i.e., substantially different patterns of change in niche breadth over time among species).

STRUCTURING OF GENETIC VARIATION
Several complementary approaches were used to test for genetic structure and characterize genetic variation in each species. Pairwise F ST values between populations were calculated using the populations program in STACKS (Catchen et al. 2013). Genetic similarity among individuals was evaluated with a principal component analysis (PCA) implemented in the R packages SNPRelate (Zheng et al. 2012) and VEGAN (Oksanen et al. 2017). STRUCTURE (Pritchard et al. 2000) was also used to measure genetic structure within each species. We utilized the admixture model and correlated allele frequencies, and other settings were set as default. For values of K ranging from 1-11, 10 independent MCMC runs using 100,000 burn-in iterations were run followed by 150,000 iterations. Independent runs were combined using CLUMPP (Jakobsson and Rosenberg 2007) the results were visualized using DISTRUCT (Rosenberg 2004). The optimal K value for each species was determined using the ΔK method (Evanno et al. 2005) implemented in STRUCTURE Harvester (Earl and vonHoldt 2012). Phylogenetic relationships were inferred for each species using SVDQuartets (Chifman and Kubatko 2014), which accommodates differences in the coalescent history of loci, to investigate if relationships between individuals and populations were similar between co-distributed species pairs. Specifically, an exhaustive search strategy for sampling all possible quartets was used with 100 replicates of nonparametric bootstrapping to generate consensus trees using 50% majorityrule, which were visualized using the R packages APE (Paradis et al. 2004) and GGTREE (Yu et al. 2017). An additional SVDQuartets phylogeny was inferred using a set of shared loci from the four focal species plus two outgroup species to confirm that the four species are closely related, which serves as a control for similar life-history traits (Wiens et al. 2010).
For each species, linear regression was used to test for correlations between pairwise F ST values and the Euclidean geographic distance separating populations was measured to test for IBD using the R function dist with method = "euclidean" (R Core Team 2016). To test for an association between environmental and genetic variation after controlling for the effects of geography (i.e., IBE; Wang and Bradburd 2014), we similarly used linear regression to test for correlations between the position of populations in a PCA of environmental variables and the residuals from the linear regression of population pairwise geographic distance and F ST values. These correlation analyses were run sequentially starting with PC1 and carried through to PC3 of the environmental variables. If an environmental PC is highly correlated with the residuals of geographic distance and F ST , or with residuals from a subsequent regression, the environmental PC explains a substantial portion of the variation not explained by geographic distance.
Phylogeographic concordance in the co-distributed species pairs was assessed using two approaches. First, we compared the strength of the association between genetic variation (i.e., pairwise genetic distance) and geography across species using Procrustes analyses (Wang et al. 2010;Knowles et al. 2016), as quantified using the procrustes and protest functions in the VEGAN R package (Oksanen et al. 2017). A Procrustes analysis quantitatively compares the similarity between genes and geography across taxa by rotating and scaling the axes of a genetic PCA onto a geographical map with the goal of maximizing similarity (i.e., the sum of squared differences between the two input datasets are minimized; Knowles et al. 2016). Because each species pair was sampled from 10 identical locations, the association between genes and geography can inform the relative amount of concordance between species pairs. The Procrustes test statistic t 0 values were tested for significance using 10,000 permutations and results were compared between co-distributed species. The absolute magnitude of each t 0 value, which represents the strength of the association between genes and geography, is not as important for assessing concordance as the relative similarity in t 0 within co-distributed species pairs. For example, if the t 0 scores for the two large-range species are much more similar (e.g., within 10% of one another) than the t 0 scores for the two small-range species, it would imply that the two large-range species are more phylogeographically concordant.
A second approach compared pairwise F ST matrices for co-distributed species; first we standardized all pairwise F ST matrices by dividing each matrix by the highest value of pairwise F ST for that species. Then we used Mantel tests implemented in the VEGAN R package (Oksanen et al. 2017) to measure the similarity of pairwise F ST matrices of the small-range codistributed species and of the large-range co-distributed species; the magnitude and significance of the Mantel test would indicate the relative amount of concordance between species pairs. As with the Procrustes analysis, because both of the species pairs were sampled from identical locations, correlation between two matrices and a lack of correlation in the other two matrices, or substantially higher (e.g., twofold) correlation between two matrices relative to the correlation of the other two, implies greater concordance in phylogeographic patterns within one species pair when compared to the other pair. The goal of the concordance analyses is not to make a strong statement about the presence or absence of concordance in one or both species pairs, but rather to ascertain if one pair is more concordant relative to the other.

ENVIRONMENTAL SPACE
The environmental space occupied by the two species with broad niche breadth was substantially larger than that of the two species with narrow niche breadth, regardless of range size (Fig. 4). PC1 was strongly (defined as r 2 > 0.5) and significantly positively correlated with six input variables: BIO1 (annual mean temperature; r 2 = 0.896, P < 0.001), BIO3 (isothermality; r 2 = 0.575, P < 0.001); BIO5 (mean temperature of the warmest month; r 2 = 0.718, P < 0.001), BIO10 (mean temperature of the warmest quarter; r 2 = 0.802, P < 0.001), BIO11 (mean temperature of the coldest quarter; r 2 = 0.824, P < 0.001), and BIO15 (precipitation seasonality; r 2 = 0.612, P < 0.001). Because five of the six environmental variables correlated with PC1 are variables associated with temperature, we refer to PC1 as an environmental axis of temperature variation. PC2 was strongly and significantly negatively correlated with BIO12 (annual precipitation; r 2 = 0.737, P < 0.001), BIO13 (precipitation of the wettest month; r 2 = 0.541, P < 0.001), BIO16 (precipitation of the wettest quarter; r 2 = 0.566, P < 0.001), and BIO19 (precipitation of the coldest quarter; r 2 = 0.560, P < 0.001), suggesting it is representative of precipitation variation. None of the soil variables were strongly and significantly correlated with any of the PC axes.
The MAXENT results indicated that environmental niche breadth was not static over time; each of the four species underwent an increase in its niche breadth from the mid-Holocene to the present (Fig. S1). Because it is impossible to co-estimate niche breadth and range size simultaneously using niche modeling, our focus was whether the relative niche breadths within species pairs were consistent over time: the niche breadth of C. bella is larger than C. chalciolepis in both the past and the present, as is the niche breadth of C. epapillosa relative to C. pelocarpa (Table S4). Comparisons of Levin's B confirmed these results and showed an increase in niche breadth for each species since the mid-Holocene. Carex bella and C. chalciolepis experienced a 43.1% and 40.0% increase in niche breadth, respectively; meanwhile C. epapillosa increased by 8.7% and C. pelocarpa by 26.9% (Table S4).

STRUCTURING OF GENETIC VARIATION
Measures of genetic variation showed significant structure in all species (Table 1; Figs. S2-S4). Genetic structure was both qualitatively and quantitatively more similar between species based on their range size, but did not show consistent similarities based on their niche breadth. Specifically, the two small-range species (C. bella and C. chalciolepis; Fig. 2) exhibited similar genetic patterns based on PCAs and STRUCTURE results, with some noted differences. In both of the small-range species, genetic similarity of individuals tracks a general east-west gradient, with the most westerly population (Tushar) dominating PC2 (Fig. 2). In C. chalciolepis, but less so in C. bella, individuals from the different mountain ranges were genetically more distinct (Fig. 2), which is reflected in a higher average pairwise F ST in C. chalciolepis (F ST = 0.26) than C. bella (F ST = 0.10; Table 1; Fig. S5). The STRUCTURE plot for C. bella (optimal K = 3) conformed with results from the genetic PCA; La Sals and Tushar each formed distinct entities while Lamphier is intermediate between La Sals and the other seven easternmost populations (Fig. S4). The plot of the optimal K ( = 4) for C. chalciolepis mirrored the genetic PCA and similarly revealed geographic clustering of populations (Fig. S4). The four North-central populations cluster together (orange) as do three southernmost populations (purple). As in C. bella, La Sals and Tushar are distinct entities; additionally, Lizard Head is split between La Sals and the three southern populations-recapitulating genetic PC1. Within the small-range species pair, the greater geographic structure is found in the species with a narrow niche breadth (i.e., C. chalciolepis).
The degree of structure in the small-range species was less than that of both the large-range species, C. pelocarpa and C. epapillosa (F ST = 0.35 and F ST = 0.35, respectively), which did not show any corresponding difference associated with niche breadth (i.e., C. pelocarpa occupies a small environmental space relative to C. epapillosa). The two large-range species showed both qualitative and quantitative similarity in genetic structure based not only on the aforementioned F ST analysis, but also the PCAs and STRUCTURE analysis ( Fig. 3; Fig. S4). For both large range species, K = 2 was optimal. In C. epapillosa, the STRUC-TURE assignment of each individual entirely to one of the two groups completely aligned with genetic PC1 (Fig. S4). Notably, certain populations (Sawtooths, Steens, Wallowas) had some individuals occurring on opposite extremes of genetic PC1 and individuals from these populations were assigned to both STRUC-TURE clusters ( Fig. 3; Fig. S4). The STRUCTURE results for C. pelocarpa corresponded to genetic PC1: the purple STRUC-TURE cluster represented populations with extreme negative values on PC1, the orange STRUCTURE cluster represented populations with extreme positive values on PC1, and three populations near zero on PC1 were depicted as split between the two STRUC-TURE clusters ( Fig. 3; Fig. S4). In contrast to C. epapillosa, no individuals from the same population of C. pelocarpa were assigned to different STRUCTURE clusters or occurred in opposite extremes of PC space. The SVDQuartets coalescent trees were largely consistent with the genetic PCAs and STRUCTURE analyses in that individuals and/or populations that were phylogenetically proximate were also close together in PC space and in STRUCTURE clusters ( Fig. 2 and 3; Fig. S2 and S4). The rooted SVDQuartets species tree recapitulated interspecific relationships from previous studies ) and confirmed the species are closely related (Fig. S3).
A significant relationship between pairwise geographical distance and genetic differentiation between populations (i.e., IBD) was observed in the small-range species (C. bella: r 2 = 0.321, P < 0.001; C. chalciolepis: r 2 = 0.355, P < 0.001; Fig. 5), but not the two large-range species (C. epapillosa: r 2 = 0.0573, P > 0.05; C. pelocarpa: r 2 = 0.0719, P > 0.05; Fig. 5). When comparing the degree of genetic differentiation associated with environmental differences, after controlling for the effects of geographic distance among individuals, each species shows a significant association (Table 2, Fig. 6). However, the two small-range species show a response to the precipitation axis (i.e., PC2), while the two large-range species respond to the temperature axis (i.e., PC1). For all species, an environmental axis explains a significant amount of the genetic variation not explained by geography; range size, not niche breadth shows similarity in patterns of IBE. Notably, the regression coefficients were also much lower in the two small-range species relative to the large-range species (Fig. 6).
All measures of phylogeographic concordance detected greater concordance in the small-range species pair. The

Figure 5. Regressions comparing pairwise geographic distance between populations (in degrees) and pairwise F ST between populations.
In the small-range species, there was a significant relationship, implying a pattern of isolation by distance. In the large-range species, there was no significant relationship between geographic and genetic distance.
Procrustes test statistic (t 0 ) values that measure the similarity of genes and geography were more similar between the two smallrange species (t 0 = 0.619, P < 0.001 for C. bella vs. t 0 = 0.669, P < 0.001 for C. chalciolepis) than they were between the two large-range species (t 0 = 0.477, P < 0.001 for C. epapillosa vs. t 0 = 0.659, P < 0.001 for C. pelocarpa). While the magnitude of the Procrustes test statistics (t 0 ) does not explicitly assess concordance, we would expect that the difference between t 0 values in co-distributed taxa decreases as concordance increases, especially in our study where the co-distributed species were sampled in identical locations. Mantel comparisons of F ST matrices of co-distributed pairs showed correlation between standardized matrices in the small-range species (Mantel statistic r = 0.597; P < 0.05) but not in the large-range species (Mantel statistic r = 0.0308; P = 0.441). The significant and strong correlation between F ST matrices in the small-range species did not implic-itly mean these species are concordant, but it demonstrated that the small-range species are more concordant than the large-range species, which showed no association between the F ST matrices.

Discussion
Although range size and niche breadth have been independently shown to structure population genetic variation within species (e.g., Lovell andMcKay 2015, Huang et al. 2017), disentangling their effects is notoriously difficult because they are often tightly linked (Brown 1984;Slatyer et al. 2013), and their relationship with genetic structure can be scale dependent (Kambach et al. 2019). With de-coupled variation along these axes (Fig. 1) in pairs of co-distributed montane sedge species (Figs. 2 & 3), we show that similarity in phylogeographic metrics covaried with range size, but not with niche breadth. Although only four species were investigated, this result was supported by similarity and dissimilarity, respectively, across multiple phylogeographic metrics, including signatures of IBD, IBE (including the dominant climatic axis-temperature versus precipitation), tests of phylogeographic concordance, and measures of genetic structure (Figs. 2,3,5,6). Given the design of this natural experiment, there could have been a range of outcomes: either niche breadth or range size could have dominated genetic structure, or neither factor could have a clear effect, or a confounding result in which one variable predominated at one spatial scale and the other variable predominated at the other scale. However, all evidence supports range size, but not niche breadth, as the primary determinant shaping phylogeographic patterns. Below we discuss the broader implications for inferences of processes in comparative phylogeography, and the spatial scale of concordance among species in montane systems.

FROM CONCORDANT GENETIC VARIATION
Although a signal of similarity predominates within codistributed species pairs, subtle differences between codistributed species reveal insights about this study system and the processes inferred to be important in a comparative phylogeographic framework. In this framework, highly concordant co-distributed species would display congruent spatial patterns of genetic variation, whereas co-distributed species with low phylogeographic concordance would show substantial differences in the geographic distribution of genetic variation (e.g., Table 2. For each species, the correlation coefficient (r 2 ) comparing genetic distance (pairwise F ST ) and geographic distance ("geography" column), residuals from the first correlation and environmental PC1 ("PC1" column), residuals from the second correlation and environmental PC2 ("PC2" column), residuals from the third correlation and environmental PC3 ("PC3" column) are displayed. In the three rightmost columns, the total proportion of genetic variation explained by the environment ("Environment" column), the total genetic variation explained by geography plus environment ("Total" column), and the ratio of the variation explained by geography relative to environment ("Geo/Env" column). For the first four columns, sampling locations clustering into different genetic groups). One model often invoked in montane comparative phylogeography studies is the 'expanding-contracting archipelago' (ECA) model (Hewitt 1996;DeChaine and Martin 2005). This model posits that during Quaternary glacial cycling, populations of montane taxa moved downslope during glacial periods, and upslope into fragmented sky islands during warmer interglacial periods (Hewitt 1996;Jackson and Overpeck 2000;DeChaine and Martin 2005). The ECA model predicts that lower elevation taxa should have greater connectivity between populations during glaciations. At our smaller spatial scale, the contrasting phylogeographic patterns in C. bella and C. chalciolepis support the ECA model because of the relative amount of genetic differentiation (Fig. 2, Table 1). A key difference between the two small-range species, which have largely concordant genetic variation (Fig. 2), is the magnitude and range of F ST values in the two species (i.e., higher variance and greater overall genetic differentiation among populations of C. chalciolepis than C. bella; Table 1, SI Fig. S5). One explanation for this pattern may relate to differences in the niche breadth of the taxa. Stronger genetic structure may result in the narrow-niche breadth species C. chalciolepis because of greater amounts of unsuitable habitat separating populations. A stronger signature of IBD in C. chalciolepis relative to C. bella is reflected by a larger regression coefficient (Fig. 5), and provides additional evidence of greater isolation experienced by C. chalciolepis during Pleistocene glacial cycling than C. bella, which by virtue of its distribution in lower elevations, would have experienced more historical connectivity among what are now isolated populations across current mountain ranges (Lorenz-Lemke et al. 2010).
The ECA model may also apply to the large-range species, but signatures of greater historical connectivity in C. epapillosa are masked by multiple refugia and sampling locations containing individuals from multiple ancestral sources-at least two distinct glacial refugia may have led to multiple historically isolated lineages within C. epapillosa. In some cases, extremely high or low pairwise F ST values in C. epapillosa are driven by populations comprised of two groups of genetically distinct individuals exhibiting no signs of admixture (Fig. S4). For example, within each of the Sawtooths, Wallowas, and Wind River populations, there are individuals that appear on both extremes of PC1 (Fig. 3) and are in separate SVDQuartets clades (SI Fig. S2, SI Fig. S3). Cryptic species within C. epapillosa, possibly caused by genetic incompatibilities due to holocentric chromosomes common in the genus (Hipp 2007), could also explain the instances of high differentiation between individuals within populations. On smaller scales within the large-range species-specifically, certain geographically-united subsets of populations with low F ST values in C. epapillosa (e.g., only the easternmost populations: Beartooths, Bighorns, Uintahs, Wind River)-there was likely greater historical connectivity in the lower elevation C. epapillosa than the higher elevation C. pelocarpa, which has less variance in pairwise F ST values (Table 1). Therefore, we argue the ECA model applies at both spatial scales, but it is masked at larger spatial scales due to the increased probability of idiosyncratic responses to historical processes at larger scales (i.e., multiple distinct refugia in this case). In contrast to many previous studies, our study enabled an examination of the degree to which species follow the ECA model at differing spatial scales, revealing an important effect of geographic scale on interpretation of the ECA model.

PHYLOGEOGRAPHY
In previous studies, when neither concordance nor discordance predominated at all spatial scales, typically patterns of largescale concordance and small-scale discordance were reported (DeChaine and Martin 2005; Barber et al. 2006;Roe et al. 2011;Chen et al. 2016;Barrow et al. 2018). These patterns were attributed to the operation of scale-dependent processes (e.g., Comes andKadereit 2003, Galbreath et al. 2010). The large and small scales reported in these studies approximately correspond to our study-"large scale" refers to a continental scale (or the marine equivalent) whereas 'small scale' represents a collection of populations that span a geographically continuous subset of the species range, although 'small scale' may still represent hundreds of thousands of km 2 . At large scales, shared historical processes frequently explain concordance (Avise 2009;Papadopoulou and Knowles 2016), whereas at smaller scales, discordance is explained by species' individualistic responses to environment changes (DeChaine and Martin 2005) or speciesspecific trait differences, such as dispersal ability (Roe et al. 2011). In our small-range species, greater concordance may be explained by shared historical processes because the range, although small relative to our large-range species, is ∼250,000 km 2 . The lack of concordance in the large range species could be caused by species-specific traits, or by chance. Across large geographic scales, species may be subject to the idiosyncracies of climate-induced distributional shifts and a tipping point may exist beyond which expanding the spatial scale will reduce the likelihood of concordance. Montane species may be readily displaced into multiple different refugia, especially if species differ in their respective niches. Alternatively, but not mutually exclusively, we may be simultaneously observing concordance due to similar effects of historical processes and discordance due to species-specific attributes.

LIMITATIONS OF THE NATURAL EXPERIMENT
We acknowledge that it is notoriously difficult to control for all the myriad processes that influence phylogeographic patterns (e.g., Dellicour et al. 2015); there are covariates (e.g., degree of population isolation) in our study that may magnify or diminish genetic signals attributed to niche breadth or range size. In addition, there may be other species-specific attributes that affect discordance on different spatial scales; niche breadth is an obvious possibility because it was explicitly incorporated. The different habitat affinities of each species appear correlated with niche breadth but there could be additional axes of variation confounding results. Habitat patchiness, common in montane species occupying heterogeneous environments, is another important axis of variation that may complicate interpretations. The two higherelevation narrow-niche species have less connectivity between populations presently and historically during glaciation than the two lower-elevation species (Fig. S1, Table S4). Therefore, in this system habitat patchiness appears to be negatively associated with niche breadth. Other studies may be impacted by an interaction effect due to an independent habitat patchiness axis because there is no universal relationship between habitat patchiness and niche breadth (Kellner et al. 2019). The experimental design in this study can tease apart the two primary axes of variation, but also represents a limitation in that it is difficult to add replicates (i.e., additional co-distributed species). We use a total of four species because the required sampling effort in the field and lab to generate sufficient genomic data to confidently infer phylogeographic patterns precluded adding more species to the present study. Future studies investigating this question in the same geographic region would be valuable for comparison with the present study.

Conclusions
Here, we provide evidence that similarity between species' range sizes influenced the structuring of genetic variation to a greater extent than similarity in environmental space occupied by species. Reviews of empirical studies have noted a positive (Slatyer et al. 2013), although scale-dependent (Kambach et al. 2019) correlation between niche breadth and range size, which suggests that their interaction may have a confounding effect. Our use of a replicated natural experiment with two key axes of variation allowed us to isolate and test the effect of the differences in range size and niche breadth. Although we cannot control for all possible covariates in a natural experiment, we found strong evidence that range size was a more important factor than niche breadth for shaping genetic patterns in montane systems. Our results suggest that the degree of concordance due to historical processes, and discordance due to species-specific attributes, can co-exist and simultaneously explain the degree of phylogeographic concordance in co-distributed taxa. The interaction of niche breadth and range size highlighted a tipping point where the scale gets so large that idiosyncratic histories may predominate over shared historical processes. The comparative framework of our study supports the growing body of literature that investigates the causes of phylogeographic concordance, or a lack thereof, in montane systems.