https://orcid.org/0000-0001-8485-5556
Abstract
Theory predicts that range expansion results in genetic diversity loss in colonizing populations. Rapid reduction of population size exacerbates negative effects of genetic drift, while sustained isolation decreases neutral variation. Amid this demographic change, natural selection can act to maintain functional diversity. Thus, characterizing neutral and functional variation is critical for disentangling the evolutionary forces that shape genetic variation in newly established populations. Coyotes (Canis latrans) provide an ideal study system for examining the genetic effects of urban colonization. Capable of thriving in environments ranging from natural to highly urbanized, this mobile carnivore recently established a breeding population in New York City (NYC), one of the most densely populated areas in the USA. In the present study, we characterized neutral and functionally linked genetic diversity on a regional scale, traversing NYC and its surrounding counties in the New York metropolitan area. We report decreased variation and significant genotypic differentiation in NYC coyotes following recent colonization of this super-urban environment. In accordance with our hypotheses, we observed evidence for a recent population bottleneck as coyotes entered NYC. Counter to our expectations, we found only minimal support for selection maintaining diversity at immune-linked loci. These findings suggested that stochastic processes, such as genetic drift, more likely drove patterns of decreased variation in super-urban coyotes. This work not only improves our understanding of NYC’s newest inhabitants, but also contributes to the growing body of knowledge surrounding urban colonization ecology. It highlights the importance of examining both neutral and functional variation when assessing the roles of drift and selection in newly established populations. When combined with similar studies across diverse systems, these insights can aid wildlife management and green design to better facilitate gene flow and maintain healthy populations of wildlife in an increasingly urban world.
Introduction
Population genetic theory predicts that range expansion results in genetic diversity loss at the expansion front (Slatkin and Excoffier 2012). Newly established populations founded by a small number of migrants can undergo genetic bottlenecks, where rapid reduction in effective population size exacerbates the negative effects of genetic drift (Nei, Maruyama, and Chakraborty 1975). If populations remain small and isolated from conspecifics, further loss of variation may result in inbreeding depression (Keller and Waller 2002; Frankham 2005; Oakley 2013). This can lower a population’s ability to respond to environmental stressors, such as disease, and threaten long-term viability (Coltman et al. 1999; Spielman et al. 2004; Smith, Acevedo-Whitehouse, and Pedersen 2009).
This predicted loss of genetic diversity is typically observed in selectively neutral markers (Lande and Barrowclough 1987; Tarr, Conant, and Fleischer 1998; Clegg et al. 2002; Wandeler et al. 2003; Fabbri et al. 2007; Evans et al. 2009; Zachos et al. 2009). Yet in the midst of demographic changes, selection may act to maintain functional variation at adaptive loci. In particular, immunogenetic diversity at hypervariable gene families, such as the major-histocompatibility complex (MHC), is often maintained by balancing selection (Bernatchez and Landry 2003; Ferrer-Admetlla et al. 2008). As these loci function in foreign antigen recognition, MHC diversity is thought to buffer populations from disease risk by facilitating detection of numerous pathogens (Spielman et al. 2004; DeCandia, Dobson, and vonHoldt 2018). This may be especially important for small, isolated populations or those colonizing novel environments with altered disease risk. Examples of balancing selection acting on immunogenetic loci include Channel Island foxes (Urocyon littoralis; Aguilar et al. 2004), Soay sheep (Ovis aries; Charbonnel and Pemberton 2005), house sparrows (Passer domesticus; Loiseau et al. 2009), water voles (Arvicola terrestris; Oliver and Piertney 2012), bobcats (Lynx rufus; Serieys et al. 2015a) and lynx (Lynx spp.; Marmesat et al. 2017).
Notably, the presence of balancing selection does not preclude drift from acting on immune loci, particularly when bottlenecks are sustained (Alcaide 2010; Ejsmond and Radwan 2011). For instance, both selection and drift have shaped immunogenetic variation in Atlantic salmon (Salmo salar; Landry and Bernatchez 2001), New Zealand robins (Petroica spp.; Miller and Lambert 2004), tuatara (Sphenodon spp.; Miller, Allendorf, and Daugherty 2010), prairie-chickens (Tympanuchus cupido; Bollmer et al. 2011) and golden snub-nosed monkeys (Rhinopithecus roxellana; Luo et al. 2012). In each of these examples, seemingly opposing forces shaped immunogenetic variation simultaneously. It is therefore critical to consider both neutral and functional variation when examining the genetic effects of colonization and the evolutionary processes that drive them.
As urbanization transforms the human–wildlife interface on a global scale, natural lands are remodeled for human habitation and resource production, creating novel environments for wildlife to (re)colonize (Bateman and Fleming 2012; Magle et al. 2012). Though the definition of ‘urban areas’ changes through time and space, cities share numerous environmental features (Johnson and Munshi-South 2017). Increased human densities, impervious surface cover and fragmentation of natural lands create a challenging matrix for wildlife to navigate (Zipperer et al. 2012). In addition to diversity losses expected from founder events, barriers to dispersal, such as freeways and buildings, may cause increased differentiation between urban colonizers and nonurban source populations (Riley et al. 2006; Johnson and Munshi-South 2017). Such phenomena have been observed at neutral microsatellite loci in blackbirds (Turdus merula) colonizing European and north African cities (Evans et al. 2009) and red foxes (Vulpes vulpes) establishing territories in Zurich, Switzerland (Wandeler et al. 2003; DeCandia et al. 2019). In each case, genetic drift shaped variation during urban colonization. However, novel selective pressures, such as altered exposure to pollutants or disease, may counteract drift to maintain diversity at adaptive loci during urban colonization events. For example, dark-eyed juncos (Junco hyemalis) colonizing the University of California, San Diego campus exhibited similar levels of immunogenetic diversity as their source population in the mountains, despite reductions in neutral variation (Whittaker et al. 2012). Similarly, immunogenetic markers examined in urban bobcats outside of Los Angeles showed evidence of balancing selection following a severe mange epizootic (Serieys et al. 2015a).
Coyotes (Canis latrans) inhabiting the New York metropolitan area present an ideal system for exploring the neutral and functional genetics of urban colonization. Having expanded their range to include every state in the continental USA (Gompper 2002; Heppenheimer et al. 2018; Hody and Kays 2018), coyotes thrive in habitats along the entire urbanization gradient from undeveloped landscape through city centers (Bateman and Fleming 2012). Beginning in the 1940s, they established territories in the highly altered landscapes of the northeast (see Fig. 3, Hody and Kays 2018). In the last two decades, this expansion included New York City (NYC), one of the most densely populated areas in the USA (Weckel et al. 2015). Examining human census data compiled in 2010, the New York metropolitan area stands out as a super-urban area containing over 19 000 000 people (U.S. Census Bureau 2017). As points of comparison, the second largest metropolitan area in the USA is Los Angeles, CA with over 12 000 000 people, and the third largest is Chicago, IL with over 9 000 000 people (U.S. Census Bureau 2017). Despite NYC’s high population density and impervious surface cover, coyotes continue expanding into the area, allowing researchers to study super-urban colonization in real time (Weckel et al. 2015; Nagy et al. 2017). For instance, Nagy et al. (2016) used camera traps to assess coyote distribution, breeding status and seasonal occupancy in NYC parks over a 3-year period. Through this footage, the authors directly observed colonization of habitat patches within NYC’s urban landscape from initial occupancy through establishment of viable breeding sites.
In the present study, we examined the genetic effects of this super-urban colonization. We surveyed nine microsatellite markers found in neutral regions of the genome and 14 microsatellite markers linked to functional MHC genes in coyotes inhabiting New York, New Jersey and Connecticut. Following their recent colonization of NYC, we hypothesized that super-urban coyotes would exhibit reduced neutral diversity compared with coyotes living outside of the city as a result of genetic drift. In contrast, we anticipated no significant differences between NYC and non-NYC coyotes at immune-linked loci. We hypothesized that balancing selection would maintain this functional variation (as seen in Aguilar et al. 2004; Whittaker et al. 2012; Serieys et al. 2015a), due to the novel selection pressures posed by the urban disease landscape. Though it is possible that directional selection may decrease immunogenetic variation by selecting for specific alleles, we anticipated no single force exerting strong enough pressure to selectively prune immunogenetic diversity this soon after colonization. Instead, we predicted that multiple, potentially competing forces, have maintained variation through balancing selection. In addition to system-specific information gained, this study contributes to our understanding of the genetic effects of urban colonization and demonstrates the value of examining neutral and adaptive loci in the broader context of urban colonization ecology.
Methods
Sampling area
Sampling locations spanned 51 613.15 km2 of land area across 27 counties in New York, New Jersey and Connecticut (Supplementary Table S1). The non-NYC coyote samples (n = 146) came from the archived tissue collection maintained by the New York State Museum (Albany, NY) and were collected between 2000 and 2010. NYC coyote samples came from scat collected at parks in the Bronx, Queens and Manhattan between 2010 and 2016 (n = 245) where coyotes had been observed with camera traps (Nagy et al. 2016). Blood (n = 3) and tissue (n = 6) samples were collected opportunistically from roadkill and translocated coyotes within NYC. Dry samples were collected whole and stored in paper bags, each containing two 5 g packs of silica desiccant. Samples that were moist upon collection were stored in plastic bags and frozen within 24 h at −20°C.
DNA extraction and species identification
We extracted DNA from tissue following the DNeasy Blood & Tissue Kit (Qiagen, Inc.) standard protocol. For scat samples, we used a QIAamp DNA Stool Mini Kit (Qiagen, Inc.) with a modified stool extraction protocol (Chaves, Dias, and Pomilla 2010). We assessed DNA concentration with Qubit™ fluorometric quantitation and standardized high concentration DNA to 5 ng/μl for downstream analyses.
To distinguish coyote scat from domestic dogs (Canis familiaris), red foxes (V. vulpes), or other species, we used polymerase chain reaction (PCR) to amplify samples with the dye-labeled (6FAM) forward primer SIDL and two different reverse primers, H16145 and H3R (De Barba et al. 2014) that target the mitochondrial DNA control region. These primers amplify two fragments in coyotes and dogs and one fragment in red foxes. As fragment lengths vary by species (coyotes amplify 115–120/360–364 bp; dogs amplify 123–128/365–368 bp; red foxes amplify —/343–347 bp), we were able to identify the species of each sample. The 14 µl PCR consisted of 7.0 μl Qiagen Multiplex PCR Master Mix, 1.4 μl Qiagen Q solution, 0.406 μl SIDL, 0.28 μl H16145, 0.140 μl H3R, 2 µL template DNA and 2.774 μl sterile water. Each primer was at a concentration of 10 µM. PCR cycling conditions included initial denaturation of 94°C/15 min, followed by 35 cycles of (94°C/30s, 46°C/90s, 72°C/60s), and a final extension of 60°C/30 min. PCR products were sent to Genewiz (Frederick, MD) for fragment analysis and assigned genotypes by visual inspection of electropherograms in Geneious v11.1.4 with the microsatellite plug-in (Kearse et al. 2012).
Neutral microsatellite genotyping
To confirm individual identity and assess neutral genetic diversity, we used the nine microsatellites designed for domestic dogs and validated for use on coyotes by Mumma et al. (2014; Supplementary Table S2). All nine dye-labeled (6FAM, NED, PET or VIC) microsatellite primers were combined into one 14-µL PCR multiplex (following Mumma et al. 2014). The multiplex consisted of 7.0 μl Qiagen Master Mix, 1.4 μl concentrated Qiagen Q Solution, 0.322 μl Cxx.119, 0.196 μl FH2670, 0.140 μl FH2611 and DBX6, 0.126 μl FH3725 and FH2001, 0.084 μl FH2054, FH2137, and DBY7, 0.056 μl FH2088 and C09.173, 2 µL template DNA and 2.676 μl sterile water. All primers were once again 10 µM. Cycling conditions included initial denaturation 94°C/15 min followed by a touchdown step of 13 cycles of (94°C/30s, 62°C/90s, 72°C/60s) decreasing 0.4°C each cycle, followed by another round of 28 cycles of (94°C/30s at, 57°C/90s, 72°C/60s), finishing with a final extension of 60°C/30 min. PCR products were sent to Genewiz (Frederick, MD) for fragment analysis. Genotypes were again assigned using Geneious v11.1.4 with the microsatellite plug-in.
Samples derived from scat were PCR amplified up to seven times. To generate consensus genotypes, we required at least two identical PCR results to verify heterozygous alleles and at least three identical PCR results to verify homozygous alleles (Prugh et al. 2005). To identify the number of loci required for individual identification from scat, we calculated the cumulative probability of identity for unrelated individuals P(ID) and for siblings P(ID)sib in Gimlet 1.3.3 (Valière 2002). We used a P(ID)sib < 0.01 as our standard for the minimum number of loci needed for individual identification (Waits, Luikart, and Taberlet 2001). We additionally used a threshold of 30% missing loci to exclude poor quality samples from downstream analyses (following vonHoldt et al. 2008; Heppenheimer et al. 2018).
MHC-linked microsatellite genotyping
To characterize immunogenetic diversity, we used 14 MHC-linked microsatellites designed for domestic dogs and successfully amplified in coyotes (Supplementary Table S3; Debenham et al. 2005; Wayne, pers. comm.). These include markers near MHC Class I (n = 3 loci), Class II (n = 10 loci) and Class III (n = 1 locus) regions.
We conducted multiplex PCRs in 10 μl reactions consisting of 5 μl Qiagen Master Mix, 1 μl primer mix, 2.1 μl deionized water (diH2O), 0.4 μl bovine serum albumin (BSA; Roche, 10 mg/ml) and 1.5 μl template DNA. For DNA derived from scat samples, template DNA was increased to 2 μl and diH2O reduced to 1.6 μl. Primer mixes consisted of 84 μl diH2O mixed with 2 μl 100 μM dye-labeled (6FAM, NED, PET or VIC) forward primer and 2 μl 100 μM reverse primer for each primer pair included in that reaction.
Cycling conditions for DNA extracted from tissue consisted of initial denaturation 95°C/5 min followed by (95°C/30s, 60°C/90s, 72°C/30s) × 28 cycles and final extension 60°C/30 min. Cycling conditions for DNA extracted from scat samples consisted of initial denaturation 94°C/15 min followed by touchdown cycling (94°C/30s, 62°C–0.4°C/90s, 72°C/60s) × 13 cycles, cycling with a consistent annealing temperature (94°C/30s, 57°C/90s, 72°C/60s) × 28 cycles and final extension 60°C/30 min.
We sent all MHC-linked PCR products to the Cornell Institute of Biotechnology (Ithaca, NY) for fragment analyses and subsequently assigned genotypes by visual inspection of electropherograms in Geneious v6.1.6 (Kearse et al. 2012). Scat-derived samples were run up to six times to ensure confidence in genotype calls, with two and three identical results required for heterozygous and homozygous calls, respectively (Prugh et al. 2005). As with neutral microsatellite loci, we used a missing data threshold of 30% to exclude poor quality samples from statistical analyses (vonHoldt et al. 2008; Heppenheimer et al. 2018).
Classifying coyote subgroups in the New York metropolitan area
We explored multiple lines of evidence to identify subgroups of coyotes for comparative analyses. These included behavioral (estimated coyote home ranges), environmental (impervious surface cover and human population density) and molecular (individual assignments to genetic clusters) variables to enable biologically relevant subgroup designation.
To explicitly consider environmental modifications, geography and coyote movement behavior, we used the National Land Cover Database (NLCD) Percent Developed Imperviousness layer in ArcMap 10.6.1 (Homer et al. 2015; ESRI 2016) to visualize the footprint of human infrastructure across our sampling area. We then plotted each sampling location with an estimated coyote home range buffer of 98 km2, or the maximum home range estimate for transient coyotes in the Chicago metropolitan area (Supplementary Fig. S1; Gehrt, Anchor, and White 2009). We traced an outline around overlapping home ranges (being careful to mind natural borders, such as the Atlantic Ocean and Great Lakes) to estimate the dispersal capabilities of coyotes in our sampling area.
We next considered human population density as a metric for urbanization, as coyotes thrive in human modified environments (e.g. using roads as corridors), while often shifting activity patterns to avoid human contact (George and Crooks 2006; Gehrt, Anchor, and White 2009; Gehrt, Riley, and Cypher 2010; Toews, Juanes, and Burton 2017). Higher human population densities, rather than increased urban infrastructure, may therefore pose a greater obstacle for coyote habitation. We retrieved census data for 2010 from the United States Census Bureau’s American Fact Finder database (https://factfinder.census.gov) for each county sampled, as the resolution of sampling locations was predominantly at the county level. We then calculated human population density as people per square kilometer land area.
To examine genetic clusters (K), we used STRUCTURE 2.3.4 (Pritchard, Stephens, and Donnelly 2000) without any sampling location information given a priori. We set a range of K = 1–10, with 10 independent runs at each K-value to estimate group assignment. For each run, we set 100 000 burn-in iterations and 500 000 Monte Carlo Markov Chain (MCMC) iterations. We then used Structure Harvester v06.94 (Earl and vonHoldt 2012) to select the model with the highest rate of change in log likelihood values (ΔK) based on the Evanno method (Evanno, Regnaut, and Goudet 2005), and aligned output from independent replicates in CLUMPP for downstream visualization (Jakobsson and Rosenberg 2007). Since coyotes typically exhibit minimal geographic structuring throughout their range (Heppenheimer et al. 2018), we used the same parameters to run STRUCTURE with location information (based on our environmental variables) specified a priori. This method enables detection of weak population structure, without identifying spurious clusters (Pritchard, Stephens, and Donnelly 2000).
Statistical analyses
We evaluated linkage disequilibrium (LD) between all pairs of loci across the entire sample set and tested for deviations from Hardy–Weinberg equilibrium (HWE) in each designated subgroup with the R package Genepop (Rousset 2008). We identified significant deviations from HWE with Fisher’s exact test, and applied a modified false discovery rate correction to adjust for multiple testing (Benjamini and Yekutieli 2001).
We then used GenAlEx v6.503 to calculate diversity statistics such as observed heterozygosity (HO), expected heterozygosity (HE), allelic richness (NA), number of private alleles (PAS) and inbreeding coefficient (FIS) for neutral and MHC-linked microsatellite loci genotyped in each subgroup (Peakall and Smouse 2012). We implemented paired t-tests in R to compare mean HO and NA values calculated per locus between subgroups (Luikart et al. 1998; R Core Team 2013). Due to the sample size sensitivity of allelic diversity metrics, we performed rarefaction analyses for allelic richness in ADZE v1.0 (Szpiech, Jakobsson, and Rosenberg 2008).
We estimated effective population size for both populations using the program LDNe (Waples and Do 2008) with the lowest allele frequency cutoff of 0.01. To test for population structure, we used Genepop to calculate the fixation index (FST) and implement the G-test for genotypic differentiation between subgroups (Weir and Cockerham 1984; Rousset 2008). We subsequently estimated gene flow by calculating the number of genetic migrants (NM) from FST values using FST ≈ 1/(1 + 4NM) (Wright 1921; Riley et al. 2006). We additionally quantified population differentiation with Dest using the R package DEMEtics (Gerlach et al. 2010). As Dest is unaffected by within-population heterozygosity (unlike FST), it has been shown to better characterize population differentiation when using highly polymorphic markers, such as microsatellites (Jost 2008; Heller and Siegismund 2009).
We next estimated the effects of genetic drift versus natural selection in driving patterns of diversity. To assess whether subgroups recently underwent a population bottleneck, we used the program BOTTLENECK to test for heterozygosity excess in neutral microsatellite loci (Cornuet and Luikart 1996; Piry, Luikart, and Cornuet 1999). As the number of unique alleles decreases faster than heterozygosity in recently bottlenecked populations, we would expect to see heterozygosity excess in populations following founder events. Here, we considered P-values up to 0.10 as potentially biologically relevant, as small sample sizes common in carnivore studies may limit our ability to detect statistically significant excess (Krausman 2017). We restricted these analyses to neutral loci, as the MHC-linked loci violated program assumptions of independence and neutrality.
We then implemented the Ewens–Watterson homozygosity test of neutrality in PyPop (Lancaster et al. 2007) to examine whether selection is acting on MHC-linked loci in any subgroups. In this test, the normalized deviate of homozygosity (FND) compares observed levels of homozygosity to those expected under neutrality, with significantly negative FND values indicating balancing selection. In these scenarios, observed homozygosity is lower than expected homozygosity under neutral conditions, which suggests selective maintenance of diversity. Positive FND values suggest directional selection, as observed homozygosity is higher than that expected under neutral conditions. This suggests maintenance of specific alleles (rather than overall diversity). We should note that demographic events, such as bottlenecks, can also cause significant deviations from neutrality. However, selection is distinguished from past demographic events by examining which loci exhibit these deviations. Bottlenecks or other demographic events should act on all loci (both neutral and functional), whereas selection should only affect relevant functional loci (or markers linked to functional genes, as is the case with our MHC-linked microsatellites). For example, if a hypothetical population is threatened by extracellular pathogens, we might expect to see evidence of selection in MHC Class II genes, which function in antigen presentation of extracellular proteins.
Results
Sample selection and microsatellite genotyping
Following DNA extraction, we discarded 56 non-NYC (tissue) and 35 NYC (scat) samples for low DNA concentrations. We discarded an additional 53 scat samples that did not amplify at either mitochondrial locus used to differentiate species.
The PID analysis in Gimlet returned a minimum of any six neutral loci required for individual identification (up to 30% missing data) using the threshold of P(ID)Sibs < 0.01. The cumulative P(ID) and P(ID)Sibs for the six loci were 1.305E-7 and 0.005 and for all nine loci were 1.334E-10 and 4.152E-4. We removed 74 genotypes presumed to be duplicates generated from multiple scat samples collected from the same individuals. Genotypes were considered duplicates if they matched at six or more loci.
Following PCR amplification of remaining individuals at nine neutral and 14 MHC-linked microsatellite loci, we discarded 77 samples with poor amplification rates (i.e. missing data ≥30% at either neutral or MHC-linked loci; vonHoldt et al. 2008; Heppenheimer et al. 2018). This produced a final dataset of 105 unique coyotes.
Classifying coyote subgroups in the New York metropolitan area
Considered together, behavioral, environmental and genetic evidence supported designation of two primary subgroups for comparative analyses: NYC (n = 21) versus non-NYC (n = 84) coyotes (Fig. 1). Regarding movement behavior, we noted that all coyotes fell within a single area of overlapping home range estimates (Fig. 1A). This suggested wide-scale movement capability throughout the sampling area. Regarding environmental variables, we observed that impervious surface cover was far more concentrated in NYC (indicated in purple) than elsewhere in the sampling area (Fig. 1A). We further noted that the three counties comprising NYC (Bronx, New York and Queens counties) exhibited uniquely high human population densities (Fig. 1B;Supplementary Table S1). As such, NYC represented a super-urban area distinct from all other sampling locations with regard to human infrastructure and presence, potentially rendering it a more difficult matrix for coyotes to navigate and successfully colonize.
Behavioral, environmental and genetic factors considered when classifying subgroups of coyotes in the New York metropolitan area. (A) Map of percent developed imperviousness from NLCD 2011, individual sampling locations and an outline of overlapping estimated home range buffers (98 km2) around each point. (B) Human population density in each US county sampled. (C–E) STRUCTURE results for (C) K = 4 with no location information specified a priori and (D) K = 3 and (E) K = 2 with location information specified a priori. Considered together, these factors supported designation of two coyote subgroups for comparative analyses: NYC and non-NYC.
Following these analyses, genetic evidence supported the existence of one large intermixing population, with restricted movement around NYC. STRUCTURE results generated without prior location information suggested an optimal K of four clusters using the Evanno method (Fig. 1C). Examination of individual assignments, however, revealed widespread admixture throughout the sampling area, with NYC coyotes (gray) and a collection of individuals sampled in Hunterdon, Franklin and a few other counties (light green) emerging as somewhat distinct clusters. We therefore ran STRUCTURE a second time with location information (NYC vs. non-NYC) specified a priori. Here, the optimal K was three clusters, with one large non-NYC group containing most coyotes sampled (dark green), a small group of non-NYC coyotes (light green) and the NYC coyotes comprising their own cluster (gray; Fig. 1D). We then visualized K = 2 (Fig. 1E), as this K-value is more consistent with the behavioral (dispersal) and environmental (impervious surface cover and human population density) factors considered previously in subgroup designation. In this analysis, non-NYC coyotes represented one cluster (dark green) with NYC coyotes representing the other (gray).
In summary, NYC coyotes consistently clustered together (Fig. 1C–E) and were sampled in areas where environmental variables associated with urbanization (such as impervious surface cover and human population density) exhibited extremely high values when compared with other locations sampled. This matched our understanding of the area’s colonization history, with NYC representing a newly established population within the larger metropolitan area. We therefore concluded that our sampling area contained one large population of highly mobile carnivores that recently expanded into a novel, super-urban area. As a result, we delimited our two subgroups as NYC and non-NYC coyotes for subsequent comparative analyses.
LD and HWE
Only one of 36 pairs of neutral loci displayed evidence of LD across all 105 coyotes (Supplementary Table S4), whereas LD tests of MHC-linked loci revealed 72 of 91 locus pairs with significant correlations (Supplementary Table S5). We anticipated this result of statistical linkage at immunogenetic markers due to their presumed physical proximity and functional similarity, as they derive from the same gene complex (Debenham et al. 2005). We therefore retained all neutral and MHC-linked loci for further analysis.
Two neutral loci in the non-NYC subgroup deviated significantly from HWE while no deviation was reported from the NYC subgroup (Supplementary Table S6). At MHC-linked loci, four loci in non-NYC and one locus in NYC coyotes exhibited significant deviation from HWE (Supplementary Table S7). In each case, heterozygote deficiency drove the observed pattern. As no loci were consistently deviated in both NYC and non-NYC coyotes, we retained all loci for downstream analysis.
Genetic diversity at neutral microsatellite markers
We found no significant difference in mean observed heterozygosity between non-NYC and NYC coyotes (HO, non-NYC = 0.745 ± 0.042, NYC = 0.783 ± 0.040; one-tailed paired t-test, t = 0.899, P = 0.803; Table 1; Supplementary Table S6). Conversely, the mean number of alleles was significantly lower in NYC coyotes (NA, non-NYC = 10.333 ± 0.898, NYC = 6.444 ± 0.626; one-tailed paired t-test, t = 5.039, P < 0.001), with 35 PAS across all loci present in the non-urban population. To circumvent sample size sensitivity of these metrics, we used rarefaction analysis to compare alleles detected with increasing sample size for non-NYC and NYC coyotes. We found that NYC coyotes carried fewer alleles than non-urban coyotes at any given sample size (Fig. 2A).
Diversity statistics for NYC and non-NYC coyotes at neutral (left) and MHC-linked (right) microsatellite loci
| Subgroup | NA | HO | HE | PAS | FIS |
|---|---|---|---|---|---|
| NYC | 6.444 ± 0.626/5.286 ± 0.384 | 0.783 ± 0.040/0.569 ± 0.047 | 0.715 ± 0.037/0.634 ± 0.032 | 0/4 | −0.103 ± 0.046/0.110 ± 0.055 |
| Non-NYC | 10.333 ± 0.898/8.214 ± 0.990 | 0.745 ± 0.042/0.667 ± 0.033 | 0.779 ± 0.029/0.724 ± 0.030 | 35/45 | 0.047 ± 0.029/0.080 ± 0.020 |
| Subgroup | NA | HO | HE | PAS | FIS |
|---|---|---|---|---|---|
| NYC | 6.444 ± 0.626/5.286 ± 0.384 | 0.783 ± 0.040/0.569 ± 0.047 | 0.715 ± 0.037/0.634 ± 0.032 | 0/4 | −0.103 ± 0.046/0.110 ± 0.055 |
| Non-NYC | 10.333 ± 0.898/8.214 ± 0.990 | 0.745 ± 0.042/0.667 ± 0.033 | 0.779 ± 0.029/0.724 ± 0.030 | 35/45 | 0.047 ± 0.029/0.080 ± 0.020 |
Metrics include mean number of alleles (NA), observed heterozygosity (HO), expected heterozygosity (HE), number of PAS and inbreeding coefficient (FIS)
Diversity statistics for NYC and non-NYC coyotes at neutral (left) and MHC-linked (right) microsatellite loci
| Subgroup | NA | HO | HE | PAS | FIS |
|---|---|---|---|---|---|
| NYC | 6.444 ± 0.626/5.286 ± 0.384 | 0.783 ± 0.040/0.569 ± 0.047 | 0.715 ± 0.037/0.634 ± 0.032 | 0/4 | −0.103 ± 0.046/0.110 ± 0.055 |
| Non-NYC | 10.333 ± 0.898/8.214 ± 0.990 | 0.745 ± 0.042/0.667 ± 0.033 | 0.779 ± 0.029/0.724 ± 0.030 | 35/45 | 0.047 ± 0.029/0.080 ± 0.020 |
| Subgroup | NA | HO | HE | PAS | FIS |
|---|---|---|---|---|---|
| NYC | 6.444 ± 0.626/5.286 ± 0.384 | 0.783 ± 0.040/0.569 ± 0.047 | 0.715 ± 0.037/0.634 ± 0.032 | 0/4 | −0.103 ± 0.046/0.110 ± 0.055 |
| Non-NYC | 10.333 ± 0.898/8.214 ± 0.990 | 0.745 ± 0.042/0.667 ± 0.033 | 0.779 ± 0.029/0.724 ± 0.030 | 35/45 | 0.047 ± 0.029/0.080 ± 0.020 |
Metrics include mean number of alleles (NA), observed heterozygosity (HO), expected heterozygosity (HE), number of PAS and inbreeding coefficient (FIS)
Rarefaction curves for cumulative alleles detected per sample size in NYC and non-NYC coyotes genotyped at (A) nine neutral and (B) 14 MHC-linked microsatellite loci. Even when sample sizes were equal, non-NYC coyotes had greater allelic diversity than NYC coyotes.
Genetic diversity at MHC-linked microsatellite markers
Across multiple diversity metrics, NYC coyotes exhibited reduced immunogenetic diversity compared with non-NYC coyotes (Table 1; Supplementary Table S7). Mean values of observed heterozygosity were significantly lower in NYC versus non-NYC coyotes (HO, non-NYC = 0.667 ± 0.033, NYC = 0.569 ± 0.047; two-tailed paired t-test, t = 3.646, P = 0.003). This trend was reflected by 12 of 14 MHC-linked microsatellite loci (Supplementary Table S7).
Reduced immunogenetic diversity was also evident in allelic richness and number of PAS. When comparing NYC versus non-NYC coyotes, the average number of alleles per locus was significantly lower in the NYC subgroup (NA, non-NYC = 8.214 ± 0.990, NYC = 5.286 ± 0.384; two-tailed paired t-test, t = 3.685, P = 0.003). The NYC population also possessed fewer unique alleles (4 PAS) than the non-NYC population (45 PAS). Rarefaction analyses further suggested that NYC coyotes possessed lower allelic diversity than non-NYC coyotes even when sample sizes were equivalent (Fig. 2B).
Population differentiation and substructure
Results from LDNe imply an effective population size of 42.6 individuals for NYC and 1884.4 individuals for the non-NYC subgroup. The G-test for genotypic differentiation returned a significant difference between the distribution of alleles at neutral loci in non-NYC and NYC subgroups (χ2 = 63.881, df = 18, P < 0.001), with evidence of genetic structure (FST=0.035). Using this value, we calculated the number of genetic migrants (NM) as 6.89 between NYC and non-NYC coyotes to estimate gene flow between super-urban coyotes and their neighbors, a value similar to the number of NYC coyotes with partial STRUCTURE assignment to the non-urban cluster (Fig. 1E). Dest calculated for neutral loci (pairwise Dest = 0.141) was higher than FST (Table 2).
Pairwise FST and Dest values calculated between NYC and non-NYC coyotes with neutral and MHC-linked loci
| Neutral | MHC-linked | |
|---|---|---|
| FST | 0.035* | 0.052* |
| Dest | 0.141* | 0.130* |
| Neutral | MHC-linked | |
|---|---|---|
| FST | 0.035* | 0.052* |
| Dest | 0.141* | 0.130* |
Asterisks indicate significant genotypic differentiation (P < 0.05) as determined by the G-test for genotypic differentiation for FST and bootstrap resampling for Dest.
Pairwise FST and Dest values calculated between NYC and non-NYC coyotes with neutral and MHC-linked loci
| Neutral | MHC-linked | |
|---|---|---|
| FST | 0.035* | 0.052* |
| Dest | 0.141* | 0.130* |
| Neutral | MHC-linked | |
|---|---|---|
| FST | 0.035* | 0.052* |
| Dest | 0.141* | 0.130* |
Asterisks indicate significant genotypic differentiation (P < 0.05) as determined by the G-test for genotypic differentiation for FST and bootstrap resampling for Dest.
Significant differentiation was also observed between non-NYC and NYC coyotes at MHC-linked microsatellite loci (χ2 > 151.674, df = 28, P < 0.001), with FST higher (pairwise FST=0.052) and number of genetic migrants (NM = 4.56) lower than those calculated at neutral loci. Once again, Dest (pairwise Dest = 0.130) was higher than the FST value calculated at immune-linked loci (Table 2).
Population bottleneck and selection tests
Tests for heterozygosity excess implemented in BOTTLENECK under the two-phase model of mutation reported greater excess in NYC coyotes (Wilcoxon test, P = 0.064; HE = 0.715, HO = 0.783) than non-NYC coyotes (Wilcoxon test, P = 0.150; HE = 0.779, HO = 0.745). Though these tests were above the commonly used significance threshold of P = 0.05, we considered values up to P = 0.10 potentially biologically relevant and worthy of consideration. These results were supported by inbreeding coefficients calculated for NYC (neutral FIS = −0.103; MHC-linked FIS = 0.110) and non-NYC (neutral FIS =0.047; MHC-linked FIS = 0.080) coyotes. Low FIS at neutral markers matched our expectation for heterozygosity excess in the recently bottlenecked population, and higher FIS at MHC-linked markers was consistent with decreased genetic diversity observed at immune-linked loci.
To examine whether selection was acting on MHC-linked loci, we implemented the Ewens–Watterson homozygosity test of neutrality in Pypop. One locus in non-NYC coyotes returned a significantly negative value (CFA12-17, FND = −1.399, P = 0.007). In the remaining MHC-linked loci, FND was not significantly positive or negative, though nearly all FND values were negative (13/14 loci for non-NYC coyotes; 11/14 loci for NYC coyotes; Fig. 3; Supplementary Table S8). Notably, four loci exhibited FND values of different signs (i.e. one positive and one negative) for NYC and non-NYC coyotes, though these were not statistically significant (Fig. 3; Supplementary Table S8).
Values for the normalized deviate of homozygosity (FND) calculated for NYC and non-NYC coyotes at MHC-linked microsatellite loci. Negative FND values suggest balancing selection, and positive FND values suggest directional selection. Significant values are denoted with an asterisk.
Discussion
Our results suggested that stochastic processes, such as random genetic drift, primarily drove patterns of reduced genetic diversity in super-urban coyotes. Among neutral loci, we found evidence of a population bottleneck in lower allelic richness in NYC versus non-NYC coyotes despite similar values of observed heterozygosity in both groups. Out BOTTLENECK results confirmed this pattern, as NYC coyotes exhibited greater heterozygosity excess (Wilcoxon test, P = 0.064; HE = 0.715, HO = 0.783) than non-NYC coyotes (Wilcoxon test, P = 0.150; HE = 0.779, HO = 0.745). We further observed inbreeding coefficients (FIS) in NYC coyotes consistent with heterozygosity excess at neutral loci.
Evidence of a population bottleneck was not surprising, as colonization of natural and urban environments is often characterized by reduced genetic diversity and small effective sizes (as seen in our NYC coyotes). For example, gray wolves (Canis lupus) colonizing the Alps exhibited lower levels of genetic variation than wolves from their source population in the Apennines (Fabbri et al. 2007). Similarly, red foxes colonizing the urban center of Zurich, Switzerland had reduced diversity and significant allelic differentiation when compared with their nonurban counterparts (Wandeler et al. 2003; DeCandia et al. 2019). This result therefore matched our expectation for reduced neutral diversity in NYC coyotes as a result of random genetic drift, though we currently see this reduction in allelic richness only.
At MHC-linked loci, we observed reduced genetic variation in NYC coyotes across multiple diversity metrics. Allelic richness, observed heterozygosity and number of PAS were significantly lower in NYC coyotes, with FIS calculated at MHC-linked loci higher in the super-urban population. This was counter to our original hypothesis that balancing selection would maintain high levels of immunogenetic variation, as seen in Channel Island foxes (Aguilar et al. 2004) and Los Angeles bobcats (Serieys et al. 2015a).
Despite diversity losses, we did observe some evidence for selection in the Ewens–Watterson homozygosity test of neutrality. Almost all FND values for NYC and non-NYC populations were negative across MHC-linked loci. This indicated maintenance of diversity (presumably through balancing selection) in both subgroups, as observed homozygosity was lower than expected homozygosity at most loci. Notably, four loci exhibited different signs (i.e. positive vs. negative) for FND values calculated in NYC versus non-NYC coyotes. This suggested that the processes shaping genetic diversity in each subgroup might differ. Further, positive FND values at three loci in NYC coyotes suggested that directional selection may be acting on certain loci, possibly contributing to the decreased variation observed at these markers.
The NYC colonization may be too recent to definitively detect selection, as these relationships were not statistically significant. Whereas foxes have inhabited the Channel Islands for thousands of generations (Aguilar et al. 2004), coyotes have only been documented in NYC since the mid-1990s (Nagy et al. 2016). Consequently, the overall trend of reduced immunogenetic diversity in super-urban coyotes may have more likely resulted from changes in population size during NYC’s recent colonization. This result is consistent with bottleneck-induced reductions in immunogenetic variation observed across other taxa (New Zealand robins, Miller & Lambert 2004; prairie-chickens, Bollmer et al. 2011; meta-analysis including amphibians, birds, fish, mammals and reptiles, Sutton et al. 2011; golden snub-nosed monkeys, Luo et al. 2012).
Selective pressure may additionally be relatively weak, as there has been little evidence of disease-induced mortality within NYC coyotes as of yet. This phenomenon may occur more broadly, as well, given relatively few cases of population-level effects of disease on coyotes (Pence et al. 1983; Pence and Windberg 1994; Gehrt, Riley, and Cypher 2010). This contrasts with Los Angeles bobcats, which experienced a severe mange outbreak associated with anticoagulant rodenticide (AR) exposure between 2002 and 2005 (Serieys et al. 2015a,b).
It is also possible that the small size of the super-urban coyote population renders drift a more likely and efficient evolutionary force than selection. Future studies may find indications of balancing or directional selection once more generations of coyotes have made NYC their home, or after the population has expanded into additional areas within the urban matrix.
Regarding population substructure, we observed significant allelic differentiation between NYC and non-NYC coyotes at both neutral and functional loci. Our FST estimates (0.035 for neutral loci and 0.052 for MHC-linked loci) were of a similar magnitude to those reported for other urban canids. For example, pairwise FST calculated for Swiss foxes ranged from 0.027 to 0.054 when comparing rural and urban populations (Wandeler et al. 2003). Similarly, coyotes separated by a freeway in California had FST values ranging from 0.030 to 0.037 (Riley et al. 2006). Further, our estimated number of genetic migrants between NYC and non-NYC groups at neutral loci (NM = 6.89) matched the Californian coyotes almost exactly (NM = 6.5 and 8.1; Riley et al. 2006). Given the higher FST at MHC-linked loci, our estimate for genetic migrants was lower (NM = 4.56). Taken together, these values suggested that NYC coyotes maintain gene flow with non-NYC coyotes, but may exhibit more restricted movement, perhaps a result of increased habitat fragmentation and difficulty navigating the urban matrix (Johnson and Munshi-South 2017). Estimates of population structure were even higher when calculating Dest (neutral Dest = 0.141; MHC-linked Dest = 0.130), which is often considered a better measure of differentiation for highly polymorphic loci. These results were further supported by the cluster assignments in STRUCTURE, as NYC coyotes repeatedly emerged as their own cluster, with admixture present with the non-NYC cluster(s).
In the short term, decreased diversity at neutral and immune-linked loci may not adversely affect NYC coyotes. Following initial colonization, population bottlenecks can be mitigated by increasing effective size and gene flow (Nei, Maruyama, and Chakraborty 1975). Further, urban environments may create more homogenous disease landscapes. Pathogen prevalence and risk of transmission may elevate with higher host density and increased contact with other urban wildlife (some of which serve as disease reservoirs), but pathogen species richness may be lower in urban environments (Gehrt, Riley, and Cypher 2010; Watts and Alexander 2012; Brearley et al. 2013; Grigione et al. 2014).
An additional threat exists in urban areas, however. Increased exposure to pollutants and pesticides may impose selective pressures or alter disease dynamics near cities. As mentioned previously, ARs have been linked to notoedric mange susceptibility in bobcats sampled near Los Angeles (Serieys et al. 2015b). In this system, both urbanization and persistent, sublethal AR exposure were associated with numerous immune anomalies that may have contributed to population-level declines during a recent mange outbreak (Fraser et al. 2018; Serieys et al. 2018). Similar mechanisms may alter immunocompetence in other urban wildlife species. Toxicant exposure may even cause mortality, as observed in coyotes exposed to ARs (Riley et al. 2003; Way et al. 2006; Gehrt, Riley, and Cypher 2010).
Outside of the realm of toxicant exposure, negative genetic effects may also arise if small population size and restricted gene flow persist. In these scenarios, inbreeding can result in further differentiation and loss of genetic variation in NYC coyotes, thereby reducing the population’s ability to respond to novel environmental pressures, such as disease (Keller and Waller 2002; Spielman et al. 2004; Frankham 2005; Oakley 2013).
To assess the relative risk of these negative effects, we recommend further monitoring of the genetic health of NYC coyotes alongside comparative analyses of pathogen and toxicant exposure in coyotes throughout the New York metropolitan area. These analyses will improve our understanding of NYC coyotes, while also informing management of this super-urban population should disease or negative genetic effects become evident and increased gene flow become necessary.
On a broader scale, our results contribute to the growing body of literature surrounding urban colonization ecology and are consistent with the genetic effects reported in other systems (Wandeler et al. 2003; Riley et al. 2006; Evans et al. 2009; Whittaker et al. 2012; Serieys et al. 2015a; DeCandia et al. 2019). They additionally highlight the importance of examining both neutral and functional genetic variation to better understand the contribution of genetic drift and natural selection in shaping diversity in newly colonized populations. As urbanization continues to grow, more species will live alongside humans in heavily modified areas. By combining insights from numerous case studies of urban colonization (such as NYC coyotes), we can anticipate potential issues (such as increased disease risk or inbreeding) before they arise. This may prove especially critical for small populations of rare or threatened species, and may ultimately necessitate corridors connecting urban green spaces with nearby nonurban areas. As urbanization transforms the human–wildlife interface worldwide, we must continue collecting data on the effects our activities have on local wildlife to ensure peaceful coexistence in perpetuity.
Acknowledgments
We thank the New York State Museum (NYSM; Albany, NY) and Gotham Coyote Project (GCP; New York City, New York) for providing samples necessary for this study.
Funding
Funding for field work and sample collection in New York City was provided by Rusticus Garden Club, Bedford Garden Club, Dorr Foundation, Louis Calder Center, Clare Boothe Luce Graduate Fellowship Program and the Mianus River Gorge Research Assistantship Program. We thank Joseph Bopp (NYSM) for his help curating samples collected in New York, New Jersey and Connecticut. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1656466.
Data availability
Microsatellite genotypes for MHC-linked and neutral markers are available through the Dryad Digital Repository: doi:10.5061/dryad.c0282c8.
Conflict of interest statement. None declared.
References
ESRI. (
R Core Team. (
U.S. Census Bureau. (
Author notes
Alexandra L. DeCandia and Carol S. Henger authors contributed equally.
