DOI: 10.5061/dryad.tx95x6b50
title: Maintenance of a narrow hybrid zone between native and introduced red foxes (Vulpes vulpes) despite conspecificity and high dispersal capabilities

description: Human-facilitated introductions of nonnative populations can lead to
 secondary contact between previously allopatric lineages, resulting in
 either homogenization of the lineages or stable hybrid zones that are
 maintained by pre-zygotic (e.g., behavioral) or post-zygotic (e.g.,
 reduced hybrid fitness) reproductive barriers. We investigated patterns of
 gene flow between the native Sacramento Valley red fox (Vulpes vulpes
 patwin) and an introduced conspecific population of captive-bred
 (fur-farm) origin in California’s Central Valley. Considering their recent
 divergence (i.e., ~50 kya), we hypothesized that pre-zygotic mechanisms
 primarily impede gene flow, rather than post-zygotic barriers.
 Additionally, some genes originating in nonnative foxes may confer higher
 fitness in the currently human-dominated landscape resulting in selective
 introgression into the native population. Genetic analysis of 682 red
 foxes (255 native, 427 nonnative) at both mitochondrial (cytB + Dloop) and
 nuclear loci (~19,000 SNPs) revealed significantly narrower cline widths
 than expected under a simulated model of unrestricted gene flow,
 consistent with the existence of pre- or post-zygotic reproductive
 barriers. We identified several loci with reduced introgression linked to
 behavioral divergence in captive bred foxes, which supports pre-zygotic
 mechanisms as a putative driver of the narrow hybrid zone. Additionally,
 several loci with elevated gene flow from the nonnative into the native
 population, were near genes associated with adaptation to human dominated
 landscapes.  Overall, this study contributes to our understanding
 of hybrid zone dynamics in vertebrates, particularly in the context of
 species introductions and landscape changes, underscoring the importance
 of considering multiple mechanisms that may be at play in maintaining
 lineages at both the species and subspecies level.
description: We obtained both scat (n = 195) and tissue (n = 490) samples as
 well as 4 hair and swab samples from red foxes in the native SVRF (n =
 293) and nonnative ranges (n = 396). Samples were collected between
 1982–2021 and many of them (n = 472) were first described in previous
 studies. Samples were obtained through non-invasive den and transect
 surveys, roadkill collection, live captures (following American Society of
 Mammalogy animal care guidelines and with University of California, Davis
 Animal Care and Use Committee approval, IACUC No.17860), and nonnative
 species removal programs in the case of the introduced red
 foxes.
description: # Maintenance of a narrow hybrid zone between native and introduced red
 foxes (Vulpes vulpes) despite conspecificity and high dispersal
 capabilities This dataset was used to generate results for the following
 manuscript which has been accepted in the Journal of Molecular Ecology
 (May 13th, 2024) ## Description of the data and file structure This
 repository contains data and scripts required to recreate analyses in the
 manuscript *Maintenance of a narrow hybrid zone between native and
 introduced red foxes (Vulpes vulpes) despite conspecificity and high
 dispersal capabilities *which has been accepted into the Journal of
 Molecular Ecology as of May 13th, 2024. The raw GBS sequencing data
 associated with this study can be found in the NCBI Sequence Read Archive
 under project number: **PRJNA1069152** : Reduced-representation
 genotyping-by-sequencing of low-elevation red foxes in California
 **Maintenance of a narrow hybrid zone between native and introduced red
 foxes despite conspecificity and high dispersal capabilities** ## **Sample
 Info** *
 [SampleList.md](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/SampleList.md) * **Table S1.** Information for red fox DNA samples from the Sacramento Valley, San Joaquin Valley, and Coastal Lowlands of California, USA, used in this study, including sample type (tissue, scat, etc.), geographic location (latitude, longitude), year collected, mitochondrial DNA haplotype and ancestry (native vs nonnative), whether samples were included in nuclear genotyping-by- sequencing (GBS) and whether they passed basic quality filtering (GBS Success), geographic range delineation, and citations for samples that were previously described in other studies. ## Raw data location The raw GBS sequencing data associated with this study can be found in the NCBI Sequence Read Archive under project number: **PRJNA1069152** : Reduced-representation genotyping-by-sequencing of low-elevation red foxes in California ## **01. Data Processing** * Demultiplex Data [demux.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/01.Data-Processing/demux.sh) * Trim reads [trim.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/01.Data-Processing/trim.sh) * Merge Duplicate Runs [merge_dup.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/01.Data-Processing/merge_dup.sh) ## 02. Alignment & Initial SNP Calling * Align to reference [align.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/02.Alignment-SNPcalling/align.sh) * Stacks reference map SNP pipeline [refmap_PL_redfox.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/02.Alignment-SNPcalling/refmap_PL_redfox.sh) * Run populations on output of above refmap pipeline [poprun_redfox.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/02.Alignment-SNPcalling/poprun_redfox.sh) ## 03. Data filtering You can check out details on the Stacks populations filters [here](http://catchenlab.life.illinois.edu/stacks/comp/populations.php) * Using data generated from the reference genome pipeline described above, we tested a variety of filtering approaches to maximimze the number of sites and individuals while removing low quality data. The pipeline here represents the approach described in the manuscript. * Remove the lowest quality individuals (&gt;75% of SNPs missing) [plinkfilter_mind0.75.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/03.SNP-filtering/plinkfilter_mind0.75.sh) * Remove sites that were called in &lt;80% of individuals [plinkfilter_geno0.2.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/03.SNP-filtering/plinkfilter_geno0.2.sh) * Again remove lower quality individuals (&gt;20% of SNPs missing) using a threshold based approach [plinkfilter_mind0.2.sh](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/03.SNP-filtering/plinkfilter_mind0.2.sh) ## 04. Sex-Biased Dispersal See [README.md](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/tree/main/04.Sex-Biased-Dispersal-Patterns) ## 05. Geographic Patterns of Ancestry See [README.md](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/tree/main/06.Hybrid-Zone-Width) * MDS Plot * Admixture (K=1-5) ## 06. Width of the Hybrid Zone ### Expected Cline Width (assuming neutral gene glow) In order to run this you will also need: * The locations file [HZAM_HZAR_location_distances.csv](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/HZAM_HZAR_location_distances.csv). 1. Simulate gene flow of maternally inherited mtDNA across the landscape using HZAM [SimulatedGeneFlow_mtDNA_Disp100%_100sims.R](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/SimulatedGeneFlow_mtDNA_Disp100%25_100sims.R) * Assumes sex-biased dispersal * Assumes 18 generations (conservative estimate) of gene flow post contact 1. Simulate gene flow of biparentally inhertied nuDNA across the landscape using HZAM [SimulatedGeneFlow_nuDNA_Disp100%_100sims.R](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/SimulatedGeneFlow_nuDNA_Disp100%25_100sims.R) * Assumes sex-biased dispersal * Assumes 18 generations (conservative estimate) of gene flow post contact 1. Estimate cline widths for both mtDNA and nuDNA using HZAR (simulated data) * [SimulatedClineWidths.R](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/SimulatedClineWidths.R) ### Observed Cline Width (using empirical data) In order to run this you will need: * The empirical cline data for mtDNA: [GeogClineData-mtDNA.csv](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/GeogClineData-mtDNA.csv) * The empirical cline data for nuDNA: [GeogClineData-nuDNA.csv](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/GeogClineData-nuDNA.csv) * Sampling Intervals for [mtDNA](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/GeogCline-SamplingIntervals-mtDNA.csv) and [nuDNA](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/GeogCline-SamplingIntervals-nuDNA.csv) 1. Generate summary tables of the empirical data for use in HZAR cline models [EmpiricalData_SummaryTables.R](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/EmpiricalData_SummaryTables.R) 2. Estimate cline widths for both mtDNA and nuDNA using HZAR (empirical data) * [ObservedClineWidths.R](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/blob/main/06.Hybrid-Zone-Width/ObservedClineWidths.R) ## 07. Genomic Cline Analysis * Bayesian Genomic Cline - See [README.md](https://github.com/squisquater/Red-Fox-Hybrid-Zone-Dynamics/tree/main/07.Genomic-Cline)