Inferring evolutionary relationships among recently diverged lineages is necessary to understand how isolating barriers produce independent lineages. Here we investigate the phylogenetic relationships between three incompletely isolated and closely related mushroom-feeding Drosophila species. These species form the Drosophila subquinaria species complex and consist of one Eurasian species (D. transversa) and two widespread North American species (D. subquinaria and D. recens) that are sympatric in central Canada. While patterns of pre- and post-mating isolation among these species are well characterized, previous work on their phylogenetic relationships is limited and conflicting. In this study, we generated a multi-locus dataset of 29 loci from across the genome sequenced in a population sample from each species, and then we inferred species relationships and patterns of introgression. We find strong statistical support that D. subquinaria is paraphyletic, showing that samples from the geographic region sympatric with D. recens are most closely related to D. recens whereas samples from the geographic region allopatric with D. recens are most closely related to D. transversa. We present several lines of evidence that both incomplete lineage sorting and gene flow are causing phylogenetic discordance. We suggest that ongoing gene flow primarily from D. recens into D. subquinaria in the sympatric part of their ranges causes phylogenetic uncertainty in the evolutionary history of these species. Our results highlight how population genetic data can be used to disentangle the sources of phylogenetic discordance among closely related species.