Urbanization and a green corridor do not impact genetic divergence in common milkweed (Asclepias syriaca L.)

Urbanization is altering landscapes globally at an unprecedented rate. While ecological differences between urban and rural environments often promote phenotypic divergence among populations, it is unclear to what degree these trait differences arise from genetic divergence as opposed to phenotypic plasticity. Furthermore, little is known about how specific landscape elements, such as green corridors, impact genetic divergence in urban environments. We tested the hypotheses that: (1) urbanization, and (2) proximity to an urban green corridor influence genetic divergence in common milkweed (Asclepias syriaca) populations for phenotypic traits. Using seeds from 52 populations along three urban-to-rural subtransects in the Greater Toronto Area, Canada, one of which followed a green corridor, we grew ~ 1000 plants in a common garden setup and measured > 20 ecologically-important traits associated with plant defense/damage, reproduction, and growth over four years. We found significant heritable variation for nine traits within common milkweed populations and weak phenotypic divergence among populations. However, neither urbanization nor an urban green corridor influenced genetic divergence in individual traits or multivariate phenotype. These findings contrast with the expanding literature demonstrating that urbanization promotes rapid evolutionary change and offer preliminary insights into the eco-evolutionary role of green corridors in urban environments.


Study system
Common milkweed is an herbaceous perennial plant native to eastern North America.Although common milkweed grows in discrete patches of one to thousands of ramets (stems), often in abandoned agricultural fields 51 , urban populations tend to be smaller and inhabit public parks, railway and transmission rights-of-way, and roadsides, as well as private lawns and gardens.Plants can reproduce vegetatively through rhizomes that generate clonal ramets, or sexually through the cross-pollination of hermaphroditic flowers, each with five pollen sac pairs collectively called pollinaria [51][52][53] .Fertilization following successful pollination by insects, such as the western honey bee (Apis mellifera), bumblebees (Bombus spp.), and Halictidae spp.(Hymenoptera) in urban areas 49,54,55 , yields follicles (i.e., fruits) filled with wind-dispersed seeds 52,56 .Common milkweed is also recognized for its conservation significance because milkweeds (Asclepias spp.) provide essential resources for the threatened Monarch butterfly (Danaus plexippus) 57 .
Multiple traits protect common milkweed against herbivores.Milkweeds contain a pressurized, milky sap called latex that physically interferes with herbivores chewing tissue by gumming their mouthparts [58][59][60] .Herbivores that can overcome this barrier must also tolerate cardenolides, a suite of secondary metabolites which can disrupt sodium-potassium pumps (Na + /K + -ATPases) required for maintaining membrane potential 61 .At least 16 herbivores have coevolved with milkweeds and developed tolerance to these defenses 61,62 , and at least nine, including the Monarch butterfly, can enhance their own toxicities by sequestering cardenolides 62,63 .Compared to growth and reproduction, heritabilities of defense traits (e.g., latex and cardenolides) are often high [64][65][66][67][68] , indicating an increased likelihood for traits in this category to evolve.

Urbanization metrics
We quantified urbanization using two methods as described in Breitbart et al. 49 .Briefly, we first measured the distance from each population to the Toronto urban center (43.6563, − 79.3809), as a proxy for the degree of urbanization.Distance from the urban center is correlated with numerous environmental factors associated with urbanization 69 , and, for this city specifically, has been shown to influence the ecology and evolution of other plant, herbivore, and pollinator communities 10,70,71 .We also calculated an urbanization score for each population with the UrbanizationScore software [72][73][74] , which ranged from − 3.56 (least urban) to 3.37 (most urban).For each population, the software downloaded a 1 km radius aerial at 100 m resolution from Google Maps (Google LLC, Mountain View, California).The amount of vegetation, buildings, and roads was quantified and used to calculate landscape-cover variables, which were finally combined with principal component analysis into an urbanization score.Urbanization metrics were highly correlated (F (1,49) = 122.952,p < 0.001, R 2 adj = 0.709) (Supplementary Fig. S1).

Common garden experiment
To test if urbanization and an urban green corridor drive genetic divergence in phenotypic traits, we conducted a common garden experiment in 2019-2022.We germinated ~ 10 full-sibling seeds from each of 5 families per population, then grew seeds in pots in a growth chamber for six weeks.Seedlings were then snipped at the base of the ramet to ensure equal baseline heights and transported to the University of Toronto's Koffler Scientific

Trait measurements
In total, we measured 27 traits from 3 functional categories: plant defense/damage, growth/development, and reproduction, as well as herbivore abundance to contextualize defense trait expression (Table 1, Fig. 2).Measuring this large number of traits with diverse ecological functions allowed us to assess how the multivariate phenotype of common milkweed was evolving in response to urban environmental gradients.
We assessed five ecologically important traits related to plant defense/damage (Supplementary Text S1).We measured leaf herbivory by chewing herbivores before and after flowering by selecting the five oldest leaves pointing east (for consistency) and visually estimating the percent leaf area removed as described in Johnson et al. 75 .We quantified stem damage by a specialist milkweed stem weevil as the sum length of the oviposition scars per plant, which strongly predicts the number of eggs deposited into the ramet 76 .Latex exudation (the amount of latex released from plant tissue) was quantified by snipping the tip (~ 0.5 cm) off the youngest fully expanded intact leaf and collecting the latex exudate on a pre-weighed 1 cm filter paper disc as described in Agrawal et al. 64 .www.nature.com/scientificreports/Discs were placed in a pre-weighed microcentrifuge tube on dry ice, then transferred to a − 80 °C freezer until weighed to the nearest 1 µg.
We generated population-level estimates of leaf cardenolide concentrations by freeze-drying the leaf used to assess latex exudation and its opposite at − 80 °C, then finely grinding 50 mg of pooled tissue per population and following the protocol described in Petschenka et al. 77 .
In mid-July and mid-August 2020, we surveyed the abundance of ten specialist herbivore species located on all aboveground parts of the plants (Fig. 2).The total number of leaf splotch mines represented the abundance of the milkweed leaf-mining fly.Eggs, caterpillars, chrysalises, and butterflies (i.e., all life stages) contributed towards Monarch butterfly abundance.In 2021, we repeated the surveys of three species: the Monarch butterfly, the milkweed leaf-mining fly, and the swamp milkweed leaf beetle (Labidomera clivicollis).
We assessed eight traits associated with plant growth and development.We measured the height of all ramets before and after flowering began as the length of the ramet from soil level to the apical meristem, a surrogate of ramet biomass 49 .Plant mortality was recorded in early September, and we counted the number of ramets per pot before and after flowering began as each single shoot emerging from the soil or a shoot connected to another shoot at the soil surface.We calculated relative growth rate as: To assess specific leaf area (SLA) and leaf dry matter content (LDMC), we collected the youngest fully expanded intact leaf and measured wet mass, dry mass, length, and width.SLA and LDMC were calculated as leaf area/dry mass (m 2 /g) and dry mass/wet mass, respectively.
We assessed nine traits associated with plant reproduction.Flowering success represented whether a plant produced any fully reflexed flowers.We recorded flowering duration as the number of days between the dates of the first fully reflexed flowers on the oldest and youngest inflorescences (i.e., flower clusters).For the three oldest  www.nature.com/scientificreports/inflorescences per plant, we recorded the number of flowers per inflorescence once at least half the flowers in an inflorescence became fully reflexed.From three flowers per inflorescence, we assessed the number of pollinaria removed and flower size.Flower size was measured as the length and width of each flower's corolla (petal), hood (nectar-containing structure), and distance from one hood tip to the opposite side of the reproductive whorl (Supplementary Fig. S3).We counted the total number of inflorescences produced throughout the entire growing season.We recorded the date of the first mature follicle per plant when it exceeded 5 cm in length as smaller follicles are the most likely to be aborted 56 , as well as the total number of mature follicles per plant.

Statistical analyses
Quantitative genetic parameters All analyses were performed using R v.4.1.1 78 .We calculated three quantitative genetic parameters to understand how urbanization could affect divergence between populations and genetic variation within populations.To estimate the variances explained by population and family for each trait, we used the glmmTMB package 79 (version 1.1.2.2) to fit the following general linear models using restricted maximum likelihood: In these models, data was restricted to the last year of measurement to minimize impacts of maternal effects.Block was treated as a fixed effect while Population and Family were random effects, with Family nested within Population.
We extracted population and family-level variances with the "VarCorr" function from the lme4 package 80 (version 1.1.27.1), and residual variances with the "sigma" function from the stats package 78 (version 4.1.1).
We calculated estimates of full-sibling broad-sense heritability (H 2 ), the ratio of genetic variance to total phenotypic variance 50 , for each phenotypic trait as: where genetic variance for full-sibs was calculated as (family-level variance)× 2 and total phenotypic variance was calculated as family-level variance + population-level variance + residual variance.We calculated Q ST , a standardized measure of the genetic differentiation of a quantitative trait among populations 81 , as: We calculated the coefficient of genetic variation (CV G ), a dimensionless measure of evolutionary potential that is closely related to the coefficient of additive genetic variation (CV A ) and is useful for comparing the magnitude of genetic variation among traits 82 , as:

Genetic variation within and between populations (Q1)
To quantify the statistical significance of heritable genetic variation for phenotypic traits within and between populations, we refitted the models described above with the "lmer" and "glmer" functions from the lme4 package, for response variables with Gaussian and non-Gaussian distributions, respectively.This step was necessary because the original models fit with glmmTMB were not compatible with the "ranova" and "PBmodcomp" functions described below.
For models with Gaussian distributions, we tested the significance of Population and Family using the "ranova" function from the lmerTest package 83 (version 3.1.3),then divided p-values by 2 because these were 1-sided tests (i.e., variance ≥ 0) 84 .We then obtained percent variance explained (PVE) for each random effect after extracting variances with lme4 and calculating: For models with non-Gaussian distributions, we tested the significance of Population and Family using the "PBmodcomp" function from the pbkrtest package 85 (version 0.5.1) with 1000 simulations, then divided p-values, which were calculated assuming that the likelihood ratio test has a χ 2 distribution, by 2 because these were 1-sided tests (i.e., variance ≥ 0) 84 .Results for ten models analyzed with 1000 simulations were identical to those analyzed with 100, so we used 100 simulations for the remaining models.We refitted models to Gaussian distributions, then extracted variances and calculated PVE for each random effect as described above.We inspected model diagnostics with the DHARMa 86 (version 0.4.3) and performance 87 (version 0.7.3) packages and transformed response variables to meet the assumptions of normality and homogeneity of variance when necessary (Supplementary Tables S1-2).Response variables associated with plant herbivory and milkweed stem weevil damage were analyzed with hurdle models 88 manually to account for an excess of zeroes as a combination of two separate models: one to evaluate if herbivory/damage was present and the other to quantify the damage.To test whether the number of phenotypic traits with heritable genetic variation within and between populations was significantly different from the expected number of these traits, we performed binomial expansion tests 89,90 with the "binom.test"function from the stats 78 package (version 4.1.1).Because we set alpha to 0.05 for the previous hypothesis tests of individual traits, we set the probability of success to 0.05.The tests were onesided (alternative = "greater") to identify scenarios wherein the proportion of "successes" (traits with heritable genetic variation) out of all "trials" (all traits) was significantly higher than the expected proportion.Since this test functions under the assumption of independent trials, and there are inherent correlations among traits (Supplementary Figures S4-S5), we acknowledge that these correlations could inflate type I error and should be interpreted with this caveat.

Genetic clines along an urban-rural gradient (Q2)
To test if urbanization caused genetic divergence among populations, we fitted linear mixed effect models with two metrics of urbanization added as fixed effects using the same packages described above.For both general and generalized linear mixed effect models, urbanization was added to models as both Distance to the City Center and Urbanization Score, separately.To test the significance of random effects (i.e., Population & Family), we repeated the remaining steps as described for Question 1 above.If urbanization contributed to heritable genetic variation within and/or between populations (i.e., p-values for Question 1 models were significant for Family and/or Population), then these p-values should have been higher in Question 2 models' results.To test the significance of fixed effects (i.e., urbanization), we adjusted the models to use maximum likelihood and then computed χ 2 and p-values from a type II sums-of-squares ANOVA with the car package 91 (version 3.0.11).We used type II sums-of-squares because this method is unbiased by the order of effects, especially for unbalanced data, unlike type I.In contrast to type III, type II has more statistical power and is based on the assumption that interactions are minimal or absent 92 .If urbanization contributed to genetic variation within and/or between populations, then some variance should have shifted from Population and/or Family to urbanization and the p-values for urbanization should have become significant.
We performed this analysis for each trait, regardless of whether it showed a significant effect of population, to complement our initial test for population divergence.For the first test, we evaluated whether the variance among populations, treated as a categorical factor, was greater than zero.However, urbanization can also be treated as a quantitative metric in a regression framework to test for clines in traits.This second analysis allowed us to potentially tease out more subtle quantitative differences among populations that are harder to detect when populations are treated as categorical factors.Thus, the latter analysis is an important complement to the first even when there is no initial evidence of population differentiation.
To test for differences associated with urbanization in overall phenotype, as opposed to differences in individual traits that could reflect specific selection pressures, we performed a multivariate phenotype analysis.This analysis incorporated all traits in the dataset and accounted for non-independence between correlated traits (Supplementary Figures S4-S5).We computed best linear unbiased predictions (BLUPs) for each model with the "ranef " function from lme4, placed these values in a response matrix, and then used the mvabund package 93 (version 4.2.1) to fit general linear models and examine how multivariate phenotype varied with urbanization (both Distance and Urbanization Score) using one-way ANOVA.

Genetic differentiation between a green corridor and urban matrix (Q3)
To evaluate how genetic variation was associated with a green corridor between urban populations ("Urban: Corridor" and "Urban: Non-Corridor"), we added Subtransect and Urbanization:Subtransect interaction terms to the linear mixed effect models described for Question 2. For both general and generalized linear mixed effect models, we tested the significance of the random effects (i.e., Population & Family) by repeating the steps as described for Question 2 above.To test the significance of the fixed effects (i.e., Urbanization, Subtransect, and their interaction), we fitted reduced models without the interaction term, ranked full and reduced models based on Akaike's information criterion (AIC) 94 , and selected the model with the lowest AIC as the best model.We adjusted the models to use maximum likelihood and then computed χ 2 and p-values from a type III sums-ofsquares ANOVA if the best model contained an interaction with p ≤ 0.1 because this method tests for main effects after testing for interactions 92 ; otherwise, we reran the analyses with type II sums-of-squares.
We also fitted models with multiple years of data and found that the effects of urbanization and a green corridor were qualitatively identical to those reported in Tables 2, 3 in 85% of cases (Supplementary Text, Supplementary Table S3).We then repeated the multivariate analyses for these models as described for Question 2. Cardenolides were not analyzed due to there being a low sample size among urban populations.
We did not perform phenotypic or genotypic selection analyses because our experiment was only designed to quantify quantitative genetic variation.We could not obtain accurate estimates of variation in lifetime fitness because common milkweed is a long-lived perennial (members of this genus can live for over two decades 51 ).Moreover, since our common garden was located in a rural area, we could only measure selection in one environment.It would have been more appropriate to measure selection if we had comparisons among urban and rural environments, but this would have required planting multiple common gardens.

Genetic variation within and between populations (Q1)
There was heritable variation for multiple phenotypic traits within populations and evidence of genetic differentiation between populations for some traits.At least one trait per category exhibited significant heritable genetic variation within populations (nine traits, total).Heritabilities and coefficients of genetic variation ranged from 0-0.588 (mean ± SE: 0.110 ± 0.031) and 0-4.456 (mean ± SE: 0.467 ± 0.207), respectively (Table 1, Supplementary Table S4).The highest statistically significant heritabilities were observed for herbivory after flowering (binary) (H 2 = 0.588, p = 0.004), pollinaria removed (H 2 = 0.544, p = 0.034), and milkweed stem weevil damage (binary) (H 2 = 0.418, p = 0.030), while nine traits, at least one per category, exhibited near-zero heritabilities.Milkweed stem weevil damage (binary) (CV G = 0.995, p = 0.030), herbivory after flowering (binary) (CV G = 0.948, p = 0.004), and Monarch butterfly abundance (CV G = 0.767, p = 0.021) exhibited the highest, statistically significant coefficients of genetic variation, while seven traits exhibited near-zero values.It was unlikely that the high frequency of heritable genetic variation within populations was due to chance (binomial expansion test: p < 0.001).Three traits associated with plant defense and reproduction exhibited statistically significant genetic divergence among populations: latex exudation (Q ST = 0.174, p = 0.016), flowering success (i.e., whether plants flowered) (Q ST = 0.492, p = 0.014), and no. of inflorescences (Q ST = 1, p = 0.035) (Table 1, Supplementary Table S4).Q ST ranged from 0 to 1 (mean ± SE: 0.241 ± 0.068) with seven traits exhibiting near-zero values.The relatively few instances of genetic divergence between populations could be due to chance (binomial expansion test: p = 0.150).Thus, these results suggest moderate heritable genetic variation within populations and little phenotypic divergence among populations, the first of which is a prerequisite for adaptation.

Genetic clines along an urban-rural gradient (Q2)
We found little evidence for genetic clines along the urbanization gradient studied.When quantified as distance to the city center, urbanization did not significantly impact phenotypic traits (Fig. 3, Supplementary Figures S6-S14, Table 2, Supplementary Tables S6-S8).When quantified as an urbanization score, we detected relationships with latex exudation (χ 2 = 4.029, p = 0.045, R 2 m = 0.037, N = 699 individuals & 51 populations) and herbivory before flowering (quantitative) (χ 2 = 6.221, p = 0.013, R 2 m = 0.010, N = 430 individuals & 50 populations).When adjusted for false discovery rates (i.e., controlling for type I error with the Benjamini-Hochberg procedure), our results did not provide strong evidence that urbanization influenced genetic divergence in phenotypic traits between populations (Supplementary Tables S9-10).Moreover, the multivariate analysis indicated that urbanization did not impact overall phenotype when urbanization was quantified as distance to the city center (F 1,40 = 14.771, p = 0.711) or urbanization score (F 1,40 = 6.006, p = 0.984).Consistent with our results that urbanization had little effect on genetic divergence, the variation explained and statistical significance of the effects of population and genetic family did not substantially change once urban metrics were included in models (Supplementary Table S5).Taken together, these multiple lines of testing reveal little support for genetic divergence along an urbanization gradient in phenotypic traits in common milkweed.

Genetic differentiation between a green corridor and urban matrix (Q3)
Proximity to a green corridor did not strongly influence genetic divergence in phenotypic traits among populations, or the presence of clines along an urbanization gradient (Fig. 4, Supplementary Figures S15-S23, Table 3, Supplementary Tables S7 & S12).In the single-trait analysis, proximity to a green corridor had a statistically significant yet very small impact on flowers per inflorescence when urbanization was quantified as distance to the city center (χ 2    ).These interactions suggest that the impact of corridors may depend on the intensity of urbanization.However, when adjusted for false discovery rates, our results did not provide strong evidence that proximity to a green corridor influenced genetic divergence in phenotypic traits between populations as none of these effects were statistically significant (Supplementary Tables S13-14).Additionally, the multivariate analysis indicated that proximity to a green corridor did not impact overall phenotypic divergence when urbanization was quantified as distance to the city center (F 1,27 = 0.346, p = 1.000) or urbanization score (F 1,27 = 0.756, p = 1.000).Consistent with our results that proximity to a green corridor had little effect on genetic divergence, the variation explained and significance of the effects of population and genetic family did not substantially change once urban metrics and proximity to a green corridor were included in models (Supplementary Table S11).These results also reveal little support that proximity to a green corridor is associated with genetic divergence among urban common milkweed populations.

Discussion
In this study, we tested the hypotheses that urbanization and a single urban green corridor drive genetic divergence in phenotypic traits among populations of common milkweed.Though we observed moderate heritable genetic variation within populations and some weak phenotypic divergence among populations, we found low support that urbanization or proximity to an urban green corridor influenced genetic divergence.These results suggest that common milkweed has not undergone rapid evolutionary change in response to urban landscapes as a whole, nor a common component of urban environments: a green corridor.Despite these results, we maintain that our study system and experimental design are well-suited for addressing our research questions and that these findings are important for developing an accurate understanding of how urbanization impacts evolution.Below, we discuss the implications of these results for understanding how species evolve in response to rapid environmental change in heterogeneous urban environments.

Genetic variation within and between populations (Q1)
Genetic divergence of phenotypic traits along an urbanization gradient necessitates that populations exhibit genetic divergence, which is more likely to occur when those populations contain heritable phenotypic variation.These criteria have been observed in previous studies of common milkweed, although not in an urban context [64][65][66][67][68]95 . Forinstance, many phenotypic traits are heritable, especially those associated with defense/damage and growth; examples include Tetraopes spp.damage, carbon:nitrogen ratio, SLA, water percentage, trichome density, ant abundance, root and shoot mass, constitutive and induced latex and cardenolides [64][65][66][67][68]95 .In addition, substantial heritable phenotypic variation within populations was found at both a continental scale (i.e., within populations sampled across the species' native North American range 96 ) and a local scale (i.e., within single genetic populations 65,66 , but see Potts and Hunter 97 ).Similarly, genetic divergence for growth-related traits was detected among common milkweed populations sampled from a > 1500 km latitudinal gradient from New Brunswick, Canada, to North Carolina, USA 68 .In our common garden experiment, we found moderate heritable genetic variation within populations and weak phenotypic divergence among populations sampled from an urbanization gradient.Our results, which indicate that these populations have met the prerequisites to evolve in response to urbanization and/or an urban green corridor, are mostly consistent with previous findings from non-urban contexts 64,67 .Our estimates of genetic variation for phenotypic traits within populations are comparable to previous studies 64,67 .For example, heritability estimates for plant height after flowering and the number of ramets (range: H 2 = 0.132-0.161)were close to estimates by Vannette et al. 67 (aboveground biomass: H 2 = 0.11) and Agrawal et al. 64 (vegetative biomass: H 2 = 0.12).We found significant and relatively high heritability estimates for herbivory after flowering (H 2 = 0.588) and milkweed stem weevil damage (H 2 = 0.418) when measured on a presence/ absence basis, but not quantitatively when measured as percent leaf area removed.In comparison, Agrawal et al. 64 found moderate heritability for herbivory when measured as the percentage of leaves with foliar damage due to chewing herbivores (H 2 = 0.284), but not for milkweed stem weevil damage when measured as the length of stem scars (H 2 = 0.037).We also found that two reproductive traits had moderate to high heritabilities (i.e., pollinaria removed: H 2 = 0.544; date of first follicle: H 2 = 0.235), suggesting a higher evolvability of these traits in response to environmental change.Overall, these results confirm the capacity for these common milkweed populations to evolve in response to ecological disturbances, such as urbanization.

Genetic clines along an urban-rural gradient (Q2)
Many taxa exhibit genetic divergence between urban and nonurban populations for various traits.In plants, this has been documented for traits associated with phenology 16,27,29,30 , size 16,29 , fecundity 16,29 , defense 16 , and competitive ability 98 .Additionally, in both native and introduced ranges, common milkweed exhibits clines for growth and leaf physiology traits that correspond with a defining feature of urban environments: temperature 96 .Yet despite surveying several suites of traits in common milkweed, we found low support for genetic divergence along an urbanization gradient.Furthermore, we detected only small effect sizes (range: 0.01-0.037)for the few traits associated with urbanization even with the large scale and replication afforded by our experimental design.Thus, multiple lines of evidence suggest the lack of such divergence at present, though we do not, and cannot, rule out the possibility of genetic divergence emerging in the future.
The relative rarity of this outcome in the existing urban evolution literature presents a valuable opportunity to explore circumstances that could prevent urbanization from influencing genetic divergence among populations.In this case, evolutionary change may have been precluded by urban ecological pressures that were possibly too small or brief, as the city of Toronto has contained ≥ 50,000 residents for only ca. 150 years 99 .The life history traits of common milkweed could also slow evolutionary change.For example, the vegetative reproduction inherent to common milkweed can lead to clonal growth and the loss of genotypic diversity within populations 100,101 , which is compounded by the species' long-lived nature.Long-distance pollen and seed flow could yield high gene flow among populations and low genetic drift within populations, as common milkweed is largely self-incompatible (i.e., self-pollinated plants only produce viable fruit at very low rates 102 ) , seeds are wind-dispersed, and flowers are frequently pollinated by insects that can travel > 1 km (e.g., Bombus spp., Apis mellifera) 103,104 .
We acknowledge that urban-nonurban environmental gradients are complex and that our data cannot fully capture the heterogeneity of the landscape mosaic.Urbanization is a multifaceted process that involves complex change along axes including, but not limited to, environmental, economic, and sociological dimensions.There are manifold ways to define "urban" vs. "rural" vs. "natural" landscapes 105 , and environmental heterogeneity functions across temporal and spatial scales 106,107 .However, our methods account for both temporal and spatial processes.Distance from the city center is highly correlated with percent impervious surface and numerous other environmental variables associated with urbanization 69 , and in Toronto reflects urban expansion in concentric zones, with each zone's environment correlated with the time since development.By including urban and rural environments as transect endpoints and densely-sampled transitional areas in between, distance from the urban center captures extensive multivariate environmental change associated with urbanization-an important feature, since we did not know which environmental factors might impact common milkweed a priori.A complementary metric is urbanization score, which quantifies the immediate land use/land cover surrounding each population, without consideration of the site's proximity to the urban core.Therefore, despite this limitation, we argue that our estimates of urbanization are credible and reliable.
This study lends valuable context to prior work demonstrating how urbanization impacts the reproductive success of common milkweed populations 49 .Specifically, as we observed relatively little phenotypic divergence among populations within the common garden, our results suggest that the previously observed in-situ phenotypic divergence was mostly consistent with phenotypic plasticity as opposed to genetic divergence.This finding underscores the importance of investigating the genetic basis of phenotypic divergence observed in cities and reiterates the role of phenotypic plasticity in shaping how populations respond to human-caused change.Future studies could investigate the specific conditions that promote phenotypic plasticity within urban environments.

Genetic differentiation between a green corridor and urban matrix (Q3)
Very little is known about how urban green corridors influence genetic divergence among plant populations.In non-urban environments, green corridors are predicted to increase gene flow among plant populations by facilitating pollen flow and seed dispersal, which is expected to decrease genetic divergence between populations [108][109][110] .Limited research in urban environments suggests that green corridors often increase gene flow among animal populations 46,[111][112][113] (but see Angold et al. 114 ), and that urban features such as railways can function as corridors among plant and animal populations 115,116 .In our study, proximity to a green corridor did not strongly influence genetic divergence in phenotypic traits among common milkweed populations.Despite finding low support for genetic divergence along an urbanization gradient in Question 2, we proceeded with this analysis because we believe that once research questions or hypotheses have been decided, it is important to address each, so as not to change the hypothesis based on the results.
As discussed above for Question 2, urbanization may not restrict gene flow among populations.If true, there would be little opportunity for the green corridor to significantly enhance gene flow and diminish any hypothetical genetic divergence in common milkweed, as proposed for butterflies in Angold et al. 114 .Conversely, our data also suggest that proximity to a green corridor does not inherently facilitate genetic divergence either.It is also possible that urban green corridors actually impact genetic divergence but that our results reflect our inclusion of a single corridor, or the specific environmental conditions associated with our chosen corridor.For example, the efficacy of the corridor could have been impeded by a highly hostile external matrix or edge effects generated by the narrow shape 117,118 .Relatedly, the corridor we sampled would have provided minimal connectivity to plants or pollen vectors if it were not perceived as a functional corridor to either group-the plants or the pollinators.The manner in which common milkweed established in this particular corridor also influences our results, as plants growing in the corridor since establishment many generations ago are more likely to exhibit genetic divergence than plants descending from relatives in nearby non-corridor environments.Regardless, www.nature.com/scientificreports/these findings contradict the expectation that habitat fragmentation impacts the evolutionary processes of urban populations and invite further research into how and when corridors impact genetic divergence of phenotypic traits in urban populations.

Conclusion
Urbanization is associated with phenotypic trait divergence for many species, yet in most cases, the genetic basis of these trait differences is unresolved.Additionally, the extent to which specific aspects of the urban landscape impact genetic divergence in phenotypic traits among plant populations is virtually unknown.Here, we show that neither urbanization nor an urban green corridor impacted genetic divergence in phenotypic traits in common milkweed, a native plant of conservation importance.These results demonstrate an example in which our measures of urbanization have not substantially contributed to genetic divergence among populations.To our knowledge, this study is also the first to investigate if urban green corridors impact genetic divergence in phenotypic traits in plant populations and ultimately suggests the absence of such effects in our study system.
To further understand how urban environments and urban green corridors impact eco-evolutionary dynamics in plant populations, future research should verify the consistency of these findings by testing similar questions in plants with diverse life history, and in a range of cities with heterogeneous landscapes.

Figure 1 .
Figure 1.Map of 52 common milkweed populations sampled along Toronto's urban-rural gradient and location of common garden experiment at the Koffler Scientific Reserve (green star).Urban: Non-Corridor populations (squares; N = 17).Urban: Corridor populations (triangles; N = 19).Rural populations (circles; N = 16).The color of the symbols indicates urbanization score, where positive values indicate a high degree of urbanization (based on the quantity of vegetation, buildings, and paved roads per 1 km 2 ).The Stamen terrain basemap shows urban and suburban areas in light gray, nonurban agricultural and forested areas in green, and Lake Ontario in blue.Map tiles by Stamen Design, under CC BY 3.0.Data by OpenStreetMap, under ODbL.
RGR =log height after flowering − log height before flowering Days between measurements

Figure 2 .
Figure 2. Schedule for measurement of traits associated with plant reproduction, growth, and defense/damage throughout 2019-2022 field seasons and the number of plants alive per year.Grey boxes indicate years of measurement.

Figure 3 .
Figure 3.The effect of urbanization on representative traits from each main category: (a) plant defense/ damage (e.g., the presence of milkweed stem weevil damage); (b) plant reproduction (e.g., date of first follicle); (c) herbivore abundance (e.g., milkweed leaf-mining fly abundance); and (d) plant growth (e.g., height before flowering) when urbanization was quantified by distance from the urban center and all populations were included.Regression lines represent predicted values with a 95% confidence interval and points which represent family-level means are shown for general and generalized linear mixed effects models.These traits were modelled using the following distributions: binomial (a), negative binomial (b, c), and Gaussian (d).

Figure 4 .
Figure 4.The effects of urbanization and proximity to a green corridor on representative traits from each main category: (a) plant defense/damage (e.g., the presence of milkweed stem weevil damage); (b) plant reproduction (e.g., date of first follicle); (c) herbivore abundance (e.g., milkweed leaf-mining fly abundance); and (d) plant growth (e.g., height before flowering) when urbanization was quantified by distance from the urban center and only urban populations were included.Regression lines represent predicted values with a 95% confidence interval and points which represent family-level means are shown for general and generalized linear mixed effects models.These traits were modelled using the following distributions: binomial (a), negative binomial (b, c), and Gaussian (d). https://doi.org/10.1038/s41598-023-47524-8

Table 2 .
Results from general and generalized linear mixed models examining the effects of urbanization on all phenotypic traits.All populations were included.Shown are maximum likelihood χ 2 and p-values obtained from type II sums-of-squares ANOVA.Though not shown, block was included as a fixed effect and often explained significant variation in the common garden.Significant p-values (p ≤ 0.05) are bolded.Relative growth rate was only analyzed for plants exhibiting positive growth between measurements.

Table 3 .
Results from general and generalized linear mixed models examining the effects of urbanization and proximity to the urban green corridor on all phenotypic traits.Urbanization was quantified via distance from the urban center and only urban populations were included.Shown are maximum likelihood χ 2 and p-values obtained from type II and III sums-of-squares ANOVA.Though not shown, block was included as a fixed effect and often explained significant variation in the common garden.Significant p-values (p ≤ 0.05) are bolded.Relative growth rate was only analyzed for plants exhibiting positive growth between measurements.