A current goal of speciation research is to identify the loci underlying reproductive barriers between species. Locating such barrier loci in empirical data is difficult due to the often complex demographic history of diverged taxa and the heterogeneity in evolutionary forces across the genome. Here we take advantage of a natural case of hybridization between two wood ant species (Formica aquilonia and F. polyctena) to identify regions of reduced long-term gene flow using demographically explicit scans of non-admixed genomes. In addition we identify candidate Bateson-Dobzhansky-Muller incompatibilities (BDMIs) through an imbalanced recombinant haplotype frequency analysis of natural F. aquilonia x F. polyctena hybrid genomes. Both approaches find barriers that are scattered across the genome. Furthermore, candidate BDMIs significantly overlap with the long-term barriers identified by gIMble, indicating that incompatibilities have persisted despite divergence with gene flow between the wood ant species. Intriguingly, BDMIs interact in a network and the number of pairwise interactions a BDMI has correlates with its long-term barrier strength: hub-like BDMIs with many pairwise interactions reduce gene flow more effectively. Finally in regards to function, long-term barriers identified by gIMble arise outside regions of both gene coding sequences (CDS) and transposable elements. In contrast, regions where long-term barriers and BDMIs co-locate are significantly associated with introns, implying a potential role of alternative splicing or gene regulation in incompatibilities, rather than CDS divergence. Overall, our results highlight the underappreciated impact of multilocus BDMIs and the need to consider network connectivity of BDMIs in future work.