Towards the automated identification of Chrysomya blow flies from wing images

The Old World screwworm fly (OWSF), Chrysomya bezziana (Diptera: Calliphoridae), is an important agent of traumatic myiasis and, as such, a major human and animal health problem. In the implementation of OWSF control operations, it is important to determine the geographical origins of such disease‐causing species in order to establish whether they derive from endemic or invading populations. Gross morphological and molecular studies have demonstrated the existence of two distinct lineages of this species, one African and the other Asian. Wing morphometry is known to be of substantial assistance in identifying the geographical origin of individuals because it provides diagnostic markers that complement molecular diagnostics. However, placement of the landmarks used in traditional geometric morphometric analysis can be time‐consuming and subject to error caused by operator subjectivity. Here we report results of an image‐based approach to geometric morphometric analysis for delivering wing‐based identifications. Our results indicate that this approach can produce identifications that are practically indistinguishable from more traditional landmark‐based results. In addition, we demonstrate that the direct analysis of digital wing images can be used to discriminate between three Chrysomya species of veterinary and forensic importance and between C. bezziana genders.


Introduction
The placement and form of wing support veins, as well as the configuration of landmark point locations defined by wing vein intersections or vertices, have often been found to be informative indicators of family, genus and species-level distinctions across many insect orders (Comstock & Needham, 1898;Johnson & Triplehorn, 2005). The configuration of support vein vertices has also proven especially well suited to geometric morphometric studies (Morgan, 1912;Rohlf, 1993;Hall et al., 2014;Quezada-Euán et al., 2015). However, these support vein vertices rarely exhibit an even distribution across the entire wing surface and the progressive reduction in the sizes of post-branching secondary, tertiary and quaternary veins often Correspondence: Norman MacLeod, Department of Earth Sciences, Natural History Museum, Cromwell Road, London SW7 5BD, U.K. Tel.: + 44 207 942 5204; E-mail: n.macleod@nhm.ac.uk makes the forms of these branches -and hence the configuration of their vertices -indistinct under microscopic inspection, as well as in photomicrographs (Fig. 1A). More importantly, though, an exclusive focus on the configuration of wing-vein vertices discards all non-vertex information about the forms of the wing margin and the geometry of the wing veins themselves. As a result, landmark configurations derived exclusively from these structures usually provide an incomplete representation of actual wing morphology. The collection of large numbers of landmark points from insect wing images can also be a labour-intensive, time-consuming and error-prone process (MacLeod, 1990). All these issues limit the ability of field workers, who are often located far from laboratories containing even basic microscopic imaging and morphometric data collection systems, to make fine distinctions between insect populations and species correctly and consistently on site. One way of addressing this problem is to move from the current reliance on individuals or taxonomic experts to make routine identifications to a system in which taxonomic identifications are automated, based on the analysis of easily accessed and relatively large morphological features.
Chrysomya bezziana represents an interesting and important case in point. This species is a significant agent of traumatic myiasis, a worldwide human and animal health problem and the cause of severe economic losses, especially in developing countries (Hall et al., 2016). Although molecular studies have demonstrated the existence of considerable gene diversity between African and Asian populations (Hall et al., 2001;Wardhana et al., 2012), monitoring efforts would be facilitated greatly if morphological markers capable of enabling these same identifications could be found . Chrysomya megacephala and Chrysomya rufifacies are other, facultative agents of traumatic myiasis in Asia (Sukontason et al., 2005) and have potential forensic importance as indicators of the minimum post-mortem interval. Anecdotal field and laboratory experience has indicated that it can sometimes be difficult to discriminate between adults of these species using visual inspection. In particular, females of C. bezziana and C. megacephala are usually separated on the basis of subtle differences in frons morphology, the unambiguous determination of which requires experience (Irish et al., 2014). As is typically the case with taxonomic data, no controlled, empirical assessment of the correctness of even expert identifications of these species has been published to date.
The quantitative analysis of insect wing form has largely focused either on the collection and analysis of gross wing dimensions (e.g. length, width) or on configurations of landmark points placed at the intersections of major support veins (Rohlf, 1993;Hall et al., 2014). However, geometric morphometrics-based procedures have recently become available for use directly on digital images of specimens (MacLeod, 2015a). Such image-based data can circumvent the need, at least in the first instance, for laborious collection, processing and analysis of landmark and/or semi-landmark data procedures while still providing access to all the data analysis tools, summaries and model-based interpretive aids that give geometric morphometrics its power to test morphology-based biological hypotheses (Adams et al., 2013;MacLeod, 2017).
Under this digital image-based approach, wing morphologies may be specified by the brightness or colour values of the individual pixels themselves. Whereas traditional geometric morphometrics contains no information about the brightness or colour of the structures represented, the image-based analogue to geometric morphometrics can take account of such varying components of organismal form, or not, as the data analysis situation demands. Most importantly, because the pixel-based representation of forms includes all information present in the digital image, there is no need to define the location of particular features of interest a priori, to limit the analysis only to the tracking of those aspects of the form that can be defined by small sets of landmarks or by outline semi-landmarks, or to limit the sample to include only those specimens that exhibit all the features used to define landmark and/or semi-landmark locations. Relaxation of these strict comparability requirements  , proximal anal vein; A-CT, anal-cubitus transverse vein; Al-An, alulal-anal indentation. Landmarks: 1, intersection between the humeral transverse vein and costal vein; 2, position of the humeral articulation or break; 3, intersection between the subcostal and costal veins; 4, intersection between the anterior branch of the anterior radius vein (AR) and the costal margin (C); 5, intersection between the first post-anterior branch of the radius vein (PA 1 R) and the costal margin (C); 6, intersection between the second post-anterior branch of the radius vein (PA 2 R) and the costal margin (C); 7, intersection between the medial vein (M) and the costal margin (C) or wing periphery; 8, intersection between the cubital vein (Cu) and wing periphery; 9, distal intersection between the anal and alulal areas (Al-Au); 10, articulation between the humeral sclerite and the radius vein (H-RA); 11, first bifurcation of the radius vein resulting in the creation of this vein's anterior (AR) and first post-anterior (PA 1 1R) branches; 12, second bifurcation of the radius vein resulting in the creation of this vein's first (PA 1 R) and second (PA 2 R) post-anterior branches; 13, anterior intersection between the radius-medial transverse vein (R-MT) and the second post-anterior radius vein (PA 2 R); 14, anterior intersection between the distal medial-cubital transverse vein (M-CT 1 ) and the medial vein (M); 15, posterior intersection between the radius-medial transverse (R-MT) vein and the medial vein (M); 16, anterior intersection between the distal medial-cubital transverse vein (M-CT 2 ) and medial vein (M); 17, maxima of curvature in the distal portion of the medial vein (M); 18, bifurcation between the cubital (Cu) and first anal (A 1 ) veins; 19, posterior intersection between the proximal medial-cubital transverse vein (M-CT 1 ) and cubital vein (Cu); 20, posterior intersection between the proximal medial-cubital transverse vein (M-CT 2 ) and cubital vein (CU); 21, posterior intersection between the anal-cubitus transverse vein (AC-T) and first anal vein (A 1 ). [Colour figure can be viewed at wileyonlinelibrary.com]. renders the image pixel-based approach a flexible -and so a useful -tool for the exploration, discovery, documentation and interpretation of morphological variation.
Although they are common in a variety of manufacturing and security contexts, applications of this new image-based approach to morphometric analysis in biological systematics have been few to date. In particular, no head-to-head comparison of the results that can be achieved by this new approach with the results generated by a traditional landmark-based geometric morphometric analysis in insects have yet been published. So that such a comparison can be made in the context of population-level identifications of a well-known and economically deleterious disease-causing species, the present study evaluated the utility of direct analyses of whole-wing morphology to support gender and population-level identification of C. bezziana and the taxonomic species-level identification of C. bezziana, C. megacephala and C. rufifacies based on a set of small samples in the manner of a pilot, or proof-of-concept, study.

Materials and methods
African C. bezziana were selected randomly from collections held at the Natal Museum (Pietermaritzburg, South Africa) and the National Collection of Insects (Pretoria, South Africa) and at the Natural History Museum (London, U.K.) . Male and female C. bezziana from Sumba, Indonesia were taken from a colony maintained by the Parasitology Department of the Indonesian Research Centre for Veterinary Science, Bogor, Indonesia, which was established from larvae collected from a case of traumatic myiasis on that island . Sumatran C. megacephala and C. rufifacies were collected as adults on Bezzilure-baited sticky traps (Sulston et al., 2014) placed at locations in Indrapuri District. Both Sumban and Sumatran specimens were preserved in 80% ethanol and maintained at − 20 ∘ C until analysis.

Methodology
Digital images (Fig. 1A) of isolated and slide-mounted wings were obtained using standard entomological photomicrography procedures . The majority of wings imaged were right wings. Images of left wings were transformed into pseudo-right wing images by mirroring and tested statistically for systematic differences with respect to right wing morphologies. No effort was made to clean the images or to select only well-preserved, complete and intact (non-torn and non-damaged) specimens.
For traditional geometric morphometric analysis, 21 vertexor termination-based landmarks were selected (Fig. 1B). As wing vein vertices and terminations are biological homologues across the species considered by this investigation, all are of the 'type 1' variety (Bookstein, 1991). Type 1 landmarks are widely considered to be the most reliable sources of evidence for documenting morphological similarities and differences in an evolutionary context. In our image-based analysis, the same digital images used to define the landmark configurations were trimmed to remove the background pixels while leaving a 10% pixel margin around each wing, scaled to ensure each wing was represented by the same number of pixels to a tolerance of 5%, down-sampled to ensure each image frame conformed to a 40 × 100-pixel matrix, and converted from an 8-bit RBG to an 8-bit greyscale colour space in order to minimize variation that might result from differences in specimen colour (Fig. 2). Once in this format, all images were aligned by their major axes and centred within the image frame. After alignment, each image was reformatted into a column vector of pixel brightness values and combined to form a data matrix.
Covariance matrices calculated from both landmark and image data were submitted to a preliminary principal component analysis (PCA) to repack the recorded morphological variation into the smallest number of statistically independent variables consistent with preservation of 95% of the original pooled sample shape variation, and projected as scores onto the PCA axes (= eigenvectors) to achieve an ordinated representation of shape similarity relations. These score data were then submitted to a canonical variates analysis (CVA) to create an ordination space that maximized between-group separations relative to within-group dispersions (MacLeod, 2015a(MacLeod, , 2015b. Combining a pre-processing PCA step with a discriminant function or machine-learning analysis is known to often result in the improvement of group separations (Anderson & Willis, 2003;MacLeod, 2015a;Marramá & Kriwet, 2017). A 1000-iteration bootstrap variant of the log-likelihood ratio test (Manly, 1997) was used to estimate the statistical significance of the resulting group separations. Image models for both CVA axes were calculated using the method described by MacLeod (2009). All image processing and data analysis operations were performed with Wolfram Mathematica™ software, written for this project by the senior author. Mathematica™ notebooks used in this investigation are available upon request.

Geographical difference test
The dataset for this test consisted of 26 African C. bezziana female wings (20 right wings, 6 left wings) and 20 Asian C. bezziana female wings (13 right wings, 7 left wings). The initial PCA of the pooled landmark and image covariance matrices showed that 18 and 28 components, respectively, were able to represent 95% of the shape variation recovered by these alternative morphology sampling schemes. Relative to the total number of variables in the original data matrices, this operation reduced the effective dimensionality of these data by 60.9% and 39.1%, respectively.

Landmark results
The first three principal components of the landmark dataset accounted for 26.6%, 12.0% and 10.9% of the sampled shape variation, respectively. Each of the other components accounted for < 10% of the total. A plot of this sample's landmark configurations projected into the space of the first two components shows that pronounced distinctions in wing shape between African and Asian specimens are present in these data (Fig. 3A).
Whereas a clear separation between areas of geographic origin is evident along PC-1, this plot indicates that an aspect of wing shape variation captured by the PC-2 axis (and most likely other PC axes) is also responsible for these group-level shape distinctions. In order to optimize this geographic discrimination and test for the presence of consistent and statistically significant group separation, the scores obtained from the projections of these landmark configurations on all 18 PC axes were submitted to a CVA (Fig. 3B).
As suggested by the PCA result, the configurations of these 21 wing landmarks contain sufficient information to allow African and Asian C. bezziana individuals to be distinguished cleanly and consistently from one another. A test of the significance of the log-likelihood ratio value calculated from these data ( = 106.1) was found to exceed all -values calculated for 1000 non-parametric, random pseudoreplicate datasets drawn from the original data (with replacement), the maximum value of which was 44.0. Based on this non-parametric bootstrap result, these results provide evidence for a statistically significant difference between African and Asian wings at a value approaching = 0.0%. With respect to the stability of this discriminant result, a jackknife (leave-one-out) analysis of the proportion of specimens identified correctly post hoc, as if they were unknown specimens, resulted in the misidentification of only a single African specimen, thereby providing a robustly estimated correctness ratio of 97.9%. Modelling of landmark shape configurations along this discriminant axis revealed that the most striking single distinction between African and Asian C. bezziana wings was a very strong inboard migration of the marginal boundary between the wing's alulal and anal regions in the vicinity of the jugal fold, which may be characteristic of Asian populations (Appendix S1 and Fig. S1, online, give further details).

Image results
A PCA of the size-standardized image covariance matrix for the pooled African-Asian C. bezziana sample showed that 28 eigenvectors were needed to represent 95% of the shape variation recorded. The first three of these represent almost 60% of the total shape variation with none of the subsequent components accounting for > 5% of the total. Figure 4A illustrates the distribution of image pixel configurations projected into the subspace formed by the first two principal components of these data.
Despite the apparent overlap between African and Asian wing image configuration fields within the PC-1 vs. PC-2 ordination space, when all 28 PC variables were included in a CVA, African and Asian wing images were discriminated into mutually exclusive fields (Fig. 4B). A non-parametric bootstrapped test of group separation along this axis based on the log-likelihood ratio found that none of 1000 pseudoreplicate datasets created randomly from these data (with replacement) produced a -index value greater than that calculated for this empirical discriminant result ( = 85.7; maximum simulated value: 66.3). Accordingly, this result is also interpreted as statistically significant at an -level of close to 0.0%. Although the raw confusion matrix for this result produced 100% correct results for the training set, a more rigorous jackknifed identification test detected a degree of instability in the single discriminant function calculated from the image pixel data, with seven of the 47 specimens (14.9%) being misidentified post hoc, as a result of being treated as unknowns. These misidentifications were sub-equally divided among the two groups, with three African specimens being identified as Asian and four Asian specimens being identified as African. Figure 5 illustrates two hypothetical end-member images calculated from the projected positions of extreme African and Asian coordinate locations along the discriminant axis and includes a pixel change difference map that summarizes the comparison of these two images with temperature-scale colour codes. Note that, for this comparison, large-to-maximal pixel changes are confined to placements of the major wing veins comprising the pteralia and moderate changes are recorded for pixels denoting wing blade veins and the alular area margin.

Species difference test
This test was conducted in 20 Asian female C. bezziana wings (13 right wings, 7 left wings), 20 Asian female C. rufifacies  wings (16 right wings, 4 left wings) and 20 Asian female C. megacephala wings (20 right wings). After tests for consistent differences between right and left wings of C. rufifacies yielded negative results, the four C. rufifacies left wings were transformed into pseudo-right wing images. Subsequent to this operation all images were processed according to the procedure described above (see Methodology).
The first two eigenvectors of the pooled image covariance matrix accounted for 42.4% of the recorded shape variation with only two subsequent axes representing > 5% of the total. For this sample, 32 eigenvectors were necessary to represent 95% of the recorded image variation, demonstrating a reduction in effective dimensionality of 46.7% in terms of the maximum number of components with positive eigenvalues and > 99% in terms of the dimensionality of the original images.
An inspection of the ordination of images in the PC-1 vs. PC-2 subspace (Fig. 6A) clearly shows that the primary source of wing shape variation for these data lies in the distinction between C. bezziana with respect to C. megacephala and C. rufifacies. It is also clear that C. bezziana exhibits the greatest variance in wing pixel configuration and C. rufifacies the least. Owing to the low dimensionality of this representation of wing shape variation, it is not clear whether these species can be distinguished from each other reliably by wing image data alone. However, when the PCA scores on all 32 eigenvector axes necessary to represent 95% of the image covariance structure were submitted to a CVA, the between-species differentiation was both clear and compelling (Fig. 6B).
A non-parametric bootstrap test of group separation along these discriminant axes found that none of the 1000 pseudoreplicate log-likelihood ratio values for randomly created datasets (with replacement) produced -values of > 120.0, whereas the -index value calculated for empirical discriminant results was 261.1. Accordingly, this result is interpreted as being very significant statistically ( ≈ 0.0%). Irrespective of the small sample sizes involved in this test, this discriminant result exhibited remarkably good stability with a jackknife performance test resulting in only five (9%) misidentifications (four C. bezziana, one C. rufifacies). Numbers below the hypothetical image models represent the discriminant axis coordinate values at which the models were calculated. In the pixel difference map values of the differences in pixel brightness values are colour-coded using a thermometer-referenced colour scale: blue, no or small change; white, moderate change; yellow and orange, strong change; red, maximum change. CV, canonical variate. Because only three species are considered, these axes represent 100% of the difference between group centroids in the linear, eigenvalue-standardized ordination space. Figure 7 illustrates hypothetical end-member images calculated from the extreme projected positions of Chrysomya species for the image CV-1 and CV-2 discriminant axes, along with pixel brightness difference maps summarizing the comparison of these images using temperature scale colour codes. Note that large-to-maximal pixel changes along CV-1 are confined largely to the variation in veins located along the entirety of the wing's leading edge, whereas those represented by CV-2 predominantly involve separate proximal and distal leading edge vein segments, in addition to subsidiary variation among the pteralia veins, the internal mid-wing support veins, the alular margin and the wing's trailing margin.

Gender difference test
This dataset consisted of 20 images of Asian female C. bezziana wings (13 right wings, 7 left wings) and 24 images of Asian male C. bezziana wings (21 right wings, 3 left wings). After a preliminary investigation failed to find any consistent shape differences between right and left wing images for either dataset, all left wing images were mirrored into right wing images and processed according to the procedure described in the Methodology section.
The first three eigenvectors of the pooled gender difference image covariance matrix accounted for 48.7% of the observed wing shape variation with no subsequent axis accounting for > 8% of the total. Twenty-seven PC axes were necessary to represent 95% of this dataset's image-based shape variation structure. This represents a reduction in effective dimensionality of 38.6% in terms of the maximum number of components with positive eigenvalues and > 99% in terms of the original images' dimensionality.
The PC-1 vs. PC-2 ordination plot for these data (Fig. 8A) demonstrates that, although there are broad sex-based differences in wing shape, these differences do not, by themselves, account for the primary structure of wing shape variation in the pooled sample along either PC-1 or PC-2. Nonetheless, as in the previous species difference dataset, despite the apparent overlap between male and female Asian wing image configurations within the PC ordination space, when the projected scores along all 27 eigenvector variables were submitted to a CVA, male and female wings were separated into remarkably distinct and mutually exclusive fields (Fig. 8B). A non-parametric bootstrapped log-likelihood ratio test of this group separation found that none of the 1000 pseudoreplicate datasets created randomly (with replacement) from these data produced -values of > 56.5, whereas the empirical -value calculated for the gender difference discriminant result was 94.3. Accordingly, this result is interpreted as statistically significant at a level close to = 0.0%.
Once again, despite the small size of this sample, this discriminant function exhibited very good stability, with a leave-one-out jackknife analysis of post hoc identification power returning correct identifications for 93.2% of training set wings treated as unknown specimens. One male specimen was misidentified, post hoc, as female and two female specimens were misidentified as male. Figure 9 illustrates the reconstructed forms of typical male and female Asian C. bezziana wing morphologies at extreme locations along the discriminant axis, along with a colour-coded image difference map that highlights regions of pronounced sexually dimorphic difference. This result indicates that the distinction between Asian male and female C. bezziana wing morphologies primarily involves the placement of veins across the leading wing margin, that areas of maximal change occur in both proximal and distal segments and, subordinately, along the trailing edge of the main wing and alular region and that females exhibit markedly broader trailing wing margins.

Discussion
In view of the small sample sizes used in this pilot investigation, we do not advocate the use of any of the discriminant functions obtained from any of the datasets analysed herein until a larger sample of Chrysomya species, populations and genders has been obtained. Nevertheless, the comparison of these results reveals a number of interesting and significant findings.

Discrimination between African and Asian C. bezziana
Direct comparison of the PCA results from image analyses of African and Asian wings (Fig. 4A) with those obtained from the landmark-based representation of the same images (Fig. 3A) is instructive. First and most notably, although the major axes of wing image variation for both datasets were aligned with the distinction between African and Asian groups, the separation between these two groups for the image dataset was not mutually exclusive as it was in the landmark dataset. Rather, one -possibly two -distinct subgroups of African C. bezziana wings occupy the wing image space along the lower reaches of PC-1 and PC-2. Moreover, the structure of wing morphological diversity is far more evident in the image-based summary of shape variation than in the landmark-based summary. A very prominent discontinuity in wing shape variation is evident in the image PC-1 vs. PC-2 subspace, as is the existence of several subsidiary African sub-clusters as well as shape outliers.
Close comparison of the image models and difference map for the image data (Fig. 5) with the PC-1 axis models included in Appendix S1 reveals an unexpectedly high degree of consistency among the ordination trends along the PC-1 axes of both analyses. Given the gross similarities between these two sets of results (Figs 3 and 4), this similarity is not unexpected, but its presence provides a strong indication that both analyses found similar trends and distinctions in their respective representations of wing shape variability and group difference, albeit distinctions that are revealed more completely and in greater detail in the case of the higher-resolution wing image data.  Numbers below the hypothetical image models represent the coordinate values at which the models were constructed. The pixel difference map is colour-coded using a thermometer-referenced colour scale: blue, no or small change; white, moderate change; yellow and orange, strong change; red, maximum change. CV, canonical variate.

Discrimination among Chrysomya species
Interestingly, the distinction between the wing images of C. bezziana and those of the other two species is not the most prominent differentiation pattern recorded by the CVA ordination (Fig. 6B). Rather, C. rufifacies was identified as exhibiting the most distinctive wing morphology in this sample, although only by a marginal amount. A more conservative -and perhaps more reasonable -interpretation of this discriminant result is that wing morphologies in all three species exhibit sub-equal patterns of uniqueness, at least as represented by these samples.
Close graphical inspection of the present CVA-based Chrysomya species results (Fig. 7) can also be used to pinpoint the features of wing morphology that differ characteristically between these species and applied to develop reliable, objective and evidence-based qualitative guides for delivering high-quality taxonomic identifications. With regard to the difference between C. rufifacies and C. bezziana/megacephala -expressed along CV-1 in the present results -the former is characterized by wings with thinner, more elliptical blades than those of the latter two species. This narrowing of the blade is associated with the distal migration of the intersections between the costal margin and the subcostal, the anterior radius and both first and second post-anterior radius veins, proximo-lateral migration of the medial cubital and anal veins, and anterior migration of the anal blade margin. As a result of these modifications, of the three species assessed in the present study, C. rufifacies wings exhibit the most elliptically eccentric or pointed tips.
With regard to the interpretation of CV-2, which primarily captures the distinction between C. megacephala and C. bezziana, wings that project to positions low along this axis (= C. megacephala) are characterized by relatively narrow subcostal and anterior radial regions with more proximally placed intersections between their costal margins, subcostal, anterior radii, and both first and second post-anterior radius veins, with more anteriorly placed second post-anterior, medial and cubital veins, and with a moderate posterior migration of the posterior wing margin in the mid-region of the blade. Conversely, wings that project to positions high along this axis (= C. bezziana) are characterized by slightly expanded subcostal and anterior radial regions with more distally placed intersections between their costal margins, subcostal, anterior radii, and both first and second post-anterior radius veins, more posteriorly placed second post-anterior, medial and cubital veins, and a pronounced posterior migration of the posterior wing margin in the distal region of the blade. Overall, these modifications impart a distinctly elliptical character to C. rufifacies wings, a moderately elliptical character to C. megacephala wings, and a distinctly quadrate character to C. bezziana wings (see also Appendix S1). The larger and more important point here is that the use of image data produces results at a level of detail that not only facilitates interpretation, but is unprecedented in comparison with traditional or geometric morphometric data.

Discrimination of gender in C. bezziana
A comparison of the models shown in Fig. 9, including the CV-1 image difference map, shows that the primary distinctions between male and female Asian C. bezziana wings fall into three categories. The first and most important of these is an increase in the space between the costal margin and both the lateral segments of the subcostal and primary anterior radius veins in females relative to the distinctly more compact arrangement of these features in males. This expansion is most obviously manifested as a migration of the costal margin in an anterior direction with particularly intensive distinctions between pixel greyscale values occurring at the anterior terminus of the humeral transverse vein, the humeral articulation or break, the anterior termini of the subcostal anterior radius, and the first and second post-anterior radius veins. Coincident with these distinctions, males also exhibit wings that are narrower distally than female wings. This difference is manifested as an anterior migration of the posterior wing margin, especially in the middle region of the wing blade, but also involves a more subtle anterior migration of the distal medial and medial-cubital transverse veins. In general, these changes can be summarized as an overall narrowing of the main wing blade in male Asian C. bezziana relative to the typical female condition, which accentuates the ellipticity of the distal wing tip. Finally, the indentation in the posterior wing margin that marks the intersection between the anal alular regions is, on the whole, far more distinct and deeper in male than in female wings. Once again, the level of morphological detail achieved by this image-based analysis is far greater than those delivered in traditional or geometric morphometric investigations (e.g. Hall et al., 2014).

Use of image analysis in Chrysomya taxon discrimination/identification
The present analyses of C. bezziana, C. megacephala and C. rufifacies wing images have shown that this appendage not only contains sufficiently well-structured morphological variation to distinguish among these three species with an accuracy of > 90%, but that similar analyses can be used to distinguish between African and Asian C. bezziana, and even between genders of this species, with comparably high accuracy levels. These identifications were obtained via the direct analysis of digital wing images with minimal image pre-processing and in wings that included a wide variety of both structural and photographic imperfections.
Despite the small sample sizes employed in this investigation, these results present compelling evidence that, in principle at least, unexpectedly detailed and marked improvements in the ability of non-expert investigators to assess and identify a range of important biological parameters can be obtained by identification systems based on digital images, and that these identification procedures can achieve levels of accuracy that match or exceed those offered by landmark-based morphometric analyses (e.g. Hall et al., 2014;Giordani et al., 2017). In the present study, all discriminant results were determined to be statistically significant using robust non-parametric tests. Moreover, all identification accuracy ratios obtained indicate that the performance of these preliminary linear discriminant functions is at least as good, if not better, than would be expected on the basis of experienced expert identifications according to the levels of correctness and consistency in identification tests conducted in other organismal groups [see Culverhouse et al. (2013) and references therein] as no similar studies have been conducted in any insect group.
Our discriminant function stability test results indicate that groups characterized by larger sample sizes exhibited more accurate and consistent identifications than groups represented by smaller sample sizes. This finding, in turn, implies that the current results represent, in a sense, worst-case scenarios with regard to the accuracy, consistency and effectiveness of the present approach to the analysis of insect wing variation. It is likely that these results will only improve if and when larger wing morphology training sets are acquired. The ability of our current image-based approach, not only to identify Chrysoma species correctly, but also to recover geographic, group and gender-based distinctions in wing morphology-aspects of variation that, with a single possible exception (Spradbery, 2002), have never been recognized by traditional qualitative entomological taxonomists or reported by previous morphometric investigators-contain important implications for other insect-morphological investigations in particular and biological-ecological investigations more generally. Not only can an image-based approach to morphological analyses have real value as an aid to biological investigation, but there are indications that an untapped wealth of important and useful information remains to be discovered in morphological data generally.
From a methodological point of view, our results also throws the inherent ambiguity of PCA into sharp focus in terms of its ability to deliver true representations of between-group distinctions. Over the years many entomological biometers have relied on PCA ordinations to provide estimates of group distinctiveness (e.g. Orlofske & Baird, 2014;Schwarzfeld & Sperling, 2014;Jaramillo et al., 2014), despite the fact that standard PCA operates only on the pooled covariance or correlation matrices of datasets irrespective of whether datasets possess subsidiary structure or whether the assessment of such structure is the point of the investigation. The differences obtained in the present work between the group-level distinctions in the subspace formed by the first few PCA eigenvectors, compared with those obtained from the same PCA data when subjected to a subsequent CVA, should be considered as a cautionary note to those considering group-level morphometric analyses. The excellent quality and stability characteristics of the CVA results obtained from our present image data also indicate that, even in cases in which sample sizes are much smaller than the dimensionality of the original data, procedures are available that will allow entomological data analysts to test specific hypotheses with a high degree of precision and to obtain statistically valid results.
Finally, it is worth noting that, while DNA barcoding could be employed, in principle, for the purposes of species, population and gender identification, a more efficient, practical, rapid and in many ways interesting approach to this problem is to identify a morphological character complex that has sufficiently well-structured patterns of variation to support automated or quasi-automated identification, especially if such identification systems can be located in the field, employ inexpensive and widely available technology (e.g. smartphones with cameras), and be used by non-expert parataxonomists or even agricultural workers. Based on the results obtained in the present investigation of African and Asian Chrysomya species, it seems clear that the morphology of Chrysomya wings may represent such a character complex.

The future of image-based identification
The primary goal of this investigation was to determine whether images of the Old World screwworm fly (genus Chrysomya) characteristics could be used directly as inputs to a geometric morphometric-style investigation of various group-level morphological distinctions. The wing was selected as the subject for analysis because this appendage is not only a fundamental aspect of the body plan in most insects, but it has been employed previously in numerous qualitative, semi-quantitative and biometric entomological investigations, including landmark-based geometric morphometric analyses (e.g. Rohlf, 1993;Giordani et al., 2017). The overall data analysis approach employed in the present study is identical to that used in standard geometric morphometric analyses (e.g. Hall et al., 2014).
The digital insect wing images employed in the current investigation captured all the biologically relevant morphological information available in the sample and used the totality of this information as the basis for analysis. That this more generalized and straightforward approach to the analysis of morphological data was able to return results fully comparable with those returned by more labour-intensive, landmark-based geometric morphometric analyses, and superior to published results of population and gender identifications based on qualitative inspection of wing images and/or actual wings by experienced taxonomists, suggests the present approach may be able to deliver reliable, accurate and fully automated identifications for these (and other) insect species in the very near future. Our approach to morphometric investigation is also consistent with the spirit, if not entirely with the letter, of recent criticisms made by Bookstein (2016) regarding the appropriateness of standard geometric morphometric data and procedures for characterizing and resolving questions pertaining to biological hypotheses. Furthermore, the productive cross-fertilization that can exist between image-based morphometrics and image-based computer vision hinted at by Bookstein (2016), and demonstrated by the present results, provides a clear signpost to the directions in which both fields might want to consider moving in the future.
A landmark analysis of the wings used in the geographic group difference test confirmed the existence of differences between African and Asian wing morphologies and returned a marginally higher group identification accuracy ratio. The ease, speed and convenience of the analysis of wing images directly, along with the high identification accuracy ratios achieved in this investigation, leads us to conclude that this strategy, if and when it is automated fully, may be sufficient to provide field taxonomists, forensic researchers, medical specialists and agricultural workers with a quick, easy and reliable means of identifying species, population origins and even genders in the field or in the laboratory without the need to routinely consult trained taxonomic specialists. If more accurate identifications are needed, these could be generated for a subset of specimens for which automated identifications are judged to have low or marginal confidence values, either by expert taxonomists or according to the results of morphometric data analyses in adequately equipped research laboratories using traditional landmark/semi-landmark data. In this context it is also important to note that the present investigation was designed only to provide data for which application of the same morphometric procedures to different datasets could be compared directly. The results detailed in this report, of course, ignore the additional performance improvements that, in principle, could be gained by moving from linear morphometric to non-linear or machine-learning approaches to the analysis of such data (MacLeod, 2007(MacLeod, , 2017.

Supporting Information
Additional supporting information may be found online in the Supporting Information section at the end of the article. DOI: 10.1111/mve.12302 Appendix S1. Supplementary materials contents. Figure S1. Directions of landmark migration between typical African (landmark configuration and displacement vector tails) and Asian (displacement vector heads) wings. See Fig. 1 for landmark names and defining criteria. The lengths of these vectors have been exaggerated (x3) in order to illustrate their directions and relative magnitudes. In addition, the landmark icons have been colour-coded to classify them into relative change categories: white, no or small change; light grey, moderate change; dark grey, strong change; black, strongest change.
previous version of this article which resulted in substantial clarification of the arguments put forward, but are in no way responsible for any remaining commission or omission errors. Their contributions are gratefully acknowledged. The authors are also grateful to the International Atomic Energy Agency (IAEA; Vienna, Austria) for funding the insect collecting in Indonesia (grants TC14842 and RC14860). All image processing and data analysis operations were performed with Wolfram Mathematica™ software written for this project by the senior author. All Mathematica™ notebooks used in this investigation are available from the corresponding author upon request.