Geometric morphometric wing analysis represents a robust tool to identify female mosquitoes (Diptera: Culicidae) in Germany

Accurate species identification is the prerequisite to assess the relevance of mosquito specimens, but is often hindered by missing or damaged morphological features. The present study analyses the applicability of wing geometric morphometrics as a low-cost and practical alternative to identify native mosquitoes in Germany. Wing pictures were collected for 502 female mosquitoes of five genera and 19 species from 80 sampling sites. The reliable species identification based on interspecific wing geometry of 18 landmarks per specimen was tested. Leave-one-out cross validation revealed an overall accuracy of 99% for the genus and 90% for the species identification. Misidentifications were mainly due to three pairings of Aedes species: Aedes annulipes vs. Aedes cantans, Aedes cinereus vs. Aedes rossicus and Aedes communis vs. Aedes punctor. Cytochrome oxidase subunit I (COI) gene region was sequenced to validate the morphological and morphometric identification. Similar to the results of the morphometric analysis, the same problematic three Aedes-pairs clustered, but most other species could be well separated. Overall, our study underpins that morphometric wing analysis is a robust tool for reliable mosquito identification, which reach the accuracy of COI barcoding.

To the authors' knowledge, only three studies applied morphometric wing analysis to native mosquito species in Europe. Intraspecific variations of mosquito wings were analysed for European population of Aedes vexans (Meigen) 22 and Cx. pipiens s.l. 23 . Börstler et al. 21 used morphometric diagnosis to discriminate female Cx. pipiens s.s. and Cx. torrentium. The here presented study gives a comprehensive evaluation of this low-budget tool for the identification of a broad range of mosquito species. The morphometric wing characteristics for 19 of the most common mosquito species of Germany were analysed to evaluate the usability for species discrimination. In addition, we sequenced the COI gene to validate morphological and morphometric identification.

Methods
Mosquito sampling. The right wings of 502 females of 19 mosquito species were analysed (Table 1). Varying environmental conditions can influence the intraspecific shape and size of mosquito wings 24,25 . Hence, specimens were collected at 80 study sites (Fig. 1), covering different landscapes of Germany to capture a broad morphological variation per species. Mosquitoes were sampled over the seasons of 2017 and 2018. Adults were sampled with two different methods: CO 2 -baited Biogents-traps (Biogents, Regensburg, Germany) or handmade aspirators modified from Vazquez-Prokopec et al. 26 . Immature stages were sampled with a standard dipper (Bioquip, CA, USA) in breeding sites and subsequently reared to adults in the laboratory. All specimens were identified by morphology 14,27 and stored at − 18 °C in a freezer until further analysis. Specimens of each species were selected from at least three different sampling sites ( Table 1). The wing pictures including coordinates of the sampling location, sampling date and sampling method for each specimen are given in the supplementary material (https ://doi.org/10.5061/dryad .zs7h4 4j5s).
Measuring wing shape. The right wing of each specimen was removed and mounted under a cover slip (15 × 15 mm) with Euparal (Carl Roth, Karlsruhe, Germany). Pictures of each wing were taken under 20× magnification with a stereomicroscope (Leica M205 C, Leica Microsystems, Wetzlar, Germany). Fiji 28 as bioscience package of ImageJ 29 was used to digitize 18 landmarks ( Figure S1). The landmark selection was in accordance with other studies analysing mosquito wing morphometry [18][19][20]30 . The wing pictures were divided among two observers (authors LE and FGS) and digitalised in random order to minimize a memory biased landmark collection between the mosquito specimens of the same species. One month later, the measurement was repeated for three specimens per species by four observers (authors LE, LJ, RL and FGS) to assess the degree of observer error and repeatability in landmark collection 31 . Genetic identification. DNA isolation was performed from the whole mosquito body, except of the right wing. Individual specimens were placed into 2 ml tubes and about 10 pieces of 2.0 mm zirconia beads (BioSpec Products, Bartlesville, USA) as well as 1 ml of cell culture medium (high-glucose Dulbecco's modified Eagle's medium; Sigma-Aldrich, St. Louis, MO, USA) were added. The homogenization was performed with a Tissuelyser II (Qiagen, Hilden, Germany) for 2 min at 30 oscillations/s and 200 μl of the homogenate were used  32 using the primers by Folmer et al. 33 . Sanger sequencing was applied for all positive amplicons (LGC Genomics, Berlin, Germany). Furthermore, morphologically identified Cx. pipiens s.l. and An. maculipennis s.l. specimens were typed to species level (Cx. pipiens pipiens form pipiens resp. Anopheles messeae) using two molecular assays 34,35 . Statistical analysis. Generalised Procrustes analysis of the two-dimensional landmark data set was performed to create aligned Procrustes coordinates and centroid size of each specimen with the R package geomorph 36 . The centroid size is defined as the square root of the sum of squared distances between centroid and each landmark 37 and can be used as a proxy for wing size 38 . Mean difference in centroid size per species were compared by an analysis of variance (ANOVA), followed by pairwise comparisons of species using t-tests with Bonferroni-adjusted p-values. During the Procrustes analysis, raw landmark coordinates were centred, scaled and rotated, so that the resulting Procrustes coordinates describe the wing shape in itself 39 . Allometric effects of the centroid size on wing shape were assessed by the "procD.allometry" function from the geomorph package 36 using 500 permutations. Wing shape variables were analysed by linear discriminant analyses (LDA) to classify genus and species, respectively (R package MASS 40 ). Subsequently, each specimen was reclassified by a leaveone-out cross validation in order to test the classification obtained by the discriminant analysis. The mean shape configuration of the landmarks was used to visualize shape differences between the genera. In addition, a neighbour joining tree (NJ) was constructed to display the Mahalanobis distances between species means in a canonical variate analysis, using the functions "CVA" in Morpho 41 and "nj", "boot.phylo" and "drawSupportOnEdge" in ape 42 . Robustness of nodes was assessed by 1000 bootstrap replicates. To assess the degree of observer error, the shape variance due to the repeated measurements per specimen was compared with the mean shape variance per species of the original data set without repeated landmark collections. Therefore, shape variances were calculated by the "morphol.disparity" function in geomorph 36 using 500 permutations. Furthermore, we assessed the fidelity of morphological characterization of our landmark data set through a Landmark Sampling Evaluation Curve (LaSEC) by using the function of Watanabe 43 . All analyses and visualisation including the R package ggplot2 44 were carried out in the R environment 45 . Sequences were pre-processed with Geneious 9.1.8 (Biomatters, Auckland, New Zealand). COI sequences were trimmed to 550 bp. To assess the phylogenetic relationship of the mosquito species identified in this study, a maximum likelihood tree was constructed using functions of the R packages seqinr 46 , ape 42 and phangorn 47,48 .
The HKY + G model was identified as the best-fit model of nucleotide substitution by the phangorn's "modelTest" function based on calculations of Akaikeʼs information criterion. Robustness of nodes was assessed by 1000 bootstrap replicates. The COI sequences generated in this study were deposited in the GenBank database under the accession numbers MT731082-MT731276.

Results
The mean centroid size differed significantly among the mosquito species (ANOVA, F = 143.3, p < 0.001). Highest centroid size was observed for Culiseta morsitans (Theobald) with almost twice the centroid size compared to the smallest species in our study: Aedes caspius (Pallas), Aedes sticticus (Meigen) and Culex territans (Walker) (Fig. 2). A small, but significant allometric effect was observed (R 2 = 3.7%, p = 0.002). However, allometric residues were not removed since size variations were considered to be important for the species identification process 21 . A pairwise cross-validated reclassification test (leave-one-out method) to compare the five genera (Aedes, Anopheles, Coquillettidia, Culiseta and Culex) revealed an accuracy greater 99% (Fig. 4). Hereby, the genus of three out of 502 specimens was misidentified. One Aedes specimen was falsely identified as Coquillettidia and two Culex specimens were misidentified as Aedes. Differences in the shape variation were strongest for Anopheles in contrast to other genera (Figs. 3 and 4). In addition, the superimposed landmark coordinates of the junctions of radius 2 and 3 as well as media 1 and 2 varied most obviously when comparing the further genera (Fig. 3). Species identification accuracy based on the pairwise reclassification test was 90%. The error rate of 10% was mainly due to the comparisons of three pairings: Aedes annulipes (Meigen) vs. Aedes cantans (Meigen), Aedes cinereus Meigen vs. Aedes rossicus Dolbeskin, Gorickaja and Mitrofanova and Aedes communis (de Geer) vs. Aedes punctor (Kirby) ( Table 2). In addition, the centroid size of these pairs did not differ significantly (Table S1). Other species were reclassified with a high reliability between 94 and 100% ( Table 2, Fig. 5). The NJ tree resulted in high bootstrap values for most species, including the three pairings with a high error rate in the pairwise reclassification test (Fig. 6). While members of the genus Culex and Anopheles clustered together, species from the genus Aedes and Culiseta clustered heterogeneously over the NJ tree. The landmark collections of three randomly selected specimens per species were repeated by four different observers. The mean shape variance per specimens due to the different observers was 0.00047. This was considerably lower than the mean shape variance per species (0.00123) in the original landmark data set without repeated measurements, thus indicating a relatively small observer bias for the analysis of species-specific shape variations. The LaSEC did not suggest oversampling of landmarks in our data set ( Figure S2). High quality, clean sequences were available for 195 mosquito specimens, covering 18 of the 19 taxa only missing the species Cx. territans (Fig. 7). In order to perform a verification of the consistency of the morphological and morphometric classification of the mosquito specimens, a maximum likelihood phylogenetic tree was constructed (Fig. 7). Most species clustered seperately. However, this was not the case for several members of the genus Aedes. The two species-pairs Ae. annulipes vs. Ae. cantans or Ae. communis vs. Ae. punctor did not cluster separately. In addition, the three sequences available for Ae. rossicus and Aedes rusticus (Rossi) either grouped with sequences of Ae. cinereus (Ae. rossicus only) or Ae. sticticus (Ae. rossicus and Ae. rusticus).

Discussion
We provided an overview of centroid size variation of the common mosquito species in Germany. As a multidirectional size, the centroid size is more sensitive to detect sample variation compared to one-dimensional measurements, e.g. wing length 49 . Moreover, it is an adequate indicator for the wing and body size 38 and thus supply useful background information for the discrimination of certain species pairs as we could show in a pairwise comparison of the centroid size per mosquito species. However, accurate species identification solely based on wing size is not feasible 50 . The analysis of the superimposed shape coordinates produced by the Generalised Procrustes analysis is by far more informative for species identification 51 .
The results of the linear discriminant analysis based on the wing shape variables confirmed the high reliability of mosquito species identification through geometric wing morphometrics. The five most common mosquito   18 , who compared the wings of the genera Aedes, Anopheles and Culex by using the same combination of landmarks as in this study. A high reclassification accuracy was also observed on the species level, except three species pairs within the genus Aedes: Ae. annulipes vs.
Ae. cantans, Ae. cinereus vs. Ae. rossicus and Ae. communis vs. Ae. punctor. Although the Mahalanobis distances of the species means derived from a canonical variate analysis clustered together in the NJ tree, an accurate differentiation of individual specimens among these three species pairs by their wings only does not seem appropriate. This is also reflected in the analysis of the COI sequences of these species. As demonstrated before, the species identification of mosquitoes using DNA barcoding of the COI gene can give ambiguous results. This in particular applies to members of the genus Aedes or Culex 52,53 . One explanation might be a close evolutionary www.nature.com/scientificreports/ relationship between the species, which is reflected in a high similarity of the wing geometry as shown for members of the genus Anopheles 19,54 . At least the females of these species pairs are morphological similar and difficult to distinguish by classical taxonomic features. In addition, these pairs are characterized by a similar breeding ecology and are often found in sympatry 14 . Further research is required for a clear distinction between these three pairs of Aedes species. This should include a multidisciplinary approach combining biological, ecological and molecular information 55 . Shape and size of mosquito wings are not only species-specific. Geographically separated populations of the same species can have significant intraspecific variation in their wing morphology 22,56,57 . Furthermore, abiotic and biotic conditions of breeding habitats can influence the wing geometry through carry-over effects from the immature to the adult stage [58][59][60][61] . Those factors have to be kept in mind when choosing mosquito specimens for the analyses of species-specific wing patterns. Specimens from a single sampling location would probably reduce shape variation per species, thus leading to an overestimation of the actual interspecific variation. Therefore, for  Table 2. Species reclassification rates in percent calculated by a cross validation test (leave-one-out method). Green cells highlight accurate reclassification. Red cells highlight incorrect reclassification. A list of species abbreviations is given in Table 1. www.nature.com/scientificreports/ each species we used specimens from at least three different sampling locations to capture a broad wing shape variance. Nevertheless, our results revealed a high reliability to distinguish mosquito species and so underpins the robustness of morphometric wing diagnosis as tool for mosquito identification, which is in line with Henry et al. 62 , who successfully discriminated specimens of Aedes aegypti (Linnaeus) and Aedes albopictus (Skuse) sampled from different parts of the world. Besides natural differences between wing from mosquito specimens of the same species, morphometric analysis can be influenced by different measurement errors 63 . Inter-observer variation of landmark collection is likely the most important source of errors. Hence, we compared the measurement by four observers revealing a small,  Table 1. Species labels are displayed at the calculated mean centroid of the first two discriminants. The labels of Ae. punctor, Ae. communis, Ae. cantans and Ae. annulipes strongly overlap.  www.nature.com/scientificreports/ but considerable observer bias. In particular, when interested in finer intraspecific shape variations this might influence the results potentially leading to misinterpretations. One reason for the observer variation is probably the presence of wing scales on the veins, which can impede the precise landmark collection 50 . The removal of the wing scales before taking pictures may reduce observer variation and should be considered in future studies, but also increase the effort of data collection. However, completely observer unbiased measurement results are impossible. Therefore, as suggested by Dujardin et al. 64 , we provide the wing pictures in a way third parties can use them for own analyses including metadata on the sampling location and date of sampling (https ://doi. org/10.5061/dryad .zs7h4 4j5s), e.g. allowing the investigation of intraspecific wing patterns in future studies. Our study supports the applicability of geometric wing morphometrics as a complementary technique for mosquito identification, which has certain advantages over other methods. Geometric morphometric is less expensive and time-consuming compared to genetic sequencing 21,65 , e.g. representing a low-cost routine alternative under semi-field conditions. In comparison with classical taxonomical identification it requires less experience, thus representing a useful tool to control own morphological identification, e.g. for inexperienced observers. Furthermore, specimens can be identified even if only one wing is preserved or if mosquitoes are stored or sampled in fluids (e.g. ethanol), which commonly result in a loss of scales relevant for accurate morphological identification of most female mosquitoes.

Ae_ann Ae_can Ae_cas Ae_cin Ae_com Ae_gen Ae_pun Ae_ross Ae_rust Ae_sƟc Ae_vex An_clav An_mess An_pb Cq_rich Cs_ann Cs_mors Cx_pip Cx_terr
In conclusion, most mosquito species chosen for this study could be distinguished reliably by geometric wing morphometrics. The accuracy obtained by geometric wing morphometrics was comparable to the results of COI barcoding, but the morphometric analysis of wings is less expensive and time-consuming to implement. Hence, wing morphometrics has a high potential to complement species identification in future entomological studies or mosquito surveillance programs in Europe.