We investigated parallel genomic divergence between urban and forest populations with two approaches. First, we examined polygenic divergence associated with urbanization by performing a local principal components analysis on outlier genomic regions. PCAs were implemented with the function snpgdsPCA in the R package “SNPRelate” (117 (link)) on each of the seven sets of outlier SNPs (urban GEA, urban-morphology). This analysis can provide insight into whether haplotypes are similarly diverging across urban–forest pairs (16 (link), 118 (link), 119 (link)). We then used a linear model to determine the effect of habitat (urban or forest), municipality (Arecibo, Mayagüez, and San Juan), and their interaction on the primary axes of genomic variation in the outlier sets (i.e., PC1 and PC2). A significant habitat effect would indicate that divergence associated with the urban environment or the trait (depending on the model) is associated with urbanization, whereas a significant municipality effect indicates regional variation driving divergence associated with the trait (e.g., as in ref. 16 (link)).
Second, we investigated parallel divergence at the allele level by examining effect sizes (eta2) of allele frequencies for all SNPs in our dataset (120 ). We used the etasquared function in the R package “rstatix.” We then compared the effect size of the habitat effect versus the interaction effect of habitat x municipality, where a stronger interaction effect suggests greater variation by region and the converse supporting parallelism. We compared effect sizes for all outlier SNPs identified in our two urbanization analyses (GEA, intersection of all three PCA) as well as the outlier SNPs identified by the intersection of the urbanization GEA and each morphology test (urban morphology SNPs). We compared effect sizes to the effect sizes of the background set of SNPs (SNPs not identified as outliers in any test).
Second, we investigated parallel divergence at the allele level by examining effect sizes (eta2) of allele frequencies for all SNPs in our dataset (120 ). We used the etasquared function in the R package “rstatix.” We then compared the effect size of the habitat effect versus the interaction effect of habitat x municipality, where a stronger interaction effect suggests greater variation by region and the converse supporting parallelism. We compared effect sizes for all outlier SNPs identified in our two urbanization analyses (GEA, intersection of all three PCA) as well as the outlier SNPs identified by the intersection of the urbanization GEA and each morphology test (urban morphology SNPs). We compared effect sizes to the effect sizes of the background set of SNPs (SNPs not identified as outliers in any test).
