Geometric morphometric analysis was combined with two different unsupervised machine learning algorithms, UMAP and HDBSCAN, to visualize morphological differences in wing shape among and within four Anopheles sibling species (An. atroparvus, An. melanoon, An. maculipennis s.s. and An. daciae sp. inq.) of the Maculipennis complex in Northern Italy. Specifically, we evaluated: (1) wing shape variation among and within species; (2) the consistencies between groups of An. maculipennis s.s. and An. daciae sp. inq. identified based on COI sequences and wing shape variability; and (3) the spatial and temporal distribution of different morphotypes. UMAP detected at least 13 main patterns of variation in wing shape among the four analyzed species and mapped intraspecific morphological variations. The relationship between the most abundant COI haplotypes of An. daciae sp. inq. and shape ordination/variation was not significant. However, morphological variation within haplotypes was reported. HDBSCAN also recognized different clusters of morphotypes within An. daciae sp. inq. (12) and An. maculipennis s.s. (4). All morphotypes shared a similar pattern of variation in the subcostal vein, in the anal vein and in the radio-medial cross-vein of the wing. On the contrary, the marginal part of the wings remained unchanged in all clusters of both species. Any spatial-temporal significant difference was observed in the frequency of the identified morphotypes. Our study demonstrated that machine learning algorithms are a useful tool combined with geometric morphometrics and suggest to deepen the analysis of inter and intra specific shape variability to evaluate evolutionary constrains related to wing functionality.