Characterizing genetic variation in natural populations is central to evolutionary biology. However, most non-model organisms lack integral genomic resources such as reference genomes and recombination maps, limiting accurate evolutionary inference. Here, we explore the consequences of reference genome bias on the inference of genetic diversity, demographic histories, and recombination rates using gray foxes (Urocyon cinereoargenteus), which, like most members of Canidae, are traditionally mapped to the dog (Canis lupus familiaris) reference genome. Whole genome sequence data from gray foxes were mapped to the gray fox reference genome and two heterospecific canid references (dog and Arctic fox; Vulpes lagopus). Our results reveal that reference bias significantly affects population genomic analyses. Mapping to the conspecific gray fox genome improved read pairing, increased detection of SNPs, especially rare variants, and reduced spurious variants. Estimates of nucleotide diversity (π) and genetic differentiation (FST) were higher using the gray fox genome. We observed that mapping to heterospecific references leads to underestimates of population sizes, distorted demographic trajectories, and more variable recombination rates. These effects are further complicated by population-specific biases, which vary in their magnitude and direction across populations, highlighting the need for tailored approaches to mitigate reference bias. Importantly, FST outlier detection also differed among references, affecting functional interpretations. Collectively, this work addresses a critical gap in the rapidly expanding field of non-model species genomics by demonstrating the importance of using conspecific genomic resources in evolutionary research and illustrating how reliance on distantly related reference genomes can distort population genetic analyses.
