wings were removed if they were heavily sclerotized and prevented the wings from laying flat. We plotted 9 homologous wing venation landmarks (following Rattanawannee et al., 2015) onto each wing image using tpsDig software version 2.31 (Rohlf, 2015); (Figure 3). Cartesian landmark data was exported as a TPS file and analyzed in R version 4.2.2.

Data analysis

We Procrustes-aligned landmark coordinates using R package ‘geomorph’ version 4.0.0 (Adams et al., 2022; Baken et al., 2021). To test for statistical differences between the two H. tripartitus populations and among the three species, we ran one-way multivariate analysis of variance (MANOVA) tests using R package ‘RRPP’ version 1.3.1 (Collyer & Adams, 2018, 2019). To visualize separation among groups, we generated density plots with discriminant analysis of principal components (DAPC) using R package ‘adegenet’ version 2.1.10 (Jombart, 2008; Jombart & Ahmed, 2011). To test the accuracy of using wing landmarks to predict an unknown bee’s species or population, we utilized DAPC cross-validation. Cross-validation also informed the number of principal components (PCs) retained in each analysis, which is non-trivial (Jombart & Collins, 2015)