Species delimitation in tree species is notoriously challenging due to shared polymorphisms among species. An integrative survey that considers multiple operational criteria is a possible solution, and we aimed to test it in a species complex of aspens in China. Genetic [four chloroplast DNA (cpDNA) fragments and 14 nuclear microsatellite loci (nSSR)] and morphological variations were collected for 76 populations and 53 populations, respectively, covering the major geographic distribution of the Populus davidiana-rotundifolia complex. Bayesian clustering, analysis of molecular variance (AMOVA), Principle Coordinate Analysis (PCoA), ecological niche modeling (ENM), and gene flow (migrants per generation), were employed to detect and test genetic clustering, morphological and habitat differentiation, and gene flow between/among putative species. The nSSR data and ENM suggested that there are two separately evolving meta-population lineages that correspond to P. davidiana (pd) and P. rotundifolia (pr). Furthermore, several lines of evidence supported a subdivision of P. davidiana into Northeastern (NEC) and Central-North (CNC) groups, yet they are still functioning as one species. CpDNA data revealed that five haplotype clades formed a pattern of [pdNEC, ((pdCNC, pr), (pdCNC, pr))], but most haplotypes are species-specific. Meanwhile, PCA based on morphology suggested a closer relationship between the CNC group (P. davidiana) and P. rontundifolia. Discrepancy of nSSR and ENM vs. cpDNA and morphology could have reflected a complex lineage divergence and convergence history. P. davidiana and P. rotundifolia can be regarded as a recently diverged species pair that experienced parapatric speciation due to ecological differentiation in the face of gene flow. Our findings highlight the importance of integrative surveys at population level, as we have undertaken, is an important approach to detect the boundary of a group of species that have experienced complex evolutionary history.