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 digitised 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 variables. Linear permutation models and redundancy analysis 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, the sexes of each subspecies and between sample locations within each subspecies range, respectively. Significant differences in CS and wing shape were observed between G. morsitans sexes, subspecies and sample locations within each subspecies range. A neighbour-joining cladogram 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 random genetic drift in G. m. centralis demes and local adaptation to environmental conditions in G. m. morsitans populations. We therefore recommend molecular studies to estimate the levels of gene flow between these populations and identify possible barriers to genetic flow.