Wing geometric morphometrics and COI barcoding of Culex pipiens subgroup in the Republic of Korea

Two members of the Culex pipiens subgroup, Culex pallens and Culex pipiens f. molestus, are known to occur in the Republic of Korea (ROK). These species exhibit morphologically similar features and are challenging to distinguish below the species level. Therefore, this study utilized wing geometric morphometrics (GM) on the right wing of the Culex pipiens subgroup, alongside sequencing of the cytochrome c oxidase subunit I (COI) region. Mosquitoes were collected from 11 locations between June and October (2020–2022) to minimize regional and seasonal variations. Additionally, Culex pipiens f. pipiens, which is not native to the ROK, was included in the analysis. Culex tritaeniorhynchus, Aedes albopictus, and Anopheles sinensis, the primary vectors in the ROK, were used as outgroups for comparison. All three taxa in the Culex pipiens subgroup could be identified with an 82.4%–97.0% accuracy using GM. However, a comparison of the COI regions of the Culex pipiens subgroup revealed no clear differences between the taxa. These data can be used for accurate identification, contributing to effective mosquito control, in addition to providing a foundation for evolutionary and ecological studies on wing shape differences.

methods are needed to compensate for the occasional limitations of molecular approaches.Additionally, GM has been actively employed in mosquito research, proving successful in various fields [24][25][26] .However, there is currently a severe lack of quantitative morphological analyses of pathogen vectors in the ROK.
Therefore, this study examined Cx. pallens and Cx.pipiens f. molestus of the Culex pipiens subgroup in the ROK using GM.Furthermore, Cx. pipiens f. pipiens of the Culex pipiens subgroup, which is absent in the ROK, was included in the analysis.Furthermore, major pathogen vectors in the ROK, including Culex tritaeniorhynchus Giles (the causative agent of Japanese encephalitis), Aedes albopictus Skuse (the causative agent of Zika), and Anopheles sinensis s.s.Wiedemann (the causative agent of Malaria), served as outgroups in the analysis [27][28][29] .Finally, sequencing of the DNA barcoding region [cytochrome c oxidase subunit I (COI)] was conducted to compare the results of the quantitative morphological analysis using GM.

Sample collection and identification
To eliminate regional and seasonal biases, specimens of the Culex pipiens subgroup were gathered from a total of 11 different sites in the ROK spanning from June to October (Fig. 1).The selected sampling sites were located around buildings in urban areas, taking into consideration the habitat preferences of the Culex pipiens subgroup BG-sentinel traps (Biogents, Regensburg, Germany) with dry ice were employed for collecting the Culex pipiens subgroup.Additionally, black light traps (BioTrap, Seoul, ROK) were deployed in a cow shed to capture Cx. tritaeniorhynchus and An.sinensis, which served as outgroups.As for Ae.albopictus, those captured alongside the Culex pipiens subgroup in Daegu were utilized.All captured mosquitoes were preserved at − 20 °C until identification.

GM analysis
Only the right wings of the identified female mosquitoes were selectively used for the GM analysis.The right wing of each mosquito species was dissected, and wing scales were removed using a brush.Afterward, the wings were mounted on microscope slides and covered with coverslips using Canada balsam (Duksan, Seoul, ROK).These prepared wing specimens were photographed under 20 × magnification using an Olympus SZ61 Stereo Microscope (Olympus Corp., Tokyo, Japan).Wing images of Cx. pipiens f. pipiens collected in Germany, a species not distributed in the ROK, were incorporated into our experiments 32,33 .
Landmark coordinates were collected using TPSdig2 (2.31) (Fig. 2), and 17 LMs were established for the GM analysis of mosquito wings 34 .Beriotto et al. 35 determined that 9, 13, and 17 LM settings produced the same total error rate in studies involving Cx. pipiens.Since this study involved three taxa of the Culex pipiens subgroup, 17 LMs were chosen for a more precise comparison.Subsequently, the TPS files for each mosquito taxon with established LMs were analyzed using various packages in R v.4.2.1 36 .
Procrustes analysis was conducted with the 'geomorph' package (v.4.0.5) to transform the difference in position and direction of each data coordinate 37 .The results were superimposed to generate Procrustes coordinates.The adequacy of the established LMs was confirmed using the 'LamBDA' package (v.0.1.1) 38.Centroid size (CS), defined as the square root of the sum of the squared distances of all landmarks on the object from their centroid, was calculated from the Procrustes coordinates to estimate wing size 39 .The allometric effect, which represents the relationship between shape and size, was also evaluated using the 'geomorph' package (1000 permutations) 40 .Mean CS values for each species were compared using analysis of variance (ANOVA), and pairwise comparisons were conducted using t-tests with the 'agricolae' (v.1.3-7)package (Bonferroni adjusted p-value) 41 .To analyze differences in wing shape between groups, linear discriminant analysis (LDA) was performed using the 'tidyverse' (v.2.0.0) and 'MASS' (v.7.3-60) packages 42,43 .The convex hull algorithm, encompassing all points for each group in the scatterplot, was also applied.To compare the wing shape of each taxon, mean shapes for each LM coordinate were calculated and compared.The Mahalanobis distance, a statistical indicator of the difference between each group, was calculated using the 'CVA' function of the 'Morpho' package (v.2.11) with 10,000 permutations 44 .Finally, a jackknife leave-one-out cross-validation was performed to determine the accuracy of the classification based on wing shape.www.nature.com/scientificreports/

COI barcoding
Genomic DNA was extracted from the whole body, except the right wing, of the species-identified mosquitoes using the Clear-S™ Quick DNA Extraction Kit (InVirusTech, Gwangju, ROK), following the manufacturer's protocol.
The mean CS differed significantly between mosquito species (F = 125.7,p < 0.0001).On average, An. sinensis exhibited the largest CS, whereas Ae. albopictus had the smallest CS (Fig. 3).Pairwise comparisons of Cx. pallens vs. Cx.pipiens f. molestus, Cx. pallens vs. Cx.pipiens f. pipiens, and Cx.pipiens f. molestus vs. Cx.pipiens f. pipiens showed no statistically significant differences in CS (p = 1) (Supplementary Table S1).The allometric effect was found to be statistically significant (R 2 = 3.3%, p < 0.0001), and it was analyzed together with the shape data, as it could provide useful information for the identification process 54,55 .
Among the outgroup species, An. sinensis and Ae.albopictus were clearly separated, as shown in the LDA scatterplot results for all species analyzed together (Fig. 4a).Cx. tritaeniorhynchus, also set as an outgroup, was generally well distinguished from mosquitoes in the Culex pipiens subgroup, with minor overlap in some specimens.Species within the genus Culex pipiens subgroup clustered together, displaying slight overlap for some species.Differences in mean wing shapes were observed in LM1-2 and LM16-17 for An.sinensis and Ae.albopictus, and mainly in LM1-2 and LM12-17 for An.sinensis and Culex mosquitoes (Fig. 4b).Significant differences in mean wing shape between Ae. albopictus and Culex mosquitoes were noted, particularly at LM2 and LM16-17.Subsequent LDA using only Culex mosquitoes revealed a clear separation of each taxon (Fig. 5a).A comparison of mean wing shapes within the genus Culex (Fig. 5b) indicated no major differences in wing shape for each taxon, with some distinctions primarily in LM16-17.The Mahalanobis distance calculated using canonical variate analysis showed that the closest taxa were Cx.pallens and Cx.pipiens f. pipiens (4.3152, p < 0.0001), followed by Cx. pallens and Cx.tritaenioryhnchus (4.4569, p < 0.0001).The most dissimilar species were An.sinensis and Cx.pipiens f. pipiens (17.2230, p < 0.0001) (Supplementary Table S2).Cross-validation (leaving-one-out method) for classification accuracy showed an overall average accuracy of 96.8% ( As shown in the ML tree, An. sinensis, Ae. albopictus, and Cx.tritaeniorhynchus, designated as outgroups, formed a monophyletic clade with strong bootstrap support.However, species within the Cx.pipiens subgroup were unresolved (Fig. 6).

Discussion
Cx. pipiens f. molestus and Cx.pallens differ in the D/V ratio of male genitalia and ommatidial number of females 56 .However employing such methods for identification is challenging for inexperienced researchers, timeconsuming, and expensive, particularly when dealing with large numbers of individuals.Molecular biological approaches serve as valuable tools to overcome the challenges of morphological identification 23,57,58 and remain the most powerful methods for accurate species identification.In response to these challenges, the analysis of mosquito wings using GM has gained acceptance and has proven to be a reliable and efficient method, complementing traditional morphological classification methods 24,32,54,55 .This study compared three taxa within the Culex pipiens subgroup (Cx.pipiens f. molestus, Cx. pipiens f. pipiens, and Cx.pallens) with Cx. tritaeniorhynchus, Ae. albopictus, and An.sinensis as outgroups.The results demonstrated that species identification using GM can be successfully applied not only at the level of evolutionarily distant genera but also in very closely related relationships.However, cross-validation results indicated that reliable classification of Cx. pipiens f. pipiens could not be achieved with high accuracy (14/17:82.4%).For example, Cx. pipiens f. pipiens and Cx.torrentium, known to co-occur sympatrically in Europe, cannot be easily differentiated via traditional morphological approaches, whereas GM was able to successfully discriminate them 59 .There are also examples of high heterogeneity in  spatially separated populations of the same species 60 .However, in the case of Cx. pipiens f. pipiens used in this study 32,33 , it was not clearly identified using GM, despite its absence in the ROK and geographical separation.Therefore, further studies with larger geographic sampling are necessary.Nevertheless, only Cx. pipiens f. molestus and Cx.pallens of the Culex pipiens subgroup are currently known to occur in the ROK.GM-based methods are thus expected to be successful in identifying various pathogen vector mosquitoes, including the Culex pipiens subgroup in the ROK.Analyzing wing size using CS as an indicator revealed no statistically significant differences in wing size among the three taxa of the Culex pipiens subgroup.Furthermore, wing size in insects, including mosquitoes, is strongly influenced by the environment 61 , suggesting that using wing size for discrimination is not suitable for classifying close relatives like the Culex pipiens subgroup.However, LDA based on wing shape variation demonstrated that wing vein structure is homologous, species-specific, and less influenced by the environment than wing size.This makes it an extremely useful morphological trait for species identification, phylogenetic analysis, and tracking evolutionary relationships 62 .
In this study, the analysis of the Culex pipiens subgroup was based on data from Cx. pallens and Cx.pipiens f. molestus from the ROK and Cx.pipiens f. pipiens from Germany 32,33 .While Cx. quinquefasciatus, another member of the Culex pipiens subgroup, was not included in this experiment, its wide distribution in tropical regions and tendency to co-occur with other subgroup members suggests the need for further experiments to determine the feasibility of accurate species identification using GM compared to other members of the Culex pipiens subgroup 8,16,63 .As suggested by Dujardin et al. 64 , the wing images used in this study will be made available through an open-access repository, allowing researchers to use them in their own studies (https:// doi.org/ 10. 5061/ dryad.gtht7 6ht6).
The sequences of the COI region of Cx. pallens and Cx.pipiens f. molestus from the ROK was compared with that of Cx. pipiens f. pipiens, but no clear differences were observed.The COI region is commonly used for accurate species delimitation in insects and identifying genetic divergence in closely related species, and it has been extensively successful 65 .In the case of the Cx.pipiens subgroup, the influence of Wolbachia infection on mitochondrial diversity and the occurrence of sympatric populations capable of hybridization may contribute to limited genetic variation in the COI region.Thus, future comparisons between populations should adopt a multi-loci approach, considering both mitochondrial and nuclear genes 12,66,67 .
In the ROK, despite the Culex pipiens subgroup being a diverse pathogen vector, its accurate identification been challenging due to difficulties in morphological classification 5,20,68 .The successful use of GM in classifying the Culex pipiens subgroup, along with Cx. tritaeniorhynchus, Ae. albopictus, and An.sinensis, suggests the potential applicability of GM-based classification to various mosquito species in the ROK.This method could benefit researchers unfamiliar with mosquito taxonomy and those facing challenges with molecular biological approaches, providing a means for accurate species identification.Furthermore, GM analysis is not confined to interspecific variation and has been applied in various studies, including sexual dimorphism and parasite detection 69,70 .Therefore, the results of this study can serve as foundational data for future ecological and evolutionary research on the Culex pipiens subgroup.The wing images utilized in this study are also expected to be valuable resources for researchers in various fields.

Figure 2 .
Figure 2. The 17 LMs set up on the right wing of a female adult mosquito for GM analysis.

Figure 3 .
Figure 3. Boxplot showing the results of centroid size (CS) comparisons for each taxon (SIN = An.sinensis, PAL = Cx.pallens, MOL = Cx.pipiens f. molestus, PIP = Cx.pipiens f. pipiens, TRI = Cx.tritaeniorhynchus, and ALB = Ae.albopictus).The box represents the first and third quartiles, and the black line in the middle represents the median.The asterisks represent outliers.

Figure 6 .
Figure 6.Maximum likelihood (ML) phylogenetic tree (GTR-I).The bootstrap value is shown next to each node (1,000 replicates).Unresolved taxa are shown in red text.Drosophila melanogaster (GenBank accession number: MT807020) was used as the outgroup.

Table 1 .
Mosquitoes collected from 11 sites in the ROK between June and October.

Table 2 .
Cross-validation (leaving-one-out method) test results for six taxa based on wing shape (overall accuracy = 96.8%).