Expected Changes of Montenegrin Climate , Impact on the Establishment and Spread of the Asian Tiger Mosquito ( Aedes albopictus ) , and Validation of the Model and Model-Based Field Sampling

Aedes albopictus has become established in many parts of Europe since its introduction at the end of the 20th century. It can vector a range of arboviruses, of which Chikungunya and Dengue are most significant for Europe. An analysis of the expected climate change and the related shift in Köppen zones for Montenegro and impact on the establishment of Ae. albopictus was conducted. Outputs of a mechanistic Aedes albopictus model were validated by 2245 presence/absence records collected from 237 different sites between 2001 and 2014. Finally, model-based sampling was designed and performed at 48 sites in 2015, in a previously unexplored northern part of Montenegro, and results were validated. The Eta Belgrade University (EBU)-Princeton Ocean Model (POM) regional climate model was used with the A2 emissions scenario for the 2001–2030 and 2071–2100 integration periods. The results point to a significant increase in suitability for the mosquito and a vertical shift to higher altitudes by the end of the century. The model showed excellent results with the area under the receiver operating characteristic curve (AUC) of 0.94. This study provides a tool for prioritizing surveillance efforts (model-based surveillance), especially when resources are limited. This is the first published analysis of Climate Change that incorporates observations from the national synoptic grid and the subsequent impact on Ae. albopictus in Montenegro.


Introduction
Aedes (Stegomyia) albopictus (Skuse), also known as the Asian tiger mosquito, is a highly invasive species native to the tropical and subtropical regions of Southeast Asia [1].It is classified as one of the world's top hundred worst invasive species [2].It has a vector capacity for a wide range of vector-borne diseases (VBD) such as Dengue, Chikungunya, Yellow fever, and Zika [3][4][5].Vector-borne diseases specific to the Aedes genus are an increasing burden worldwide and expect to become an even greater problem by the end of the century with the expected changes in climate which will allow the mosquito to spread to higher latitudes.
It was first introduced in Europe through international transport routes for used tires and 'lucky bamboo' [6].The drought-resistant eggs can travel for long periods of time and survive until they reach a suitable environment for development.The first record of the Asian tiger mosquito in Europe was in Albania 1979 [7].Since the later introduction and establishment in Italy [8], it has been reported in 31 European countries (Figure 1) including Albania, Austria, Belgium, Bosnia and Herzegovina, Bulgaria, Croatia, Czech Republic, France, Germany, Gibraltar, Greece, Georgia, Hungary, Italy, Malta, Monaco, Montenegro, the Netherlands, Northern Macedonia, Portugal, Romania, Russia, San Marino, Serbia, Slovakia, Slovenia, Spain, Switzerland, Turkey, United Kingdom, and the Vatican.
The pressure exerted by climate change (CC) could lead to a dramatic increase in the spread as well as a drastic increase in the invasion speed of Ae. albopictus with current limit occurring due to overwintering conditions [9][10][11][12].An important factor to take into consideration for the future spread of the vector is the existence of natural barriers such as mountain ranges as well as the importance of road transport through main traffic routes in the future spread of the vector [13,14].A depiction of this is the presence of Ae. albopictus in Switzerland, where the mosquito is firmly established in the Canton of Ticino, but only a few isolated populations within the country are present north of the Swiss Alps, even though the area is climatically suitable for the establishment and overwintering of the vector [15].Moreover, Ae. albopictus has already been reported in the south of Germany since 2015, with road traffic as the most likely introduction route [16][17][18][19][20].
Atmosphere 2018, 9, x FOR PEER REVIEW 2 of 18 borne diseases (VBD) such as Dengue, Chikungunya, Yellow fever, and Zika [3][4][5].Vector-borne diseases specific to the Aedes genus are an increasing burden worldwide and expect to become an even greater problem by the end of the century with the expected changes in climate which will allow the mosquito to spread to higher latitudes.It was first introduced in Europe through international transport routes for used tires and 'lucky bamboo' [6].The drought-resistant eggs can travel for long periods of time and survive until they reach a suitable environment for development.The first record of the Asian tiger mosquito in Europe was in Albania 1979 [7].Since the later introduction and establishment in Italy [8], it has been reported in 31 European countries (Figure 1) including Albania, Austria, Belgium, Bosnia and Herzegovina, Bulgaria, Croatia, Czech Republic, France, Germany, Gibraltar, Greece, Georgia, Hungary, Italy, Malta, Monaco, Montenegro, the Netherlands, Northern Macedonia, Portugal, Romania, Russia, San Marino, Serbia, Slovakia, Slovenia, Spain, Switzerland, Turkey, United Kingdom, and the Vatican.
The pressure exerted by climate change (CC) could lead to a dramatic increase in the spread as well as a drastic increase in the invasion speed of Ae. albopictus with current limit occurring due to overwintering conditions [9][10][11][12].An important factor to take into consideration for the future spread of the vector is the existence of natural barriers such as mountain ranges as well as the importance of road transport through main traffic routes in the future spread of the vector [13,14].A depiction of this is the presence of Ae. albopictus in Switzerland, where the mosquito is firmly established in the Canton of Ticino, but only a few isolated populations within the country are present north of the Swiss Alps, even though the area is climatically suitable for the establishment and overwintering of the vector [15].Moreover, Ae. albopictus has already been reported in the south of Germany since 2015, with road traffic as the most likely introduction route [16][17][18][19][20].  [21].Established (red)refers to an established population (evidence of reproduction and overwintering) of the species has been observed in at least one municipality within the administrative unit; Introduced (yellow)-the species has been detected (but without confirmed establishment) within the administrative unit; Absent (green)-field surveys or studies on mosquitoes were conducted but the species has not been detected within the administrative unit; No data (dark grey)-no sampling has been performed and no data on the species is available within the administrative unit; and Unknown (light grey)-it is unknown whether there are field studies on this species within the administrative unit.[21].Established (red)-refers to an established population (evidence of reproduction and overwintering) of the species has been observed in at least one municipality within the administrative unit; Introduced (yellow)-the species has been detected (but without confirmed establishment) within the administrative unit; Absent (green)-field surveys or studies on mosquitoes were conducted but the species has not been detected within the administrative unit; No data (dark grey)-no sampling has been performed and no data on the species is available within the administrative unit; and Unknown (light grey)-it is unknown whether there are field studies on this species within the administrative unit.
The establishment of Aedes albopictus is influenced by a wide array of abiotic factors such as mean annual temperature [9,22], total annual precipitation [23], mean January temperature [9,24,25], mean summer temperature (June-July-August), frequency of precipitation [3,26,27], relative humidity [27] and land use [28][29][30]; biotic factors like NDVI [28,30]; and various socio-economic factors such as human population density [6,29,30], night-time light [6,29], water-keeping habits, rate of urbanization, housing quality, etc. [31,32].Temperature and precipitation have been identified as the climatic elements that have the most significant effect on the suitability for the establishment and inter-annual activity of the vector.On the other hand, the mosquito is known to survive in arid conditions (total annual precipitation ~200 mm) [30] if there is a sufficient density of artificial containers which can provide a suitable breeding site.For this reason, as well as host-availability, human population density is a crucial variable in the establishment of Ae. albopictus.The main biological feature that has supported the spread of Ae. albopictus towards central Europe is its ability to lay diapausing eggs which can survive relatively cold winters.Thus, the mean temperature of the coldest month is usually the most significant variable limiting the northward expansion of the mosquito with a direct influence on diapausing eggs [25,33,34].The mean annual temperature has the most significant impact on the inter-annual population density.Night-time light was identified as the most important variable for predicting the occurrence of Ae. albopictus in a study by Ducheyne et al. [29] using the Gini impurity criterion to assess variable importance in a Random Forest modeling framework.
The first specimen of Ae. albopictus in Montenegro was registered in the capital city Podgorica on 21 August 2001.During the period from 2002-2011, the species rapidly spread southwards and inhabited most of the human dwellings at the Adriatic coast and around Skadar Lake [35,36].The meteorological and climatological conditions in Montenegro are mainly affected by the fact that it is a Mediterranean country with a very complex orography.Elements of the global atmospheric circulation affecting the whole Mediterranean basin are the descending part of the Hadley circulation on the South and mid-latitude circulation in the North [37].Observed and projected CC identified the Mediterranean region as one of the most critical CC hotspots [38].We suppose that the upcoming CC will significantly affect suitability for the appearance of the tiger mosquito.
Studies that examine the projected climate parameters and their influence on the spread of Ae. albopictus point to a pronounced northward spread of the vector and reduced suitability in the current hotspots in southern Europe due to warmer and drier summers, as well as a vertical shift to higher altitudes [39][40][41][42].There is already extensive work on mapping the global risk for Ae.albopictus together with studies focusing on Europe.However, these do not capture the observed data that was collected in Montenegro since 2014.Furthermore, these models fail to capture the regional heterogeneity caused by a very complex orography and a characteristic spatial heterogeneity of human population density and would not support surveillance and monitoring on the country level due to differences in resolution and scale.
In this study, we aim to show the expected shift of Köppen climate zones in Montenegro for two projected periods (2001-2030 and 2071-2100) and to examine the expected impact on the vector population.Firstly, an analysis of CC based on the projections with the Eta Belgrade University (EBU)-Princeton Ocean Model (POM) regional climate model [43] and observations from the national grid of synoptic weather stations for the reference period 1981-2010 is presented.Secondly, a mosquito suitability model-based on the Multi-Criteria Decision Analysis (MCDA) framework with temperature and precipitation parameters is applied and examined in the scope of expected CC.Finally, the model is verified against observed presence/absence data for Aedes albopictus.
First of all, the model was validated against field observations (2001-2014) of Ae. albopictus presence/absence in Montenegro.Then, the model output was used to choose areas in the north of the country where samplings were never performed.Lastly, we revalidate the model using the new findings of Ae. albopictus presence/absence obtained in 2015.

Location and Current Climate of Montenegro
Montenegro is a country located in Southeastern Europe, between 41 • and 44 • N and 18 • and 21 • E. It is a mostly mountainous state covering a relatively small area (13,812 km 2 ), but with a huge terrain configuration diversity.
The specific position of Montenegro, from the aspect of weather and climatic conditions, stems from the fact that it is located in the zone of intensive thermal asymmetries, between the cold regions of northern Europe and the hot regions of northern Africa.The intensive exchange of the warm air masses (which move towards the north) and the cold air masses (which move toward the south) takes place above Montenegro.Very often there are collisions and mixing of air masses with extremely different physical and meteorological characteristics.Orography, mountains and valleys, their orientation, the proximity of the sea and the meridian provision of the Adriatic Sea, the vicinity of a large water surface such as the Mediterranean Sea and the vicinity of a large land surface on the north have a significant influence on weather and climate in Montenegro.
The position of Europe and Montenegro is such that large action systems such as the Genovian cyclone, the Adriatic cyclone, the Icelandic depression, the Black Sea depression, the Azores anticyclone, the Siberian anticyclone, the Central European anticyclone, the Arctic cold front from the north, and the warm tropical front from the south have a substantial impact on weather and climate.
From the mosquito distribution point of view, of particular importance is the spatial heterogeneity of the human population distribution.According to the Statistical office data for 2018 Montenegro, the population consists of 629,219 inhabitants, while 31% is situated in two major cities Podgorica and Nikšić.In the coastal region, the number of citizens in towns varies between 6000 and 20,000 but the number of tourists, which is increasing each year (2 million in 2017), significantly changes this picture particularly in June (130,000), July (160,000) and August (180,000).
Climate characteristics of Montenegro are derived from the climatological database of the Institute for Hydrometeorology and Seismology of Montenegro for the reference period 1981-2010 for synoptic stations listed in Table 1.The complete set of data used in the study includes daily values for mean air temperature ( • C) and daily total precipitation (mm).According to the Köppen classification [44], the climate of Montenegro is dominated by the C-mild temperate/mesothermal and D-continental/microthermal climate types.Climate zones and their characteristics can be shortly described as follows: • Csa-mild temperate climate with extremely hot summer and exerted dry period during the summer.The mean temperature of the coldest month is between −3 • C and 18 • C while the mean temperature of the warmest month is above 22 The climate of Montenegro is dominantly affected by the cyclone of Genoa and the Siberian anticyclone, combined with a large variety of orographic effects.The Central and Northern part of Montenegro, together with mountain regions on the eastern border, exert characteristics of the continental and mountain climate, without a dry period and with an apparent influence of air coming from the Mediterranean Sea.
It particularly affects precipitation regimes and the mean temperatures of the coldest and the warmest month.It is worth mentioning that in the mountains above the Kotor bay, the rainiest area in Europe is situated with total annual precipitation of 4600 mm.Part of the coastal and central region, the southern and south-eastern parts of the country are characterized by hot and dry climate particularly in the summer.

Field Experiment
From 2001 to 2014 Ae. albopictus eggs and larvae were sampled with different sampling intensity.The differences of the sampling efforts were due to the scope of the different projects realized during the period.All sites were chosen in urban settlements (Aedes albopictus is an urban mosquito species) by an experienced entomologist.Sampling has been performed at settlements distributed over the entire length of the coastline of Montenegro (293 km) and in the central part of the country.All settlements were sampled during one or more year within this period.In total, 237 different sites were sampled, and 2245 collections of eggs and larvae were accomplished.In the survey program, plastic ovitraps were used [45].The ovitraps (black plastic containers, dimensions: bottom diameter = 10 cm, top diameter = 12 cm, height = 13 cm, total volume = 1.0 Liter, water volume = 750-800 mL, masonite strip for egg collection 15 × 2 cm) have been positioned from July to September (in the period of the year when Ae. albopictus is most active) and checked every 7-14 days.The ovitraps were positioned on the ground in vegetated and shaded sites, with a free space above of at least 1 m [46].The samples were labeled and brought to the laboratory for identification.The number of eggs oviposited on masonite strips was counted in a stereomicroscope.After counting, the strips were placed in glass pots with water collected from the traps for egg hatching and the development of adults.Adult mosquitoes were identified morphologically using a diagnostic key [47].In addition, during ovitrap positioning and checking, adult Ae. albopictus specimens were sampled by Human Landing Collection (HLC) [48] using manual aspirators.
The sampling protocol in 2015 was the same as before but designed to provide additional validation and reveal the presence/absence of Ae. albopictus in the northern part of Montenegro, where no samplings were performed before (this, predominantly mountainous part of the country was considered unsuitable for the establishment of Aedes albopictus).Areas and habitations with different suitability levels predicted by our model were selected.Within them, sites where the main roads coming from the highly infested Montenegrin coastal areas intersect different suitability areas were chosen for sampling/trapping locations.In 2015, 240 ovitrap collections at 48 sites distributed over northern Montenegro were performed.In addition, 60 sites (406 ovitrap collections) in the southern part were sampled.During 2016, 192 ovitrap collections at 64 sites, 4 larval dipping session, 11 HLC, 25 gravid trap nights, and 35 nights of sampling by dry ice-baited EVS traps were performed in northern Montenegro.

Global Circulation Model and Regional Climate Model Description
To study the impact of CC on the suitability of the establishment of the Asian tiger mosquito in Montenegro, we analyzed (a) expected changes in temperature and precipitation regimes in order to make a better assessment of the local changes of climate and (b) the shift of climate zones according to the Köppen climate classification.The observed and expected impact of CC on temperature and precipitation regimes is assessed comparing the observed climatology (1971-2000 and 1981-2010) with the expected climate conditions (2001-2030 and 2071-2100).
Climate simulations used in the study are the result of the Eta Belgrade University (EBU)-Princeton Ocean Model (POM) model runs for the A2 scenario over the 2001-2030 and 2071-2100 integration periods [43] for 7 selected sites (Table 1).The EBU-POM is a two-way coupled regional climate model (RCM).The atmospheric part is the Eta/National Centers for Environmental Prediction (NCEP) limited area model (resolution 0.25 • × 0.25 • on 32 vertical levels; centered at 41.5 • N, 15 • E, with boundaries at ±19.9 • W-E and ±13.0 • S-N; defined on the Arakawa E-grid), while the oceanic part is the POM (resolution 0.20 • × 0.20 • on 21 vertical levels).The driving global circulation model (GCM) is the ECHAM5 (ECMWF Hamburg Atmospheric Model 5) model [49] coupled with the Max Planck Institute Ocean Model (MPI-OM) [50].More details about model integrations and performed bias correction for the Serbia and Montenegro regions can be found in papers Mihailović et al. [51] and Petrić et al. [42].The POM model was set over the Mediterranean Sea without the Black Sea; for other open seas, the sea surface temperature from the GCM was used as a bottom boundary condition.To reduce the systematic model error (model bias) in the key climate variables, the statistical method of bias correction [52] was applied to surface air temperatures and daily precipitation using the modeled and observed daily time series of variables over the period .Current climate scenarios are used in the study as a continuation of work on Climate change national communication for Montenegro which is based on the same climate simulations.

Vector Distribution Model
A Mechanistic Multi-Criteria Decision Analysis model was used for modeling the climatic suitability for the establishment of Ae. albopictus.This model employs sigma fuzzy membership functions to generate a single suitability score on a scale of 0-100.The membership function is a smooth, continuous curve that defines how each value of the input variable is mapped to a membership value corresponding to suitability.The model combines five climatological parameters based on standard climate normals for the observed and projected values for temperature and precipitation.The normals were computed as period averages for 1981-2010, 2001-2030, and 2071-2100.The model is based on extending previous work [27,29,40,42,53] on the abiotic factors and governing thresholds which affect the life cycle of Ae. albopictus.
Based on the local characteristics of the study area, we selected the following climate normals as input for the process based MCDA model: (1) mean annual temperature (T a ), (2) total annual precipitation (H a ), (3) mean summer temperature (June-July-August) (T jja ) (4) mean January temperature (T Jan ) and (4) precipitation frequency expressed as number of days within a year with rain >1 mm (H ν ).In addition to this, we use night-time light data as a proxy for human population density and an elevation raster as a corollary variable for mean annual temperature.
Since the empirical abiotic thresholds introduce sharp transitions in the simulated climatic suitability for the considered climactic parameters, which does not correspond to the reality of the process, sigma functions are used so each value of the parameter is mapped to a membership value using a sigma fuzzy membership function for the suitability on a scale of 0-100.
T c and S c in Table 2 are the critical threshold and saturation values which were used to calculate the sigma functions.The threshold value is generally defined as the value at which no further decrease in the output variable occurs despite a further decrease in the input variable.While the saturation value is the value at which no further increase in the output variable occurs despite a further increase in the input variable.For T jja , a symmetric sigmoid function was constructed with maximum suitability corresponding to values between 16 • C and 32 • C. The threshold values defined by Proestos et al. [27] and Petrić et al. [42] were used (Table 2).The final suitability is calculated using a weighted geometric mean with equal weights.The output for all 7 points in the meteorological surface network (Table 1) is interpolated using a Cokriging spatial interpolation algorithm with an elevation layer from the SRTM30 dataset, aggregated to a 30 arc second resolution (Diva-GIS, http://www.diva-gis.org/gdata)to create a density grid for the simulated climatic suitability for the establishment of Ae. albopictus.Finally, as a proxy for population density, we use cloud-free composites of NLT (https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html).The NLT raster is scaled to 0-100 and the weighted value is applied to the MCDA output for each pixel using the Raster Calculator tool in the QGIS software [54].
The performance of the model was evaluated by calculating the receiver operated characterizing (ROC) and the area under the ROC curve (AUC) for observed presence-absence data for Ae.albopictus for 2015.The AUC can be understood as the probability of a classifier ranking a randomly chosen positive event (the presence of mosquito) higher than a randomly chosen absence.Useful predictive models have an AUC of about 0.7, excellent models would have an AUC > 0.9, while a random method simulation would have an AUC < 0.5.

Current and Future Shift of Climate Conditions in Montenegro
The recent and future changes in the precipitation and temperature regimes for Montenegro are presented in this section.

Expected Changes in Precipitation Regime as a Result of CC
As all Mediterranean countries, Montenegro strongly depends on precipitation during the cold part of the year.The observed variation of annual precipitation (Figure 2) indicates a slight decrease of precipitation during spring and summer but a more pronounced increase during winter time for 2001-2030.Particularly interesting is the observed (19.9%) and predicted (20.3%) increase of precipitation in June which is expected to be typical for mountain regions (KO, PLJ, ZA), while in coastal areas (BA, HN, UL), a slight decrease in precipitation in summer is expected, which can additionally enhance the already hot and dry weather in summer (Table 3).During the winter (DJF) and spring (MAM) periods in the region of the Adriatic coast, an increase of precipitation is expected, while in the mountains, a slight decrease of precipitation can appear particularly during winter.Since high amounts of precipitation are typical for this region, expected changes will not significantly affect the overall precipitation regime.It is expected that the ending of the century will be highlighted with a significant drop in annual precipitation, particularly from spring to autumn (April-September) (Figure 2).Still, individual differences can be noticed between the coastal and the high mountain region.
While the summer decrease is particularly exerted in the coastal area (Table 4), an increase in precipitation is expected during winter.This increase is in accordance with observed mean latent heat flux trends for 1958-2011 period in the Montenegro region [55].During the winter (DJF) and spring (MAM) periods in the region of the Adriatic coast, an increase of precipitation is expected, while in the mountains, a slight decrease of precipitation can appear particularly during winter.Since high amounts of precipitation are typical for this region, expected changes will not significantly affect the overall precipitation regime.It is expected that the ending of the century will be highlighted with a significant drop in annual precipitation, particularly from spring to autumn (April-September) (Figure 2).Still, individual differences can be noticed between the coastal and the high mountain region.
While the summer decrease is particularly exerted in the coastal area (Table 4), an increase in precipitation is expected during winter.This increase is in accordance with observed mean latent heat flux trends for 1958-2011 period in the Montenegro region [55].The expected overall decrease of precipitation in Montenegro is a result of the changes in general atmospheric and oceanic circulation which will, in the same way, affect the whole Mediterranean area through an increase in the mean surface pressure and frequency of anticyclonic circulation as well as the weakening of the local Mediterranean storm track [56][57][58][59].However, the expected increase of precipitation in mountain areas (Table 3) is more related to local orographic effects and the expected increase of the frequency of extreme weather events.Namely, the current and the projected increase of sea surface temperature [60], particularly in the warm period of the year, will produce an increase of the air specific humidity.Since warmer air can hold more water (1 • C in temperature increase, raises the water holding capacity by 7%), local winds blowing from the sea to inland will transfer more humidity to the mountain region close (1-5 km) to the coast.Due to high slopes, warm moist air will be forced to rise, producing intensive cumulus convection and a high intensity of precipitation.It will affect the mean precipitation per rainy day and the number of days with excessive precipitation on high attitudes.

Expected Changes in Temperature Regime as a Result of CC
Variations of the temperature normals presented on Figure 3 indicate that the current climatology is warmer than the previous one (with an exception in February) with a particular temperature increase during spring (0.8 • C) and summer months (1.0 • C)-the most critical periods for mosquito development and activity [10].
According to expected temperature changes, the beginning of the XXI century should bring a warmer winter and slightly colder autumn temperatures (Table 5).However, according to the scenarios, these differences will become much more pronounced towards the end of the century leading to a warmer winter and a much warmer summer period (Table 6).Obtained results are completely in agreement with findings presented by Ciardini et al. [37] for the Mediterranean region.Brunet et al. [61] and Toreti et al. [62] pointed out that the increase of temperature measured in the coastal areas of the Mediterranean is larger than at the global scale and with remarkable seasonal and sub-regional differences.Our findings indicate that in addition to coastal areas, inland locations in Montenegro (KO, PL, PO) experience a significant increase of annual and (all) seasonal temperatures by the end of the century (Table 6).Data presented in Table 7 can be used in favor of the above sentiments, indicating higher temperature deviations for 2001-2010 with respect to the 1981-2010 climatology than it was projected by the climate scenarios.4.During the 2001-2030 period, it is expected that Csb will be replaced with Csa, while Cfbx and Dfbx, will be replaced with Cfb (Figure 4).By the end of the century, there is an obvious expansion of Csa and it is expected that Cfb will be completely replaced by Csax".4.During the 2001-2030 period, it is expected that Csb will be replaced with Csa, while Cfbx and Dfbx, will be replaced with Cfb (Figure 4).By the end of the century, there is an obvious expansion of Csa and it is expected that Cfb will be completely replaced by Csax".

Expected Changes in the Climatic Suitability of Ae. Albopictus
The current (1981-2010) and expected distribution of Ae. albopictus for the periods 2001-2030 and 2071-2100 are shown in Figure 5.The species is mainly abundant across the coast and in the valley of Ribnica and Morača, and the Zeta plain on the southeast.This spatial distribution is preserved through all three integration periods with a pronounced shift to higher altitudes by the end of the century resulting in a high suitability for the whole country.Opposing contributions from an expected increase in winter and decrease in summer temperatures for the northern regions result in little change in suitability for the first half of the XXI century.On the other hand, there is a significant increase in suitability as a result of the expected change in both the seasonal and mean annual temperature expected by the end of the century.

Expected Changes in the Climatic Suitability of Ae. Albopictus
The current (1981-2010) and expected distribution of Ae. albopictus for the periods 2001-2030 and 2071-2100 are shown in Figure 5.The species is mainly abundant across the coast and in the valley of Ribnica and Morača, and the Zeta plain on the southeast.This spatial distribution is preserved through all three integration periods with a pronounced shift to higher altitudes by the end of the century resulting in a high suitability for the whole country.Opposing contributions from an expected increase in winter and decrease in summer temperatures for the northern regions result in little change in suitability for the first half of the XXI century.On the other hand, there is a significant increase in suitability as a result of the expected change in both the seasonal and mean annual temperature expected by the end of the century.

Discussion
The observed and projected CC identified the Mediterranean basin as one of the most important CC hotspots [38].The current (1981-2010) climatology is warmer than the previous one  with a particular temperature increase during spring and summer months, the most important periods for mosquito development and activity [10].Warmer winter and slightly colder summer temperatures are observed.However, according to the scenarios, these differences will become much more pronounced towards the end of the century leading to a warmer winter and much warmer The contribution of this change is most significant for the northern highlands.While for the coastal areas where the annual temperature was high CC did not induce insignificant changes in suitability.Similarly, the expected drop in precipitation does not significantly affect the suitability since the annual precipitation is still above the threshold for egg survival. The

Discussion
The observed and projected CC identified the Mediterranean basin as one of the most important CC hotspots [38].The current (1981-2010) climatology is warmer than the previous one  with a particular temperature increase during spring and summer months, the most important periods for mosquito development and activity [10].Warmer winter and slightly colder summer temperatures are observed.However, according to the scenarios, these differences will become much more pronounced towards the end of the century leading to a warmer winter and much warmer

Discussion
The observed and projected CC identified the Mediterranean basin as one of the most important CC hotspots [38].The current (1981-2010) climatology is warmer than the previous one  with a particular temperature increase during spring and summer months, the most important periods for mosquito development and activity [10].Warmer winter and slightly colder summer temperatures are observed.However, according to the scenarios, these differences will become much more pronounced towards the end of the century leading to a warmer winter and much warmer summer temperatures.An overall decrease in annual precipitation is expected by end of the century with an exception for mountain regions in which this is mitigated by local orography and an expected increase in sea surface temperature which will affect mean precipitation per rainy day and the number of days with excessive precipitation.Since high amounts of precipitation are typical for this region, the expected changes will not significantly affect the overall suitability for the mosquito, however, the extreme shift in the mean January and annual temperature will have a significant impact on the spatial spread of the vector in Montenegro.The projected shift of the Köppen climate zones derived from climate scenarios used in this study is in accordance with results presented in European study carried out for slightly different integration periods [63].This indicates a less diverse climate in the whole country, dominated by a C (temperate) climate type with an extension of the hot-dry summer region and a possible secondary maximum of precipitation in the mountains.At the end of the XXI century, further unification of climate is expected with a complete change of climate with Csb leading to a hotter, drier, and overall more extreme climate.
We anticipated that the upcoming change of climate will significantly increase the suitability for the establishment and spread of Ae. albopictus in Montenegro.Our MCDA model is based on the EBU-POM regional climate model for Montenegro.The model captures the regional heterogeneity and was used to support the planning of the monitoring/surveillance programs on the country level.The model was validated against field observations of Ae. albopictus presence/absence in Montenegro, furthermore, it was used to plan model-based monitoring in the yet unexplored northern part of the country and validated against the findings.
When the model is compared to the output generated by Caminade et al [40], Kramer et al. [6], Proestos et al. [27], and Fischer et al. [41], we see that the model differs in the northern and central parts of the country while the areas of high suitability along the coast are simulated similarly.The MCDA model by Caminade et al. [40] fails to capture the suitable region along the Zeta river in the Podgorica-Nikšić stretch.Furthermore, by employing a more rigorous threshold of 400 mm for annual precipitation, their model simulates a slight drop in suitability for the 2030-2050 period following a projected drop in precipitation.We expect this not to play a significant role for the densely populated coast where human water-keeping habits can provide ample breeding opportunity for the mosquito.The suitability simulated for the same integration periods for Serbia [42] shows a complementary trend on the border with Montenegro, with a very low suitability for the 2001-2030 period, and the border becoming increasingly suitable by the end of the century allowing for an easier spread of the vector through road traffic between Serbia and Montenegro.
The overall model's predictive performance was high, and the second validation proved the model usefulness for planning and performing the targeted surveillance in northern Montenegro in 2015.A model-based focusing of the field sampling efforts might be very useful in situations where only limited resources are available for the surveillance of Ae. albopictus.A good example of this is the surveillance, based on a similar MCDA model [42], that was conducted in 2018 in the Vojvodina Province (northern Serbia), which revealed the presence of Ae. albopictus in 15 out of 16 targeted localities [64].
The mean temperature of the coldest month (indicating overwintering conditions) and the mean annual temperature are identified as the most significant variables limiting the spread and population density of Ae. albopictus [25,33,34].An expected decrease in suitability in the current hotspots in southern Europe due to warmer and drier summers is indicated by several models [39][40][41].Even though the CC scenarios indicate that these differences will become much more pronounced for Montenegro towards the end of the century, leading to warmer winters and much warmer summers, the suitability of current hotspots of Ae. albopictus in Montenegro will be preserved through all three integration periods.The model projections show a pronounced shift to higher altitudes.In Montenegro, the current vertical limit for climatic suitability for Ae.albopictus based on the simulated results is 1000 m.a.s.l.This is in accordance with the results presented for Serbia, where the current vertical limit is estimated to be 1042 m.a.s.l.[42].We expect the mosquito to be present in altitudes up to 1200 m.a.s.l. by the end of the century.This will result in a high suitability for the whole region.
In Asia, Ae. albopictus is better adapted to peridomestic settings, with vegetation as its preferred resting sites than domestic areas dominated by Ae. agypti [65].In Europe, it invades both domestic and peridomestic environments since highly populated urban areas provide abundant breeding sites with low predator pressure [66]."Container type" breeding sites in populated areas, preferred by larvae of Ae. albopictus, are much less dependent on total annual precipitation levels than breeding sites supporting the development of other mosquito species [67,68].In Spain, the species has already been observed in dry areas (200 mm) which do not meet the conservative threshold for annual precipitation (≥500 mm/year) and should not support establishment [30,69].In our model, the expected drop in annual precipitation (particularly from April to September) expected for the Montenegrin climate at the end of the century did not influence the spatial distribution of Ae. albopictus.
It seems that the existence of high mountain ranges (1000-2000 m.a.s.l.) with only a few main traffic routes slowed down but did not prevent the northward spreading of Ae. albopictus in Montenegro.Aedes albopictus is a generalist that readily adapts to diverse environmental conditions in both tropical and temperate regions [70].Spread and population density of this species is also influenced by a wide array of abiotic and biotic variables such as transportation intensity, land use, NDVI, night-time light, temperature, total annual precipitation, and relative humidity [6,9,[22][23][24][25]27,29]; and various socio-economic factors such as water-keeping habits, rate of urbanization, housing quality, etc. [31,32,71].It is advisable to observe the mosquito population and environmental parameters for a particular local population for better understanding the vectorial capacity and planning the control measures [72].The same could be applied in downscaling the vector model, taking into consideration the local particularity of the biotic, abiotic, and the socio-economic factors, e.g., using models which capture the regional heterogeneity and provide fine resolution maps for model-based monitoring/surveillance.New studies identify the human presence (host and container availability) as the most important factor [29,30].
From a public health perspective, Aedes albopictus is a very important invasive species because of its role in the transmission cycles of pathogens such as Dengue and Chikungunya, with over 20 million Dengue cases occurring annually [1].There is a significant threat of infected air passengers arriving in Europe to a country where the mosquito vector is established [73].For this reason, screening for infectious diseases should be performed at international airports with special attention being dedicated to areas that are highly suitable for the mosquito population.
Field measurements following the transition of the vector to higher altitudes at a regional level might offer an excellent opportunity to examine the temperature thresholds related to the overwintering and seasonal activity of the species.
Further care should be taken when selecting the thresholds for the environmental variables that define the climatic suitability for the species (e.g., scale and resolution, laboratory measurements vs field observation, climate vs weather parameters).

Conclusions
The overall predictive performance of the Model is high and the second validation proved the model usefulness for planning and performing the targeted surveillance.Model-based prioritizing of the field sampling efforts should be used in situations where only limited resources are available for Ae.albopictus surveillance.

Figure 1 .
Figure 1.The current distribution of Aedes albopictus in Europe in June 2018[21].Established (red)refers to an established population (evidence of reproduction and overwintering) of the species has been observed in at least one municipality within the administrative unit; Introduced (yellow)-the species has been detected (but without confirmed establishment) within the administrative unit; Absent (green)-field surveys or studies on mosquitoes were conducted but the species has not been detected within the administrative unit; No data (dark grey)-no sampling has been performed and no data on the species is available within the administrative unit; and Unknown (light grey)-it is unknown whether there are field studies on this species within the administrative unit.

Figure 1 .
Figure 1.The current distribution of Aedes albopictus in Europe in June 2018[21].Established (red)-refers to an established population (evidence of reproduction and overwintering) of the species has been observed in at least one municipality within the administrative unit; Introduced (yellow)-the species has been detected (but without confirmed establishment) within the administrative unit; Absent (green)-field surveys or studies on mosquitoes were conducted but the species has not been detected within the administrative unit; No data (dark grey)-no sampling has been performed and no data on the species is available within the administrative unit; and Unknown (light grey)-it is unknown whether there are field studies on this species within the administrative unit.

Figure
Figure The deviation of the normal monthly temperature for current (1981-2010) and future (2001-2030, 2071-2100) climate in respect to past climate conditions (1971-2000).The deviation for 1981-2010 is shown in blue; for 2001-2030 in orange; and for 2071-2100 in gray.
Climate zones over Montenegro according to the Köppen classification obtained from the simulation by the EBU-POM model for the periods 1981-2010, 2001-2030, and 2071-2100 are shown in Figure Climate zones over Montenegro according to the Köppen classification obtained from the simulation by the EBU-POM model for the periods 1981-2010, 2001-2030, and 2071-2100 are shown in Figure
Figure S2: Observed catches of Aedes albopictus for 2015, presences are represented by red dots, absences are represented by blue dots.

Table 1 . The geographical position and climate characteristics (1981-2010) of the selected locations (*-Stations included in the MCDA mosquito model; † -Stations included in the Köppen climate analysis). Station Name with Short Code Latitude and Longitude Altitude (m) Annual Temperature ( • C) Annual Precipitation (mm)
Csb-mild temperate climate with warm summers and a drought period which commonly occurs during the summer.Mean temperature of the warmest month is below 22 • C, while the mean temperature of at least four months is above 10 • C; • Csbx-this is a specific climate type which differs from the previous ones because of the presences of a secondary precipitation maximum (x"); • Cfb-climate type with same temperature characteristics as Csb but with a uniform distribution of precipitation over the year; • C; • • Cfbx"-same as Cfb but with present secondary precipitation maximum (x"); • Dfbx"-boreal climate without an exerted dry period over the year.The mean temperature of the coldest month is below −3 • C, while the mean temperature of the warmest month is above 10 • C. The mean temperature of the warmest month is below 22 • C, the mean temperature of the least warm four months is above 10 • C. Primary precipitation maximum appearing in autumn and secondary in spring.

Table 2 .
The climatic thresholds limiting the establishment of Ae. albopictus.

Table 4 .
The relative deviation (%) of average annual precipitation for the 2071-2100 climate in respect to referent climatology 1981-2010.

Table 6 .
The deviation of the normal temperature for 2071-2100 with respect to the referent climatology of 1981-2010.

Table 5 .
The deviation of the normal temperature for 2001-2030 with respect to the referent climatology of 1981-2010.

Table 6 .
The deviation of the normal temperature for 2071-2100 with respect to the referent climatology of 1981-2010.