Understanding Phlebotomus perniciosus abundance in south-east Spain: assessing the role of environmental and anthropic factors

Leishmaniosis is associated with Phlebotomus sand fly vector density, but our knowledge of the environmental framework that regulates highly overdispersed vector abundance distributions is limited. We used a standardized sampling procedure in the bioclimatically diverse Murcia Region in Spain and multilevel regression models for count data to estimate P. perniciosus abundance in relation to environmental and anthropic factors. Twenty-five dog and sheep premises were sampled for sand flies using adhesive and light-attraction traps, from late May to early October 2015. Temperature, relative humidity and other animal- and premise-related data recorded on site and other environmental data were extracted from digital databases using a geographical information system. The relationship between sand fly abundance and explanatory variables was analysed using binomial regression models. The total number of sand flies captured, mostly with light-attraction traps, was 3,644 specimens, including 80% P. perniciosus, the main L. infantum vector in Spain. Abundance varied between and within zones and was positively associated with increasing altitude from 0 to 900 m above sea level, except from 500 to 700 m where it was low. Populations peaked in July and especially during a 3-day heat wave when relative humidity and wind speed plummeted. Regression models indicated that climate and not land use or soil characteristics have the greatest impact on this species density on a large geographical scale. In contrast, micro-environmental factors such as animal building characteristics and husbandry practices affect sand fly population size on a smaller scale. A standardised sampling procedure and statistical analysis for highly overdispersed distributions allow reliable estimation of P. perniciosus abundance and identification of environmental drivers. While climatic variables have the greatest impact at macro-environmental scale, anthropic factors may be determinant at a micro-geographical scale. These finding may be used to elaborate predictive distribution maps useful for vector and pathogen control programs.


Background
Phlebotomine sand flies (Diptera: Psychodidae) are haematophagous insects that transmit Leishmania spp., protozoan parasites endemic in tropical and temperate zones, including the Mediterranean subregion [1]. Among over 800 sand fly species worldwide, 12 have been identified in Spain [2]. These include Phlebotomus perniciosus and P. ariasi, vectors of Leishmania infantum responsible for zoonotic visceral leishmaniosis in western Mediterranean countries; and P. papatasi and P. sergenti, vectors of L. major and L. tropica, respectively, that cause cutaneous leishmaniosis in Northern Africa and the Middle East.
The risk of L. infantum infection in endemic areas is geographically variable, depending on sand fly density [3,4]. Unlike mosquitoes, sand flies breed on terrestrial sites protected from desiccation and with organic matter for larvae to feed on such as animal burrows and shelters, abandoned buildings, caves and stone walls [5]. Efforts to collect immature sand fly stages from the natural environment are very unproductive, and the precise microhabitats for sand fly breeding are poorly characterised [6]. Hence, the great majority of sand fly distribution studies monitor adult stages only. Their activity is typically seasonal; they can be found over a broad altitudinal range and temperature is considered the main artifice of sand fly phenology patterns in Mediterranean countries [4]. The role of other climatic and environmental variables on sand fly abundance is still inconclusive due to the wide variety of natural habitats in which sand flies are found, and the complex interconnections between the multiple factors affecting sand fly biological cycles. Species may have preferential macrohabitats, and in western Europe P. perniciosus is widespread while P. ariasi prevails in cooler, more humid regions [2,7,8]. Locally, sand fly presence and abundance may vary depending on climate, orientation, predominant vegetation, soil types, proximity to livestock and other factors [6,[9][10][11][12][13][14]. Accurate mapping of sand fly densities is further constrained by the wide variability in study designs, sand fly collection methods and statistical methods used to analyse distributions. They are commonly collected using light-attraction and/or adhesive interception traps, and they may lead to significantly different sand fly density estimations [15,16]. Complex data statistical analysis is required for quantitative longitudinal study designs that recognise the strong spatial and temporal aggregation, so typical of sand fly populations [4].
The Murcia Region in southeastern Spain is endemic for canine and human leishmaniosis caused by L. infantum. Recent studies have shown that asymptomatic infection is widespread in rural areas and prevalence is associated with specific environmental factors [17,18]. There are no investigations of the spatial distribution of sand flies in most of the areas covered by the previous leishmaniosis studies. Surveys performed in areas close to Murcia City in the 1980s identified eight sand fly species and P. perniciosus was the most frequent vector, active between March and October, with peaks in July and September [19]. The Murcia Region is geographically and bioclimatically diverse due to its relatively large size (11,300 km 2 ), altitudinal gradient (0-2,000 m above sea level) and distance range to the Mediterranean Sea (0-200 km). This makes it an ideal place for a quantitative investigation of environmental factors driving sand fly abundance on a large geographical scale. With this objective in mind, the present longitudinal study used a standardised sampling procedure to estimate sand fly abundance in rural areas in the main five bioclimatic zones in Murcia. Mixed generalised linear models for count data were then used to analyse overdispersed sand fly distributions in relation to macro/micro environmental and human-driven factors.

Study area and design
The Murcia Region, with a permanent population of 1,470,000, has an agricultural and tourism based economy. It has a typical semi-arid Mediterranean climate, with long, dry summers and an average annual rainfall of 350 mm, which is commonly delivered over a few intense precipitation events. Sand fly abundance was monitored in 25 animal premises including 5 premises, 3 sheep sheds and 2 dog kennels, in each of the main five geographical zones that are traditionally recognized in the Murcia Region: N (north), S (south), C (central), W (west) and SE (southeast) (Fig. 1, Table 1). Premises were selected by local veterinarians based on owners' willingness to participate and the prevalence of Leishmania infection in the animals was unknown. Premises in each zone were situated within a narrow altitude range and at a similar distance to the sea (Fig. 1 Sand fly trapping devices included a miniature CDC (Centers for Disease Control and prevention) light attraction trap (J. W. Hock Company) and eight A4 castor oil impregnated white paper sheets (interception adhesive sticky traps) in each premise. The total number of light trap-days placed in the study were 200 (1 trap × 25 premises × 8 days) and similarly, the number of sticky trap-days was 1,600 (8 traps × 25 premises x 8 days). Traps were always placed in the same spot on each visit inside the building except for sticky traps in dog kennels where four were placed inside the dog house and the other half on in the open-air part outside the dog house. In all premises, sites selected for sticky traps were considered representative of different premise microhabitats, and they included wall surfaces and holes, fences and open windows. After each 24 h sampling, sticky traps and collection cups from light traps were gathered and taken back to the laboratory within a few hours.

Sand fly counting, sexing and morphological speciation
Collection cups were kept at -20°C for at least 2 h and sand flies were then counted, sexed and stored in absolute ethanol until speciated. Sticky traps were kept in the fridge and submitted to same procedure regularly and were completed by December 2015. Morphological identification was performed using entomological keys [20][21][22]. Briefly, the males were identified according to morphological features in the aedagus, stylopodite and coxopodite and the females were identified based on the pharynx, cibarium and spermatheca.
Climatic data from meteorological stations included daily averages for the study period and monthly averages for the 2006-2015 time series of the following variables: absolute maximum, maximum, mean, minimum and absolute minimum RH and temperature, maximum and total rainfall, and maximum and mean wind speed. Furthermore, annual, May-October (adult sand fly activity period in Murcia), November-April (period of null or low adults sand fly activity) mean values were calculated for climatic variables and used as independent variables in some of the multivariable regression models described below.
Data relating to the premises, animal management, structural features of the buildings, the frequency of use of disinfectants in the building and insecticides on the animals and position of the traps, were collected by inspecting and taking measurements of the building and interviewing the owner.

Statistical analysis
The distribution of sand flies and environmental variables and the association between the presence/absence and sand fly counts in positive traps and other variables were analysed. Yates-corrected chi-squared test and the non-parametric Kruskal-Wallis test were used to compare proportions and medians, respectively, and the correlation between numerical variables was assessed using Spearman's rank coefficient test. Multilevel negative binomial regression models were then developed to examine the independent contribution of environmental factors to sand fly abundance considering the correlation between repeated sand fly counts over time in study premises [24]. Two types of multilevel models were developed according to the data used: (i) Type I model used temperature and RH data from thermohygrometers, building characteristics and environmental data (other than temperature, RH and precipitation) from the buffer area around the premises, and (ii) Type II model used GIS-derived environmental data from the buffer area around the premises only. The later were developed to identify variables that could be used to generate a sand fly density map for Murcia Region and to compare outputs with Leishmania prevalence models. Environmental variables were used as fixed explanatory variables. They were fitted as categorical variables in the Type I models and as continuous variables in the Type II models. Random variation in sand fly counts between premises was considered both at the intercept and in the slope over time [24]. Briefly, this allows for variation between premises in the relationship between explanatory variables and the response (intercept variation), and for this variation to be different for each premise over time (slope variation).
A step-wise model building approach [25] was used beginning with a model including climatic variables. Other environmental variables significantly associated and showing a positive or negative trend with the outcome in the bivariate analysis, were subsequently added to the model. They included building characteristics, land use and soil and ground taxonomy variables. Due to the high correlation between environmental variables (for example between building age and type of wall material or altitude and temperature), several models including only variables significantly associated with sand fly counts were considered. Among them, the one with the lowest Akaike information criteria (AIC) were deemed the most parsimonious [25]. Parameter estimates were exponentiated to calculate incidence rate ratios. Significance was taken for alpha = 5% (P <0.05). R (http://cran.r-project.org/) program was used for all the statistical analysis.
CDC and sticky traps provided a similar distribution of sand flies (Table 1 and Additional file 1: Table S1). The percentage of positive CDC traps and sand fly abundance in CDC traps were highest in the W and S zones and lowest in the SE (Table 1). Abundance in sticky traps was also greatest in the W and lowest in the SE, while the percentage of positive sticky traps was highest in the S and lowest in the SE (Additional file 1: Table S1). In addition to differences between zones, sand fly abundance also varied significantly between premises in the same zone, particularly in CDC traps from the S and the N ( Table 1). Sand fly abundance was positively associated with altitude except that it was lower in the 536-705 m compared to 282-352 m altitude ranges (P < 0.05). These two ranges corresponded mostly, to premises located in the N and S of the region, respectively (Table 1) The species distribution in the 3,586 (98%) sand flies speciated is shown in Table 2. P. perniciosus represented 80% of all sand flies followed by P. papatasi (10%), P. sergenti (5%), S. minuta (4%) and P. chabaudi, P. longicuspis, P. ariasi and P. alexandri (less than 1%) ( Table 2). Species distribution varied according to trap type, zone and animal species premises. Phlebotomus perniciosus was relatively more abundant in CDC compared to sticky traps, in sheep than in dog premises and less abundant in the C zone compared to other zones (Table 2). Instead, the relative abundance of P. papatasi was greatest in sticky traps and C and N zones, P. sergenti in sticky traps, dog premises and W zone, and S. minuta in dog premises and in C ( Table 2). The overall sex ratio was similar for all species (Table 2). However, in CDC traps P. papatasi, P. sergenti and P. longicuspis females were substantially more abundant than males while in sticky traps, 76% of all sand flies were males (not shown).
Seasonal distribution of P. perniciosus in CDC traps and bivariate relationship with indoor temperature and relative humidity The spatio-temporal distribution of male and female P. perniciosus was similar (P > 0.05), and Fig. 2 shows the seasonal abundance for both sexes in CDC traps, together with the mean indoor-recorded temperature and RH, on the day when traps were collected. Abundance for this and other major species was highest in July and peaked sharply in all zones, during the second week of this month, and few sand flies were collected in May (week 1) and October (weeks 7 and 8) (Fig. 2, Table 2). The peak in P. perniciosus abundance in the week 2 of July coincided with the lowest recorded mean RH and highest mean temperature in the study (Fig. 2). This was associated with a similarly drastic change in the weather regionally, particularly the days when W, S and N were sampled; the mean maximum RH and temperature and mean wind speed in these zones were 89%, 33°C and 2.1 m/s on July 3rd and 53%, 40°C and 1.4 m/s on July 7th (Fig. 3).  The strong, negative association between P. perniciosus abundance and RH is further reflected in Table 3 showing the relationship between these two variables in the overall study data set. Moreover, altitude was negatively correlated with temperature (r = -0.21) and RH (r = -0.37) (P < 0.05).
Bivariate relationship between P. perniciosus in CDC traps, trap positioning, building characteristics and animal species and husbandry The proportion of positive traps and sand fly abundance in positive traps was greatest in older, small and not frequently disinfected buildings, with unplastered stone   or brick walls and traditional ceilings made of cane and plaster, wood or bricks (Table 4). Sand fly abundance was also numerically greater in poorly ventilated buildings with straw or soil bedding. Neither the proportion of positive traps or sand fly abundance in positive traps was associated with animal species (sheep or dogs), average number of animals or animal density in the building or the use of insecticidal treatments on the animals (Table 4). Although sand fly abundance was associated with trap distance to the wall, the relationship did not follow a density positive or negative trend (Table 4).
Bivariate relationship between P. perniciosus abundance in CDC traps and external temperature, relative humidity, rainfall and wind speed A summary table of climatic variables recorded at weather stations, showing a negative or positive association with the proportion of P. perniciosus positive traps and/or its abundance in positive traps is presented as supplementary material (Additional file 1: Table S2). Both the proportion of positive traps and abundance in positive traps were consistently negatively associated to May-October mean RH and maximum wind speed (Additional file 1: Table S2). Similar consistently negative associations were observed between the proportion of positive traps and the maximum annual and maximum and mean November-April wind speed and maximum annual rainfall, and between abundance in positive traps and May-October absolute maximum temperature. Other climatic variables were negatively associated with either the proportion of P. perniciosus positive traps or abundance, but the relationship was less consistent (Additional file 1: Table S2). In contrast, maximum RH in November-April was positively associated with sand fly abundance (Additional file 1: Table S2).
Bivariate relationship between P. perniciosus abundance in CDC traps and land use, soil and ground types After excluding significant associations between the percentage of P. perniciosus positive traps/abundance and land uses and soil and ground types present in comparatively small amounts (for example land used as urban fabric and coniferous forests), and those in which neither a consistent positive or negative trend with sand fly abundance was observed, results may be summarized as follows. The percentage of positive traps and abundance was greater in areas with moderate or large amounts of non-irrigated arable land, sparse vegetation and sandy grounds compared to areas with little or no amounts of these land types (Additional file 1: Table S3). The proportion of positive traps was positively associated with fluvisols grounds and abundance was negatively associated with coluvial soils and positively associated with petrocalcic xerosols (Additional file 1: Table S3).

Multivariable relationship between P. perniciosus sand fly abundance and environmental variables
Phlebotomus perniciosus count data adequately fitted a negative binomial distribution (P > 0.05). The most parsimonious Type I model included a variable combining indoor temperature and RH and external mean maximum wind speed between May and October, and building age ( Table 5). Incidence rate ratios (IRR) were greatest for lowest RH and highest temperature and decreased with increasing RH, reaching their lowest value for RH > 60% and lowest temperature < 22°C. Moreover, IRR increased with decreasing wind speed and increasing building age   (Table 5). It was not possible to include random effects in this model as this led to model convergence failure. Among Type II models, the one with the lowest AIC indicated that sand fly abundance was negatively associated with precipitation, maximum temperature and maximum wind speed in May to October. The model revealed wide variation between premises in the sand fly count baseline (intercept) although this variability remained constant during the study (slope) ( Table 6).

Discussion
In a recent study investigating the presence/absence of P. perniciosus in southern Spain, the probability of finding sand flies increased with altitude up to 769-1,153 m, reflecting the positive association between the sand fly presence and temperature in this altitude range [26]. Sand fly abundance in the present study was similarly lowest in coastal areas and highest in the 844-849 m altitude range (W zone). However, it was significantly lower at 536-794 m (N zone), indicating that altitude or temperature alone, are inadequate predictors of sand fly abundance. RH, which was strongly, negatively associated with sand fly abundance, was similarly low in the N, W and S (265-352 m) zones but sand fly counts were much greater in W and S than in the N. Models greatly improved when maximum wind speed was fitted because wind exposure was highest and most variable in the N and SE (44-83 m) zones. Sand flies are poor fliers [15], and the wind may prevent them from entering buildings and probably generates drafts inside animal buildings, discouraging adult sand fly activity there.
Climate was responsible for the observed seasonality and the marked fluctuations in sand fly abundance over a short period. The huge increase in the second week of July coincided with a "heat wave" characterised by a sharp increase in temperature and a drop-in RH and wind speed. As far as we are aware, there are no previous reports of similar increases in sand fly abundance following heat waves typical of Mediterranean summers. Notwithstanding, Branco et al. [11] reported highest sand fly density in central Portugal associated to highest average monthly temperature, lowest RH and absence of strong wind. A similar relationship between temperature and RH and the abundance of the sand fly Lutzomyia shannoni was reported in Florida in the USA [27]. Rainfall and RH were also negatively associated with sand fly activity in other Mediterranean regions [10,28,29]. Although high RH is required by sand fly instars to develop and adult sand flies are very sensitive to desiccation [30], low RH favours adult activity, possibly during short spells in search of food.
Animal building age was also strongly associated with sand fly abundance. It indirectly accounted for several factors that impact on sand fly survival. Sand flies are very sensitive to disinfectants and insecticides, but they were not frequently used except in the most modern dog premises. Besides, old buildings with stone walls and accumulated organic matter are considered ideal for sand flies to breed and rest. They were also poorly ventilated, and carbon dioxide (CO 2 ) is a strong attractant for blood-searching females [16].
The role of land use, soil and ground types on sand fly abundance remains unclear. Many such variables were associated with sand fly abundance in the univariate analysis, but in most cases, there was no evidence of a consistently positive or negative trend. Exceptions were increasing sparsely vegetated and non-irrigated arable land and petrocalcic xerosol ground and decreasing coluvial soils associated with greater sand fly counts. However, none of these variables was retained in the final multivariable model. This may not be surprising given the strong correlation between environmental variables. The wide variety of environments in which P. perniciosus can thrive suggests that its density on a large geographical scale depends more on climatic conditions than on specific terrains and land uses. This conclusion, however, may not be extended to other regions and species [31]. Moreover, multilevel models revealed considerable unexplained variation between study premises in the same zone, so clearly, microhabitat factors not accounted for in this study can have a profound effect on sand fly density.
The strong correlation between outdoor and indoor climatic variables allowed using the former to model sand fly abundance and may be used to generate and validate sand fly abundance density maps, and identify areas that require further studies of vector and pathogen distribution. In previous epidemiological studies on canine and human leishmaniosis in Murcia Region, seroprevalence in dogs was highest in the S, lowest in the N and variable in the SE [18]. Similarly, human PCR prevalence was highest in the S zone, lowest in the N zone and C and variable in the SE zone [17]. Murcia coastal SE is climatically variable, and this could be associated with a higher sand fly and leishmaniosis spatial overdispersion. Leishmaniosis foci associated to P. perniciosus have been reported in coastal areas in Italy [32]. Further entomological and epidemiological studies are needed to in Murcia's SE zone, as well as in the C and W zones where information on Leishmania prevalence is presently incomplete.
The study focused mainly on CDC light trap captures after observing that sand fly distributions in sticky traps were similar but had comparatively few sand flies. P. perniciosus was the most abundant species in both trap types. Light traps are particularly suited for sand flies with strong phototropism such as P. perniciosus females [32]. In contrast, sticky traps sample sand flies by interception providing an unbiased estimate of insect activity in a place [15]. Sergentomyia minuta feeds on reptiles and are not strongly phototropic and was the dominant species in most studies in Spain using sticky traps (reviewed by Galvez et al. [10]). However, the number species identified in the present study was the same and their relative abundance in light traps similar, to that reported in south-east Spain 30 years ago [33][34][35][36]. According to the later author, less common sand fly species have narrower preferential bioclimatic conditions and among them, P. papatasi, P. sergenti and P. alexandri favour arid zones [35]. While no definite conclusions can be drawn from the present study in this respect, P. papatasi was relatively more abundant in the most arid C zone, but P. sergenti was most common in the least arid W zone.
In summary, this study confirms the presence of sand flies in the Murcia Region including the two main L. infantum vectors, P. perniciosus and P. ariasi, and provides a quantitative analysis of their spatial distribution in relation to environmental variables. Sand fly abundance is heterogeneously distributed, strongly depending on temperature, RH, rainfall, wind speed and microenvironmental factors. These findings may be extrapolated to other Mediterranean regions to improve our understanding of P. perniciosus and L. infantum infection dynamics. Moreover, the methods used in this study may be a model to perform standardised and optimized abundance studies on sand flies.

Conclusions
Phlebotomus perniciosus is the predominant sand fly species in the countryside in the Murcia Region, and its abundance is spatially and temporally heterogeneous. Climate, including relative humidity, temperature and wind speed and not land use or soil characteristics, have the greatest impact on sand fly density on a large geographical scale. Microenvironmental factors such as animal building characteristics and husbandry practices can significantly affect sand fly counts on a small geographical scale. These finding may be used to developing predictive vector and pathogen distribution maps.

Additional file
Additional file 1: Table S1. Percentage of "sticky" traps with sandflies (positive traps) and sand fly abundance in positive traps in dog kennels and sheep flocks in Murcia Region, southeast Spain in 2015. Table S2. Percentage of CDC traps with P. perniciosus and abundance in positive traps according to climatic variables from meteorological stations. A study of sand fly abundance in dog kennels and sheep flocks in Murcia Region, southeast Spain, in 2015. Table S3