Urban colonization through multiple genetic lenses: The city‐fox phenomenon revisited

Abstract Urbanization is driving environmental change on a global scale, creating novel environments for wildlife to colonize. Through a combination of stochastic and selective processes, urbanization is also driving evolutionary change. For instance, difficulty in traversing human‐modified landscapes may isolate newly established populations from rural sources, while novel selective pressures, such as altered disease risk, toxicant exposure, and light pollution, may further diverge populations through local adaptation. Assessing the evolutionary consequences of urban colonization and the processes underlying them is a principle aim of urban evolutionary ecology. In the present study, we revisited the genetic effects of urbanization on red foxes (Vulpes vulpes) that colonized Zurich, Switzerland. Through use of genome‐wide single nucleotide polymorphisms and microsatellite markers linked to the major histocompatibility complex (MHC), we expanded upon a previous neutral microsatellite study to assess population structure, characterize patterns of genetic diversity, and detect outliers associated with urbanization. Our results indicated the presence of one large evolutionary cluster, with substructure evident between geographic sampling areas. In urban foxes, we observed patterns of neutral and functional diversity consistent with founder events and reported increased differentiation between populations separated by natural and anthropogenic barriers. We additionally reported evidence of selection acting on MHC‐linked markers and identified outlier loci with putative gene functions related to energy metabolism, behavior, and immunity. We concluded that demographic processes primarily drove patterns of diversity, with outlier tests providing preliminary evidence of possible urban adaptation. This study contributes to our overall understanding of urban colonization ecology and emphasizes the value of combining datasets when examining evolutionary change in an increasingly urban world.

Urban colonization further constitutes a founder event, whereby decreased genetic variation may follow reduced effective population size (Greenbaum, Templeton, Zarmi, & Bar-David, 2014;Nei, Maruyama, & Chakraborty, 1975). This phenomenon is observed during colonization of natural (Fabbri et al., 2007) and urban (DeCandia et al., 2019;Evans et al., 2009) areas, but may be exacerbated by the difficulty of navigating the urban matrix. If effective population size remains small through colonization and establishment, the resultant diversity loss may threaten long-term viability due to increased risk of inbreeding (Frankham, 2005(Frankham, , 2008Oakley, 2013;Spielman, Brook, Briscoe, & Frankham, 2004). Deleterious alleles may accumulate and decrease the fitness of urban colonizers. This diversity loss is especially problematic for threatened or sensitive species colonizing environments characterized by novel stressors, such as cities, where local adaptation may be necessary for persistence.
Despite decreases in effective population size, selection can maintain diversity during urban colonization and increase the frequency of adaptive alleles in urban populations (Donihue & Lambert, 2015). For example, immunogenetic variation in the major histocompatibility complex (MHC) gene family remained stable in urban colonizing dark-eyed juncos (Junco hyemalis) amid losses of neutral variation (Whittaker, Dapper, Peterson, Atwell, & Ketterson, 2012).
In some cases, urban areas may even serve as "theaters for evolution," where novel selective pressures drive rapid adaptation (Littleford-Colquhoun, Clemente, Whiting, Ortiz-Barrientos, & Frère, 2017). For example, signatures of urban adaptation have been reported for genes associated with lipid and carbohydrate metabolism (Harris & Munshi-South, 2017), harm avoidance behavior (Mueller, Partecke, Hatchwell, Gaston, & Evans, 2013), and toxicant exposure (Reid et al., 2016). Thus, adaptation is a potent force that can influence divergence in rural-urban conspecifics during initial colonization and longer-term persistence. Teasing its effects apart from drift through examination of genome-wide SNPs and functionally linked markers will help determine how selection versus stochastic events shape wildlife populations across diverse rural-urban gradients.
One of the best documented urban wildlife colonizers has been the red fox (Vulpes vulpes; Figure 1), a generalist species that thrives across rural-urban gradients around the globe (Coman, Robinson, & F I G U R E 1 Red foxes (Vulpes vulpes) have successfully colonized urban areas in Europe since the 1930s. Photo credit: © L. Hamelbeck-Galle/stadtwildtiere.at Beaumont, 1991;Harris, 1981;Marks & Bloomfield, 1999). Red fox urbanization is most pronounced in Europe, where in Switzerland, city foxes dramatically increased in density during the mid-1980s following successful eradication of rabies (Gloor, Bontadina, Hegglin, Deplazes, & Breitenmoser, 2001). Though prophylactic culling continued until 1995, release from rabies enabled fox populations to expand in rural and urban environments, with a large breeding population of more than 10 adult foxes per km 2 detected in Switzerland's largest city, Zurich (Gloor et al., 2001). Mortality records and resident surveys reported a fourfold increase in fox sightings between 1985 and 2000, as family groups established territories throughout the Zurich metropolitan area (Gloor, 2002;Gloor et al., 2001).
Further analysis of movement patterns revealed that foxes settled in rural or urban areas remained almost exclusively within one habitat type, despite proximity to the rural-urban border (Gloor, 2002). This restricted movement suggested the existence of subpopulations within Zurich's urban and adjacent rural areas.
To explore the genetic consequences of these observations, Wandeler, Funk, Largiadèr, Gloor, and Breitenmoser (2003) used 11 neutral microsatellite loci to report on patterns of genetic diversity five to seven generations after initial colonization (with an assumed generation time of two to three years). As predicted in founder events, urban foxes exhibited decreased variation and significant genotypic differentiation when compared to rural foxes, a result supported by previous radio-tracking data (Gloor, 2002). Further, pairwise F ST was highest (F ST = 0.068 ± 0.020) and immigration rates lowest between the two urban subpopulations, suggesting that Lake Zurich and the Limmat River in the Zurich city center presented substantial isolating barriers for urban foxes. The authors concluded that fox colonization of Zurich likely occurred through two founder events (one east and one west of Lake Zurich, the Limmat River, and Zurich's city center), and neutral demographic processes subsequently shaped patterns of genetic diversity.
This work was among the first to document genetic patterns of urban colonization. However, through use of neutral microsatellites, it could only assess neutral genetic structure without considering potential adaptive differences. For example, Zurich foxes may have different selection pressures associated with anthropogenic food availability (Contesse, Hegglin, Gloor, Bontadina, & Deplazes, 2004), which could influence selection on metabolic pathways (Harris & Munshi-South, 2017). City foxes may also experience different disease pressures than rural foxes, as seen with the zoonotic cestode Echinococcus multilocularis transmission cycle within Zurich city foxes (Hofer et al., 2000;Otero-Abad, Rüegg, Hegglin, Deplazes, & Torgerson, 2017). These variable selection pressures, combined with shifts in behavior to avoid human interactions (Gloor et al., 2001;Sih et al., 2011;Soulsbury, Baker, Iossa, & Harris, 2010), could facilitate strong pressures on city foxes not present in rural conspecifics and shape genetic differentiation between the subpopulations.
In the present study, we revisited the question of how urbanization influences population divergence and evolutionary change by re-assessing the city-fox phenomenon in Zurich. Through examination of genome-wide SNPs and MHC-linked microsatellites, we addressed the following questions with a subset of the original Wandeler et al. (2003)  3. Can we identify evidence of selection at MHC-linked markers or outlier SNPs associated with urban colonization?
Given the recent timing of this urban colonization and the generalist life history of red foxes, we hypothesized that fine-scale population structure would separate adjacent rural and urban foxes. Though foxes are capable of dispersing long distances, the localized radio-tracking data reported by Gloor (2002) suggested minimal movement across the rural-urban border. We anticipated further structure between subpopulations east and west of Lake Zurich, the Limmat River, and the Zurich city center, following the results of Wandeler et al. (2003).
Regarding patterns of diversity, we predicted that urban foxes would exhibit decreased variation consistent with founder events when compared to adjacent rural populations. We likewise anticipated a smaller effective population size of urban-dwelling foxes. Lastly, though we predicted that all marker types would exhibit decreased diversity in urban foxes, we anticipated identification of outlier SNPs associated with potentially adaptive urban phenotypes (such as carbohydrate metabolism or exploratory behavior). We additionally expected to find evidence of balancing selection at MHC Class I and II markers (which function in parasite recognition) and directional selection at MHC class III markers (which function in protein secretion). Whereas the neutral markers examined in Wandeler et al. (2003) were limited to demographic inference, the functionally linked and genomic datasets presented herein allowed for consideration of both stochastic demographic events and local adaptation as important evolutionary forces shaping genetic diversity during urban colonization.

| Study system
We re-evaluated urban and rural fox samples collected in the greater Zurich area from spring 1996 to autumn 1998 (Princeton IACUC #2056A-16). Sample selection is detailed in Wandeler et al. (2003), but briefly, we included samples from two urban environments (U east and U west ) that were separated by Lake Zurich, the Limmat River, and the dense Zurich city center (Gloor, 2002).
Three rural areas (R east , R west , and R north ) were selected from within a 20 km radius of Zurich's city center to represent the surrounding rural fox population. We only included individuals believed to be resident foxes based on sampling location, season, sex, and age to avoid erroneous inclusion of dispersing juveniles or males traveling through multiple territories during mating season. We found the methods described by Wandeler et al. (2003) sufficient to establish residency, especially given the relatively small territory sizes (mean = 0.28-0.40 km 2 ) of foxes sampled around Zurich and the restricted movement observed between habitat types (Contesse et al., 2004;Gloor, 2002;Soulsbury et al., 2010).
Due to sample availability and quality requirements, we analyzed a subset of the 128 fox samples used in the original Wandeler et al. (2003) study. Our final samples sizes were 50 foxes (31 rural, 19 urban) for genome-wide SNP analyses and 100 foxes (57 rural, 43 urban) for MHC-linked microsatellite analyses (see Tables 1 and 2 for the number of foxes sampled in each sampling area).

| DNA extraction and restriction-associated DNA sequencing
Genomic DNA was previously extracted with a standard phenol chloroform protocol (Bruford, Hanotte, Brookfield, & Burke, 1992), a QIAamp tissue extraction kit (Qiagen, Inc.), or a Chelex-protocol (Goossens, Waits, & Taberlet, 1998), depending on tissue type (see Wandeler et al., 2003 for details). We quantified DNA with PicoGreen assays and screened for quality with 1% agarose gels, only including samples that had high molecular weight. We then standardized samples to 5 ng/μl.
To generate genome-wide SNPs, we conducted restriction-associated DNA sequencing (RADseq) on all samples that passed DNA quality standards, using a dual-adapter/barcode-ligated library preparation method as described by Ali et al. (2016). DNA was digested with restriction enzyme sbfI, individually barcoded with ligated biotinylated adapters, pooled, and sheared to 400 bp with a Covaris LE220. We used a streptavidin bead-binding assay Invitrogen) to select fragments with the ligated barcode and eluted selected DNA in 55 μl low TE. Libraries were then prepared following standard recommendations for NEBNext Ultra ll DNA Library Prep Kit, where we used Agencourt AMPure XP magnetic beads for library purification and size selection. We standardized final libraries to 10 nM and used paired-end DNA sequencing (2X150nt) on an Illumina HiSeq 2500.
We processed raw sequencing reads with a custom Perl script (sbfI_ flip_trim_150821.pl, see Supporting information) which put all forward and reverse reads with a restriction enzyme cutsite into a single file.
We then used STACKS v1.42 (Catchen, Hohenlohe, Bassham, Amores, & Cresko, 2013) to continue data processing. We demultiplexed samples with process_radtags, where reads were discarded if they had more than two barcode errors or a quality score below 90% based on a sliding window assessment of sequence quality. We next filtered reads for PCR duplicates using clone_filter default parameters and mapped reads to the reference dog CanFam3.1 assembly (Lindblad-Toh et al., 2005) using STAMPY v1.0.21 (Lunter and Goodson 2011). To call genomewide SNPs, we used the ref_map pipeline with -m (minimum number of raw reads needed to create a stack) set to the default parameter 3.
SNPs were further filtered with populations, where we retained only the first SNP per locus (--write-single-snp flag) and removed loci with missing data for more than 10% of foxes (-r 0.90). After SNP genotyping, we additionally removed individual foxes with missing data for >10% genotyped loci. Lastly, we filtered for statistical linkage disequilibrium (LD) in PLINK with --indep-pairwise 50 5 0.3 flag. We evaluated this final, filtered SNP dataset for outliers and aberrant genotypes by implementing principal component analysis with flashPCA to confirm final sample selection (Abraham & Inouye, 2014).

| MHC-linked microsatellite genotyping
We selected twelve MHC-linked microsatellites designed for domestic dogs (C. familiaris) after successful cross-species amplification in TA B L E 1 Mean diversity statistics calculated with 10,149 SNP loci for each subpopulation sampled red foxes (Debenham et al., 2005; R. Wayne, personal communication; Supporting information Table S1). These included microsatellites linked to MHC Class I (n = 2 loci), Class II (n = 8 loci), and Class III (n = 2 loci) genes. Multiplex PCRs were conducted in 10 μl reactions containing 5 μl Qiagen Multiplex PCR Master Mix (Qiagen, Inc.), 1 μl primer mix, 2.1 μl deionized water (diH2O), 0.4 μl bovine serum albumin (Roche, 10 mg/ml), and 1.5 μl template DNA (standardized to 5 ng/μl with diH2O). In each multiplex, forward primers were marked with one of four dye labels (6FAM, NED, PET, and VIC), and primer mixes consisted of 84 μl diH2O mixed with 2 μl 100 μM dye-labeled forward primer and 2 μl 100 μM reverse primer for each primer pair included. Cycling conditions consisted of initial denaturation 95ºC for 5 min followed by 28 cycles of 95ºC for 30 s, 60ºC for 90 s, and 72ºC for 30 s and a final extension 60ºC for 30 min.
We diluted 1 μl of PCR product in 20 μl diH2O and thoroughly mixed 2 μl of diluted product with 9.5 μl of a formamide denatur- We used the R package Genepop (Rousset, 2008) to evaluate LD between all pairs of loci and test for deviations from Hardy-Weinberg equilibrium (HWE) in each subpopulation across all loci and at each locus. Following Wandeler et al. (2003), we implemented the complete enumeration method for loci with fewer than five alleles and the Markov chain method for loci with more than four alleles. We then used Fisher's exact test to identify significant deviations from HWE at the subpopulation level. We determined significance levels using a modified false discovery rate (FDR) with a 5% threshold to account for multiple testing (Benjamini & Yekutieli, 2001).

| Population structure
As traditional methods often fail to detect fine-scale structure in urban populations (Combs, Puckett et al., 2018), we combined several related analyses (described below) to evaluate population admixture, differentiation, and possible colonization routes in Zurich foxes sampled in five areas (i.e., R east , R west , R north , U east , and U west ) in and around Zurich. Briefly, we used ( frameworks. This increased confidence that observed substructure patterns resulted from true signal rather than statistical relicts of any singular method. We first used the admixture model in STRUCTURE 2.3.4 with the SNP dataset to assess genetic cluster assignments based on shared ancestry (Pritchard, Stephens, & Donnelly, 2000). A more traditional approach, STRUCTURE's admixture model often serves as a starting point when assessing population structure. We tested a range of clusters (K) spanning 1-10 (or 1-2n, where n is the number of distinct sampling areas), with 10 independent runs at each K value. Each run had 10,000 burn-in length and 25,000 Monte Carlo Markov Chain repetitions. We determined the optimal K value based on evaluation of log probabilities. Given the geographic proximity of our five sampling areas, we ran this model with sampling location given a priori using LOCPRIOR. This parameter setting aids detection of weak structure without forcing erroneous evolutionary clusters (Pritchard et al., 2000). We next completed a PCA using the SNP dataset through flash-PCA to explore clustering patterns in PC space without additional input information (Abraham & Inouye, 2014). Rather than using shared ancestry to assign population clusters (as in STRUCTURE analyses), PCA-based approaches rely on ordination. This allowed us to consider the variables underlying emergent clusters without including prior assumptions about sampling location. To later incorporate geographic information, we used Spearman's rank correlation tests to compare PC1 with the Swiss coordinate system's X-and Y-coordinates (following Combs, Puckett et al., 2018). We further explored the impact of sampling area on structure by using DAPC in the R package adegenet (Jombart, Devillard, & Balloux, 2010) to visualize differences between our five subpopulations assigned a priori. This approach differs from PCA (which describes total variance in allele frequencies) by maximizing between-group differences while minimizing within-group variance, thus enabling DAPC to detect subtle differences between populations. We visualized the primary discriminate functions in adegenet with scatter and constructed a simple neighbor-joining (NJ) tree with our genetic distance matrix comparing all foxes in the R package ape (Paradis, Claude, & Strimmer, 2004;Saitou & Nei, 1987).
We lastly calculated pairwise F ST to examine genotypic differentiation between subpopulations in STACKS using the p-value-corrected AMOVA F ST method for RADseq data (Catchen et al., 2013).
For MHC-linked microsatellites, we calculated F ST and implemented the G test for genotypic differentiation in Genepop (Rousset, 2008;Weir & Cockerham, 1984). To further estimate allele sharing between subpopulations, we calculated private allelic richness for pairwise combinations of subpopulations in ADZE v1.0 using the SNP dataset (Szpiech, Jakobsson, & Rosenberg, 2008).
The genetic clusters identified through these analyses (considered alongside life history, movement, and sampling data previously collected for Zurich foxes; Gloor, 2002;Soulsbury et al., 2010;Wandeler et al., 2003) provided the focal groups examined in subsequent calculations of genetic diversity statistics and tests for genetic outliers.

| Genetic diversity statistics
We estimated the effective population size (N e ) with our genomewide SNP dataset using the linkage disequilibrium method implemented in NeEstimator v2.1, although these estimates were restricted to "urban" and "rural" foxes (rather than individual subpopulations) due to small sample sizes (Do et al. 2014). We used an allele frequency threshold of 0.01 to minimize bias caused by rare alleles. Though we could not predict census sizes from N e estimates (Frankham, 1995 Due to sample size differences in both datasets, we used a rarefaction approach to estimate allelic richness detected in each subpopulation with increasing sample size in ADZE v1.0 (Szpiech et al., 2008). For RADseq data, we set the maximum standardized sample size to 100 and the missing data tolerance to 25% (Szpiech et al., 2008), in order to maximize rarefied sample size (g) while retaining the majority of loci. For microsatellite data (where more foxes were genotyped at fewer loci), we set the maximum standardized sample size to 200 and the missing data tolerance to one, in order to retain all loci for analysis.
To explore potential effects of inbreeding, we used the SNP dataset to identify runs of homozygosity (ROH) in each subpopulation with the R package detectRUNS (Biscarini, Cozzi, Gaspa, & Marras, 2018).
Following Marras et al. (2015), we used the consecutive runs method for ROH detection and required that ROHs contain a minimum of 15 SNPs, a maximum inter-SNP distance of 1 Mb, and no heterozygous or missing genotypes. We set the minimum ROH length to 100 Kb due to the high incidence of short ROHs detected in mammals (Ceballos, Joshi, Clark, Ramsay, & Wilson, 2018;Frazer et al., 2007).

| Identifying outlier loci
To control for relicts of demographic history, we adopted a similar multifaceted approach to identify outlier loci associated with urbanization. We performed pairwise subpopulation comparisons in pcadapt to detect outliers in clusters identified through previous structure analyses (Luu, Bazin, & Blum, 2017). We then performed an overall comparison of all urban and rural foxes in GEMMA while controlling for sampling location and relatedness (Zhou & Stephens, 2012. Given the current debate surrounding the use of We next used the software package GEMMA to identify outlier loci that may be associated with urban fox colonization (Zhou & Stephens, 2012. GEMMA fits a univariate linear mixed model that tests each SNP for association with the specified phenotype (i.e., urban or rural, based on individual sampling location) while accounting for population stratification and designated covariates. We calculated a centered relatedness matrix from the SNP dataset using the --gk 1 flag in GEMMA and included sex and subpopulation as covariates. We applied a modified FDR to adjust the 5% significance threshold to correct for multiple testing (Benjamini & Yekutieli, 2001
This dataset consisted of 57 rural and 43 urban foxes. Mean successful amplification rate for the nine polymorphic loci included in this dataset was 97.67%, with a minimum of 94.00% for C2_BF_2 and maximum of 100.00% for CFA12-21 and CFA12-19.
Significant LD was found in eleven pairs of MHC-linked loci across all individuals after applying modified FDR correction (Benjamini & Yekutieli, 2001; Supporting information Table S2). As these loci derive from the same gene family, we anticipated statistical linkage at these markers due to their presumed physical proximity or functional similarity (Debenham et al., 2005). We therefore retained all loci for analysis.
No subpopulation reported significant deviation from HWE across all loci (Fisher's exact test; Supporting information Table S3).
We additionally tested each locus within each subpopulation (exact HW test) and observed no significant deviations after correcting for multiple testing (modified FDR; Benjamini & Yekutieli, 2001

| Population structure
We observed minimal differences between results using our full SNP dataset (n = 10,149 loci) and a subdivided dataset of intergenic sites (n = 5,439 loci; Supporting information Appendix S2). We therefore retained all SNPs for downstream analyses to maximize informational content.
Results from the STRUCTURE admixture model with LOCPRIOR conditions returned an optimal K of one evolutionary cluster based on the mean log probability of K (Supporting information Figure S1).
These results suggested shared ancestry between all foxes sampled within and around Zurich and positioned all subsequent analyses to focus on finer-scale substructure.

Within this evolutionary cluster, the overall Mantel test and
Mantel correlogram supported patterns of spatial genetic structure (Supporting information Figure S2) To further explore this relationship, we implemented DAPC with 10 retained PCs and a priori subpopulation assignments to R east , R west , R north, U east , and U west (Figure 3). The first discriminate function showed considerable overlap between sampling areas, with all five clusters evident nonetheless (Figure 3a). Retention of two discriminate functions (Figure 3b) closely matched PCA results, with less overlap between subpopulations and tighter clustering within subpopulations. The divide between sampling areas east (R east and U east ) and west (R west and U west ) of Lake Zurich, the Limmat River, and Zurich's city center further emerged in this analysis, with large branches of the NJ tree broadly corresponding to this geographic subdivision across the lake and river (Figure 3c). Once again, R north was situated between the east-west clusters in the DAPC plots, with assignment to different branches of the NJ tree.
We next considered differentiation between the five clusters observed in our PCA and DAPC plots. At MHC-linked microsatellite loci, significant differences were reported for most pairs of populations (G test) after correcting for multiple testing, with the exceptions of R north + R east , R north + R west , and R east + U east (Supporting information   Table S4). This result suggested little immunogenetic differentiation between rural populations and adjacent rural and urban populations on the east side of the Limmat River.  Tables S4 and S5).
We calculated F ST of a similar magnitude between urban and rural subpopulations separated by Lake Zurich, the Limmat River, and Zurich's city center, with the lowest values reported between rural subpopulations and adjacent rural-urban subpopulations.
As such, private allelic richness for pairwise combinations of subpopulations was highest for adjacent rural-urban subpopulations (R west + U west = 0.0360; R east + U east = 0.0324) and lowest for two pairings separated by the center city barriers (U west + U east = 0.0189; R west + U east = 0.0181; Supporting information Table S6).
Taken together, these results (considered alongside radio-tracking and home range data; Gloor, 2002;Gloor et al., 2001;Soulsbury et al., 2010) supported analysis of five subpopulations in subsequent calculations of genetic diversity and outlier detection tests.
Though comprising one evolutionary lineage, foxes sampled within and around Zurich exhibited geographic substructure confirmed by diverse yet complementary analyses.
F I G U R E 2 Principal components calculated for 50 foxes across all 10,149 SNPs. When plotted against the Swiss Y-coordinate (northing), PC1 recapitulates the geographic sampling area (inset; adapted from Wandeler et al. 2003), thus mirroring the Swiss X-coordinate (easting).

| Genetic diversity statistics
Estimates for N e using the LD method were 41.0 for urban (U east + U west ) foxes and 5,915.5 for rural (R east + R west + R north ) foxes.
Though N e estimates should not be mistaken for projected census sizes, they do suggest far greater numbers of foxes outside of the city limits.
Genetic diversity statistics for the SNP dataset are reported in Table 1. Both urban subpopulations exhibited reduced allelic diversity when compared to rural subpopulations with similar sample sizes, as measured by the number of private alleles. To compare allelic diversity across all areas, we used a rarefaction approach implemented in ADZE v1.0 to calculate allelic richness at standardized sample sizes (Szpiech et al., 2008). Here, U east exhibited the lowest allelic richness of all subpopulations, with the remaining four subpopulations clustered together (Supporting information Figure S4A).
Heterozygosity deviated from this trend of reduced urban diversity (Table 1) We lastly observed a strong positive correlation between the sum of ROH and the number of ROH calculated per fox (r 2 = 0.918), with urban foxes (particularly U east ) more densely concentrated in the upper right quadrant of the plot than rural foxes (Figure 4c). This trend is reflected by mean values calculated for each subpopulation.
As we expect newly colonized populations to exhibit more ROHs than larger source populations, these analyses aided our understanding of fox colonization history in and around Zurich.

| Identifying outlier loci
Following our population structure results, we analyzed foxes east (R east + U east ) and west (R west + U west ) of Lake Zurich, the Limmat River, and Zurich's city center separately in pcadapt. We excluded R north from this analysis due to its intermediate position between R east and R west in the PCA and DAPC analyses. For both runs of pcadapt, we regressed the SNP dataset (n = 10,149 loci) against the first two PCs, since K = 1-2 captured the relevant population structure in this and previous (DAPC) analyses. Using a significance threshold of q < 0.05, we detected 514 outliers in east subpopulations and 508 outliers in west subpopulations, with 42 loci appearing as outliers in both analyses. These overlapping SNPs included 17 annotated as intergenic, 27 as intronic, three as exonic, and two near promoters (note that some SNPs have multiple annotations). We queried genic outliers significant in both analyses (Supporting information Table S7) in the Ensembl, OMIM, and GeneCards databases, and observed numerous functions related to energy metabolism, drug tolerance, and immune processes. Ensembl VEP annotations included "low impact" (i.e., minimal effect on protein), "moderate impact" (i.e., may alter protein effectiveness), and "modifier" (i.e., no known effects on protein).
We next examined SNPs associated with urban versus rural habitats in the full fox dataset (n = 50 individuals), using GEMMA to fit a univariate linear mixed model with sex and sampling area as covariates. We further controlled for population structure by including a relatedness matrix of all foxes sampled. We analyzed all 10,149 SNPs and adjusted our likelihood ratio test (lrt) significance threshold using the modified FDR (adjusted p = 0.005). This analysis produced 91 SNPs significantly associated with the urban phenotype, with 48 annotated as intergenic, 42 as intronic, 11 as exonic, and three near promoters (again, note that some SNPs have multiple annotations).
Although no SNPs overlapped with outliers identified in pcadapt, possible gene functions once again related predominantly to immunity and metabolism (among other cellular processes), and VEP analysis included "low impact," "moderate impact," and "modifier" annotations.
To identify outlier loci in the MHC-linked microsatellite dataset, we implemented the Ewens-Watterson homozygosity test of neutrality in PyPop (Supporting information Figure S5;

| D ISCUSS I ON
In the present study, we revisited the city-fox phenomenon in Zurich, Switzerland, to examine the genetic effects of urban colonization and the evolutionary processes that shaped them.
Through concurrent analysis of genome-wide SNPs and MHClinked microsatellites, we addressed the following questions: (1) Do we observe subpopulation structure between urban and rural foxes in Zurich? (2) What patterns of genetic diversity characterize urban fox colonization? (3) Can we identify evidence of selection at MHC-linked markers or outlier SNPs associated with urban colonization?
Considered together, our results suggested the presence of one evolutionary cluster subdivided into smaller groups clustered by geographic sampling area. We observed significant differentiation across natural (i.e., Lake Zurich and the Limmat River) and anthropogenic (i.e., urban infrastructure separating adjacent rural and urban foxes) barriers that were consistent across multiple data-types and complementary analytical frameworks. These results matched previous radio-tracking data that suggested smaller home range sizes and localized movement of foxes inhabiting rural and urban areas around Zurich (Gloor, 2002 populations. As rare alleles typically decline faster than genomewide heterozygosity, there often exists a temporal lag between these diversity metrics before populations reach a new equilibrium (Cornuet & Luikart, 1996). Given the recent nature of this colonization event and subtleties of within-population clustering, we anticipate no long-term negative effects of these declines.
These results should be interpreted with caution, due to the small sample size, minimal overlap of specific outlier sites, and ongoing debate on the use of RADseq data for detecting genomic outliers (Catchen et al., 2017;Lowry et al., 2016;McKinney et al., 2017 (Gilroy, Phillips, Richardson, & Oosterhout, 2017).
Through examination of genome-wide SNPs and MHC-linked microsatellites, this study expanded upon the results reported by Wandeler et al. (2003) to provide a more comprehensive understanding of the genetic effects of urban fox colonization and the evolutionary processes that shaped them. We reported patterns of genetic variation consistent with founder events and increased differentiation consistent with natural and anthropogenic barriers to dispersal. We additionally identified outlier loci with putative gene functions related to urban-associated processes, such as energy metabolism, behavior, and immunity. In addition to system-specific information gained, this study contributes toward our overall understanding of the genetics of urban evolution, an exciting frontier within urban evolutionary biology (Johnson & Munshi-South, 2017;Santangelo et al., 2018). It further emphasizes the value of combining datasets to parse the roles of stochastic and adaptive processes underlying evolutionary change in an increasingly urban world.

ACK N OWLED G M ENTS
This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1656466.

CO N FLI C T O F I NTE R E S T
None Declared.

AUTH O R CO NTR I B UTI O N S
ALD, KEB, and BMvH designed the study; GC, PW, and CD contributed samples; ALD, KEB, and CC conducted the laboratory work; ALD and KEB processed the raw data; ALD analyzed the data; ALD, KEB, EH, and BMvH provided analytical assistance; ALD, KEB, and BMvH prepared the manuscript; all authors contributed to and approved the final manuscript.

DATA ACCE SS I B I LIT Y
RADseq clone filtered and demultiplexed fastq files and CanFam3.