Wing morphometrics as a tool in species identification of forensically important blow flies of Thailand

Correct species identification of blow flies is a crucial step for understanding their biology, which can be used not only for designing fly control programs, but also to determine the minimum time since death. Identification techniques are usually based on morphological and molecular characters. However, the use of classical morphology requires experienced entomologists for correct identification; while molecular techniques rely on a sound laboratory expertise and remain ambiguous for certain taxa. Landmark-based geometric morphometric analysis of insect wings has been extensively applied in species identification. However, few wing morphometric analyses of blow fly species have been published. We applied a landmark-based geometric morphometric analysis of wings for species identification of 12 medically and forensically important blow fly species of Thailand. Nineteen landmarks of each right wing of 372 specimens were digitised. Variation in wing size and wing shape was analysed and evaluated for allometric effects. The latter confirmed the influence of size on the shape differences between species and sexes. Wing shape variation among genera and species were analysed using canonical variates analysis followed by a cross-validation test. Wing size was not suitable for species discrimination, whereas wing shape can be a useful tool to separate taxa on both, genus and species level depending on the analysed taxa. It appeared to be highly reliable, especially for classifying Chrysomya species, but less robust for a species discrimination in the genera Lucilia and Hemipyrellia. Allometry did not affect species separation but had an impact on sexual shape dimorphism. A landmark-based geometric morphometric analysis of wings is a useful additional method for species discrimination. It is a simple, reliable and inexpensive method, but it can be time-consuming locating the landmarks for a large scale study and requires non-damaged wings for analysis.


Background
Blow flies are considered to be of medical importance worldwide. Adults are mechanical vectors of several pathogens in humans, i.e. viruses, bacteria, protozoan cyst, helminth eggs, and fungi [1][2][3][4]. Their larvae are myiasis-producing agents in living humans and vertebrate animals, particularly the genera Cochliomyia, Chrysomya, Lucilia and Calliphora [5][6][7]. In addition, blow flies are forensically important insects as immature stages feed on human corpses and can be used in forensic investigations [8][9][10]. Forensic entomology is the analysis of insect evidence for forensic and legal purposes and is most frequently used for the estimation of the minimum time since death (PMI min ) [11].
Correct species identification of the blow fly is a crucial step in understanding its biology for not only designing fly control programs but also determining the PMI min precisely. Misidentification may impact the effectiveness of fly control strategies and bias the calculation of developmental times, eventually leading to an incorrect PMI min . Several identification techniques based on morphological [12,13] and molecular characters exist [14,15]. However, the use of classical morphology such as bristles on the body or male genitalia are very difficult to apply for non-experts and DNA identification can still remain ambiguous as, for example, some forensically important fly species are not, or insufficiently, represented in the reference libraries [16], many of the existing sequences in those online libraries just represent "dark taxa", i.e. they are not identified to species level [17], and DNA barcodes failed to distinguish among certain closely related species [14,[18][19][20][21]. This problem is even more serious in regions, where the accurate knowledge of relevant species is a challenge.
Besides classical morphology or DNA identification, the use of morphometrics has shown to be a valuable tool for interspecific discrimination. Morphometrics is defined as the quantitative studies of biological size and shape, shape variation, and its covariation with other biotic or abiotic factors [22] and can be a valuable tool for interspecific discrimination. In recent years, a landmarkbased geometric morphometric analysis of insect wings has been extensively applied in entomology, particularly in taxonomy [23][24][25][26] and geographic variation of species [27][28][29] due to its simplicity, low costs and high reliability. Several orders have been studied such as the Diptera [23,30], Hymenoptera [31], Coleoptera [26] and Odonata [32].
However, just a few wing morphometric analyses of blow fly species have been published [24,28,33], and the focus of such analysis for identification was mainly on medical and not in forensic entomology [28,34,35]. Therefore, the aim of the study is to test a landmarkbased geometric morphometric analysis of wings for species identification of medically and forensically important blow flies of Thailand, one of the global hotspots of biodiversity.

Specimen collection
Adult blow flies used in this study were collected during 2013-2014 in several locations of Thailand (Table 1, Fig. 1). A total of 372 blow flies were captured by using semi-automatic funnel trap invented by Kom Sukontason (Additional file 1), or by using a hand-held fly net, sweeping it over a bait of 1-day-old beef offal (~300 g). Specimens were either killed by ethyl acetate in the field or kept alive in the same fly net used to catch the flies and then transported back to the laboratory within 1 h and frozen at -20°C for 2 h. Specimens were pinned and identified based on the taxonomic key of Kurahashi and Bunchu [13]. Sampled specimens belong to 12 species including Chrysomya megacephala, Chrysomya chani, Chrysomya pinguis, Chrysomya rufifacies, Chrysomya villeneuvi, Chrysomya nigripes, Lucilia cuprina, Lucilia papuensis, Lucilia porphyrina, Lucilia sinensis, Hemipyrellia ligurriens and Hemipyrellia pulchra (Table 1), therefore including the most relevant taxa from medical and forensic point of views.

Slide preparation
The right wing of each fly was removed with fine forceps. A drop of Permount™ Mounting Medium was added on a microscope slide, and then the wing was placed onto the drop and covered with the cover slip. Before placing the wing onto the drop, the wing was submerged in xylene to facilitate the montage and avoid bubbles. All mounted wing slides were kept as thin as possible using the minimum mounting medium to maximise wing flattening and then dried at room temperature for a week. Each wing was photographed using a digital camera attached to a stereomicroscope at 1.5× magnification. Images were used to build tps files by using the TpsUtil V. 1.64 software [36] to minimise a possible bias when digitising the landmark locations. Altogether 19 landmarks as used by Hall et al. [28] (Fig. 2) were digitised using TpsDig2 V.2.20 software [37]. Each wing was digitised twice to reduce the measurement error [38].

Geometric morphometric analysis
The tps files were subjected to the MorphoJ software [39], and then the raw landmark coordinates of all specimens were aligned and superimposed using Procrustes Fit function to remove variation due to differences in scale, position and orientation from the coordinates. The centroid size (the square root of the sum of the squared distances between the centre of the configuration of landmarks and each landmark) [40] and Procrustes coordinates obtained from landmark data were used for further statistical analyses. For establishing a possible measurement error, the Procrustes coordinates of each specimen were averaged after a generalised Procrustes analysis in MorphoJ. Centroid size was also averaged for each specimen.

Size variation
Wing size was estimated by its centroid size [40]. Wing size difference among species was analysed using Kruskal-Wallis H-test, followed by Mann-Whitney U-test with Bonferroni correction applied for the significance level (0.05). Statistical analysis was performed in SPSS V.17.0 software for Windows (SPSS Inc., Chicago, Illinois, USA).

Shape variation
Wing shape variation was analysed using MorphoJ software [39]. Canonical variate analysis (CVA) was used to determine the most important feature as a possible discriminator between groups (genera or species). The statistical significance of pairwise differences in mean shapes was analysed using permutation tests (10,000 rounds) with Mahalanobis distances and Procrustes distances. Additionally, a cross-validation test in discriminant function analysis (DFA) was used to assess the accuracy of classification based on Mahalanobis distances in a permutation test with 10,000 rounds using MorphoJ software [39].

Sexual dimorphism
Sexual dimorphism consisted of sexual size dimorphism (SSD) and sexual shape dimorphism (SShD). Differences in size between sexes for each species were tested by Mann-Whitney U-test. Statistical analysis was performed in SPSS V.17.0 software for Windows (SPSS Inc., Chicago, Illinois, USA) at a significance level of 0.05. For  wing shape dimorphism, DFA was performed, and shape differences between males and females of each species were estimated based on Mahalanobis distances in permutation test with 10,000 rounds using MorphoJ software [39]. In addition, cross-validation test was performed to assess the accuracy of the classification.

Phenetic relationships of wing shape among blow fly species
To examine phenetic relationships among 12 blow fly species based on wing morphology, UPGMA (unweighted pairgroup method with arithmetic averages) was performed using PAST V.3.09 software [41]. The UPGMA dendrogram was constructed based on Mahalanobis distances obtained by pairwise comparison of analysed species from CVA.

Allometric effects
Allometry tries to describe how the characteristics of creatures change with size. For example, wing size sometimes affects wing shape variation (allometry) [34,42]. To estimate such an allometric effect, the regression of Procrustes distance (dependent variable) on centroid size (independent variable) was analysed among species and within each species separately. Moreover, the sex-dependent effect of size on shape was analysed by a multivariate regression of shape, pooled within sex, on centroid size using MorphoJ software [39]. We also evaluated the effect of removing allometry on species and sex discrimination. The residuals from the regression of Procrustes coordinates on centroid size from the previous analyses were used for assessing the differences in shape without the size effect (allometry-free variables). The residuals from regression were subjected to a cross-validation test in DFA based on Mahalanobis distances in a permutation test with 10,000 rounds using MorphoJ software [39].
cuprina were significantly different from the other ten species (Mann-Whitney U-test, P < 0.0008). In this regard, the difference among species could not be explained by their wing size difference. For the effect of sexes on size, most species (7/11) showed no significant difference between males and females (Mann-Whitney U-test, P > 0.05), except for Ch. rufifacies, L. cuprina, L. sinensis and He. ligurriens (Mann-Whitney U-test, P < 0.05), which males were smaller than females (Fig. 4).

Shape variation
The canonical variate analysis (CVA) was used to maximise variation between groups and minimise intraspecific variation. The shape difference in each genus or species on the shape space were scattered on the first two canonical variate axes (CV1 and CV2). The results from CVA were clearly discriminated in both genus and species. At the genus level, the first two canonical variates (Fig. 5) accounted together for 100% of the total variation (CV1 = 93.46%, CV2 = 6.54%), and showed that specimens clustered into distinct groups belonging to the same genus, and successfully placed all three genera in their respective subfamily. The scatter plot from CV1 and CV2 (Fig. 5) shows that each genus can be clearly separated from each other. Mahalanobis distances obtained by pairwise comparisons of all three genera revealed highly significant differences (permutation 10,000 rounds in MorphoJ: P < 0.0001), ranging from 4.9309 (Lucilia and Hemipyrellia) to 11.3103 (Chrysomya and Hemipyrellia) ( Table 2). Procrustes distances also showed highly significant differences between genera (permutation 10,000 rounds in MorphoJ: P < 0.0001), ranging from 0.0336 (Lucilia and Hemipyrellia) to 0.0741 (Chrysomya and Hemipyrellia) ( Table 2). Visualised shape changes along CV1 axis were found with landmarks 3, 10, 19, 2, 4, 5, 6 and 9, whereas shape changes along CV2 axis were most clear using landmarks 10, 4, 3, 7, 5, 11 and 12 (Fig. 5). The crossvalidation test showed a high percentage of correctly classified specimens with 94.3% (Hemipyrellia), 97.7% (Lucilia), and 100.0% (Chrysomya) ( Table 3).

Sexual shape dimorphism
The DFA in wing shape between males and females of most species revealed highly significant differences (permutation 10,000 rounds in MorphoJ: P < 0.0001, P < 0.01, and P < 0.05), except for Ch. pinguis and Ch. nigripes (permutation test with 10,000 rounds in MorphoJ: P > 0.05). Moreover, the percentage of correctly classified specimens after a cross-validation test ranged from 50% (Ch. pinguis) to 100.0% (Ch. chani) ( Table 5).

Allometric effects
Regression of the Procrustes coordinates on centroid size among species showed a highly significant difference (permutation test with 10,000 rounds in MorphoJ: P < 0.0001), allometry explained 2.3% of total shape variation. The relationship between shape and size within each species showed that wing shape variation was significantly correlated to size in most species (permutation test with 10,000 rounds in MorphoJ: P < 0.05), except for Ch. nigripes, Ch. villeneuvi, L. papuensis, L. porphyrina and He. pulchra (permutation test with 10,000 rounds in MorphoJ: P > 0.05) ( Table 6). Although the regression of shape variation on size was significant, the percentage of variation in wing shape explained by size changes was relatively low. Additionally, allometry causes significant differences between sexes in Ch. chani, Ch. megacephala, and Ch. villeneuvi (Table 6).
After removing the effect of size on shape variation, a cross-validation test showed a high percentage of correctly classified specimens at the generic level with 94.3% (Hemipyrellia), 97.7% (Lucilia), and 100.0% (Chrysomya). Based on the species level, the cross-validation test also showed a high percentage of correctly classified specimens in most species (>75%), except for L. sinensis (62.5%) and He. pulchra (0%) ( Table 7). Additionally, the results of shape variation between sexes after removing the allometric effects, showed wing shape between males and females were not significantly different in most species (permutation 10,000 rounds in MorphoJ: P > 0.05), apart from Ch. chani, Ch. megacephala, Ch. rufifacies and Ch. villeneuvi (permutation 10,000 rounds in MorphoJ: P < 0.0001 and P < 0.01). The percentage of correctly classified specimens between males and females of each species after a cross-validation test ranged from 35.3% (Ch. nigripes) to 100.0% (Ch. chani) ( Table 7).

Discussion
Wing size is known to be easily affected by environmental factors [27,34] and our results clearly show that size cannot be used to separate blow fly species. Only Ch. nigripes and L. cuprina were clearly separated from the other ten species included in the present study by wing size alone. Moreover, the majority of species do not show significant differences between males and females, except for Ch. rufifacies, L. cuprina, L. sinensis and He. ligurriens, where males were smaller than females. In contrast, wing shape showed to be a stable character compared to size [24,25,30] and very informative on the phylogenetic and evolutionary relationship of organisms [43][44][45]. Therefore, it was not surprising to see that our CVA results proved that wing shape could be used to separate medically and forensically relevant blow flies of Thailand, not only at the genus level but also at the species level. However, the latter depends on the genus to which the species belongs. Wing shapes of Chrysomya spp. were clearly distinct from Lucilia and Hemipyrellia species. But while the percentage of correctly classified specimens from the cross-validation test was very high within Chrysomya (>90.6%), wing shape largely overlapped within Lucilia and Hemipyrellia spp., leading to a much lower percentage of correct assignment (33.3-87.5%). It is not surprising that wing shape of most of Lucilia species overlaps. The same pattern is seen using morphology; distinguishing among Lucilia species is very difficult because most of them look alike [46], and about molecular data, several studies have shown that Lucilia species have low interspecific variation among closely related species [20,47,48]. In this regard, using a landmark-based characterization of wing morphology is a reliable technique for classifying Chrysomya spp., but a much less precise technique to separate Lucilia spp. and Hemipyrellia spp. The large

Males Females
He. pulchra

(1/3) --
Statistically significant differences between males and females based on Mahalanobis distances are denoted with asterisks (permutation 10,000 rounds in MorphoJ: ***P < 0.0001; **P < 0.01; *P < 0.05). Hemipyrellia pulchra has only females, thus it could not be used for classifying between sexes overlapping among species within Luciliinae shows that wing shape between Lucilia spp. and Hemipyrellia spp. is very similar to each other. Therefore, using wing morphometric analysis for species identification within Luciliinae should be done carefully and should be performed in combination with additional morphological methods for accurate species identification. Nevertheless, this study demonstrates that a landmark-based analysis of wing morphometry can be a good tool for identification of Thailand blow fly species, as it was already shown in previous studies on South American taxa Ch. megacephala and Ch. albiceps [33], and Cochliomyia hominivorax and Cochliomyia macellaria [24]. Our cluster analysis using UPGMA dendrogram based on the wing morphology of the 12 blow fly species clearly placed all species into their respective subfamily, either Chrysomyinae (Chrysomya spp.) or Luciliinae (Lucilia spp., Hemipyrellia spp.). The phenotypic relationships between species of Chrysomya detected here are in accordance with their molecular phylogenetic tree [49,50]. In the molecular analysis, Ch. rufifacies and Ch. villeneuvi belong to a different clade than other taxa of Table 6 Percentage of predicted indicates the amount of size-related shape variation of wings in each blow fly species and between sexes of each species  As for the Luciliinae, the phenotypic relationships between Lucilia spp. and Hemipyrellia spp. detected in this study are congruent with molecular studies. Hemipyrellia spp. formed a clade within the Lucilia spp. that provided strong support for the synonymy of Hemipyrellia and Lucilia [51]. Such a result suggests that wing morphology could detect some phylogenetic signal in Lucilia and Hemipyrellia. Thus, a landmark-based morphometric analysis of wings could be used as a valuable tool in taxonomy and systematics. In comparison with molecular techniques, a landmark-based analysis of wing morphology is simple, reliable and inexpensive, and just requires non-damaged wings for analysis.
Our allometric analysis suggests that wing size explained part of the variation in wing shape. However, the percentage of total variation in wing shape explained by changes in size was very low (2.3%). Thus, allometric effects seem not to be the main factors for shape variation among species. We also found intraspecific allometric effects in most analysed taxa, indicating that size-related shape changes varied among individuals within the same species. In addition, there were also significant allometric differences between sexes of Ch. chani, Ch. megacephala and Ch. villeneuvi. Due to its significant impact, it was important to perform the shape analysis with allomery-free variables. When removing the effect of size on shape variation, the percentage of correctly classified specimens among genera remained the same. The percentage of correctly classified specimens among species, however, decreased slightly. These results show that the removal of allometric effects does not improve species separation and that allometry is not an important factor in wing shape variation among species. Moreover, sexual shape dimorphism was often found in most species when allometric effects were included, but less relevant after removing the allometric effects. In general, male wings were narrower when compared to female wings. Similar results have been reported in blow flies, Co. hominivorax and Co. macellaria, which sexual shape dimorphism showed that male wings were narrower than female wings in both species [24]. The study of Hall et al. [28] showed significant sexual shape dimorphism in Chrysomya bezziana. This suggests that allometry is an important factor of sexual shape dimorphism in wings, which is commonly found in other insects [42]. Therefore, the estimation of the allometric effects is a necessary step in any study of phenotypic variation.
Due to a small number of specimens for L. sinensis and He. pulchra in this study, further studies including more specimens of these two species are recommended to increase the reliability of wing shape for species discrimination. Although wing landmark-based analysis can be a time-consuming process (e.g. in locating the landmarks for a large-scale study) this technique is simple and high reliable. The reliability of wing morphometric analysis depends on (i) wing preparation, i.e. wings should always be processed in the same way, with flattened slide-mounted wings providing the most accurate method of wing measurement [28]; (ii) morphometric analysis, e.g. using the same photographic equipment with the same conditions, operated by the same person to produce the data, the user should have some skills in collecting landmark coordinates, and digitization should be repeated at least once to reduce the measurement error by averaging the repeated digitizations [35,38].

Conclusions
Our results demonstrate that a landmark-based analysis of wings could be used to separate medically and forensically relevant blow flies of Thailand at both genus and species levels, even though it is performed with and without the effects of allometry. Using wing landmarks was a highly reliable method for classifying Chrysomya species, but less reliable for species discrimination of Lucilia and Hemipyrellia. Allometry did not affect species separation but had an impact on sexual shape dimorphism. Therefore, an estimation of possible allometric effects is a necessary step in any study of phenotypic variation by morphometrics methods. The use of wing morphometric analysis could be an alternative method used for both species and sex discrimination. In addition, the congruence between wing morphometric analysis in the present study and molecular phylogenetic tree from the previous studies, suggest that wing morphology is a valuable tool in taxonomy and systematics.