Using Center for Disease Control (CDC) light traps, sandflies were collected in 2011 in seven different geographical locations across Thailand (Fig. 1). Isolated islands and mainland location sites were selected in this study. Three isolated islands, including Lang Ga Jiew (LGJ), Ka (KA), and Mapraw (MAP) Islands, were selected for sandfly collection; these three islands are located in the Gulf of Thailand, in Chumphon province, southern Thailand. LGJ Island is a small offshore island that covers an area of 0.063 km² and is located approximately 8 km away from the shore. MAP and KA islands are located approximately 1 km away from the shore. Four other mainland location sites were selected: a small sample from an urban location in Nonthaburi province (NTB), which is located in central Thailand, and other from limestone caves of Kanchanaburi province (KCB, western Thailand), Ratchaburi province (RCB, central Thailand) and mountainous areas of Kamphaeng Phet province (KPP, northern Thailand).

Collected female sandflies were stored in 95% ethyl alcohol and mounted in Hoyer medium on glass microscopic slides. Morphological identification was based mainly on the keys of Lewis [9,10]. Subsequently, P. stantoni and S. hodgsoni were photographed using a Nikon DS-Ri1 SIGHT digital camera connected to stereo-microscope Nikon AZ 100 M (Nikon Corp., Tokyo, Japan) with a size scale provided on the picture. A total of 157 right wing pictures from P. stantoni and S. hodgsoni across different locations were processed for landmark-based GM analysis (Table 1).

Twelve landmarks were digitized (Fig. 2). When connected after a Procrustes superimposition onto the consensus configuration, they appear as polygons allowing an easier visual comparison of mean wing shape between sandfly species and populations. The Procrustes superposition was conducted according to the Generalized Procrustes Superimposition (GPA) procedure [11,12], generating both size (centroid size [CS]) and shape (partial warps [PW]) variables.

CS is defined as the square root of the sum of the squared distances between the center of the configuration of landmarks and each individual landmark [13]. Size comparisons were performed using a non-parametric permutation test (1000 cycles), and comparisons were illustrated by quantile boxes.

The shape variables (partial warps scores [PW]) were processed by standard multivariate analysis to statistically compare groups. The principal components (PCA) of the PW (i.e., the relative warps [RW]) were computed to visualize the morphospace for sandfly species from different locations. The discriminant analysis (DA) used as input variables the RW, or a few first of them according to the sample sizes, and was illustrated by the corresponding factor maps. For these DA, the two localities with very low sample sizes were removed: NTB (6 individuals) and KA (7 individuals). Validated reclassification scores were estimated after a DA on the two taxa, after the DA on the geographic populations within each species. In this validated reclassification, each individual was allocated to its closest group (Mahalanobis distance) without being used to help determine a group center [14]. The statistical significance (P ≤ 0.05) of shape variation was analyzed using non-parametric methods (1000 runs) applied to the pairwise Mahalanobis distances. A neighbor-joining tree was generated based on Procrustes distances between sandflies from different locations.

The allometric effect on the discrimination between collecting sites was explored by regressing the discriminant factors on centroid size and computing the coefficient of determination.

CLIC version 97 [15,16], freely available at http://mome-clic.com, was used to conduct the digitization of the images and subsequent multivariate analyses on the coordinate data. The following modules were used: COO for landmarks collection, MOG to generate centroid size and shape variables, as well as to compute Procrustes distances and VAR to compare the size between groups. The Procrustes distances were used to compute the neighbor-joining tree using the PHYLIP neighbor module [17], and the resulting tree was visualized using the NJPLOT software [18]. The module PAD of CLIC was used to perform discriminant analyzes between groups, as well as to compute Mahalanobis distances and the corresponding statistical significance.

Mean wing CS comparison of P. stantoni and S. hodgsoni from different geographical locations (Table 1) showed size variation (Fig. 3, Table 2). Mean wing CS was significantly (P < 0.05) larger in P. stantoni than in S. hodgsoni (1.22 mm and 0.91 mm, respectively). P. stantoni from mainland in Kamphaeng Phet province had the largest wings (mean, 1.28 mm), significantly different from P. stantoni from Lang Ga Jiew Island and an urban location in Nonthaburi province (P < 0.05). However, the mean wing CS of P. stantoni from caves on Lang Ga Jiew Island and Nonthaburi province were not significantly different (P > 0.05). The wing CS of S. hodgsoni from a mainland cave in Ratchaburi province was the largest (mean, 1.02 mm), whereas S. hodgsoni from Lang Ga Jiew Island had the smallest wing CS (0.89 mm). Almost all populations of S. hodgsoni were significantly different (P < 0.05) from each other. However, there was no significant (P > 0.05) difference of wing size between S. hodgsoni from Ka Island and Lang Ga Jiew Island, Ka Island and Kanchanaburi province, or Mapraw Island and Kanchanaburi province.

After superimposition of the mean landmark configurations onto the consensus one, P. stantoni and S. hodgsoni wing shapes were clearly different (Fig. 4). Polygons of mean wing shape of within species samples are shown in Fig. 5 (P. stantoni) and Fig. 6 (S. hodgsoni). Consistent intraspecific variation was visible in both sandfly species.

In agreement with the polygon superimpositions, the RW–based morphospace showed no overlapping of wing shapes between P. stantoni and S. hodgsoni, as well as between some geographic localities within S. hodgsoni (Fig. 7). Within each species, significant geographic differences (P < 0.05) could be found (Table 3). P. stantoni from Lang Ga Jiew Island and Kamphaeng Phet province were significantly different (P < 0.05), as were P. stantoni from Nonthaburi province and Kamphaeng Phet province (P < 0.05).

The factor map of the DA clearly showed the large shape divergence between taxa, as well as some non-overlapping groups within both species: P. stantoni was almost completely subdivided into island and mainland samples (Fig. 8), while in S. hodgsoni the trend was similar, showing divergence between mainland (KCB, RCB) and island samples (LGJ, MAP), and, within these areas, between sampling sites (KCB versus RCB, and LGJ versus MAP) (Fig. 9). Statistical significance was frequently found, except with localities NTB and KA, represented by only a few individuals (Table 3). For S. hodgsoni the contribution of size to shape-based discrimination ranged from 0.099 (with DF1) to 0.165 (with DF2), and for P. stantoni (a single DF since two populations were compared), it was 0.117.

The validated reclassification scores (Table 4) confirmed the complete separation of the two taxa, providing an average correct assignment of 79% between two samples in P. stantoni (scores ranging from 71 to 88%), and of 65% between four samples in S. hodgsoni (scores ranging from 63 to 71%).

Based on Procrustes distances, the neighbor-joining tree first separated the two genera and then, within species, continental and islands locations (Fig. 10).