The cactophilic fly Drosophila mojavensis exhibits considerable intraspecific genetic structure across allopatric geographic regions and shows associations with different host cactus species across its range. The divergence between these populations has been studied for more than 60 years, yet their exact historical relationships have not been resolved. We analysed sequence data from 15 intronic X-linked loci across populations from Baja California, mainland Sonora-Arizona and Mojave Desert regions under an isolation-with-migration model to assess multiple scenarios of divergence. We also compared the results with a pre-existing sequence dataset of 8 autosomal loci. We derived a population tree with Baja California placed at its base and link their isolation to Pleistocene climatic oscillations. Our estimates suggest the Baja California population diverged from an ancestral Mojave Desert/mainland Sonora-Arizona group around 230-270 Kya, while the split between the Mojave Desert and mainland Sonora-Arizona populations occurred one glacial cycle later, 117-135 Kya years ago. Although we found these three populations to be effectively allopatric, model ranking could not rule out the possibility of a low-level of gene flow between two of them. Finally, the Mojave Desert population showed a small effective population size, consistent with a historical population bottleneck. We show that model-based inference from multiple loci can provide accurate information on the historical relationships of closely related groups allowing us to set into historical context a classic system of incipient ecological speciation.