Disentangling the processes underlying geographic and environmental patterns of biodiversity challenges biologists as such patterns emerge from eco-evolutionary processes confounded by spatial autocorrelation among sample units. The herbivorous insect, Belonocnema treatae (Hymenoptera: Cynipidae), exhibits regional specialization on three plant species whose geographic distributions range from sympatry through allopatry across the southern USA. Using range-wide sampling spanning the geographic ranges of the three host plants and genotyping-by-sequencing of 1,217 individuals, we tested whether this insect herbivore exhibited host-plant-associated genomic differentiation while controlling for spatial autocorrelation among the 58 sample sites. Population genomic structure based on 40,699 SNPs was evaluated using the hierarchical Bayesian model ENTROPY to assign individuals to genetic clusters and estimate admixture proportions. To control for spatial autocorrelation, distance-based Moran’s eigenvector mapping was used to construct regression variables summarizing spatial structure inherent among sample sites. Distance based redundancy analysis (dbRDA) incorporating the spatial variables was then applied to partition host-plant-associated differentiation (HAD) from spatial autocorrelation. By combining ENTROPY and dbRDA to analyze SNP data we unveiled a complex mosaic of highly structured differentiation within and among gall former populations finding evidence that geography, HAD and spatial autocorrelation all play significant roles in explaining patterns of genomic differentiation in B. treatae. While dbRDA confirmed host association as a significant predictor of patterns of genomic variation, spatial autocorrelation among sites explained the largest proportion of variation. Our results demonstrate the value of combining dbRDA with hierarchical structural analyses to partition spatial/environmental patterns of genomic variation.