Intraspecific independent evolution of floral spur length in response to local flower visitor size in Japanese Aquilegia in different mountain regions

Abstract Geographic differences in floral traits may reflect geographic differences in effective pollinator assemblages. Independent local adaptation to pollinator assemblages in multiple regions would be expected to cause parallel floral trait evolution, although sufficient evidence for this is still lacking. Knowing the intraspecific evolutionary history of floral traits will reveal events that occur in the early stages of trait diversification. In this study, we investigated the relationship between flower spur length and pollinator size in 16 populations of Aquilegia buergeriana var. buergeriana distributed in four mountain regions in the Japanese Alps. We also examined the genetic relationship between yellow‐ and red‐flowered individuals, to see if color differences caused genetic differentiation by pollinator isolation. Genetic relationships among 16 populations were analyzed based on genome‐wide single‐nucleotide polymorphisms. Even among populations within the same mountain region, pollinator size varied widely, and the average spur length of A. buergeriana var. buergeriana in each population was strongly related to the average visitor size of that population. Genetic relatedness between populations was not related to the similarity of spur length between populations; rather, it was related to the geographic proximity of populations in each mountain region. Our results indicate that spur length in each population evolved independently of the population genetic structure but in parallel in response to local flower visitor size in different mountain regions. Further, yellow‐ and red‐flowered individuals of A. buergeriana var. buergeriana were not genetically differentiated. Unlike other Aquilegia species in Europe and America visited by hummingbirds and hawkmoths, the Japanese Aquilegia species is consistently visited by bumblebees. As a result, genetic isolation by flower color may not have occurred.

Local adaptation of floral traits to pollinators may have occurred across multiple regions, but there is little evidence as to whether variation in floral traits has occurred independently among regional populations (but see Anderson et al., 2014;Toji et al., 2021). As in textbook examples of ecological speciation (Nosil, 2012), one useful approach to understand the interaction between trait diversification and speciation in angiosperms is to combine a field analysis of local plant evolutionary adaptations with a population genetic analysis that examines genetic relationships among populations. This approach will allow us to unravel the historical patterns of floral trait evolution. Local adaptation of floral traits to pollinators may lead to speciation via the establishment of prezygotic reproductive isolation (Grant-Stebbins model;Grant & Grant;Anderson et al., 2014;Johnson & Anderson, 2010;Stebbins, 1970) because one possible result of specialization of a trait to a particular pollinator is a lack of pollinator sharing among plant populations (Anderson & Johnson, 2008;Herrera et al., 2006;Newman et al., 2015). About 25% of plant diversification events may be associated with pollinator shifts (Van der Niet & Johnson, 2012); thus, combined analyses of local adaptation of floral traits and population genetics can shed light on the mechanisms of plant diversity (Thompson, 2005;Thompson et al., 2013). Furthermore, regional pollinator difference is the topic that is closely related to speciation and morphological diversification, yet as mentioned above, we have very limited knowledge about historical patterns of floral trait evolution. In this study, we hypothesized that the Aquilegia plants in each mountain region have their respective common ancestor, which would be confirmed by population genetic analysis, and their flower size independently evolved in each population as an adaptation to the local pollinator size. The purpose of this study is to provide a further perspective to understand convergent evolution related to regional pollinator differences.
In genus Aquilegia (Ranunculaceae), adaptive radiation to different pollinators (bumblebees, hummingbirds, and hawkmoths) has occurred. Mainly, flower color, spur length, flower orientation, and pistil length have evolved to adapt to each pollinator (Fulton & Hodges, 1999;Hodges et al., 2004). Moreover, molecular phylogenetic evidence also indicates that pollinator shifts have led to morphological diversification and speciation within this genus (Whittall & Hodges, 2007). According to Whittall and Hodges (2007), a more ancestral floral state of Aquilegia is purple, downward facing, shortspurred flowers, which are pollinated by bumblebees. From plants with this floral state, taxa with red, downward facing flowers with protruding stamens and intermediate length spurs, which are pollinated by hummingbirds, were derived. Then, taxa with white and yellow long-spurred, lateral and upward facing flowers, which are pollinated by hawkmoths, were derived from those taxa. Their results reveal an interesting pattern of species-level diversification as a consequence of pollinator shifts, although evidence for flower trait diversification at the earlier stages of speciation is lacking. Aquilegia is a great model for pollinator-driven speciation, but to observe early stages of speciation, it is useful to investigate the pattern of evolutionary morphological diversification within a single species (Sobel & Streisfeld, 2015).
In this study, we focused on evolutionary processes leading to spur length and flower color differentiation in Aquilegia buergeriana var. buergeriana. In this species, geographic variation in spur length has previously been observed in six populations in Utsukushigahara and Norikura regions, but the relationship between spur length and flower visitors in these populations is not known . Yellow-flowered individuals are dominant in this species, and bumblebees seem to be the main flower visitors. In some populations, red-flowered individuals occur orthotopically with yellowflowered individuals, but the genetic relationship between red-and yellow-flowered individuals is unknown. Differences in flower color in Aquilegia can lead to genetic isolation even between neighboring or sympatric populations and is likely to be important in speciation (Hopkins & Rausher, 2012;Schemske & Bradshaw, 1999). Our study system is the best one to describe the early genetic isolation of the plant because there is variation not only in spur length but also in flower color.
Here, we first investigated the correspondence between variation in floral spur length and flower-visiting insect size in 16 Aquilegia populations in four mountain regions. The results showed a morphological correlation between spur length and average visitor size in each population, even within the same mountain region; spur lengths were shorter in populations visited by smaller flower visitors, and spur lengths were longer in populations visited by larger flower visitors. Next, we identified genome-wide single-nucleotide polymorphisms (SNPs) by the MIG-seq (multiplexed inter-simple sequence repeat genotyping by sequencing) method (Suyama & Matsuki, 2015) to clarify the genetic relationships among the populations. These results showed that genetic relationships tended to be clustered by mountain region and, therefore, that spur length evolved in parallel in each mountain region. Individuals with different flower colors were not differentiated in the genome-wide SNPs analysis, however. This result suggests that pollinator isolation by flower color has not occurred in these populations. Instead, the red flower color is maintained in various populations in which most individuals have yellow flowers.

| Plant species and study site
Aquilegia buergeriana var. buergeriana f. flavescens is a perennial, protandrous herb that grows along forest edges in mountain area, endemic to Japan except for Okinawa (Hayashi, 2009). The spur and sepals of its flowers are pale yellow (yellow-flowered individual) or reddish brown (red-flowered individuals) (Figure 1a,b).
Flowers of both colors face downward. In the study area, in the central Japanese Alps, yellow-flowered individuals are more common. Japanese Aquilegia species are mainly visited by bumblebees Itagaki & Sakai, 2006;Tamura & Shimizu, 1999), even though, in general, yellow-flowered Aquilegia are pollinated by hawkmoths (Hodges et al., 2004). Unlike most Aquilegia with yellow flowers, however, the yellow flowers of Japanese A. buergeriana do not have protruding anthers and pistils and are not visited by hawkmoths (Toji's personal observation by camera trap for about 3 days in UT-1640 population).
We studied Aquilegia populations in four mountain regions (Utsukushigahara, Norikura, Ontake, and Iizuna) of the central Japanese Alps (Figure 1d; Table S1).  width did not differ among the five petals of each flower. We considered spur length to be the most important trait because of its relation to visitor size. Therefore, in subsequent analyses we focused on spur length. Since the normality of the data was not supported in some populations, we performed a non-parametric multiplecomparison Steel-Dwass test to compare average spur length between populations.
We also examined spatial autocorrelation (i.e., whether the variation in spur length could be explained by physical distance) by using the "moran.test" function in the "spdep" package in the R Software Environment ver. 4.0.2 (R Core Team, 2013) to run Moran's I test.
This analysis used the average spur length and the latitude and longitude of each population.

| Flower visitor assemblages and size variation
To

| Factors influencing local spur length
The Before conducting the LMM analysis, the variance inflation factor (VIF) statistic was calculated to check for correlation (multicollinearity) between variables, using VIF = 0.5 as the threshold (Neter et al., 1996). For all variables, VIF was less than 0.25, confirming that no multicollinearity existed. Next, we conducted a likelihood ratio test using the parametric bootstrap method (Hoel et al., 1971)  Then, we adopted the model with the lowest Akaike information criterion (AIC). We excluded the nine ON-1000 population, which was not visited by bumblebees, from this analysis, treating it as a missing value.

| Genetic structure of A. buergeriana var. buergeriana populations
Leaf samples were obtained from 6-21 individuals in each study population (a total of 232), and DNA was extracted by the CTAB method (Doyle & Doyle, 1990). We performed MIG-seq (Suyama & Matsuki, 2015) to detect genome-wide SNPs. A MIG-seq library was prepared following the protocol of Suyama et al. (2022). A first PCR was performed to amplify inter-simple sequence repeat regions using MIG-seq primer set 1 (Suyama & Matsuki, 2015), and then a second PCR was performed on the purified/equalized first PCR product to add the index sequences necessary for sequencing on Low-quality reads and extremely short reads containing adaptor sequences were removed by using trimmomatic 0.39 (Bolger et al., 2014). De novo SNP discovery was performed by using the Stacks 2.41 software pipeline (Catchen et al., 2013;Rochette et al., 2019).
For de novo SNP discovery, we used the following parameters: minimum depth of coverage required to create a stack (m) = 3, maximum distance between stacks (M) = 2, and maximum mismatches between loci when building the catalog (n) = 2. Three different filtering criteria were applied for quality control of the SNP data. First, SNPs that were retained by 80% or more samples were included in the SNP dataset. Second, SNPs with a minor allele frequency of less than 0.05 were removed. Third, loci containing SNPs with extremely high observed heterozygosity (Ho ≥ 0.6) were removed. Fourth, after performing a Hardy-Weinberg equilibrium test on each population, we excluded loci where allele frequencies deviated from the Hardy-Weinberg equilibrium at p < .01 in three or more populations.
The following population genetic statistics of SNP sites for each population were calculated with the Stacks populations module: expected heterozygosity He, observed heterozygosity Ho, nucleotide diversity π, and inbreeding coefficient F IS (Hartl & Clark, 1997). The population genetic structure was examined by PCA using PLINK 1.9 (Purcell et al., 2007). In addition, a Bayesian clustering analysis was performed with STRUCTURE software version 2.3.4 (Falush et al., 2003;Pritchard et al., 2000). For the STRUCTURE analysis, simulations were performed with 100k burn-in iterations and 100k Markov chain Monte Carlo iterations. The number of genetic clusters (K) was calculated 10 times for each possible K value from 1 to 10, and the appropriate number of clusters was estimated based on the ΔK value (Evanno et al., 2005). Then, to examine the genetic structure within each mountain region in more detail, we performed additional STRUCTURE analyses. First, SNP re-detection was performed in each of three mountain regions, the Utsukushigahara, Norikura+Ontake, and Iizuna regions, based on results of the initial analysis. The population structure was obtained based on all samples, with the above filtering criteria used in SNP detection.
Second, 10 independent STRUCTURE analysis runs were performed for each mountain region with 100,000 burn-in steps and an additional 100,000 steps with the admixture model; log-likelihood values were estimated for each possible K value (K = 1-10), and the appropriate number of clusters was estimated based on the ΔK.

| Isolation by distance and isolation by phenotype
We investigated whether the genetic structure of A. buergeriana var.
buergeriana reflects geographic distance or spur length differences.
In general, populations separated by greater distances are more genetically differentiated than populations close together (Wright, 1943). On the other hand, if populations with similar traits are also genetically similar, then we can expect to find a correlation between differences in traits between populations and the degree of genetic differentiation. We used GenoDive software version 3.0 (Meirmans, 2020) to calculate the genetic isolation (F ST ) between populations.
The geographic distance between populations was calculated from the latitude and longitude of the populations, and the difference in the average spur length of each population was used as the trait difference. We calculated the relationship between pairwise F ST or F ST / (1 -F ST ) and geographic distance between populations, as well as the relationship between pairwise F ST or F ST /(1 -F ST ) and trait difference between populations, following methods in Rousset (1997) and Noutsos et al. (2014). The relationship between genetic isolation and geographic or trait distance was tested by Mantel tests using the R package "ade4" with 10,000 Monte Carlo permutations.

| Spur length and flower visitor size
The average spur length of each population of A. buergeriana var.

| Factors influencing local spur length
The LMM model with the lowest AIC value was the model that included only average visitor size (only bumblebees) as an explanatory variable (Table 1; Table S3). A very strong linear relationship was found between average spur length and average visitor size (only bumblebees) (p < .0001, Figure 2).

| Isolation by distance and isolation by phenotype
A significant relationship between geographic distance and genetic isolation (F ST , F ST /(1 − F ST )) was detected (Table 2). On the other hand, differences in average spur length between populations were not related to genetic isolation.

| Intraspecific independent evolution of spur length among mountain regions
The spur length of A. buergeriana var. buergeriana was correlated with the average flower visitor size (only bumblebees) of the population; spur length varied with the average visitor size even among spatially close populations in the same mountain region (Figures 2 and 4a).
The genetic results obtained by PCA and the STRUCTURE analysis of genome-wide SNPs suggest that populations within each mountain region are more closely related to each other than to populations in other mountain regions (Figures 3 and 4; Figure S4). Genetic isolation was proportional to geographical distance and did not reflect spur length differences (Table 2). These results suggest that genetic origin is derived from each mountain region, and spur length of A. buergeriana var. buergeriana evolved independently in each mountain region.
In the Iizuna region, populations at different altitudes seem to belong to different genetic clusters ( Figure 4d); furthermore, flowers in lower-altitude populations were visited by B. diversus and those in higher-altitude populations were visited by B. consobrinus (Table   S1). These results suggest that genetic differentiation may occur between higher-and lower-altitude populations because of a lack of pollinator sharing. Further studies are needed to determine whether gene flow by pollination is hindered between populations at higher and lower altitudes in the Iizuna region.
In anole lizards, leg length has evolved independently on different islands to suit local habitats (Losos, 2010), and in stickleback

TA B L E 1
The GLM model that best explained variation in average spur length among populations of Aquilegia buergeriana var. buergeriana fishes, the evolution of marine to freshwater forms (sticklebacks that move between rivers and the sea) occurred independently in different marine and freshwater locations in various regions of the world (Jones et al., 2012). We propose that plant species distributed across a wide geographic range with site-specific, different-sized pollinators constitute another model suitable for testing independent adaptive radiation. We have demonstrated that spur length in an Aquilegia species may have evolved independently among mountain regions by using a population genetics approach to compare traits among mountain regions. Independent evolution in different mountain regions has recently been examined in various model systems: for example, the independent evolution of upland and short-winged forms of scorpionfly Panorpodes (Panorpodidae) (Suzuki et al., 2019), the independent evolution of Potentilla matsumurae (Rosaceae) in fellfield and snowbed environments on different mountains in Japan (Hirao et al., 2019), and the independent evolution of alpine morphology in Antirrhinum species (Antirrhineae) (Durán-Castillo et al., 2021). Furthermore, we recently presented a case in which we used microsatellite markers to show the independent adaptation of floral tube size in Lamium album var. barbatum (Lamiaceae), associated with flower visitor size, in the Utsukushigahara and Norikura regions of the Japanese Alps (Toji et al., 2021). These examples show that comparisons between mountain regions can be used to study independent trait evolution in various organisms, and similar patterns might be found in many places around the world.

| Flower color does not contribute to genetic isolation
Although red-flowered individuals were observed in some popula- wide SNP data has shown that they are genetically distinct (Sobel & Streisfeld, 2015). In the hybrid zone between the two ecotypes, the MaMyb2 gene, which is involved in the synthesis of flower pigments, is geographically maintained despite neutral gene flow occurred (Sobel & Streisfeld, 2015;Stankowski & Streisfeld, 2015;Streisfeld & Kohn, 2005). Gene flow between yellow and red flower M. aurantiacus ecotypes in the early stages of speciation seems to be limited mainly by differences in pollinator preference (Sobel & Streisfeld, 2015). Similarly, gene flow between two closely related Aquilegia species: hummingbird-pollinated, red-flowered A. formosa and hawkmoth-pollinated, yellow-flowered A. pubescens are also limited by pollinator isolation when the two species are distributed parapatrically (Fulton & Hodges, 1999;Noutsos et al., 2014).
Why yellow-and red-flowered individuals in A. buergeriana var.
buergeriana may not become genetically isolated? In the central Nagano region, where this study was conducted, bumblebees appear to be abundant and many flowers depend on bumblebees for pollination (e.g., Egawa & Itino, 2020), whereas potential pollinators such as birds and butterflies that prefer red flowers are scarce. In another Japanese mountain region (the Taisetsu mountains), flowers are dominantly visited by bees and flies at the community level (Mizunaga & Kudo, 2017). A recent review has reported that Lepidoptera account for less than 10% of insect visitors to flowers in many parts of Asia, whereas bees and flies account for more than half (Funamoto, 2019

| CON CLUS IONS
Two main conclusions follow from our results that (1)

F I G U R E 4
Bayesian genetic clustering analysis STRUCTURE results based on SNPs data. In each case, the appropriate value of K was determined from ΔK (see Figure S3). Population numbers, corresponding to those in Figure  independently in different mountain regions, and (2) the few redflowered phenotypes that occur within the species may not lead to genetic differentiation. First, given that the independent evolution of floral size in different mountain regions has also recently been reported in L. album var. barbatum (Toji et al., 2021), the independent evolution of floral size among mountain regions may be a generalized event that occurs commonly in different taxa. The approach used here to test for independent evolution among mountain regions is applicable to any taxon and a variety of traits. In particular, morphological analyses combined with MIG-seq (Suyama & Matsuki, 2015), which can be used to obtain genome-wide SNPs from non-model organisms, constitute a powerful method for elucidating patterns of morphological and genetic diversification within species. Second, we found no relationship between flower color and the degree of genetic differentiation, despite the fact that pollinator isolation caused by differences in flower color has been reported in two closely related species of Aquilegia (Fulton & Hodges, 1999;Noutsos et al., 2014). We infer that in the mountainous region of Japan, where bumblebees are locally abundant large pollinators, shifts to different pollinator taxa are unlikely to occur, and the polymorphism in A. buergeriana var. buergeriana flower color is likely maintained by random genetic drift. Thus, our results are an important exception to diversification in genus Aquilegia, which is well known to have occurred by both flower color and pollinator shifts (Whittall & Hodges, 2007).

ACK N OWLED G EM ENTS
We thank S. Duhon for English editing. We thank the Chubu District Forest Office (Forestry Agency), the Chubu Regional Office for Nature Conservation (Ministry of the Environment), and the Matsumoto Regional Office (Nagano Prefectural Government) for permission to work in the area.

CO N FLI C T O F I NTE R E S T
There are no conflicts of interest.