Understanding the relative roles of geography and ecology in driving speciation, population divergence and maintenance of species cohesion is of great interest to molecular ecology. Closely related species that are parapatricly distributed in mountainous areas provide an ideal model to evaluate these key issues, especially when genomic data are analyzed within a spatially and ecologically explicit context. Here we used three closely related species of Primula that occur in the Himalayas, the Hengduan Mountains and Northeast Qinghai-Tibet Plateau (QTP) to examine the effects of geography and ecology on interspecific divergence and maintenance of species cohesion. We used genomic data for 770 samples of the three species using restriction site-associated DNA (RAD) sequencing and combined approximate Bayesian computation (ABC) modeling, Bayesian generalized linear mixed modeling (GLMM) and species distribution modeling (SDM). The three species are clearly delimited by the RADseq data. Further ABC modeling indicates that the three species originated in the Himalayas and diverged from each other following the uplifts of the Hengduan Mountains and the Northern QTP during the Pliocene. After a long period of divergence, the three species came into secondary contact triggered by past climatic changes but with no significant introgression. The three species display complex and different drivers of genomic variation, which provides further insights into the effects of geographical and ecological factors on maintaining species cohesion. Our findings highlight the significance of combining the use of population genomics with environmental data when evaluating the effects of geography and ecology on interspecific divergence and maintenance of closely related species