An important consequence of the discontinuous distribution of insect populations within their geographic range is phenotypic divergence. Detection of this divergence can be challenging when it occurs through subtle shifts in morphological traits with complex geometries, such as insect wing venation. Here, we used landmark-based wing geometric morphometrics to investigate the population-level phenotypic variation of the two subspecies of Glossina morsitans, G. m. centralis Machado and G. m. morsitans Westwood that occur in Zambia. Twelve homologous landmarks digitized on the right wings of 720 specimens collected from four and five sites (80 per site with 1:1 sex ratio) within the G. m. centralis and G. m. morsitans range respectively, were subjected to generalised Procrustes analysis to obtain wing centroid size (CS) and wing shape variable. Linear permutation models were then used to compare CS and wing shape between male and female G. morsitans, the two subspecies G. m. centralis and G. m. morsitans, and between sample locations within each subspecies range. Significant differences in CS and wing shape were observed between G. morsitans sexes, subspecies and sample locations within each subspecies range. Neighbour-joining trees derived from the analysis of Procrustes distances showed that tsetse within each subspecies range were highly divergent. We conclude that G. morsitans populations in Zambia exhibit significant population-level variation in fly size and wing shape which suggests high levels of population structuring. The main drivers of this structuring could be adaptation to local climatic condition, especially temperature and random genetic drift. We therefore recommend molecular studies to estimate the levels of gene flow between these populations and determine their levels of genetic isolation.