Ecology and morphological variations in wings of Phlebotomus ariasi (Diptera: Psychodidae) in the region of Roquedur (Gard, France): a geometric morphometrics approach

Phlebotomus ariasi Tonnoir, 1921, is the predominant sand fly species in the Cevennes region and a proven vector of Leishmania infantum, which is the main pathogen of visceral and canine leishmaniasis in the south of France. Even if this species is widely present in Western Mediterranean countries, its biology and ecology remain poorly known. The main goals of this work are to investigate the phenotypic variation of P. ariasi at a local scale in a region characterized by climatic and environmental fluctuations, and to determine if slope and altitude could affect the sand fly phenotypes. Sand flies were captured along a 14 km-long transect in 2011 from May to October. At the same time, environmental data such as altitude and slope were also collected. Morphological analysis of P. ariasi wings was performed by a geometric morphometrics approach. We found morphological variation among local populations of P. ariasi. Strong shape and size variations were observed in the course of the season (particularly in June and July) for both genders. During June, we highlighted differences in wing phenotypes according to altitude for both sexes and to slope and station for females. The phenotypic variations observed in P. ariasi along the studied transect indicated these populations are subjected to environmental pressures. Nevertheless, it seems that sand flies are more sensitive to extrinsic factors in June and July, suggesting a phenotypic plasticity.


Background
Phlebotomus ariasi Tonnoir, 1921 is the predominant sand fly species in the Cevennes region [1] and one of the two proven vectors of leishmaniasis caused by Leishmania infantum in the south of France, with Phlebotomus perniciosus Newstead, 1911 [2]. The resting sites of P. ariasi are found in houses, animal sheds, caves and holes in walls ("barbacanes"), near roads or in villages. This species is active during dusk and night, present during summer in temperate regions, abundant in periurban and rural environments and it often lives close to human and domestic animal populations [3]. This species is widely present in Western Mediterranean countries [4][5][6].
Currently, biology and ecology of P. ariasi remain poorly known. During the past 10 years, the risk of emergence or re-emergence of leishmaniasis has increased in France [7]. The current expansion of the distribution of this disease underlines the need to increase knowledge on this vector. However, only few studies have been performed to understand the population structure of this species and they were based on the analysis of cuticular hydrocarbons [8], isoenzymes [9], random amplified polymorphic DNA [10], sequences of the mitochondrial cytochrome b gene [5] and microsatellite markers [11]. Until now, no study using a geometric morphometrics approach has been done to investigate phenotypic variations of P. ariasi in France. Significant morphometric wing shape variations were found in other sand fly species within populations originating from various Mediterranean regions, forming clusters or latitudinal clines in other species, such as P. papatasi [12,13] and P. sergenti [14,15]. These studies highlighted the existence of local phenotypic variations linked to environmental factors.
The main goals of this work were to investigate phenotypic variations of P. ariasi populations collected at a local scale in a region characterized by climatic and environmental fluctuations, and to determine if slope and altitude could affect sand fly phenotypes.

Study area
The field study was performed in the south of France, on the hill "le massif de l'Oiselette" located between two valleys: "Hérault" (Ganges, Hérault) and "Arre" (Le Vigan, Gard). Sand flies were collected along a 14 km transect between Saint Julien de la Nef and le Vigan, including Roquedur-le-Haut (at 601 m above sea level) (Fig. 1). Twenty stations were selected along this transect (Table 1 and Fig. 1). This region is subject to the Mediterranean sub-humid climate [16] and characterized by the presence of "Garrigue" species such as Quercus ilex and Quercus pubescens.
This geographical area is known to be an endemic area for human and canine leishmaniasis caused by L. infantum [1]. Various domestic animals that are potential hosts for sand flies were present throughout the transect, such as chicken, sheep, ducks, geese, horses, rabbits, cats and dogs. Furthermore, some stations were located in semi-rural areas and thus various wild animals could be also hosts for sand flies [17].

Sand fly collection and identification
Sand fly collections were performed monthly between May and November 2011 with CDC miniature light traps (John W. Hock Co. FL, USA) and flight interception traps (20 × 20 cm white paper cover with castor oil) Fig. 1 Map of the study area showing the main sites sampled. Red circles indicate the sampling localities [18]. Along this transect, 20 stations were selected. In 14 sampling sites (Table 1), one or two light traps were set up for two nights (inside and/or outside of houses, animal barns, etc.), operating between 18:00 pm and 08:00 am. A total of 105 light traps were set during 210 nights of trapping. In 17 stations (Table 1), a total of 3,589 sticky traps were used and a mean of 189 sticky traps by station were placed in various biotopes, inside and around human dwelling and animal housing, close to the vegetation and crevices in the walls. The sticky traps were settled up in the sampling sites for 2 consecutive days.
At the end of the sampling period, collected specimens were transferred individually into 1.5 ml Eppendorf tubes with 96% ethanol and labeled accordingly. Prior to mounting, heads, genitalia and wings of the sand flies were removed. The heads and genitalia were cleared in Marc-André solution (chloral hydrate/acetic acid) and mounted in chloral gum [3]. Specimen identification was individually verified based on the morphology of the pharynges and/or the male genitalia or female spermathecae, using the keys of Abonnenc [3], Lewis [6] and Killick-Kendrick et al. [19]. From the identifications, we selected the wings of the P. ariasi specimens, which were in the majority of our sample.

Wing preparation
A total of 374 specimens of P. ariasi were used for the geometric morphometrics analysis (186 males and 188 females) ( Table 2). Rohlf et al. [20] suggest using only one side for the pair of organs or limbs to avoid asymmetry bias between the two sides. In this study, only right wings of specimens were used. They were stained using the method previously described in Prudhomme et al. [13] and then, mounted in Euparal on labeled slides. The wing slides were photographed using a Leica Z16 APOA stereoscopic zoom dissection microscope with DFC 425 digital camera system, digitized, and archived.

Morphometric analysis
Pictures were first entered into tps-Util 1.60 [21]. Then, 16 landmarks were used for the analysis following the method of Rohlf et al. [20] with tpsDIG2 2.18 software [22] (Fig. 2). Landmarks are located at the intersections of wing veins with the wing margin and at the intersections of cross veins with major veins (Fig. 2). The morphometric analyses as well as graphical outputs were performed using various modules of the CLIC package [23]. The centroid sizes were analyzed as a size estimator using a nonparametric Wilcoxon-Mann-Whitney test or Kruskal-Wallis test followed by a post-hoc test using Mann-Whitney tests with Bonferroni correction, using the statistical package R, version 3.1.2 [24]. Centroid size is the square root of the sum of squared distances of a set of landmarks from their centroid, i.e. the square root of the sum of the variances of the landmarks about that centroid in x-and y-directions [25]. The landmark configurations were scaled, translated, and rotated against the consensus configuration by the GLS Procrustes superimposition method [25][26][27][28] in order to produce shape variables (partial warps, PW). The principal components (based on the partial warps) [25] were used to compare population samples. To assess the degree of similarity between populations, pairwise Mahalanobis distances between populations were calculated using CLIC software [23] and tested by nonparametric permutation tests with 1,000 iterations each. These distances were also used to perform a simple reclassification test for each individual. The percentage of correctly assigned individuals to the corresponding group was assessed. Finally, residual allometry (contribution of size to wing shape) was estimated by multivariate regression of PW on size.

Sexual dimorphism
The wing shape showed statistically significant differences between males and females (Fig. 3). Mahalanobis distances were significantly different between both sexes (Adjusted P-value < 0.0001). Moreover, simple reclassification scores to the sexes were equal to 96% for males and 93% for females, supporting the observed differentiation between the two groups. Centroid sizes were used as a measure of the overall wing size differences among populations. The wing size linked to the gender was found to be significantly different (Kruskal-Wallis test : χ 2 = 253.4024, df = 1, P < 0.0001), with females displaying larger wings than males (Fig. 4). The contribution of size to wing shape differentiation was 72% (Fig. 5).  Since a phenotypic difference between males and females was observed, the subsequent analyses were performed separately for each sex.

Differentiation by month
Mahalanobis distances, used to investigate wing shape variations, were significantly different between the months; between June-July for females (Ajusted P-value < 0.0001) and between June-July and July-August for males (Ajusted P-value < 0.05 after Bonferroni correction; eight components, 85.53% of total shape variance) (Fig. 6). For females, it was not possible to test shape differentiation for the other months due to the small number of specimens (Table 2). Indeed, for the discriminant analyses, the sampling size is a limiting factor; the low number of principal components does not allow explaining the variance of the dataset. Moreover, simple reclassification scores to the month groups were on average equal to 73.6% for males (58-88%) and 81% for females (80-82%), supporting the differentiation observed between samples.
These results showed a strong phenotypic structuring between months, in particular between June and July for both sexes. Since a phenotypic difference between these two last months was observed, the following analyses performed to determine the possible effects of slope, altitude and station were realized for females and males in June and July separately. It was not possible to test for the other months due to a too small number of specimens (Table 2).

Differentiation by slope
As detailed above to determine the possible effect of slope on P. ariasi phenotype, we performed analysis on wings considering each sex and the two months, June and July, separately (see Table 2).
Mahalanobis distances to study wing shape were significantly different between the slopes (Ajusted P-value < 0.01) only for females in June (Fig. 9). Moreover, simple reclassification scores to the slope groups were on average equal to 95% (80-90%), supporting the observed differentiation between samples. However, no significant difference was observed between the slopes for the females in July (Ajusted P-value > 0.05) and for the males for any month (Ajusted P-values > 0.05; 13 components, 96.5% of total shape variance) (data not shown).
The contribution of size to wing shape differentiation was 0% for females for both months and 0 and 32% for males in June and July, respectively (Fig. 10).

Differentiation by altitude
The effect of altitude on P. ariasi phenotype was tested for females and males in June and July separately (see Table 2). In June, Mahalanobis distances used to study wing shape were significantly different between the altitudinal groups 1 and 2 for females (Ajusted P-value < 0.01667) and between group 3 and all the other groups for males (Ajusted P-value < 0.00833, 11 components, 93.2% of total shape variance) (Fig. 11). Moreover, simple reclassification scores to the altitudinal groups were on average equal to 73% (65-80%) for females and 53% (43-75%) for males, supporting the observed differentiation between samples. However, in July, no significant difference was observed between the groups for females and males (Ajusted P-value > 0.01667, 19 components, 99.2% of total shape variance and Ajusted Pvalue > 0.01667, 4 components, 71.3% of total shape variance, respectively) (data not shown). To realize the analyses with males in July, the group 4 had to be removed because of a low number of specimens (n = 2).
Concerning wing size, in June, a significant effect of altitude on centroid size was observed for females (χ 2 = 7.0282, df = 2, P = 0.0298), the post-hoc test Fig. 7 Mean, standard deviation and error of centroid wing sizes by month. a Females. b Males showed significant differences between groups 1 and 2 (P = 0.022). For males, the difference was strongly significant between groups 1 and 3 (χ 2 = 12.4331, df = 3, P = 0.006) with a P-value of 0.0076 for the post-hoc test (Fig. 12). Conversely, no significant effect of altitude was found on the centroid sizes in July for both sexes (χ 2 = 4.9144, df = 2, P = 0.0857; χ 2 = 0.4848, df = 3, P = 0.9222, for females and males respectively) (data not shown).
The contribution of size to wing shape differentiation was 7 and 0% for females and 48 and 1% for males in June (Fig. 13) and 37 and 14% for females and 36 and 0% for males in July (data not shown).

Phenotypic differentiation by station
In order to observe a possible effect of stations on wing phenotype, several analyses were realized by sex in June and July. Due to low number of individuals in certain stations, we performed a serial of discriminant analyses considering first all the stations and then excluding stations with a low number of individuals in order to increase the number of principal components and thus the variance explained (Table 3).
For females in June, Mahalonobis distances were not significantly different between stations for the first analysis based on the 3 principal components included (data not shown), but significant between ST19 and ST20 for the second analysis based on 8 principal components (Table 4 and Fig. 14). In July, even if Canonical variate analysis (CVA) showed some grouping by station (data not shown), Mahalonobis distances were not significantly different between stations for all analyses (Table 4). For males, no significant difference was observed between stations for June and July (data not shown). The results of the shape analyses and the contribution of size to wing shape differentiation are synthesized in Table 4.

Sexual dimorphism
Our results showed significant differences in wing shape and size between males and females; females have larger wings. Sexual dimorphism is widely observed in many insects such as Mansonia [29] and Aedes mosquitoes [30], Anopheles, Culex and Ochlerotatus mosquitoes [31], Drosophila melanogaster [32] and D. subobscura  Fig. 9 Distribution of female individuals along the first discriminant factor of shape analysis according to slopes. This distribution was based on the partial warps in June. Abbreviations: SE, southeast slope; NW, northwest slope [33]. Thin-plate spline deformation analyses showed that the deformation was mostly present on the medial part of the wings. The medial of the wings of flying insects is known to play an essential role in flying capability [34]. Because of blood meal necessity, the female sand flies might need a greater flying capacity compared to males. This could explain the wing differences between them. A previous study based on a mark-release-recapture method showed that P. ariasi females present a dispersion ability of more than 1 km [35]. It would be Fig. 10 Regression of first discriminant factor of shape analysis on centroid size from females (a, b) and males (c, d). Horizontal axis: centroid size of the wing; Vertical axis: first discriminant factor, representing 100% of the total discrimination. This regression was based on the partial warps in June (a, c) and July (b, d). Regression line is shown. Signs indicate each individual. Abbreviations: SE, southeast slope; NW, northwest slope Fig. 11 Distribution of the individuals along the first two discriminant factors of shape analysis by altitude groups for females (a) and males (b). This distribution was based on the partial warps in June. Horizontal axis: discriminant factor 1; Vertical axis: discriminant factor 2. Altitude groups: 1 (200-300 m), 2 (300-400 m), 3 (400-500 m), 4 (> 500) interesting to measure the dispersion capacity of males in order to test this hypothesis.

Wing phenotype differentiation by month
We observed a difference of wing shape between June and July for both sexes and a difference of size only for females (smaller in July). It was not possible to make any conclusion for the other months (May, August and September) due to the low number of individuals captured. These differences observed between June and July may be due to environmental variations as it was observed in previous study [36]. The individuals captured in June and July have certainly emerged at the end of May and June, respectively. Major temperature and relative humidity differences were observed between these 2 months. Moreover, the temperature fluctuations between night and day are greater in May than in June or July. Even if the exact conditions are not well known, the climatic parameters influence the entry into diapause and the exit and also the larvae development of insects [37]. Phlebotomus ariasi, as most of northern insect species, presents a period of diapause (at the forth-instar larval stage) in order to survive over the winter and postpone the reproduction until favorable conditions [38]. The climatic differences between months, at the end and/or at the beginning of the sand fly season may be responsible of differences in the development of the larvae before and/or after the diapause [39]. Indeed, the differences in larvae development could impact the diapausing stage and/or adult phenotype.

Environmental factors and phenotypic variation in June
The data analyses revealed shape and size differences by altitude for males and females in June suggesting differential environmental and climatic pressures according to altitude. Previous studies highlighted the impact of altitude and ---  temperature (two correlated parameters) on sand fly biology [12,40,41]. The size of adults is largely influenced by temperature, with larger adults found at low temperatures and smaller adults at high temperatures [12]. It is known that larvae development is influenced by temperature conditions [12,42]; a high temperature will produce a rapid development and thus small individuals. Assuming that wing size reflects the size of specimens, centroid size should predict larger individuals at lower temperatures and thus at higher altitudes [43]. However, the size differences observed by altitudinal groups did not show significantly larger wings at higher altitudes. This correlation may not be observed in this study because of too small altitudinal differences and thus temperature fluctuations between sampling sites. We also found a significant spatial differentiation in June for females by slope and station. The absence of spatial differentiation observed in males, conversely to females could be a consequence of stronger sex-specific selection pressures such as the host availability for blood meal in the sampled stations. The strong sexual dimorphism observed also support the existence of sexspecific selection pressures (see above). Nevertheless, the reasons of differential shape and size between sexes have not been well explored [44].

Conclusions
In conclusion, this study showed phenotypic variations among local populations of P. ariasi of different types: sexual dimorphism, shape and size variations between months (June and July in particular) and different wing phenotypes according to altitude for both sexes and by slope and station for females.
The degree of phenotypic variation observed in P. ariasi populations seems to reflect the local environmental fluctuations in terms of climatic but also biotic characteristics to which these populations are subjected. These data underline the plasticity and the capacity of adaptation of these insects even at a sympatric level.
These wing variations may have an effect on the biology of P. ariasi in terms of dispersion, fitness and transmission of pathogens such as Leishmania. Indeed, some effects were recently demonstrated for other insects such as relationships between wing shape and reproductive mode of Lysiphlebus fabarum group (Hymenoptera) [45] or enhancement of flight performance by genetic manipulation of wing shape in Drosophila [46]. Further studies are necessary to investigate the impact of wing variations on sand fly biological and ecological traits.

Acknowledgments
We are particularly grateful to the EDENext PhDB group, Professor Rioux for the interesting and useful discussions and the habitants from the study area for their kindness and availability. We are grateful to Elisabetta Andermarcher for editing the English.

Funding
Mr. Lacoste (Inkerman foundation), IRD (Institut de Recherche pour le Développement) and CNRS provided financial support. This study was funded by EU grant FP7-261504 EDENext and is catalogued by the EDENext Steering Committee as EDENext457 (http://www.edenext.eu). The contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission.

Availability of data and material
The data supporting the conclusions of this article are included within the article.
Authors' contributions JP was the main investigator of the study, he was involved in sand fly capture and identification, wing preparation for geometric morphometrics, analysis and interpretation of data, drafting of the manuscript. CC and CT contributed to sand fly capture and identification, and wing preparation for geometric morphometrics. NR contributed to sand fly capture, dissection and identification. MH, and BV contributed to sand flies capture. JPD contributed to wing geometric morphometrics analysis. BA was involved in conception and design of the study, sand fly capture, interpretation of data and manuscript revising. DS and ALB were involved in conception and design of the study, sand fly capture, interpretation of data, critical revising of the manuscript and acquisition of funding. All authors read and approved the final manuscript.