Trait variation is the basis of evolutionary biology and a fundamental component on which natural selection acts, along with inheritance, differential fitness and competition (Geber & Griffen, 2003). At a metapopulation scale, climatic and geologic pressures contribute to genetic and morphological diversity (Anderson & Ferree, 2010;Kahilainen, van Nouhuys, Schulz & Saastamoinen, 2018). Within a population, however, there is longstanding evidence that individuals of the same species undergo similar pressures of selection which often ultimately leads to the homogenization of phenotypes within that population (Benkman & Parchman, 2013;Jones & Reithel, 2001;Smithson & Macnair, 1996;Stebbins, 1970). As a result, many studies on trait variation have focused on interspecific or among-population variation with less focus on cases of within-population phenotypic diversity (i.e. polymorphism) (Sapir, Gallagher & Senden, 2021).
However, polymorphic species provide a unique opportunity to investigate morphological diversification and genetic differentiation at both local and metapopulation scales, and provide opportunities to characterize the processes that deviate from classic evolutionary theories of homogenizing selection; those that result in broad standing phenotypic variation within species and within populations such as balancing selection, temporally varying fitness, heterozygote advantage or frequency-dependent selection (for a review, see Sapir et al., 2021).
Key to understanding the processes that maintain genetic and morphological diversity across spatial scales is the characterization of population structure (Twyford, Wong & Friedman, 2020) and functional traits (Assis, Patton, Hubbe & Marroig, 2016). One model that reconciles local adaptation in the face of gene flow is a ring species complex. Ring species provide an opportunity to explore the evolution of geographically distributed traits through the generation and maintenance of morphological gradation across genetically interconnected subpopulations (Pereira & Wake, 2015). Under a ring species model, subpopulations are expected to have accumulated mutations, and experience genetic drift and local adaptation, with differences often but not necessarily manifested into different phenotypes (Kuchta & Wake, 2016). In his early descriptions of the then-presumed ring species Ensatina eschscholtzii, Stebbins' defined the main requirements as having a founding population that diversified around poor habitat producing a ring-shaped distribution, leaving contiguous populations with the possibility of gene flow that undergo local adaptations, causing divergence in traits associated with reproduction such that the farthest extending populations can no longer interbreed when they come into contact (Stebbins, 1957). Ring species are defined by this historical geographical expansion, but patterns of phenotypic differentiation vary among documented ring species. Investigations of ring species complexes have uncovered examples of genetic divergence among close relatives with strong phenotypic similarity. In the angiosperm Euphorbia tithymaloides, Cacho and Baum leveraged molecular evidence to reveal that two subspecies considered to be closest relatives based on morphology were actually the most genetically distant, being the terminal ends of the ring having undergone subsequent convergence in phenotype (Cacho & Baum, 2012). The convergent evolution appears to be associated with the local ecology, specifically pollination biology (Cacho, Monteverde-Suárez & McIntyre, 2019). Irwin et al. propose that similar environmental pressures are responsible for the parallel evolution of song length and body size in the ring species of greenish warblers, Phylloscopus trochiloides, while other traits such as plumage patterning and song structure have diverged (Irwin, 2000;Irwin, Alcaide, Delmore, Irwin & Owens, 2016;Irwin & Irwin, 2002;Irwin, Thimgan & Irwin, 2008).
The California endemic Calochortus venustus Douglas ex Benth.
(Liliaceae), commonly known as the butterfly mariposa lily, displays hypervariable floral traits: the petal colour of individual plants ranges from white, yellow, purple, pink, to deep red; a distal petal spot is often but not always present; and anther colour varies from whiteyellow or pink-lilac (Baldwin et al., 2012) (Figure 1a). Such traits would normally justify the delineation of different species, however morphologically, C. venustus is consistently identified by having a square nectary densely covered in yellow trichomes and a dark red blotch above the nectary. In addition, no other Calochortus species has a distal petal spot (A.I.H. personal observations). These characters differentiate C. venustus from its sister species, C. catalinae, and from a close relative, C. argillosus (Patterson & Givnish, 2004). This perennial herbaceous lily is distributed around the Central Valley of California, inhabiting diverse habitats such as valley grasslands, foothill woodlands and yellow pine forests across an elevational range of 300-2700 m (Baldwin et al., 2012). While the highly polymorphic nature of this species has led some to propose taxonomic subgroupings as varieties, none of these have been formally accepted. The distribution of petal phenotypes shows variation both within and between populations, and one colour may occur more frequently in one locality while a few kilometres away another variant may predominate (Figure 1b). A previous study conducted a quantitative analysis of petal colour, form and shape and found no infraspecific groupings based on phenotype (R.K. Zebell, Master's Thesis 1993, personal communication). at both the metapopulation and regional scales with the goal of disentangling hypervariable morphological diversity from genetic differentiation. To address these hypotheses, we specifically ask the following questions:
Populations were surveyed along the Coast Ranges and Transverse Ranges in 2018 to approximate phenotypic ratios for adaxial petal pigmentation. These surveys were completed again in 2019 with the addition of populations along the Sierra Nevada. Data from both years were summarized across five broad regions (Figure 1a). Floral tissue was collected and preserved in silica for 178 individuals encompassing the full geographical distribution and floral phenotypic diversity (Figure 1b). Geographical coordinates were marked with a Garmin etrex30, and floral tissue from the same plant was photographed. Voucher specimens were collected to be deposited into the L.H. Bailey Hortorium (BH) at Cornell University.
Field photographs were used to acquire individual phenotypes for five floral traits: (1) pigmentation of adaxial petal was binned into seven states: white, yellow, pink, purple, red, half pink/half yellow and half pink/half white; (2) abaxial petal was scored as same or different from adaxial petal pigmentation; (3) spots on adaxial petals were scored as present or absent; (4) anther colour was scored as pink/purple or yellow/white; and (5) filament colour was scored as pigmented (pink) or unpigmented (white or slight pink tinge). We summarized character states within six regions by counting the number of individuals with each state for each trait in the R package 'ggplot2' v3.3.2 (Wickham, 2016). To calculate statistical significance of the frequency of each character state within each population, we ran a test for equal proportions using prop.test in R. To evaluate how traits are correlated across morphological space, we ran a principal components analysis (PCA) with the function PCAmix in the R package 'PCAmixdata' v3.1 (Chavent, Kuentz-Simonet, Labenne & Saracco, 2017).
DNA was extracted from tissue preserved in silica following a modified SDS protocol (Edwards, Johnstone & Thompson, 1991;Konieczny & Ausubel, 1993). A total of 190 DNA samples including nine individuals of Calochortus superbus, C. catalinae and C. argillosus were sent to the University of Wisconsin Biotechnology centre for library preparation for single-end 100 bp restrictionsite-associated DNA sequencing (RAD-Seq) (Andrews, Good, Miller, Luikart & Hohenlohe, 2016). In short, libraries were prepared with the apek1 restriction enzyme and multiplexed. Pooled samples were sequenced on one lane of Illumina NovaSeq 6000 S1 100SR at QB3 Genomics at UC Berkeley (Berkeley, California, USA).
Accession 1873 was chosen as a reference to build a draft genome for downstream analyses. DNA was sent to the QB3 sequencing core (Berkeley, California, USA) for Illumina sequencing using 150 PE NovaSeq S4 to generate approximately 200 GB of data. Additionally, five Oxford Nanopore flow cells were sequenced using the Rapid Sequencing library kit (SQK-RAD004). The raw sequence data were cleaned and processed as described in the Supplementary Materials and Methods.
Multiple assembly methods were performed including (1) a hybrid approach with both Illumina and error-correct Nanopore with ABYSS (Simpson et al., 2009), (2) long-read only with FLYE v2.7.1-b1590 (Kolmogorov, Yuan, Lin & Pevzner, 2019), (3) long-read only with Nanopore reads and the super contigs generated in MASURCA v3.4.0 (Zimin et al., 2013(Zimin et al., , 2017) ) and (4) a hybrid approach with MASURCA.
For a full description of assembly methods, see Supplementary Materials and Methods. A reference-based approach was used to minimize biases that arise with bioinformatic processing of RAD-Seq (Shafer et al., 2017); therefore, the MaSuRCA assembly was used since it is the most complete with the fewest number of contigs.
Raw sequence reads were demultiplexed following the STACKS v2.5 (Catchen, Hohenlohe, Bassham, Amores & Cresko, 2013) pipeline and quality checked with FASTQC v0.11.8 (S. Andrews, 2010). Reads were mapped to the reference genome with BWA-MEM v0.7.12 (Li, 2013). Single nucleotide polymorphisms (SNPs) and genotypes were called using ref_map.pl in STACKS without prior geographical information. We removed outgroups, four poorly sequenced ingroup samples, and noise from minor allele frequency (MAF), genotype and low SNP counts in VCFTOOLS v4.2 (Danecek et al., 2011) as follows: --mac 3, --max-missing 0.7, --maf 0.05. We identified linkage disequilibrium (LD) by calculating r 2 using POPLDDECAY v3.40 (C. Zhang, Dong, Xu, He & Yang, 2019), and filtered SNPs with an r 2 higher than 0.15 with 5 kb steps in a 10 kb window with PLINK v2.00a1LM (Chang et al., 2015). Different filtering regimes of MAF and read depth can potentially bias results; therefore, we compared PCA results for MAF 0.02-0.1 in 0.02 intervals and read depth 10X (Benjelloun et al., 2019;Linck & Battey, 2019;Malomane et al., 2018).
A PCA of individual genotypes was used to infer population clusters in the R packages 'SNPRelate' v1.22.0 and 'gdsfmt' v1.24.1 (Price et al., 2006;Zheng et al., 2012). Sparse non-negative matrix factorization (sNMF), a model-based algorithm implemented in the R package 'LEA' v3.0.0 (Frichot & François, 2015), was used to infer least-squares estimates of ancestry proportions (Frichot, Mathieu, Trouillon, Bouchard & François, 2014). We ran the algorithm for K = 1-10 with 10 repetitions for each K. The most likely ancestral K value was determined by cross entropy criterion which estimates assignment error based on the prediction of masked genotypes (Janes et al., 2017).
Populations were discretized following PCA results. Individual heterozygosity was calculated from variant sites in VCFTOOLS (Danecek et al., 2011) and summarized in the R package 'ggplot2' (Wickham, 2016). Population-level heterozygosity, inbreeding coefficient (F IS ) and pi were computed from variant sites in STACKS.
We also calculated global heterozygosity using both variant and invariant sites in ANGSD v0.935 (Korneliussen, Albrechtsen & Nielsen, 2014) which uses the site frequency spectrum for each individual to directly estimate heterozygosity. Pairwise population F ST values were calculated from variant sites in the R package 'SNPRelate' (Zheng et al., 2012). We explored spatial genetic differentiation at two scales by testing isolation by distance (IBD), or genetic distance by geographical distance between individuals, both across the metapopulation and within individual mountain ranges, namely the Coast Ranges, Transverse Ranges and Sierra Nevada Ranges using a Mantel test with 10,000 repetitions in the R package 'adegenet' v2.1.3 (Jombart, 2008).
Nine outgroups representing three species were incorporated into the same dataset used for population genetic analyses. A maximum likelihood tree was inferred in RAXML v8.2.11 (Stamatakis, 2014) under a GTRCAT model, 1000 bootstrap iterations, 25 distinct rate categories. In addition, a python script from https:// github.com/btmar tin72 1/raxml_ascbias was used to count and remove invariant sites. Invariant site counts 5937, 6770, 6538, 5896 for A, C, G, T, respectively, were used to apply a Stamatakis ascertainment correction to the model using 75,272 variant sites.
The resulting bipartitions tree was visualized in FIGTREE v1.4.4 (Rambaut, 2012). We also explored clustering and network relationships using NeighborNet in SPLITSTREE4 v4.16.2 (Huson & Bryant, 2006) which does not assume strictly bifurcating hierarchical relationships between samples.
Inference of ancestral range evolution was completed using the R package 'BioGeoBE ARS' v1.1.2 (Mat zke, 2013) which uses a maximum likelihood approach to calculate the probabilities of biogeographical events at each node in a phylogenetic tree. An ultrametric tree without outgroups was generated from the RAXML tree using the chronos function in R package 'ape' v5.5 (Paradis & Schliep, 2019). Four clock models were fitted to the ML tree: strict, discrete, correlated and relaxed; the strict clock model had the lowest phi information criterion and was selected to produce the ultrametric tree (Paradis, 2013).
We pruned the tree to include one accession from each clade within each region (n = 10) to avoid sample size bias. We used the Dispersal Extinction Cladogenesis (DEC) model (Ree & Smith, 2008) with the addition of a jump dispersal parameter, J, (Matzke, 2013), although use of this parameter is controversial (Klaus & Matzke, 2020;Ree & Sanmartín, 2018). This model explores dispersal-driven, founderevent range evolution consistent with a stepping-stone model of expansion and concordant with results from our other analyses.
We inferred haplotype networks to inform our understanding of haplotype evolution among individuals. We selected contigs containing more than one SNP to phase genotypes in BEAGLE v5.2 (Browning, Zhou & Browning, 2018). A second dataset was created using SNPs with at most 10% missing data (--max-missing 0.9) in VCFTOOLS. Haplotype networks were constructed under the TCS model and visualized in POPART (Leigh & Bryant, 2015).
To visualize areas of historically high and low migration, we implemented Estimating Effective Migration Surfaces (EEMS) (Petkova, Novembre & Stephens, 2016). We ran three MCMC chains with the following parameters: 350 demes, 50,000,000 MCMC iterations and 10,000,000 burn-in iterations. The posterior distributions were averaged to plot the final values in R package 'rEEMSplots' v0.0.0.9000.
We characterized phenotypic polymorphisms and associated geographical patterns in C. venustus. A strikingly wide range of petal pigmentation was observed across a spectrum of colours ranging from white, through shades of yellow, pink, and purple, to deep red (Figure 1a). Southern populations were found to be highly polymorphic, while populations gradually decreased in numbers of pigmented individuals towards the northern latitudes on both sides of the Central Valley (Figure 1b). Looking across altitudinal transect s, white phenot ypes were found to be distributed throughout the full range while pigmented phenotypes emerge mostly in higher elevations (Figure 2). The altitudinal distribution of phenotypes across six populations is illustrated in Figure S1. Petal spot presence and abaxial petal pigmentation (Figure 3a) show an eastwest differentiation with polymorphism focused in southern and higher elevation populations (Figure 3d, e, respectively). Variability in the colour of the anthers and filaments of the stamen show no regional differences except for the southern Sierra region which had a significantly lower proportion of pink filaments (Figure S2).
The PCA of morphological floral traits resulted in three suites of traits that commonly occur: (1) white adaxial petals with a distal petal spot, yellow-to-white anthers and white filaments; (2) pink adaxial petals with different abaxial pigmentation and pink filaments; and (3) red and purple adaxial petals without a distal petal spot and pink-to-purple anthers (Figure S3). These colour combinations may reflect selection on coloration patterns and/or correlated expression of genes associated with these pigmentation patterns.
The different reference assembly approaches yielded drastically different qualities of assemblies in terms of total bases assembled, number and size of contigs, and BUSCO genes found (see Supplementary Methods). MaSuRCA resulted in the most complete assembly of 2.2 GB in 275,791 contigs with a N50 of 11,884 bp. The polished MaSuRCA assembly also yielded the highest proportion BUSCO genes (Figure S4). Mapping RAD-Seq reads to the reference genome and calling SNPs resulted in 2,540,973 loci. Four samples and two outgroups were removed due to poor mapping to the reference genome. After filtering and LD pruning, a total of 94,744 SNPs out of an original 9,385,103 were kept for 174 individuals of C. venustus. Read depth per individual ranged from 2.81 to 8.97X with an average of 7.07X coverage (Table S1) and standard deviation of 0.912. The PCA explained 3.31%, 2.09%, 1.60% and 1.31% of variation in axes 1-4, respectively. A plot of axes 1 and 2 illustrates a 'Y' fork divergence between Coast Range and Sierra Nevada populations (Figure 4). These PCs capture geography, with axis 1 representative of longitude and axis 2 of latitude. Differentiation between the two northern populations is reflected in axis 3, in addition to three individuals from the central Sierra Nevada showing distinctly unique genotypes (Figure S5a). Axis 4 illustrates divergence within the southern Sierra Nevada and Transverse Range populations (Figure S5b). We found no evidence of phenotypic cohesion for adaxial pigmentation phenotype providing evidence for gene flow between phenotypes in polymorphic populations (Figure S6). Six individuals that exhibited 'marginal picotee', in which margins of the petals are a different colour from the base (Figure 1a top left), clustered separately from the rest of the Transverse Ranges individuals (Figure S6). These six individuals were geographically clustered between two polymorphic populations; therefore, we do not know whether this genetic signal is reflective of geographical pattern or the marginal picotee phenotype. The PCAs under various MAF and read depth filtering stringencies reflect very similar clustering patterns as our principal filtering scheme (Figure S7).
Based on sNMF analysis, genotypes were best explained by a model of K = 5 ancestral populations (Figure S8). The ancestry plot at K = 5 corroborates that population structure is geographically cohesive with varying levels of admixture between neighbouring populations (Figure 5c). Ancestry assignment models for K = 2-4 show in sequential order: differentiation of the Coast Ranges, increasing divergence of individuals in the Sierra Nevada and a trend of increasing genetic differentiation with increasing distance between individuals of the Transverse Ranges and the Sierra Nevada (Figure S9).
Mean population heterozygosity calculated from both invariant and variant sites was similar for all populations except the north Coast Range (Table 1), and average individual heterozygosity was 0.0187 (Table S1). Population-level summary statistics calculated from variant sites (Table 1) show similar values for mean nucleotide diversity (pi) and mean expected heterozygosity across populations.
Values using only variant sites show the same patterns of heterozygosity as those calculated from both variant and invariant sites.
A regression between the two datasets with populations shows an adjusted r 2 = 0.961 (pvalue = 3.65e-4) and with individuals the adjusted r 2 = 0.841 (pvalue = 2.2e-16). Pairwise F ST analysis reveals overall low genetic differentiation (average 0.05068), that is, there is high gene flow between populations, and a trend of increasing genetic differentiation with increasing geographical distance. Tests for isolation by distance (IBD) across individuals revealed clusters with statistically significant but low IBD for the metapopulation. This led us to test for IBD separately within each major mountain range (Table 3), revealing the strongest regional trend in the Coast Ranges.
The RAxML tree illustrated with ancestral range reconstructions shows that C. venustus is a monophyletic group that likely originated in the southern Coast Ranges (Figure 5b). Populations expanded in parallel from south to north along major mountain ranges surround- Coast Range and (2) the Transverse Ranges up to where it connects to the south Sierra Nevada population (Figure S13).
Calochortus venustus as a species displays striking variation in floral pigmentation and patterning both within and between populations, making it one of the most polymorphic systems to be described to date. Polymorphism alone does not suggest that polymorphic traits are products of selection; however, geographical and altitudinal patterns of petal pigmentation and petal spots provide preliminary evidence of and potential for local adaptation (Figure 3b-e). Notably, these patterns could also be the result of neutral processes such as random genetic drift of alleles underlying morphological traits or bottlenecks during range expansion, and reciprocal transplant experiments are necessary to definitively test between neutral and adaptive hypotheses.
Investigating population dynamics at various spatial scales in C. venustus is essential to understanding the evolutionary processes that may drive spatial phenotypic hypervariability despite ongoing gene flow.
There is strong geographical cohesion irrespective of floral phenotype across spatial scales and gene flow between all colour variations, affirming that this wide phenotypic spectrum is indeed functioning as a single polymorphic species. We found both east-west and north-south genetic divergence, and some morphological traits follow the pattern of east-west differentiation while others mirror the south-north divergence. The strongest genetic differentiation is between northern populations at terminal ends of the ring-like range expansion. Phylogeographical and range evolution analyses reveal a south-to-north bidirectional expansion around the California Central Valley with colonization from the Transverse Ranges to the southern foothills of the Sierra Nevada.
Elucidation of this evolutionary history paired with phenotypic surveys unveiled parallel evolution towards the white phenotype in the two northern populations. Our results show populations remain interconnected and are structured broadly by the landscape and at regional scales by isolation by distance, or limited dispersal.
Altogether, this suggests C. venustus has an evolutionary history like that of ring species, and we discuss our results with the recommendation that others leverage morphological and genetic data across different spatial scales to tease apart population genetic processes from phenotypic patterns, as some may follow population structure while others may not.
Populations exhibit varying degrees of petal colour polymorphism in a ring-like distribution (Figure 1) suggesting the progression of sequential change as populations expanded north bidirectionally along the Central Valley. We find within-population-level polymorphism in C. venustus and define polymorphic traits as those which persist within populations, differentiated from plasticity by the lack of developmental or environmental influences. The occurrence of five polymorphic floral traits presents hundreds of combinations of variable floral displays, yet we recovered character states that commonly occur together (Figure S3) and may function similar to floral syndromes (Fenster, Armbruster, Wilson, Dudash & Thomson, 2004). On a background of multiple polymorphic traits, these frequently observed combinations of phenotypes could be due to covariation of reproductive traits under directional selection (Assis et al., 2016).
Remarkably, some individuals exhibited marginal picotee, although this had not been previously described in C. venustus but has been documented extensively in other plant lineages (Saito et al., 2006;Yuan, Rebocho, Sagawa, Stanley & Bradshaw, 2016). While P. L. Fiedler previously described the presence of distal petal spots (Baldwin et al., 2012), we found this trait varied geographically such that Coast Range populations often showed distal petal spots while eastern populations did not (Figure 3c, d). While we continue to directly assess the ecological functions of various C. venustus floral traits, evidence from numerous other systems shows that petal spots and patterns often serve in pollinator attraction (de Jager, Willis-Jones, Critchley & Glover, 2017;Johnson & Midgley, 1997) suggesting that pollination biology may be different between western and eastern populations. This spatial phenotypic pattern coincides with previous findings by R.K. Zebell (unpublished Figure 6, personal communication). We also observed differences across developmental time in petal, anther and filament pigmentation. Thus, careful consideration of developmental timing is required to minimize errors with phenotyping to distinguish between standing polymorphism and phenotypic changes or plasticity throughout plant development (Fusco & Minelli, 2010;Wennersten & Forsman, 2012).
The altitudinal trend in which deeply pigmented petals are found at higher altitudes could be explained by the potential for anthocyanin accumulation to protect various plant tissues from UV radiation (Landi, Tattini & Gould, 2015;Steyn, Wand, Holcroft & Jacobs, 2002). As plants are subject to higher UV exposure at higher elevations, pigmentation could potentially act in both pollinator attraction and protection of cells and tissues associated with reproduction (Weiss, 1991). Lessons on genetic bifurcation in space from applying ring species concepts
Our results show that populations have differentiated from one another across space in a ring-like distribution, regardless of floral phenotype. This is similar to how we think about speciation through time, with strong geographical cohesion shaped by limitations to dispersal, while interconnected habitat, in this case along mountain ranges, provides the foundation for continued gene flow (Figure 4 and Figure S13).
Several notable trends that emerge from this combination of data as seen through the lens of incipient speciation are discussed below.
The appearance of clear genetic clusters gave a first indication that genetic structure was based on geography rather than floral phenotypic traits. We predicted that if genetic differentiation did not directly underlie morphotypes, there would be within-population genotype assemblage for each phenotype, yet we discovered that distribution was still the primary factor of population structure. Evidence of high geographical cohesion, or that gene flow was stronger within populations than between them, is in line with the morphological consistencies found by R.K. Zebell: 'different color forms from the same site were more similar to each other in petal and gland shape than either was to the geographically nearest corresponding color form' (unpublished p. 97, personal communication). This indicates that while overall colour patterns are polymorphic within populations, there are more subtle phenotypic characters like petal and gland shape that correlate with underlying genetic structure.
Under a ring species model, the metapopulation remains interconnected as a chain of populations while genetic differentiation occurs spatially through isolation by distance with the strongest genetic divergence at terminal ends of the ring (Cacho & Baum, 2012). 2) and low number of private alleles in these populations (Table 1). Similarly, the north Coast Ranges and north Sierra Nevada individuals show the highest F ST values and least admixture as they have fewer neighbouring populations from which to directly exchange genes, yet there may be indirect gene flow from neighbouring populations and inherited alleles from historical migrations. This evolutionary history is consistent with traditional ring species concepts.
F ST revealed high gene flow between populations with low values among other California herbaceous plants (Wilson, Hipkins, Rey-Vizgirdas & Kaye, 2005;Wolf, Howe & Hamrick, 2000). Despite high gene flow, population structure was shaped by isolation by distance, and deviations from IBD aligned with major geographical landmarks.
We found IBD at a regional scale which aligns with the steppingstone model of spatial divergence in ring species and is corroborated by our phylogeographical results, haplotype network, ancestral area reconstruction and Fst values. Metapopulation IBD was low due to the nature of the distribution, that is, the two northern populations are genetically the most differentiated, but the distance between them is shorter than their distance to southern populations. Samples from the Coast Range traverse about 384 km and samples from the Sierra Nevada ca. 310 km, yet IBD is much stronger in the Coast Range region (Table 3) perhaps due to its relatively flat terrain compared to the Sierra Nevada. The variable strength of IBD across each major mountain range highlights that genetics like geography are heterogeneous across the landscape (Storfer et al., 2007). As revealed by EEMS, population splits align with breaks between and within major mountain ranges concordant with the PCA and phylogenetic analysis, although the underlying mechanisms remain to be tested. The divide between the northern and central Coast Range NCR CCR TR SSN CSN NSN NCR 0 -----CCR 0.0315 0 ----TR 0.0522 0.0288 0 ---SSN 0.0580 0.0384 0.0235 0 --CSN 0.0775 0.0604 0.0476 0.0331 0 -NSN 0.0888 0.0747 0.0616 0.0496 0.0345 0 Note: CCR, central Coast Ranges; CSN, central Sierra Nevada; NCR, north Coast Ranges; NSN, north Sierra Nevada; SSN, south Sierra Nevada; TR, Transverse Ranges.
TABLE 2 Pairwise population F ST values calculated from variant sites in Calochortus venustus located throughout California based on distance. Although gene flow is limited to nearest neighbours, populations remain interconnected along the arrangement of mountain ranges which provide the habitat connectivity necessary to maintain gene flow. Studies on Quercus spp. throughout California have found that populations remain interconnected (Sork et al., 2010), which is likely facilitated by continuity of woodland and riparian habitats as well as the dispersal mechanism of wind pollination, yet another study found that the Transverse Ranges limit gene flow (Ashley, Abraham, Backs & Koenig, 2015). We found genetic cohesion throughout the Transverse Ranges despite their complex terrain and wide altitudinal gradients. This pattern has also been found in non-plant organisms such as the flightless ground beetle (Chatzimanolis & Caterino, 2007). Patterns of regional genetic diversity contrast with morphological findings indicating that an integrative approach is needed to investigate adaptive potential. The north Coast Range population shows high values of genetic diversity unique from other populations in stark contrast to the population's low phenotypic diversity, namely predominantly white distal-spotted flowers. The negative F IS value indicates random mating spanning over 190 km, which in conjunction with the exceedingly high number of private alleles, higher levels of heterozygosity (Table 1) and short branch lengths on the maximum likelihood tree (Figure 5b) is consistent with rapid population expansion (Maruyama & Fuerst, 1984;Wright, 1949).
The Transverse Ranges host high morphological diversity while genetic diversity is not particularly high. Calochortus venustus' striking floral phenotypes in the Transverse Ranges may be driven in part by the same mechanisms that led this region to be a biodiversity hotspot for southern California (Baldwin, 2014;Thorne, Viers, Price & Stoms, 2009). Studies have concluded that habitat heterogeneity is of major importance, noting that climate, especially precipitation, and complex terrain and soil diversity are among the top driving factors (Baldwin, 2014;Raven & Axelrod, 1978;Thornhill, 2017). Contrasting patterns of genetic and morphological diversity in the north Coast Range and Transverse Range populations suggest these two populations may be of interest for their ecological and adaptive potential (Angert, Bradshaw & Schemske, 2008;Fraser & Bernatchez, 2001;Vandergast et al., 2008).
Summary statistics across all populations provide evidence that this system is not selfing; however, low levels of inbreeding may stem from their bulbous habit (Vallejo-Marín, Dorken & Barrett, 2010) as bulbs may reproduce vegetatively (i.e. produce clones) through the production of bulbils and bulblets.
The phylogenetic reconstruction with ancestral range reconstruction reveals historical migration from southern populations that bifurcated along two major parallel mountain ranges towards north (Figure 5b). Phenotypic data mapped onto the phylogeny indicate that there is parallel evolution towards white petals at both terminal northern populations in which the populations are nearly monomorphic for the white phenotype (Figures 1b, 5. The Transverse Ranges and southern Sierra Nevada region were found to be a historic corridor for northeastward expansion (Figure S13), likely through a single colonization into the Sierra Nevada with subsequent expansion and fragmentation supported by ancestral range reconstruction (Figure 5 and Figure S11) and decreasing genetic diversity towards northern latitudes (Table 1). Although we cannot date our tree due to lack of fossils or a well-resolved Calochortus phylogeny, the estimated age of the genus (Patterson & Givnish, 2004) tells us this is a relatively young lineage; thus, we can approximate that the recent uplift of major mountain ranges over the last 20 MYA in part shaped population structure in this species. While C. venustus is not a classic ring species, that is, the terminal populations are not in sympatry, they are the most similar in petal colour and therefore could potentially share pollinators and other agents of selection; yet as the most genetically divergent, they may be the most apt to be reproductively isolated if they were to reach sympatry.
Our results support C. venustus is an incipient ring species. We prescribe to the concept that phenotypes are not as important as evolutionary history when defining ring species (Kuchta & Wake, 2016), Note: n = total number of individuals tested at the metapopulation scale or within each mountain range, respectively.