Does anyone have any guidance or suggestions for variant calling RADseq and WGS datasets to generate a single population genomic dataset?

Both datasets represent population-level samples for the same species. Our current workflow uses a highly complete/contiguous reference genome and BWA to call variants. However, even after stringently filtering the variants for depth and missing data neighbor-joining trees and other population genetic analyses show significant genetic distance between the WGS and RADseq datasets (see image).

Neighbor joining tree generated using bitwise distances: