Potential distribution of dominant malaria vector species in tropical region under climate change scenarios

Risk assessment regarding the distribution of malaria vectors and environmental variables underpinning their distribution under changing climates is crucial towards malaria control and eradication. On this basis, we used Maximum Entropy (MaxEnt) Model to estimate the potential future distribution of major transmitters of malaria in Nigeria—Anopheles gambiae sensu lato and its siblings: Anopheles gambiae sensu stricto, and Anopheles arabiensis under low and high emissions scenarios. In the model, we used mosquito occurrence data sampled from 1900 to 2010 alongside land use and terrain variables, and bioclimatic variables for baseline climate 1960–1990 and future climates of 2050s (2041–2060) and 2070s (2061–2080) that follow RCP2.6 and RCP8.5 scenarios. The Anopheles gambiae species are projected to experience large shift in potential range and population with increased distribution density, higher under high emissions scenario (RCP8.5) and 2070s than low emission scenario (RCP2.6) and 2050s. Anopheles gambiae sensu stricto and Anopheles arabiensis are projected to have highest invasion with 47–70% and 10–14% percentage increase, respectively in Sahel and Sudan savannas within northern states in 2041–2080 under RCP8.5. Highest prevalence is predicted for Humid forest and Derived savanna in southern and North Central states in 2041–2080; 91–96% and 97–99% for Anopheles gambiae sensu stricto, and 67–71% and 72–75% for Anopheles arabiensis under RCP2.6 and RCP8.5, respectively. The higher magnitude of change in species prevalence predicted for the later part of the 21st century under high emission scenario, driven mainly by increasing and fluctuating temperature, alongside longer seasonal tropical rainfall accompanied by drier phases and inherent influence of rapid land use change, may lead to more significant increase in malaria burden when compared with other periods and scenarios during the century; especially in Humid forest, Derived savanna, Sahel and Sudan savannas.


Introduction
With nearly half of the world's population at risk in 2017 [1], malaria continues as one of the highest killer infectious diseases in the world after lower respiratory tract infections, diarrhoeal a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 and proliferation of infectious disease vectors / disease transmission [9,19,26,27]. Over the course of the 21st century, several models have been used in predicting the impact of climate change across a variety of sectors and regions for risk analysis and adaptation strategies [28][29][30][31]. With respect to malaria risk assessment, several studies have been carried out on the impact of changing climate scenarios on the occurrence and distribution of malaria vectors across the world. Most of these studies focused on Anopheles gambiae sensu stricto (s.s.) and Anopheles arabiensis which are major siblings of the species complex-Anopheles gambiae sensu lato (s.l.) [12,32,[33][34][35][36][37]. An. gambiae s.s. and An. arabiensis are dominant most efficient vector species of human malaria in the Afrotropical Region [38][39][40][41]. Tonnang et al. [35] used CLIMEX simulation model to predict possible shifts in An. gambiae s.s. and An. arabiensis boundaries southward and eastward of Africa rather than jumps into quite different climatic environments under climate change scenarios. This was previously observed by Peterson [34] who applied ecological niche model, Genetic Algorithm for Rule-set Prediction (GARP) to predict potential distributional shift from west to east and west to south of Africa for An. gambiae s.s. and An. arabiensis, respectively [12]. Drake and Beier [37] used Low Bias Bagging for One-Class-Classification (LOBAG-OC) to project extensive reductions in An. arabiensis Anopheles species sampling points reprinted for illustrative purposes only from Okorie et al. [48] under a CC BY 4.0 license, with permission from PLOS ONE [38]. https://doi.org/10.1371/journal.pone.0218523.g001 An. gambiae species under low and high emissions scenarios in Nigeria habitat (as a result of climate change) in Western and Central Africa, and parts of some Eastern, North-eastern and Southern African countries. Also, under three climate scenarios, Ren et al. [32] used maximum entropy algorithm (MaxEnt) to predict prospective changes (increases) in future distribution of the four dominant vectors in China, different in each climate change scenario.
Previous studies demonstrated the distribution of Anopheles species under changing climate scenarios more in terms of general spatial distribution across bioclimatic domains within areas of interest, with little or no consideration of the influence of ecological subregions [12,32,[33][34][35][36][37]. In this study we estimate future spatial distribution of dominant malaria vectors and relative climate impacts under low and high emissions scenarios that may affect malaria risk pattern under different ecosystem in tropical region, using a case study of Nigeria. This study builds on an earlier study [38] by assessing prospective changes in the distribution of An. gambiae sensu lato and its two major siblings: An. gambiae sensu stricto, and An. arabiensis across bioclimatic domains within ecological zones and administrative regions in the country, in relation to elemental conditions under both low and high emissions scenarios in 2041-2060 and 2061-2080, respectively. The choice of low and high emissions scenarios may have implications for policies towards the type and level of preparedness and intervention needed to reduce malaria burden caused by malaria vectors within states and regions along ecological gradients, in a world with and without climate policies [15,28,29]. MaxEnt, a well-established and successful method for species distribution modelling [42] was used based on its reliable usage and high performance among other modelling methods [43] for estimating disease vectors distribution and disease transmissions [32,[44][45][46][47].

Study area and mosquito occurrence data
Georeferenced An. gambiae species data (Fig 1) sampled between 1900 and 2010 across Nigeria (study area) was obtained from Nigeria Anopheles vector database (https://doi.org/10.1371/ journal.pone.0028347.s001) [48]. Nigeria is located in the tropical region within Latitudes 4 o and 14 o north of the Equator and between Longitudes 2 o 2' and 14 o 30' east of the Greenwich Meridian (Fig 1) [20,38]. It has 36 regional states with Federal Capital Territory, Abuja within six broad geopolitical regions (Fig 1). Species distribution is distinctively defined by ecological zones on the basis of combinations of soil, landform and climatic characteristics in the country [38,49,50]. Climates are tropical at the coastal south within Humid forest and Derived savanna, sub-tropical further inland within Derived and Guinea savannas, semi-arid in the far north within Sudan and Sahel savannas, and temperate within Mid Altitude zone of Jos and Mambilla plateaus (Fig 1) [38,51]. Mean annual temperature ranges from 26˚C to 33˚C [20,23,50], contingent on climate zone. The north records rainfall between 500 mm and 750 mm annually, while the south records between 1,200 mm and above 4000 mm [23,38,51]. Mean elevation within the south is about 150 m above sea level, 600-700 m in the north, and 1,500-2,100 m within Mid Altitude zone (Fig 1) [20,38,52]. The country has over 190 million people [53] with about 90% at risk of malaria [54], estimated to be approximately 411 million and 794 million people in 2050 and 2100 [53,55] with 95% and above 97% of the populations projected to be at the risk of malaria, respectively [8]. Its urban population is expected to increase from 51.9% in 2019 to 72% in 2050 [55], with deforestation rate of about 3.5% per year [56].

Environmental variables
To model impact of climate change under low and high emissions scenarios on malaria vector species distribution in Nigeria, nineteen bioclimatic variables with about 1 km 2 spatial resolution were obtained for 1960-2080 from WorldClim (http://www.worldclim.org)-Global Climate Data, the global climate models (GCM)-community climate system model version 4 (CCSM4) based on RCP2.6 and RCP8.5 [57]. The bioclimatic variables were: annual mean temperature (bio1), mean diurnal range (bio2), isothermality (bio3), temperature seasonality (bio4), maximum temperature of warmest month (bio5), minimum temperature of coldest month (bio6), temperature annual range (bio7), mean temperature of wettest quarter (bio8), mean temperature of driest quarter (bio9), mean temperature of warmest quarter (bio10), mean temperature of coldest quarter (bio11), annual precipitation (bio12), precipitation of wettest month (bio13), precipitation of driest month (bio14), precipitation seasonality (coefficient of variation) (bio15), precipitation of wettest quarter (bio16), precipitation of driest quarter (bio17), precipitation of warmest quarter (bio18), and precipitation of coldest quarter (bio19). The climate data was based on three climate periods: the baseline time period 1960-1990 which observed data interpolations is referred to as current conditions [57][58][59], and future climate periods 2041-2060 (2050s) and 2061-2080 (2070s) [57]. The 1960The -1990 Climate Normals used in this study from WorldClim [57] serves as an implicit predictor of the conditions characteristic of the future up to 2005 when it was created. It also serves as a stable benchmark against which changes in climate observations in 2041-2060 and 2061-2080 are compared [59]. Also incorporated in the model were land use land cover (lulc) data obtained from U.S. Geological Survey data release [60], and Digital Elevation Model (DEM) obtained from the Consultative Group on International Agricultural Research-Consortium for Spatial Information (CGIAR-CSI) [38,61].

Modelling procedures and data analysis
Using Maximum entropy algorithm (MaxEnt) model version 3.3.3k [62], model operation was performed following the procedures in our previous study [38]. As a general-purpose machine learning technique of ecological niche modelling, MaxEnt estimates potential species distribution using species presence-only records and environmental layers [42,63]. In order to predict future occurrence and distribution of the Anopheles species, projected bioclimatic / environmental variables representing 2041-2060 (2050s) and 2061-2080 (2070) were added to projection directory in MaxEnt. Using all available data without having an independent dataset under sub-sample replicated run type, occurrence data for each species was split twenty-one times into training (75%) and testing (25%) subsets; performed to test the predictiveness / performance of the MaxEnt model specified by Area Under the Receiver Operating Characteristic (ROC) Curve (AUC) [38]. The AUC measures the ability of the model to discriminate between sites where a species is present (y = 1) against where it is absent (y = 0), given as a plot of sensitivity against specificity [38,[64][65][66]. The AUC was considered more in terms of model's predictiveness, which according to Lobo et al. [67], the value of AUC tells of the degree to which a species is restricted along the range of predictor conditions in the study area, in order that presences can be told apart from absences [43]. To guard against bias in datasets, following large ranges of the documented species [64] relative to the study area (especially in the North Eastern part of the country) (Fig 1) [38], a bias layer was created to provide MaxEnt with background samples (pseudo-absences) [43,[68][69][70]; defining locations with documented Anopheles species as 1, and "no data" in the grid of unsampled pixels [71]. More environmentally distant pseudo-absences from the presences (as applied in this study) increases the rate of well-predicted absences and the AUC scores [67]. Over-prediction / under-prediction of the relationships by the model was controlled by increasing maximum iterations from 500 (default) to 5000, allowing the model to have adequate time for convergence; and regularization was left at 1 (default) to minimise model over-fitting [38,43,65].
The combined environmental variables predicted the probability of species occurrence by producing a point-wise mean (model images), created by MaxEnt model based on maximum entropy of optimal conditions of included environmental variables that match the threshold value within species occurrence records [38,65]. These images were classified in ArcMap for distributions of the studied mosquito species, and reclassified into suitable and unsuitable habitats using 10 percentile training presence logistic threshold provided by MaxEnt. Suitable habitat was defined to include 90% of the data used to develop the model, to take care of some errors inherent in the datasets [38,65]. The mean distribution density of each Anopheles species was determined from zonal statistics in ArcGIS, across ecological and geopolitical zones, and in each state [38]. Prevalence (distribution density in percent) was derived from the mean distribution density according to the equation, p x = (y x /y max )100; where p x represents species prevalence in year x, y x represents predicted mean distribution density of species in year x, and y max represents maximum mean distribution density of species estimated by zonal statistics. Mean distribution density less than or equal to one (equivalent to �50% prevalence) defines the probability of the species not occurring, greater than one (equivalent to >50% prevalence) defines species presence, and equal to two (equivalent to 100% prevalence) represents maximum prevalence of species within a zone or state [38]. Percentage change (y) of species between baseline and future climates was derived from the equation y = ((y x -y 0 )/y 0 )100, where y 0 represents species prevalence under baseline climate, y x represents species prevalence in year x. The contribution of environmental variables to the delineation of suitable environments for the Anopheles species were examined using jackknife test of variable importance and analysis of percent contribution of each variable [38,45]. In Jackknife test the training gained by each variable is estimated as if the model was run with one variable, and then compare it to the training gain involving all variables [38,65]. Logical assessment of variables contributions was achieved by evaluating estimates of relative contributions of the environmental variables alongside jackknife plots produced for training gain, test gain and AUC (S1 File) [38,64].

Potential future distributions of An. gambiae species under RCP2.6 and RCP8.5
Noticeable changes have been projected to occur in potential distribution of the dominant malaria vector species complex, An. gambiae s.l. and its two major siblings under the chosen climate change scenarios within bioclimatic and ecological domains in the tropical country-Nigeria. The 'business as usual' high emissions scenario (RCP8.5) without large investments threatens to increase the future geographic range, distribution density and prevalence of An. gambiae species more than the mitigation (low emissions) scenario (RCP2.6). An. gambiae s.l., a species complex of eight reproductively isolated species [40,41] is projected to increase across all bioclimatic domains in the country with total percentage increase of 26.53% and 32.20% under RCP8.5, and 15.20% and 14.89% under RCP2.6 in 2050s and 2070s, respectively (Fig 2). This is expected to translate into large range expansions and high prevalence with increased distribution density within Humid forest ( Fig 3). The Southern Guinea savanna from Oyo state through Kwara state to Niger state (Fig 1) which appears less suitable for An. gambiae s.l. under current climates [38], is projected to experience large invasion with increased distribution density especially under RCP8.5 in 2050s (by 43.14%) and more increasingly so in 2070s (by 51.97%); increasing by 28.44% under RCP2.6 in 2050s but declines to 25.09% in 2070s (Table 1 (Table 1). Also, highlands within Ondo, Ekiti, Edo, Kogi, Enugu and Anambra states are projected to be less suitable for An. gambiae s.l. between 2041 and 2060; but are projected to experience An. gambiae s.l. invasion between 2061 and 2080 under RCP8.5 (Fig 3; S1 Fig). By administrative regions (geopolitical zones), the South East and South West are projected to hit species maximum prevalence in 2070s under RCP8.5 (Table 2). Other regions are estimated to experience large increased mean distribution density in 2050s and larger in 2070s under RCP8.5. Prevalence is projected to be larger in 2050s in the three northern regions than in 2070s under RCP2.6 ( Table 2). According to the model projections, fourteen states will hit maximum prevalence of An. gambiae s.l. in both 2050s and 2070s under high emission scenario (Fig 4). These states are: Abia, Akwa Ibom, Bayelsa, Benue, Ebonyi, Imo, Jigawa, Kwara, Lagos, Ogun, Osun, Oyo, Rivers and Sokoto (Fig 4). Other states are projected to experience a sequential increase in mean distribution density of An. gambiae s.l., with observed anomaly in Kaduna, Katsina and Kano states, An. gambiae species under low and high emissions scenarios in Nigeria where either current climates record higher density than 2050s/2070s, or 2050s record higher density than 2070s (Fig 4).
The mean distribution density of An. gambiae s.s. is projected to increase from current climates across all bioclimatic domains in the country by 10.48% and 14.49% in 2050s and 2070s under high emissions scenario, and by 5.03% and 5.95% under low emission scenario, respectively (Fig 2). This is expected to lead to large geographic range expansion and increased prevalence (percentage increase) within Derived savanna ( Fig 5; S2 Fig). An. gambiae s.s. is expected to maintain its maximum prevalence even beyond 2080 in all states within South Western and South Eastern regions (Tables 1 and 2). Like An. gambiae s.l. under high emissions scenario, An. gambiae s.s. is projected to experience shift in composition and geographic range in Niger state; as boundaries (areas) along Zamfara, Kaduna, Bauchi and Plateau states within Northern and Southern Guinea savannas, and Mid Altitude seem less suitable for An. gambiae s.s. in 2050s, and less so in 2070s (Fig 5; S2 Fig).  (Table 1). This will lead to species invasion and increased distribution density in the less suitable North-Eastern landmass, especially within Lake Chad areas of Borno state (Table 1; Fig 5; S2 Fig). All other states are estimated to experience sequential increases in mean distribution density of An. gambiae s.s. in 2050s and 2070s, respectively. Observed anomaly is estimated to occur in Kaduna, Sokoto and Zamfara states under RCP8.5, just as it occurs in many other states under RCP2.6, either from higher distribution density under current climates than 2050s/2070s, or in 2050s than 2070s (Figs 4 and 5).
The mean distribution density of An. arabiensis is projected to increase across all bioclimatic domains in the country by 4.43% in 2050s and 8.16% in 2070s under high emissions scenario, and by 2.03% and 2.72% under low emissions scenario, respectively (Fig 2). This projected increase may lead to large geographic expansion range and increased prevalence of An. arabiensis within ecological zones; Derived savanna ( (Fig 6; S3 Fig). This species shift results in prevalence anomaly in Kaduna state under both RCP2.6 and RCP8.5 (Fig 4). Anomaly is also observed in Akwa Ibom, Ekiti, Kebbi and Rivers states under RCP8.5; and in Kano, Kebbi, Lagos and Ogun states under RCP2.6 (Fig 4). Mid Altitude zone of Jos and Mambilla plateaus, and highlands along Cameron boundary are projected to experience more presences of An. arabiensis in 2050s under RCP8.5 with percentage increase of 1.87%, but record less presences in 2070s with percentage decrease of 1.78% (Table 1; Fig 6; S3 Fig). Similarly, the An. arabiensis is projected to decrease in population within the Mid Altitude zone under RCP2.6 by 1.78% and 1.68% in 2050s and 2070s, respectively (Table 1; Fig 6; S3 Fig). The Fresh water and Mangrove swamp forests (within Humid forest) of the Niger Delta region are projected to remain less suitable for An. arabiensis under both RCP2.6 and RCP8.5 in 2050s and 2070s Especially in Bayelsa state where environmental variables least favour the occurrence of An. arabiensis (Fig 4). By regions, the model predicted that An. arabiensis will continue its highest prevalence in South West, followed by South East, South South, North Central, North West; and lowest in North East (Table 2). Generally, large increase in population is projected for the studied An. gambiae species under RCP8.5 between 2050s and 2070s. With less than one percentage increase in prevalence, An. gambiae species under low and high emissions scenarios in Nigeria a decline or a likely state of dynamic equilibrium in vectors population is expected to be the case under RCP2.6 (Fig 2). Also, model projections show that sympatric species-An. gambiae s.s. and An. arabiensis [46,56] do not always occur together, nor with other members of the An. gambiae complex such as An. coluzzii, An. melas, etc. (Figs 2 and 4; S1-S3 Figs). An. gambiae s.s. is estimated as the most prevalent species of the An. gambiae complex. It is predicted to occur in many areas least suitable for other members of the complex, resulting in its higher prevalence than the An. gambiae s.l. In this study, An. gambiae s.l. is predicted as a complex species in areas suitable for one or more of its siblings other than only An. gambiae s.s. and/or An. arabiensis (Figs 2 and 4).

Contributions of environmental variables to future distributions of An. gambiae species under RCP2.6 and RCP8.5
MaxEnt predicted that annual mean temperature (bio1), 33˚C which may rise by about 4.9˚C on average under RCP8.5 and 1.4˚C under RCP2.6 [8,20,24], will greatly regulate (in isolation) the total energy inputs in the ecosystem towards high prevalence of An. gambiae s.l. between 2041 and 2080 in Nigeria (S1 File). Cold temperature anomalies throughout the year (bio6) is also projected to greatly influence year-round widespread of An. gambiae s.l. between 2041 and 2080, especially under RCP8.5 when combined with precipitation of driest quarter (bio17) (December-February, projected to decline by 18 mm, 15 mm and 10 mm in December, January and February, respectively, in the Humid forest and Derived savanna) [20] (S1 File). This is also estimated to be complemented by mean temperature of the wettest quarter (bio8) (June-August), temperature annual range (bio7) of more than 28˚C-36˚C by 3˚C-5˚C [20], mean temperature of the warmest quarter (bio10) (March-May), mean temperature of driest quarter (bio9), maximum temperature of warmest month (bio5) (about 35˚C in February in the south and above 40˚C in April in the north) [20,38], changes in land use land cover (lulc), and terrain surface (DEM) (S1 File). However, the occurrence and seasonal distribution pattern of An. gambiae s.l. is expected to be altered if precipitation of coldest quarter (bio19) (June-August, about 211 mm in the arid north and above 2000 mm in the coastal south) is omitted, because bio19 contains the most information absent in other variables.
Mean diurnal range (bio2) between 7˚C and above 16˚C [20,72] is projected to greatly regulate (in isolation) the distribution and prevalence of An. gambiae s.s. between 2041 and 2080 under both low and high emissions scenarios (S1 File). This is estimated to be complemented by precipitation (bio17) and mean temperature (bio9) of driest quarter, and precipitation of coldest (bio19) and warmest (bio18) quarters when used alongside all other environmental variables including terrain surface (DEM) (S1 File). Changes in land use and land cover dynamics (lulc) is estimated as the variable with the most information that is not present in the other variables, regarding the distribution and prevalence of An. gambiae s.s. as well as An. arabiensis in 2050s and 2070s (S1 File). Precipitation of driest quarter (bio17) is projected to highly influence (in isolation) the widespread and seasonal distributions of An. arabiensis between 2041 and 2080 under both low and high emissions scenarios (S1 File). Precipitation of coldest quarter (bio19), precipitation seasonality (bio15), mean diurnal range (bio2), mean temperature of driest quarter (bio9), precipitation of warmest quarter (bio18), precipitation of driest month (bio14), and temperature annual range (bio7) are expected to combine with all other environmental variables to greatly influence year-round occurrence and distribution of An. arabiensis in 2050s and 2070s, respectively (S1 File).

Discussion
In relation to climate drivers, the potential future distribution of An. gambiae s.l. and its two major siblings: An. gambiae s.s. and An. arabiensis in Nigeria is highly worrisome. Climate change under both high and low emissions scenarios is expected to cause large geographic range expansions with increased distribution density, species shifts and invasions in areas previously too cool for the Anopheles species population [34,35]. Higher prevalence of the malaria vectors projected for high emission scenario corroborates with the prediction of ecological models about shift in the distribution of world biomes in response to changes in climate system associated with increased greenhouse gases [38,73]. Confirming earlier claim that areas with low presence of the studied Anopheles gambiae species will possibly experience high prevalence with species migration and invasion, in response to changes in environmental variables fundamental to their distribution patterns [38]. This is evidenced in the scenario that the less suitable areas of Sudan, Sahel, Northern and Southern Guinea savannas within Gombe, Bauchi, Yobe, Adamawa and Borno states in North Eastern region [38], are becoming more suitable for the occurrence and distribution of these Anopheles species under the changing climate regimes. This is expressed in longer rainy season and warmer nights that promote vectors proliferation in a continuous high emissions 21st Century's world [20,25]. The general decline or near dynamic equilibrium in vector species population between 2050s and 2070s under the mitigation scenario-RCP2.6 is assumed to follow indirect effect of mitigation scenario, where increase in global mean temperature is limited to 2˚C as regional mean temperature within the tropical country is limited to about 1.4˚C [8,15,16]. Nevertheless, low environmental suitability has been projected to continue in some areas for An. gambiae s.l., An. gambiae s.s., and An. arabiensis based on topography and other ecological gradients [74,75]. This is of course highly noticeable in some high elevation areas of Kaduna, Zamfara, Kebbi, Niger, Bauchi, Gombe and Borno states within Northern and Southern Guinea savannas, and parts of Plateau, Adamawa, Bauchi, and Taraba states within Mid Altitude zone [38]. The projected scenario also supports the assertion by Siteti et al. [75], that, topographic and human settlement patterns affect the spatial distribution of malarial mosquitoes [38,76]. It is also apparent from this study that An. gambiae s.s. and An. arabiensis are likely to be more widespread with higher relative probability of occurrence in lowlands and areas with human settlements [77]. An. arabiensis exhibits low probability of occurrence in Fresh water and Mangrove swamp forests [12,33,35] along humid Atlantic coast of the Humid forest within South South region, either based on landforms [73,74], soil or climatic characteristics [38,49,50]; corroborating with Ayala et al. [77] who classified these zones unsuitable for An. arabiensis in Cameroon. However, An. gambiae s.s. is expected to continue its highest prevalence amongst other species, in line with Okwa et al. [78] that, An. gambiae s.s. is the most efficient and most widespread within the gambiae complex [38]. According to Olayemi et al. [79], the epidemiological success of An. gambiae s.s. is largely dependent on its highly dynamic ecological behaviour, which it evolved over a long time to take advantage of certain tropical weather conditions that encourage mosquito proliferation and human / vector contact.
The projected An. gambiae s.s. and An. arabiensis expand largely in range and shift in sympatric coexistence under changing future climates [12,46,56], in agreement with Peterson [34] whose future projections on suitable areas for these vectors in Africa match the suitability/ occurrence patterns in bioclimatic domains within Nigeria. The distribution of these highly anthropophilic members of Anopheles gambiae complex appears relative to ecological zones within the tropical region [40]. The prediction of their dominant occurrence in humid and sub-humid (Derived and Guinea savannas) zones of Nigeria in West Africa (in this study) corresponds to the prediction of their environmental suitability in sub-humid zone in East Africa, from west of Mount Kilimanjaro to the coastal plain of northern Tanzania [47]. It similarly corresponds to their wide distribution in humid and sub-humid domains in Madagascar within the southeast coast of Africa [80]. Their distribution within these ecological zones in Nigeria also corresponds closely to expert-based predictions of Lindsay et al. [40], Levine et al. [46] and Moffett et al. [44], who presented climate suitability zones for these species. While the projected distribution of An. gambiae s.s. agrees with Tonnang et al. [12] regarding the region (Nigeria) in a similar study on Africa as Peterson [34], especially under their chosen climate scenario 1 (similar to RCP8.5), the proposed distribution of An. arabiensis does not. Tonnang et al. [12,35] projected An. arabiensis to reduce in range and population within Humid forest and Derived savanna in Nigeria, and concentrate more in Sahel, Sudan and Guinea savannas. Also, Drake and Beier [37] projected An. arabiensis to lose suitable habitats in Nigeria within Africa in 2050, instead of gaining. Although our projections agree with Tonnang et al. [12] regarding increased prevalence of An. arabiensis in Sahel, Sudan and Guinea savannas resulting from species invasion, they do not agree with reduction in suitable areas within Humid forest and Derived savanna, modelled to experience pervasive distribution of An. arabiensis in 2050s and more in 2070s [38,44]. Discrepancies of this sort, apart from the type of climate scenario used, may sometimes be attributed to differences in data sources, the general lack of the vector's distributional data within area of interest, and the type of model used together with computational approaches [44,67,80,81]. Nevertheless, our model predictions of these Anopheles species' fundamental niche characteristically conform to the realised geographical niche in Nigeria based on confirmed occurrence records [48]. MaxEnt model was able to use the included environmental variables to distinguish between presences and potentially unsampled locations [43,67] of the studied Anopheles species with average AUC of 0.7 (S1 File). This value of AUC score associated with the widespread of the species confirms the assertion of Lobo et al. [67] that, low AUC values denote the true generalist nature of the species distribution in a widespread species, which the probability of presence increases steadily with predictor values [38].
The north-eastward species shift in both RCP2.6 and RCP8.5 aligns with Peterson [34] who predicted west to east and west to south potential distributional shift of An. gambiae s.s. and An. arabiensis within Africa, respectively. Our model result also agrees with Tonnang et al. [35] who predicted that shifts in An. gambiae s.s. and An. arabiensis boundaries may occur southward and eastward of Africa rather than jumps into quite different climatic environments under climate change scenarios. The eastward shift, though tending more northward than southward is expected to be the case in Nigeria under both low and high emissions scenarios between 2041 and 2080. This does not overlook potential southward shift based on the north to south climate variation [50] and the location of the country within the tropical region. According to predictions by Moffett et al. [44], the abundance of An. gambiae and An. arabiensis is highest in West Africa (where Nigeria is located), as human population density and urbanisation are highly critical in malaria risk [34,44]. Rapid population expansion and increased deforestation through uncontrolled urbanisation and agricultural intensification [44] converge with energy consumption (dominated by fossil fuels) and poor or inadequate public health / healthcare infrastructure, alongside increased temperatures with high humidity, more recurrent temperature fluctuations from mean diurnal temperature range and longer high seasonal rainfall [33,34,40,82], to impel future increasing prevalence of the modelled malaria vectors in ecological subregions. Densely populated and urbanised Lagos state [83,84] within the humid forest of South West, and many other potentially populated / urbanised states with increased anthropogenic activities including Abia, Akwa Ibom, Anambra, Enugu, Imo, Ogun, Ondo, Osun, Oyo and Rivers states in the humid south, and Kano and Sokoto states within Sudan savanna in the North West are thus likely to be more vulnerable [20,76].
As the country will generally become warmer with high humidity and altered rainfall patterns, particularly under high emissions scenario, seasonal malaria transmission is expected to be higher and more unbalanced across the country [38,84,85], just as vector species prevalence increases and varies within ecolgical and regional zones (subregions) based on varying tropical climate conditions. This is because Anopheles mosquitoes thrive in regions with warm temperatures, humid conditions, and high rainfall [38,47,86,87] regulated by changes in background climate that alter the abundance, range (both latitude and altitude), distribution, and behaviour of malarial mosquitoes, and the life cycle of the malarial parasite, such that patterns of malaria changes [9,34,38,88,89]. The high magnitude of change in the species prevalence due to future changes in tropical climates within the tropical subregions, especially under high emissions scenario, lends credence to the proposition that, a little increase in temperature such as half degree centigrade increase can translate into a 30% to 100% increase in mosquito abundance [26,27]. This will most likely exacerbate malaria prevalence [90] by causing increases in the population at risk in areas where malaria presence is static in the future [91]-in line with World Health Organisation projections of more than 400 million people at risk of malaria by 2070 in Nigeria, under both high and low emissions scenarios [8,9,88].

Conclusions
Anthropogenic climate change and land use dynamics are expected to increase the population and potential range of malaria vectors and exacerbate malaria burden. This study explored the effect of changing climate under low and high emissions scenarios on future composition and potential distribution of dominant malaria vector species in tropical subregions and administrative regions within the tropical country-Nigeria; using data of vectors known locations and environmental drivers of their occurrence and distribution, following algorithms defined in species distribution model-MaxEnt. Results show that An. gambiae s.l. and its two major siblings: An. gambiae s.s., and An. arabiensis will experience increased prevalence between 2041 and 2060, higher in 2061-2080 under high emissions scenario than low emission scenario in both 2041-2060 and 2061-2080; driven mainly by increasing and fluctuating temperature, longer seasonal tropical rainfall accompanied by drier phases, high humidity in dry season from precipitation during warm months, and inherent influence of rapid land use change. Humid forest and Derived savanna within all southern and most North Central states, with highest species invasion in Sahel and Sudan savannas, particularly in north-eastern states are likely to be most impacted. The predicted variability in species composition and distribution based on diversity in climate conditions in the tropical subregions across the country, may define the future spatial epidemiology of these vectors and possible malaria risk pattern. The impending high magnitude of change in prevalence of these dominant malaria vector species may lead to high malaria burden in 2050s, and higher in 2070s; especially under high emissions scenario, mostly dependent on past, current and future anthropogenic emissions and natural climate variability. Hence, there is need for more studies aimed at spatial estimation of future malaria prevalence at a finer scale with respect to vectors dynamics and distribution within the tropical subregions; and in relation to incremental differentials in climate and land use dynamics for the development of more robust adaptation strategies under global change.