Evidence of Wing Shape Sexual Dimorphism in Aedes (Stegomyia) albopictus in Mallorca, Spain

The Asian tiger mosquito Aedes albopictus (Skuse, 1894) is a highly invasive species widely distributed on the Spanish Mediterranean coast and the Balearic archipelago. Most studies involving this species in Spain have been focused on surveillance and control methods. However, micro-evolutionary studies for Ae. albopictus in Spain have been traditionally neglected. Morphological diversity could be the result of long-term evolutionary diversification in responses to selective pressures such as temperature, precipitation, food availability, predation, or competition that may influence flight activity, host-seeking, and blood-feeding behavior. Wing geometric morphometric have been used not only to study micro- and macro-evolution in mosquitoes but also in studies of population structuring and sexual dimorphism. Therefore, the main goal of this study was to investigate the wing shape patterns of Ae. albopictus populations to unveil sexual dimorphism that could provide information about their ecology and behavior. Mosquito eggs were collected using oviposition traps at the main campus of the University of the Balearic Islands (Palma de Mallorca, Spain) and reared under laboratory conditions. In order to study wing shape variation patterns in Ae. albopictus males and females, the left wing of each adult mosquito was removed and analyzed based on 18 landmarks. Our results indicated strong levels of sexual dimorphism between Ae. albopictus males and females. Furthermore, according to the cross-validated reclassification test, males were correctly distinguished from females with an accuracy of 84% and females from males 75%. We observed a significant sexual dimorphism in the wing shape patterns of Ae. albopictus when considering different seasonal patterns (spring vs. autumn). Our results suggested that selective pressures may affect males differently to females. Host-seeking, blood-feeding, and oviposition behavior of females may act as a major driver for wing shape sexual dimorphism. These results should be considered for the development of more effective and targeted mosquito control strategies.


INTRODUCTION
The Asian tiger mosquito Aedes albopictus (Skuse, 1894) is a highly invasive species currently distributed throughout the Mediterranean area including the Balearic archipelago. This mosquito species is native to south-west tropical and subtropical forests in Asia and it has spread through all the continents except Antarctica (Kraemer et al., 2015;Cunze et al., 2016). Due to its ecological plasticity, Ae. albopictus is well adapted to exploit the resources available in urban and peri-urban areas (Bonizzoni et al., 2013).
The introduction of Ae. albopictus to new geographic regions has been facilitated by globalization processes such as transport of commodities and by anthropogenic modifications of natural environments generating suitable breeding sites (Schaffner et al., 2013;Wilke et al., 2020). Worldwide trade of tires and lucky bamboo (Dracaena spp.) have been pointed out as the main pathways for its spreading (Eritja et al., 2017;Sánchez et al., 2017;Wilke et al., 2020) and recent studies have also shown that cars and commercial flights are major drivers in the expansion of Ae. albopictus to new areas (Eritja et al., 2017;Ibañez-Justicia et al., 2017).
Aedes albopictus was first introduced in Europe in Albania in 1979. Since then, it has spread throughout the continent (Gratz, 2004) and it has been recorded in at least 31 European countries Petrić et al., 2018). The most recent recordings include the south of England (Medlock et al., 2017), the north of Portugal (Osório et al., 2018), the region of Bohemia in the Czech Republic (Rettich and Kulma, 2018), and the southern Italian islands of Lampedusa, Linosa and Pantelleria (Di Luca et al., 2017).
Aedes albopictus alongside with Aedes aegypti (L.) is considered the main vector of dengue and chikungunya viruses causing major disease outbreaks in the Southwestern Indian Ocean, India, and Central Africa (Paupy et al., 2009). This species is also a competent vector for other arboviruses such as Zika, yellow fever, and filarial worms (Paupy et al., 2009;Grard et al., 2014;Gutiérrez-López et al., 2019). In Europe, Ae. albopictus is currently responsible for local outbreaks of chikungunya (Rezza et al., 2007) and dengue (Gould et al., 2010;Lazzarini et al., 2020) due to virus introduction by travelers arriving from endemic areas (Emmanouil et al., 2020). Similarly, local transmission cases of the dengue virus were reported from Italy, Croatia, and France since 2010 (La Ruche et al., 2010;Gjenero-Margan et al., 2011;ECDC, 2020). More recently and for the first time, six dengue autochthonous cases were reported in 2018 and one in 2019 in Spain (ECDC, 2019;Monge et al., 2020) in areas where Ae. albopictus is considered the primary vector.
In the Balearics, Ae. albopictus is considered established but increasing its distribution area every year (Tavecchia et al., 2017;European Centre for Disease Prevention and Control and European Food Safety Authority, 2020). This species breeds preferably in small containers, mostly human-made, that are present in urban and peri-urban areas. Aedes albopictus is known for being anthropophilic and a nuisance species. However, despite its current importance as vector and nuisance species, little information is available about the main drivers responsible for modulating its ecology, behavior, and genetics. Consequently, identifying the possible microevolutionary patterns that may affect the behavior of Ae. albopictus (Suesdek, 2019) is crucial. Environmental factors can affect mosquito bionomy traits such as longevity, body size, reproduction, larval development, and fecundity, among others (Chandrasegaran et al., 2020). Those factors vary across the mosquito spatial distribution and may change its behavior related to its flight activity, host-seeking, or feeding frequency (Bara et al., 2015;Carvajal et al., 2016;Chandrasegaran et al., 2020).
Wing geometric morphometrics (WGM) is a practical tool for describing phenotypic variation in organisms representing individuality by the relative position of morphological landmarks that define the shape of the morphological trait studied (Klingenberg, 2010(Klingenberg, , 2011(Klingenberg, , 2016Sánchez et al., 2017;Multini et al., 2019). In insects, the wing is the preferred structure for morphometric analyses due to its bi-dimensional composition reducing digitizing error (Dujardin, 2008).
This methodology has been proven effective to study micro-and macro-evolution in mosquitoes and also in studies of population structuring and sexual dimorphism (Dujardin, 2011;Louise et al., 2015;Virginio et al., 2015;Wilke et al., 2016;Lorenz et al., 2017;Multini et al., 2019). The population structuring may be affected by exogenous or endogenous pressures that could affect males and females differentially (Carvajal et al., 2016;. Morphological diversity could be the reflection of biological differences such as shape, that result from long-term evolutionary diversification (Zelditch et al., 2012). This dissimilarity in shape would mean different functional roles or responses to selective pressures (Zelditch et al., 2012).
According to Lorenz and Suesdek (2020), the wing is an important structure for sexual signaling and flight, but it is still unknown if variation in wing patterns have some ecological or behavioral role that directly influences the flight or mating. Hence, the study of sexual dimorphism patterns can provide answers to questions about mating behavior, genetic drift, and how populations react to selective pressures, including host-seeking (Chandrasegaran et al., 2020). Comprehensive knowledge of the morphometrical characteristics of the vector may help in describing the population diversity, morpho-ecological traits and could then be useful for success in control campaigns (Jirakanjanakit et al., 2005;Sendaydiego et al., 2013;Chaiphongpachara et al., 2019).
The goal of this study is to investigate the sexual dimorphism based on wing shape variations in Ae. albopictus populations from the Balearic Islands and to explore if Ae. albopictus male and female populations are driven differently by selective pressures. In addition, the study also assessed the wing shape variation between different seasons.

Mosquito Sampling
The sampled area was located in the municipality of Palma (Mallorca, Spain), the most urbanized and populated city on the island (Figure 1). Mallorca has a Mediterranean climate represented by mild wet winters and hot dry summers with mean annual precipitation of about 500 mm and a mean annual temperature of 17 • C (Gelabert et al., 2003).
A total of 72 black plastic container ovitraps (Ø 9 cm, h = 15 cm, 950 mL) were placed to collect Ae. albopictus eggs at the main campus in an experimental plot of the University of the Balearic Islands (UTM 31 S 469347.93 m E 4387387.76 m N). This suburban area is characterized by having residential houses, usually provided with gardens and human-made small containers (plant pots, vases, etc.) which are favorable for Ae. albopictus development. Ovitraps were filled with 500 mL of tap water and 2.5 g of hay as an attraction lure. Wood tongue depressors were placed inside ovitraps as oviposition substrate (ECDC, 2012). Thirteen samplings were performed during autumn 2017 and spring -summer 2018, with six samplings between October 2017 and December 2017, and seven between April 2018 and July 2018. From 21st September to 20th December was considered as Autumn, from 21st December to 20th March was Winter, from 21st March to 20th June was Spring and from 21st June to 20th September was Summer.
Collected eggs were hatched in the laboratory under controlled conditions (27 ± 1 • C, 70% relative humidity, 12 h photoperiod), flooding them into tap water in containers until the emergence of adults. The emerged adults were maintained ad libitum in 10% sucrose solution, morphologically identified, sexed, and preserved at −20 • C when dead.
In addition to ovitraps, five mosquito adults' samplings were performed using one BG-Sentinel trap (BG-1 Sentinel, Biogents, Germany) baited with BG-lure (Biogents, Germany) and CO 2 (2 Kg of dry ice per day). Adult traps were set for 24 h per sampling between 23/10/2017 and 07/11/2017. By conducting a large sampling effort and different collection methods we greatly reduced the probability of collecting sibling specimens. For the analysis, a total of 50 adults from autumn were collected with the BG-Sentinel trap. The remaining mosquitoes (191) were reared from the egg batches collected by the oviposition traps deployed during spring and summer. To further reduce the probability of analyzing sibling specimens, when the sampling sufficiency was reached (more than 30 mosquitoes per population) adult mosquitoes were chosen randomly for the analyses.

Geometric Data Acquisition
A total of 241 left wings, 123 from females and 118 from males emerged in the laboratory and collected from the field in different seasons ( Table 1) were dissected and mounted between a slide and a coverslip (without mounting medium) fixing the coverslip with adhesive tape Beriotto et al., 2021). Wings were photographed at 40x magnification using a Zeiss Scope A.1 stereoscopic microscope attached to a camera AxioCamICc1 (Zeiss, Germany) by one of the authors (JLM) to minimize error measurements (Fruciano, 2016). The images were processed using ZEN 2.3 lite program (Blue edition). Eighteen landmarks (LMs) were selected to study the geometric morphometry of Ae. albopictus wings according to Bookstein (1997); Vidal et al. (2011), and . The coordinates of the 18 landmarks represented by vein intersections (Figure 2) were obtained using tpsDig2 2.30 (Rohlf, 2005;Sendaydiego et al., 2013). These landmarks were selected for being homologous and found in all representatives of the Culicidae family (Lorenz et al., 2017;Beriotto et al., 2021).

Morphometric Analysis
The wing geometric morphometric analysis was conducted using MorphoJ 1.06d software (Oracle Corporation) that provides a wide range of tests for two-or three-dimensional landmark data (Klingenberg, 2011). Classifiers about sex and season information were generated for each specimen to examine morphometric differences. Firstly, allometric effect of the wing size on the wing shape was estimated using multivariate regression analysis of the Procrustes coordinates on the centroid size with a permutation test of 10,000 randomizations, and, finally, removed from the subsequent analyses (Klingenberg, 2016). Due to the different sample methods, we focused on wing shape variation patterns disregarding size variations, coordinates, and orientation of the landmarks (Klingenberg, 2010;Börstler et al., 2014;Virginio et al., 2015). After that, to extract shape information, landmarks were subjected to the generalized least-squares Procrustes superimposition algorithm aligning by the principal axes in order to standardize the shape of wings and determine the wing shape coordinates (partial warps) for sex (Virginio et al., 2015;Carvajal et al., 2016) and then for season classifiers . After that, a covariance matrix was performed to assess the linear relationship between the independent and the dependent variables. Principal Component Analysis (PCA) was carried out using a covariance matrix for all populations. Then, a Canonical Variate Analysis (CVA) with a permutation of 10,000 randomizations was carried out to test the pairwise distances between males and females, and between seasons (Klingenberg, 2010;Carvajal et al., 2016). In the visualization of the CVA, every possible shape corresponds to a shape point in the morphospace and every direction in the morphospace corresponds to a specific shape change (Klingenberg, 2010). Subsequently, thin-plate splines were obtained by the regression analysis of CVA scores against wing shape variation (Bookstein, 1989;Wilke et al., 2016). Each specimen was then reclassified using the cross-validated classification test of the discriminant analysis based on Mahalanobis distances to determine the degree of wing shape dissimilarity between sex and seasons Wilke et al., 2016). Finally, the Mahalanobis distances resulted from the CVA analysis with season classifier were used to group clusters defined by unweighted pair group method with arithmetic mean (UPGMA, Sneath and Sokal, 1973) to assess the similarity of the wing shape patterns (Carvajal et al., 2016). UPGMA was achieved using PAST software.

RESULTS
The allometry test predicted 0.81% (P = 0.0790) of the wing shape variation. The allometric effect was not considered as it only explained a negligible proportion of the variance. The Procrustes superimposition analysis using the classifier sex resulted in significant differences in the wing shape patterns of males and females. In the wireframe, it was observed that male wings were narrow and lengthened, whereas females were wider and shorter. The CVA showed significant differences in the wing shape patterns between males and females ( Figure 3A) and, in the principal-partial warps from the PCA analysis, it was observed that the most variable anatomical landmarks Frontiers in Ecology and Evolution | www.frontiersin.org FIGURE 3 | Comparison between male and female specimens. (A) CVA represented in the morphospace to compare the wing shape variation patterns between males (blue) and females (red). X-axis: First canonical variate; Y -axis: Frequency of each wing shape. (B) Wireframe graph from the PCA displaying the mean value of wing shape variation in the male (blue) and female (red) populations.
Female specimens yielded high values in the cross-validated classification test from the discriminant function, yielding 94.31% of accuracy when compared to male specimens ( Table 2). Conversely, males also yielded high cross-validated classification values, resulting in 93.22% accuracy when compared to females (Figure 3 and Table 2).
The CVA considering seasonality showed significant differences in wing shape between seasons of the year for both sexes (Figure 4). Discriminant analysis revealed that the dissimilarity in wing shape separately for males and females over the seasons was statistically significant (P < 0.05) except between these groups: Females-Spring vs. Females-Summer (P = 0.2275), Males-Autumn vs. Males-Spring (P = 0.1542), Males-Autumn vs. Males-Summer (P = 0.2900) and Males-Spring vs. Males-Summer (P = 0.1151). From the 482 total comparisons carried out with the 241 wings, 350 of them resulted in the cross-validated classification values based on Mahalanobis distances yielding scores higher than 70%, indicating significant differences in wing shape ( Table 2). In general, Ae. albopictus collected during the summer yielded lower reclassification values (38-79%) suggesting a higher overlap in wing shape patterns between males and females in this specific season.
Constructed UPGMA tree based on Mahalanobis distances highlighted the segregation between Ae. albopictus sex and seasons, showing a Cophenic correlation of 0.9767 (Figure 5). Male and female groups were segregated into two separate clusters. In males, spring and summer samples were clustered together, whereas Autumn was clustered apart from the others. Results indicate that the Autumn male population has a higher dissimilarity when compared to the other populations. Conversely, Autumn and Spring were clustered together in female population, and summer was clustered apart.

DISCUSSION
Since the introduction of Ae. albopictus in the Balearic Islands in 2012 (Miquel et al., 2013), it has increased its abundance (Tavecchia et al., 2017) and has become a public health concern in many places of the Mediterranean coast of Spain . In Europe, emerging and re-emerging outbreaks of diseases transmitted by mosquito vectors are increasing (Martinet et al., 2019;Emmanouil et al., 2020). For that reason, it is essential to know the bioecology of vector species as well as the preferred breeding sites, host preference, and dispersal capacity among other characteristics.
Our results showed significant wing sexual dimorphism between male and female Ae. albopictus specimens collected either as eggs or adults in a suburban area of Palma city. Male   wings were narrower and lengthened, whereas females were wider and shorter. In our study, the highest differences observed in the wireframe were found at the central contraction zones of the wing (LM 2, 10, 17, and 18) most likely due to the robust and mechanical stability of the wings. These wing morphometrics differences could be a consequence of different selective pressures such as sex-specific behaviors (i.e., host-seeking or mating) Chandrasegaran et al., 2020). In fact, shape variation in sexual dimorphism studies is mainly concentrated in the landmarks located between the middle and distal regions of the wing (Lorenz et al., 2017). Similar results of sexual dimorphism based on wing shape using WGM were observed in Ae. aegypti species collected in urban areas in the Philippines (Sendaydiego et al., 2013;Carvajal et al., 2016). Similarly, such sex-specific variation in wing shape patterns has been also confirmed within ten mosquito species belonging to the genus Culex, Aedes, and Anopheles in southeastern Brazil (Virginio et al., 2015). Furthermore, our results indicated that the wing pattern variation was associated with seasonality, finding that shape overlapping occurred on mosquitoes of both sexes collected in spring. On the contrary, female mosquitoes collected during autumn showed lower levels of wing shape pattern similarity in comparison to summer and spring. Our findings indicate that different conditions (extrinsic or intrinsic) of each season may significantly affect Ae. albopictus wing shape in male and female specimens and possibly other morphological and physiological traits (Virginio et al., 2015).  and Louise et al. (2015) used the same methodology to demonstrate that short periods of time (less than a year) are enough to result in detectable wing shape changes (Suesdek, 2019).
Our results suggest that microevolutionary changes affected the wing shape patterns of the Ae. albopictus populations analyzed in this study. Wing shape pattern evolution may express flight activity changing behavior, intrinsic to the hostseeking process, host preference, vectorial capacity, and habitat suitability, among others (Morales Vargas et al., 2013;Suesdek, 2019;Chandrasegaran et al., 2020).
The evolution of populations within species and how this process impacts vector biology and epidemiology is not the only consequence of biotic factors (Suesdek, 2019;Chandrasegaran et al., 2020). Abiotic local factors, such as temperature, humidity, landscape, and presence of breeding sites could influence microevolution (Chandrasegaran et al., 2020).
There are several ecological differences between males and females of Culicidae that may explain the difference in wing shape of both sexes (Virginio et al., 2015;Christe et al., 2019). Usually, males remain near breeding sites, feeding on sugar sources, such as flowers, or resting in vegetation waiting for females to mate. In addition, females must seek hosts to obtain blood to perform oogenesis and look for aquatic habitats for oviposition. These ecological selective pressures could be responsible for the morphological variations between sex. Other studies such as de Oliveira  and Christe et al. (2019) also demonstrated that males and females of Aedes fluviatilis (Lutz, 1904) react to selective pressures such as anthropogenic changes in the environment differently. Further, Nasci (1986) studied the wing length between host-seeking and non-hostseeking females of Ae. aegypti and found differences. These results concur with our results, indicating that wing shape patterns are sex-specific and reflect the mechanisms employed by each sex to cope with selective pressures. Landscape heterogeneity was not measured in this study; however, wing shape variations may be attributed to host and habitat availability causing local population changes (Carvajal et al., 2016;Christe et al., 2019). In fact, urbanization can act as a driver of microevolution events that can be observed in the wing shape of Culex quinquefasciatus in São Paulo (Brazil) . Multini et al. (2019) studied the population structure of Anopheles cruzii (Dyar and Knab, 1908) populations for 3 years from three locations with different urbanization levels (urban, peri-urban and sylvatic), finding different wing morphometrics patterns between years and collection sites. Their results suggested that urban disturbances may affect mosquito biology that is reflected with a phenotypic change in the wing. In this sense, evidence of habitat-related morphometric variation was found in urban parks in Sao Paulo, Brazil, where De Carvalho et al. (2017) observed heterogeneity in the wing shape of Culex nigripalpus (Theobald, 1901) females. In the case of Ae. albopictus, habitat disturbances may lead to the creation of new breeding sites (mostly human-made) compared to natural breeding sites, allowing the rapid spread of the population increasing the contact between mosquito vectors and human hosts. However, this important evolutionary force remains mostly unexplored (Medlock et al., 2015).
To the best of our knowledge, this is the first geometric morphometric work with Ae. albopictus conducted in Spain. We have observed a significant sexual dimorphism variation in the wing shape patterns of Ae. albopictus. Furthermore, we have detected microevolutionary changes in wing shape between female populations from autumn in comparison to spring and summer. However, its impact on Ae. albopictus ecology and behavior remain unknown. Nevertheless, we can only assume that the morphometric results show us morphological differences between groups being a limitation when compared with the entire population and address functional and control strategies of the species (Richtsmeier et al., 2002). Further research increasing the number of samples per season or even include different populations must be considered to address these functional implications.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JL-M designed and performed the experiment. MM designed and directed the experiment. CB assisted with the wing mounting. AB contributed to the design of the figures and the performance of the analysis. JL-M, AB, CB, and MM wrote the manuscript. All authors agreed to be accountable for the content of the work.

ACKNOWLEDGMENTS
We would like the acknowledge the students Toni Sureda and Tania Navarro for their contribution to the field and laboratory work. Also, we would like to thank the editor and the reviewers for their time and comments which have greatly improved the manuscript.