Reconstructing the phylogeny of the pine family (Pinaceae) has historically been problematic because analyses based on morphological characters and genetic analyses using small numbers of sequences have produced incongruent results. Such discrepancies may have resulted from continuous gene flow or genome mixing events among certain pine species. Here, we reconstructed a species-level tree of the Pinaceae using genome-wide data from seven pine species. A multi-locus analysis of SNPs from homologous fragments indicated that Pseudotsuga diverged from Pinus and Picea 167.26–149.19 million years ago. Analysis using the isolation-with-migration model identified significant signals of two-way gene flow between Pinus taeda and Pinus elliottii (P<0.001), while the isolation-with-initial-migration model suggested that introgression did not continue after this divergence. Gene flow has only ever occurred in the initial period of species divergence. Further analyses highlighted complex historical genome mixing and hybridization events over the evolution of the pine family. Thus, the divergence and radiation of the seven pines was a complex process accompanied by introgression and historical admixture.