Species diversity and spatial distribution of CL/VL vectors: evidence for variation in gene expression of immunogenic SP15 and LeIF proteins due to bioclimatic characteristics and physiological status of sand flies, Iran

Main approaches of this investigation were climate regionalization to recognize the spatial distribution of cutaneous/visceral leishmaniasis (CL/VL) vectors as risk-maps using ArcGIS modeling system, evaluation of species biodiversity, assessing bio-climate (BC) effect on expression plasticity of genes possessing vaccine properties isolated from wild-collected sand flies. The Inverse Distance Weighting (IDW) interpolation method was used to obtain accurate geography map using ArcGIS10.3.1 in closely-related distances. Species diversity was calculated based on Shannon diversity index using Rv.3.5.3. Expression fold change of SP15 and LeIF genes was evaluated using cDNA synthesis and RT-qPCR analysis.


Abstract
Background Main approaches of this investigation were climate regionalization to recognize the spatial distribution of cutaneous/visceral leishmaniasis (CL/VL) vectors as risk-maps using ArcGIS modeling system, evaluation of species biodiversity, assessing bio-climate (BC) effect on expression plasticity of genes possessing vaccine properties isolated from wild-collected sand flies.

Methods
The Inverse Distance Weighting (IDW) interpolation method was used to obtain accurate geography map using ArcGIS10.3.1 in closely-related distances. Species diversity was calculated based on Shannon diversity index using Rv.3.5.3. Expression fold change of SP15 and LeIF genes was evaluated using cDNA synthesis and RT-qPCR analysis.

Results
Three BC zone were identified in northeast of Iran. Phlebotomus papatasi were abundantly captured in all BC regions and the frequency was more in plains areas of mountainous BC as hot spots of CL. VL vectors were more prevalent in spatial cluster of Mediterranean BC. Semi-arid BC was identified as a major contributing factor to up-regulate SP15 salivary gene expression ( P =0.0050, P <0.05), and Mediterranean BC had considerable effect on up-regulation of LeIF-Leishmania gene in gravid and semi-gravid P. papatasi population ( P =0.0109, P <0.05).

Conclusions
The diversity and spatial distribution of CL/VL vectors associated with BC regionalization obtained in our research provide epidemiological risk maps and establish more effectively control measures against leishmaniasis. Oscillations in gene expression indicate that each gene has its own features, which are profoundly affected by bioclimatic characteristics and physiological status of sand flies. It is essential to consider BC factors affecting regulatory regions of environmentally responsive loci for genes used in vaccine designing.

Background
Leishmaniasis is an epidemic-prone infectious and neglected tropical disease (NTD) caused by obligate intra-macrophage protozoa of the genus Leishmania. The Leishmania parasites are transmitted by the bite of infected female sand flies to humans, which resulting in an extensive range of clinical manifestations mainly from self-healing cutaneous leishmaniasis (CL, the most common) to progressive fatal visceral leishmaniasis (VL) [1,2]. Geographical spread of the disease is extended throughout widespread territories and countries (98 countries) in all over the world, except the Australia, the Pacific Islands, and Antarctica [3]. More than 95% of new CL cases occur in six countries including Afghanistan, Algeria, Brazil, Colombia, Iran, Iraq and the Syrian Arab Republic [2]. Iran is one of the most important foci of leishmaniasis in all around the world, which is one of the major public health problems in more than half of the provinces [4]. Iran's geographic location provides an appropriate climate and landscape habitat for different species of rodents and dogs as main reservoirs and sand flies as principle vectors [5].
In northeast of Iran, CL and VL were reported from all counties and five out of eight counties of Northern Khorasan Province, respectively [6]. It is essential to investigate the ecology, biodiversity, adaptation to habitats, and distribution of sand flies, according to environmental variables involved in transmission, therefore control measures and spatial modeling of infectious diseases can be obtained using the Geographic Information System (GIS)-based survey [5,7,8]. Biodiversity includes intraspecific (genetic polymorphism), interspecific, and ecosystem diversity. Evaluation of biodiversity measurements can be achieved based on spatial scale [9] as followings: alpha diversity (α) (the number of species within a community like a natural habitat), beta diversity (β) (interspecific differences with which we can rapidly compare the change rate of diversity at different habitats), and gamma (γ) or regional diversity (diversity among ecosystems). Calculating the diversity indices and comparing their results is a useful method for recognizing the ecological status of sand flies. Hence, determining biodiversity indices of sand flies in different topographic locations and the study of correlation between climatic elements and incidence of leishmaniasis can helps us developing an effective prevention approach against vector-borne diseases [10,11,12]. Despite the worldwide distribution of leishmaniasis and also its significant disease burden, currently there are no commercially available human leishmaniasis vaccines delivering the necessary level of protection.
Due to the lack of efficacious vaccines to treat leishmaniasis, studying the vectors is an important component of the effort to control leishmaniasis [13,14]. Prior to the use of antigenic proteins specific to the genus Leishmania or exploited from sand fly saliva as a vaccine target, the genetic and expression changes of salivary and Leishmania proteins firstly needs to be thoroughly assessed in natural-field habitats.
In our recent ongoing research, bioinformatics analyses demonstrated that the vaccine comprising the combination of truncated forms of salivary protein from P. papatasi (SP15) and Leishmania eukaryotic initiation factor (LeIF) namely SaLeish, has an appropriate traits to increase both innate and specific cellular immune responses with acceptable population coverage in different endemic areas of the world [15]. The protein PpSP15 was the first identified salivary protein, and could be used in a sand fly salivary protein-based vaccine strategy as saliva-driven protective immunity against vector-transmitted leishmaniasis [16,17]. Indeed, the immunization obtained by SP15 from P.
papatasi and its orthologs in P. duboscqi (PdSP15) was sufficient for protection through specific delayed-type hypersensitivity (DTH) response with a Th1 profile in Non Human Primates (NHP) and humans [18,19]. It has been indicated that SP15 of P. papatasi salivary has the ability of stimulating the dendritic cells (DC), the potent antigen presenting cells, and therefore initiating the salivamediated immunity conferring the induction of rapid sand flies saliva-specific Th1 immune response [16,20]. Another protein selected to be examined was LeIF, which was obtained from natural collected-sand fly of Northern Khorasan province. LeIF is homologous and highly conserved among Leishmania spp., which can be used as a natural adjuvant for antigen-based vaccine strategy and results in inducing the early production of LeIF-specific Th1 cytokines, IL-12 and IFNγ in PBMC of both leishmaniasis patients and normal individuals [21,22]. The LeIF protein has previously been used as leishmaniasis prophylactic vaccine (Leish-111f) that elicited an increase in CD4 + T cells and thus indicating a predominant Th1-type immune response [23]. There are few studies considering the ecoepidemiological effects according to bio-climate (BC) regionalization on the expression of sand fly salivary proteins as well as their Leishmania within at peak activity of sand flies in leishmaniasis endemic regions. Besides, a number of studies have merely investigated the effect of environmental response on the evolution of gene expression, regardless of the precise determination of the locations affecting the expression of genes of target associated with GIS modeling system [24].
The aims of this investigation were 1) determining the spatial distribution of potential and principle vectors of CL and VL based on BC regions in ArcGIS tool as hazardous maps in northeast of Iran, 2) to regionalize the ecological habitats of sand flies as bio-climate zones according to the environmental factors using integrative approaches, species diversity assessment and richness, 3) similarity analyzing in sand fly populations collected in different ecotopes based on Shannon-Wiener index and Jaccard coefficient using R, 4) evaluating different gene expression of immunogenic proteins (SP15 and LeIF) in wild-caught P. papatasi salivary gland and from natural Leishmania-sand fly in relation to the effects of bio-climate factors and physiological statuses of sand flies appropriate for protein-based vaccine design.

Specifications of study area
Integrative methods were used to evaluate eco-epidemiological parameters such as bioclimatic change, landscape fragmentation, and ecological niche considering vector domiciliation and their distribution in different locations of Northern Khorasan foci. The Northern Khorasan province is located in northeast of Iran and from the epipaleolithic period onwards was substantially inhabited by nomadic hunter-gathers. Geographic location of studied area (GPS coordinates) is between latitudes of 36˚34 -38˚17 Ń and longitudes of 55˚52 -58˚20 É, and has an area of about 28,434 km 2 including eight counties: Bojnord (the capital city of province), Shirvan, Esfarayen, Maneh-o Samalghan, Jajarm, Farouj, Garmeh, Raz-o Jargalan (Fig. 1). This province is bordered Turkmenistan from north, with Khorasan-Razavi province from south and east, with Golestan province from west, and with Semnan province from the southwest (Fig. 1). Northern Khorasan district is generally regionalized into two fractions as: Mountainous areas, and flat regions with low elevation. The Kope-Dagh to the north and Aladagh, extending of the Alborz Mountains to the south, are two main ranges with fertile plains between the mountains, which determine the favorable agricultural conditions and pastoralism.
Regarding various weather conditions, this area has a thriving culture and vegetation. Regions of Bojnord, Shirvan, Faruj, and Esfarayen are cold mountainous whereas western districts of Maneh-o Samalghan and Raz-o Jargalan counties are temperate and forested. Parts of lush green forest of Golestan National Park (one of the oldest national parks in the world) are located in Northern Khorasan Province. However, this province has also desert areas such as Jajarm, where there is a wide range of temperature oscillations (cold winters, and relatively hot and dry weather in summer) (Additional file 1: Figure S1a, b, c). Furthermore, significant rainfall amount occurs at different intensities in the province. Due to vast mountainous areas, the Northern Khorasan zone is mainly known to have a cold to temperate mountainous region. The provincial climate is affected by longitude and latitude, and altitude. Also, three basic types of air masses bring variety of weather conditions to the Northern Khorasan area, including (1) Western humid air masses that enters the province in early autumn and lasts until late spring, and the peak of its activity during winter season is the main cause of rainfall in the province (Additional file 1: Figure S1b), (2) Siberian cold air masses, which affects the province from early autumn to early spring causing a significant temperature reduction (Additional file 1: Figure S1c), and (3) a warm dry air mass, which enters to the province from the south during summer, and increased the aridity and temperature (Additional file 1: Figure S1d).

Model variables and climatic regionalization
The De Martonne aridity index which is used to determine the regional climate based on annual temperature and rainfall is calculated as follows: I = P/T + 10 where (I) stands for De Martonne aridity index, P stands for annual rainfall in millimeters, and T stands for the annual temperature in ˚C ( Fig. 1, Additional file 1: Figure S1). Also, climatic regionalization was carried out as earlier described [5] by the use of Principal Component Analysis (PCA) as multivariate statistical method to minimize sources within the variation of the group. Then, Clustering Integration Method (CIM) was applied to precisely improve the classification accuracy. The southern cities have higher temperatures compared to northern regions, due to their different altitudes from sea level, entrance of various air masses, and different longitude and latitudes (Additional file 1: Figure S1). The rainfall levels were also varied in different regions of the province. The maximum and minimum rainfall is observed in spring (March) and summer (June), respectively. Due to highly variability of sand fly species, even in closely distance from its sampled location, we utilized the Inverse Distance Weighting (IDW) deterministic interpolation using ArcGIS 10.3.1 to specify the collected sample points in terms of sand fly dispersion and targeted protein expression rates (SP15 and LeIF) in local scale levels. IDW is a flexible spatial interpolation method, which estimates unknown values along with specifying the search distance, closest points, power setting and, barriers [25]. To obtain the geography map in GIS, ordinary IDW parameters related to the accuracy of the values at each point was calculated with respect to distance: where Z0 is the estimation value of variable z in point I, zi is the sample value in point I, di is the distance of sample point to estimated point, N is the coefficient that determines weigh based on a distance, n is the total number of predictions for each validation case. The Jaccardʼs similarity coefficient (JSC) is the simplest data, which compares the similarity for two sets of data between different locations of communities. This index always gives a value with a range from the lowest (0%) to the highest similarity (100%) between two populations [26].

Site selection
This cross-sectional study was performed at the peak of sand fly activity from mid-June to the late of    Collecting the Sand flies Considering geographical locations (altitude, latitude, and longitude) and weather conditions of the counties ( Fig. 1, Additional file 1: Figure S1a, b, c, d) as well as during the peak seasons of sand flies activity (according to the incidence of leishmaniasis in each of these regions), trapping was carried out using 120 sticky oil papers, four CDC light traps and, 70 funnel traps. The catching time was from five to 10 days before sunset until sunrise at each location. Dead sandflies were removed from castor oil papers using fine brush and, were put into sterile tubes containing 96% alcohol, and then in order to identify the species, they were kept at -20 ˚C for future morphological and/or molecular experiments. Blood-fed, semi-gravid, and gravid female sand flies were collected alive by the use of aspirator, funnel traps and CDC light traps, and after that were immediately moved into suspended net cages using mouth aspirator. To provide paramount protection of captured sand flies, a piece of 30% sugar-soaked cotton pad was placed on the top panel of holding cages, and were then covered with plastic bags containing a moist sponge to avoid from losing humidity and/or excessive heat.
Field-caught sand flies were first euthanized in water, dissected, and were then identified to species based on morphological characteristics of the head and terminal genitalia under compound microscope (400×) using a systematic key [27,28].

RNA extraction from salivary gland and Leishmania in wild-caught sand flies
Expression pattern of SP15 and LeIF proteins isolated from salivary glands and Leishmania protozoans of P. papatasi were incriminated in different locations of studied areas during the peak of sand fly activity from June to mid-September. Physiological status of sand flies was considered to determine SP15 gene expression level. Therefore, due to fresh-fed sand flies indicated a significant up-regulation of SP15 compared to un-fed ones [29], we selected only blood-fed sand flies to investigate the pattern of expression for SP15 transcript in different eco-epidemiologic conditions of Northern Khorasan province. After accomplishing definitive identification of female P. papatasi, salivary glands in addition to their head accessories were dissected using sterile forceps, and then were immediately re-   LF1153R, 5′-TCGATCTGCGTGTGGTAGTG-3′) with a 159 bp fragment for the LeIF cDNA. α-tubulin gene was used as a housekeeping load control mRNA to the qPCR reactions primers Tub-P24F and Tub-P24R in P. papatasi sand fly [30]. All PCR experiments were followed by melting curves to check the amplicon specificity.
Different expression ratios for salivary protein and Leishmania genes (SP15 and LeIF) were calculated and displayed as fold changes over a control, using the 2 −ΔΔCT method [31]. The fold changes were calculated by the expression 2 −ΔΔCT , where ΔΔC T = ΔC T (sample) -ΔC T (calibrator), ΔC T = ΔC T (sample) -ΔC T (alpha tubulin gene), C T = cycle at the statistically significant increase in the emission intensity over the background. The calibrator was represented by the average expression (mean ΔC T ) of the newly emergent un-fed nulliparous laboratory-reared female sand flies [30]. Moreover, negative controls were considered without the templates. Fold changes were calculated for each sample in comparison with the calibrator.

Data analysis
Statistical analyses were carried out using the software GraphPad Prism v. 5 and computed a single P value from the sum of these discrepancies. Species biodiversity used in this investigation is affected by two factors as: (1) number of species and (2) their abundance of distribution [26]. Species diversity was calculated based on Shannon-Wiener index, species richness and also evenness index using the Ade4 package of R v.3.5.3 software to estimate the biodiversity of sand fly species in the study areas. Species diversity is a measure of diversity within an ecological community, which encompasses both species richness and evenness of species abundance. Species richness is the number of species in a specific area or environment, whereas, species evenness refers to how close each species are in number in the similar location. Evenness is species abundance distribution and the more even this distribution is the greater the species evenness will be [32].
Computation of these indices is shown in Table 3. Also, estimation by similarity was calculated between different communities using Jaccardʼs similarity coefficient within different localities collection (Table 4).    Figure S2).
We analyzed the expression of SP15, protein of salivary glands, from blood-fed of uninfected female in 10 individual P. papatasi sand flies from each location at each period of the season (early, mid, or late) (Figs. 6, 7, 8, Table 2, Additional file 4: Table S2), and also the LeIF expression was quantified from Leishmania major isolated from infected female of P. papatasi sand fly with fresh blood-fed, gravid, and semi gravid status (Figs. 9, 10, 11, Table 2, Additional file 5 Table S3). Four principal climatic elements of studied region (wind, annual precipitation, temperature, and humidity) were considered to determine bioclimatic classification and distribution of sand fly species of Northern-Khorasan province (Figs. 1, Additional file 1: Figure S1). Among aforementioned factors, annual temperature, humidity and rainfall have the greatest effects on the sand flies lives. After evaluating the impact of principle climatic factors on climatic variables, three regions were classified as homogeneous climate zones (Table 1, Additional file 2: Table S1 (stations sheet), Additional file 3: Figure S2). The regions are as follows: (1) The Mediterranean with spring rains area in northwestern including Razo-Jargalan, and Maneh-o Samalqan (Table 1, Additional file 1: Figure S1, Additional file 2: Table S1). (2) The cold-mountainous area in southwestern including northern parts of Garmeh, in central including Bojnord, and in eastern including Shirvan, Faruj, and Esfarayen. (3) The low rain fall and warm-semi dry area including Jajarm, southern and western parts of Garmeh and Esfarayen, respectively (Additional file 1: Figure S1).  Table S1). In the genus Phlebotomus, P. papatasi were in all climate regions (n = 800 specimens: 68.55%). Regarding, P. papatasi was the most abundant one in spatial cluster of mountainous regions particularly in Baba Aman, Quch Qaleh, and Paygho where human dwellers are almost in all ecotopes, ranging from 11.25-18% (Fig. 2, Table 1), and the least abundant one of them was observed in Gholaman and Darkesh of Mediterranean region, and the farthest eastern part of mountainous regions (Fig. 2). P. alexandri, proven vector of VL and suspected vector of CL, was the second-most widespread species caught in mountainous climate zone (Fig. 3, Table 1, Additional file 2: Table S1 (xy sheet)). From 153 P. alexandri, 60.79% (n = 93 specimens) were caught from rodent burrows' ecotope in sylvatic habitats, while 39.21% P. alexandri (n = 60 specimens) were found in peridomestic animal shelters. Paraphlebotomus species were found only in spatial cluster of mountainous and Mediterranean bio-climate (BC) zones (Fig. 3, Table 1). Sergentomyia genus was prevalent almost in all three bio-climate areas (Fig. 4, Table 1). However, Sergentomyia was more prevalent only in spatial cluster of Mediterranean and mountainous BC zones especially more for Gholaman, Qareh Bashloo and, Quch Qalleh-Olya with the mean of range between 10.26-16.33% (Fig. 4 Table   S4).

SP15 and LeIF expressions
Quantitative data for the expression profiles of SP15 and LeIF of salivary and Leishmania genes were summarized in Table 2, Additional file 7: Table S5, in terms of the physiological status (blood-fed, gravid and semi-gravid) and bio-climate classification (Mediterranean, mountainous and semi-arid) for female P. papatasi sand flies. Fold change ratios were ranged from the least as 0.12 in June-early to the highest as 24.13 in September-late for SP15 in three Bio-Climate zones (BCz) (Figs. 6,7,8).
Expression fold of SP15 was ranging from 0.12 in Mediterranean BCz (Gholaman) to 6.41 in semi-arid BCz (Garmeh-Houmeh) in June-early (Fig. 6). One way ANOVA (Kruskal-Wallis test, KWt = 8.548) of SP15 fold expression in June-early indicated that the medians were significant (P = 0.0139, P < 0.05) in all three BC zones (Fig. 6, Table 2, Additional file 7: Table S5). Although SP15 fold changes had a marginal increase from Mediterranean BCz to mountainous BCz in June-early, it had also a remarkable rising expression in semi-arid BCz at three times of collection (Fig. 12a, b, c). Mann Whitney U test (MWt) was performed to compare SP15 expression fold change ratios between each two BCz in early (June), mid (August), and late (September). Regarding, the obtained data were analyzed and tabulated in Table 2, Additional file 7: Table S5. Expression fold ratios of SP15 were ranged from 0.31 in Mediterranean BCz (Darkesh) to 7.9 in semi-arid BCz (holy tomb of Emamzade-Ebrahim, Sankhast County) in August-mid (Fig. 7, Additional file 4: Table S2). The Kruskal-Wallis test was KWt = 11.19 for SP15 fold expression in Mid (August) and, the median was significant (P = 0.0037, P < 0.05) in three BC zones (Fig. 12b, Table 2, Additional file 7: Table S5). The expression profiles of SP15 salivary gene were ranged from 0.48 in mountainous BCz (Malkesh) to 24.13 in semi-desert BCz (Kouran district) in September-late (Fig. 8, Additional file 4: Table S2). The Kruskal-Wallis test was KWt = 10.60 for SP15 fold changes in late (September), and the median was significant (P = 0.0050, P < 0.05) in all BC zones (Fig. 12c, Table 2, Additional file 7: Table S5). The expression profiles of salivary SP15 gene in three regionalized bio-climate locations are as followings: Semi-desert > Mountainous > Mediterranean ( Fig. 12a, b, c, Table 2, Additional file 7: Table S5).
The effect of physiological status of blood meal (blood-fed, gravid, and semi-gravid) on expression of  Table 2).
Fold changes ratios for LeIF expression was ranged from the least as 0.03 for blood-fed status in semidesert BC zone (Sankhast) (Figs. 9, 12d, Additional file 5: Table S3) to 4.82 as the highest expression fold for gravid status in Mediterranean BC zone (Gholaman) (Figs. 11, 12f, Additional file 5: Table S3).
The Kruskal-Wallis test was KWt = 9.046 for LeIF fold expression in gravid status of P. papatasi, the median was significant (P = 0.0109, P < 0.05) in three BC zones (Figs. 11, 12f, Table S5). The expression profiles of LeIF from natural Leishmania-field-caught P. papatasi sand flies were as followings: gravid of A, B and, C > semi-gravid of A, B and, C > blood-fed of A and C in three regionalized BC zones (Fig. 12d, e, f, Table 2, Additional file 7: Table S5). In line with recent reports in Iran [33], current investigation not only indicated higher diversity for specific species of sand flies in mountainous region only within Sergentomyia and Paraphlebotomus species, but also has confirmed actually more diversities and richness of vectors of visceral leishmaniasis only in semi-desert and Mediterranean BC communities (Table 3, Additional file 6: Table   S4). Vectors of visceral leishmaniasis indicated no considerable diversity in mountainous region except Dartum and Mahnan counties. Of note, in contrast with the recently published data [34], our findings did not show any similarity for the evenness indices in P. papatasi population in Northern Khorasan (Table 3). In fact, the evenness were ranged from E = 0.204 in mountainous BC zone for

Discussion
Paraphlebotomus species up to E = 0.721 in all three BC ecotopes for vectors of visceral leishmaniasis and Sergentomyia species (Fig. 4, Table 3, Additional file 6: Table S4). The lack of diversity among P.
papatasi population and visceral vectors can be caused by the extension of Alborz Mountain range (Aladagh Mountain range) as an important ecological barrier in front of sand fly distribution (Table 3).
P. alexandri were the most abundant species among Paraphlebotomus population, mainly collected from mountainous BC areas at high altitudes (Fig. 3, Table 1, Additional file 1: Figure S1b, Additional file 2: Table S1 (xy sheet)). In fact, in agreement with previous findings in southwest of Iran [5], P.
alexandri were caught from sylvatic habitats (Qare-Bashloo, Dartum, Gerivan and, Mahnan) near rodents' ecotope and animal shelters in the vicinity of human dwellers (Fig. 3, Additional file 2: Table   S1 (xy sheet)). Also, P. sergenti was found mostly in rural areas (Paygho, Kohene-Oqaz, Zartanlu) of mountainous BC zone (Fig. 3, Additional file 2: Table S1 (xy sheet)). Our findings were consistent with previous results in Iran [5] which indicated that population of these species tend to be colonized more in places with human inhabitants, and their peridomestic animals even in remote mountainous areas, due to anthropophilic behavior of P. alexandri and P. sergenti, (Fig. 3, Table 2, Additional file 7: Table   S5). One of the important reasons for the higher frequency of visceral leishmaniasis vectors in districts of Gholaman and Darkesh compared to other mountainous BC regions can be caused by the existence and development of dams along with the boundary lines (Fig. 5, Additional file 2: Table S1 (xy sheet)). It is notable that each bio-climate regionalized has its environmental effects on the regulation of SP15 salivary proteins of sand flies associated with gradually increasing levels of SP15 gene expression fold reaching the highest levels of expression late in the season in semi-desert area ( Fig. 12a, b, c), when the environment is drier and sugar sources are scarce (Figs. 6, 7, 8, Table 2, Additional file 7: Table S5) [30]. Also, in addition to blood-feeding of P. papatasi on mammalian hosts, feeding with various plants species as a sugar-source diet has various either favorable or adverse effects on life cycle [35,36] and gonotrophic cycles [37,38] Table S5). Interestingly, our geospatial of BC analyses indicate that up regulation of SP15 gene expression in semi-arid area confirmed the conceivable correlation of SP15 salivary protein expression with bio-climate variation considering vegetation cover for availability of sugar sources as principle factor [39]. In fact, the expression levels of P. papatasi salivary gland genes (SP15) in flies from a dryer habitat (Kouran) (Fig. 8, Additional file 4: Table S2) are higher than the expression levels in those flies from an irrigated area (Darkesh) late in the summer season, when drought may affect the sugar content of plants [38]. Accordingly, correlation between the highest expression of SP44 gene of salivary protein of P. papatasi was also confirmed with arid-environment and the rarity of sugar sources in September late [30]. Moreover, aging has been proved as another physiological factor affecting the gene expression profiles; however, it is not considered as significant as genotype or sex determination in insects [40]. Also, similar to what has been observed by other researchers [29], the direct relationship was corroborated between the up-regulation of SP15  Table 2, Additional file 7: Table S5). Nevertheless, testing Leishmania gene expression (LeIF) in gravid and semi-gravid flies showed also similar results in Mediterranean and mountainous BC regions (Fig. 12e, f, Table 2, Additional file 7: Table S5). However, along with other factors analyzed in this study, gene expression fluctuations for LeIF fold change from June (early) to September (late) can be also caused by the existence of some specific plants in three regionalized BC zones (seim-arid, Mediterranean and, mountainous) with the different quality of sugar meals which is attractive for sugar-feeding behavior of P. papatasi [38]. Also, expression of genes involved in digestion of blood-meal or sugar meal in blood-fed, semi-gravid, and gravid of sand fly are regulated according to the type of nutrient acquisition [41] which play important roles in Leishmania gene expression associated with the distinct environmental BC zones in natural habitats of principle vectors. In contrast with the expression change of salivary protein SP15, the expression fold of natural Leishmania-field-caught sand fly gene (LeIF) was the lowest in the semi-desert BC region in all three physiological statuses of P. papatasi sand flies (gravid, semi-gravid, and, blood-fed) (Fig. 12d, e, f).
Level of LeIF expression was significantly higher in gravid and semi-gravid sand flies compared to blood-fed sand flies. In fact, in line with recently published data [42]; our research indicated less expression fold of LeIF in fresh-blood fed sandflies, and this can be resulted from toxic products of fresh blood-meal digestion. Also, we assumed that according to pruzinova et al. data, as midgut binding is a key factor of Leishmania parasites' survival in P. papatasi, therefore more availability of sugar sources in irrigated areas (Mediterranean BC zone) can be responsible for higher expression of Leishmania promastigote genes and consequently promastigote survival in association with increasing of the salivary protein content [43] to exacerbate the disease in natural geographic habitats.
In the future, it is suggested that the gene regulatory regions of salivary proteins possessing vaccine specificity should be identified against leishmaniasis. Consistent with the results of this study, natural environmental effects which are effectively triggering a significant increase in gene expression fold of specific proteins should be considered with the effects of arid/semi-arid habitat especially for salivary gland protein (SP15). Undoubtedly, the relevance of molecular biology to the mimicry of natural effect of living environment on disease agents is an underlying assumption to assist indoors lab studies using potential applications of gene drive strategy for principle vectors to serve a desirable function of gene at high enough levels, and consequently to control leishmaniasis or vector-borne diseases.

Supplementary Information
Additional file 1: Figure S1. Northern Khorasan province was incriminated based on significant factors a) Height, b) Annual rainfall, c) Temperature, d) Humidity using deterministic spatial interpolation method in GIS.
Additional file 2: Table S1. Distribution of various sand fly species based on species diversity collected from three climatic zones of Northern Khorasan province, north east of Iran (2017-2018).
Additional file 3: Figure S2. Bioclimatic regionalization of Northern Khorasan province generated by the IDW method in ArcGIS.

Ethics approval and consent to participate
Not applicable.

Consent to publish
Not applicable.

Availability of data and materials
All essential data generated or analyzed during this study are presented in this published article and supplementary section. The raw datasets are available from the corresponding author upon reasonable request.

Competing Interests
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding
Funding was available through the studentship granted from Pasteur Institute of Iran to Ali Bordbar, grateful to Ms. Sahar Ebrahimi for helping with the field work, and her assistance in sand fly dissection, morphological identifications, and laboratory experiments. We are gratefully thanked to Mr. Gholamreza Ebrahimi for helping with the field work that generously provided accommodation and transfer facilities. We appreciate indigenous people in rural areas of Northern Khorasan province that sincerely allowed us to collect sand flies from interior properties, stables and from their agricultural farms.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.