Urbanization is rapidly changing the face of the planet. As cities develop, natural habitats 64 experience drastic environmental changes, from increased temperatures and pollution to greater 65 impervious surface and habitat fragmentation [1] . Evidence is accumulating to support the 66 hypothesis that the environmental features associated with urbanization drive phenotypic 67 differences between populations in urban and nonurban habitats [2] . For example, urban Anolis 68 lizards have evolved longer limbs and more toe lamellae to improve sprint speed on the smooth 69 artificial surfaces common in cities [3, 4] . Increased impervious surfaces often lead to warmer air 70 temperatures in cities [5] , which has driven increases in thermal tolerance of urban acorn ant [6] 71 and Daphnia populations [7] . These studies, and others (see table S1 in [8] ) suggest that many 72 urban populations are adapting to the novel environments created by humans. 73 Despite accumulating evidence of evolution in response to urban-driven environmental 74 change, most studies focusing on phenotypic divergence associated with urbanization have 75 examined just one or a few traits at a time [4, [9] [10] [11] [12] . However, theoretical and empirical work in 76 other systems suggest that selection can drive the evolution of multivariate phenotypic clines 77 along environmental gradients [13, 14] . Multivariate phenotypic analyses incorporate multiple, 78 often correlated, traits to explain overarching shifts in phenotypes across environments.
Although we can make informed predictions about how particular traits might respond to 80 urbanization from work in other systems, understanding how cities are driving phenotypic 81 evolution requires quantifying divergence in many traits that simultaneously influence fitness. It Here, we investigate multivariate phenotypic evolution along an urbanization gradient in 85 Toronto, Canada. Traits involved in plant reproduction are particularly likely to show 86 multivariate, genetically-based phenotypic associations with urbanization between urban and 87 nonurban populations due to the direct effect of these traits on fitness. In animal-pollinated 88 plants, urbanization might impact fitness through changes to their pollinator communities. Both 89 the diversity and abundance of pollinators are known to change along urbanization gradients 90 [15, 16] , with positive [17] [18] [19] to negative [20] effects on pollinator visitation. Variation in 91 pollination along urbanization gradients might drive changes in the extent of pollen limitation 92 experienced by outcrossing plants [21, 22] , which can influence the strength and opportunity for 93 pollinator-mediated selection on plant reproductive traits [23] [24] [25] . Consequently, we expect that (Table 1) [26, 27] .
Field observations 110 We quantified variation in pollination among our study populations to identify possible 111 mechanisms that might underlie phenotypic divergence along the urbanization gradient. We 
We additionally collected twenty ripe infructescences (i.e., group of fruit) from the same 117 populations in which we recorded pollinator observations. We counted the number of flowers 118 and seeds in each infructescence and used the number of seeds per flower as a metric of pollen 119 limitation (fewer seeds per flower = more pollen limitation). We generated data for the same 120 metric from plants in the common garden. These data provide information about in situ variation 121 in pollen limitation (field-collected inflorescences) and evolved differences in plants' abilities to 122 set seed from pollen (common garden data).
Data analysis 125 We investigated whether urbanization was associated with multivariate phenotypic 126 divergence using canonical redundancy analysis (i.e., RDA, [29] ). The RDA regressed a matrix be correlated with cyanogenesis frequency (rpearson= 0.63) in the same species and region [10] . 132 The RDA tests whether distance to the urban center explains more phenotypic variation 133 than expected by chance. In principle, we could observe a significant RDA even if just one trait 134 was associated with distance, but our interest lies in understanding whether multiple traits are 135 evolving together along the transect (i.e., including traits that do not vary significantly on their 136 own). Therefore, we used the canonical coefficients from the RDA that describe the individual 137 contribution of phenotypic traits to the first constrained axis of the RDA (RDA1 [i.e., distance]) 138 to calculate a multivariate phenotype score for each individual. This score, referred to as 139 clinemax, is the multivariate quantitative trait that shows the strongest association with distance to 140 the urban center [13] . Regressing clinemax against distance to the urban center enables us to better 141 visualize the shape of the multivariate cline and to quantify how multivariate phenotypes are 142 changing along our urbanization gradient. Urbanization was associated with multivariate trait divergence and the evolution of a 154 multivariate phenotypic cline in Toronto, ON ( Fig. S2 and Fig. 1b) . Distance to the urban core 155 explained 2.7% of the total variation in the multivariate phenotypic composition of populations 156 (RDA, F1, 202 = 5.51, P < 0.001, Fig. S2, Fig. S3 ). The composite trait showing the strongest 157 association with distance to the urban core (i.e., clinemax) showed a positive cline along the 158 urbanization gradient (β = 0.01, t206 = 9.32, P < 0.001, r 2 = 0.29, Fig. 1b ) that was at least as 159 strong as any trait when analyzed individually ( Fig. S4 and S5, see also standardized beta urban populations, whereas honeybee visitation was greatest in nonurban populations, and sweat 175 bees showed little change in visitation across the urbanization gradient (Fig. 2a) . 176 While common garden infructescences were on average more thoroughly pollinated than 177 those in natural populations (i.e., more seeds per flower; source F1, 50 = 30.04, P < 0.001, 178 meanfield = 1.1, meangarden = 2.5, Fig. 2b ), the number of seeds per flower was highest in urban based on their loadings onto RDA1 [13] , distance explained 29% of the variation in the 191 multivariate phenotype clinemax. Urban populations had later germination and flowering (Fig. 192 S4a and S4b), greater vegetative biomass (Fig. S4c) , larger banner petals (Fig. S4d) , thinner 193 stolons (Fig. S4e) , and lower HCN frequencies (Fig. S4f) . Many of the traits (three of six) most 194 strongly associated with distance to the urban core were reproductive traits, supporting the 195 prediction that these traits are likely to show evolutionary divergence between urban and 196 nonurban habitats, although some notable traits (e.g. reproductive biomass, number of flowers, 197 and number of inflorescences) showed no association with urbanization.
Urban evolution in white clover 200 Correlational data suggests that lower cyanogenesis in urban populations has evolved due to 201 selection for tolerance to freezing damage [10] . Recent data suggest a cost to producing the two 