Current and future potential distribution of Culex (Melanoconion) (Diptera: Culicidae) of public health interest in the Neotropics

Abstract

Anthropogenic activities are altering ecosystem stability and climate worldwide, which is disturbing and shifting arbovirus vector distributions. Although the overall geographic range of some epidemiologically important species is recognized, the spatiotemporal variation for other species in the context of climate change remains poorly understood. Here we predict the current potential distribution of 9 species of Culex (Melanoconion) based on an ecological niche modeling (ENM) approach and assess spatiotemporal variation in future climate change in the Neotropics. The most important environmental predictors were the mean temperature of the warmest season (27 °C), precipitation during the driest month (50 mm), and precipitation during the warmest season (>200 mm). The best current model for each species was transferred to the future general circulation model IPSL-CM6A-LR, using 2 shared socioeconomic pathway scenarios (ssp1-2.6, ssp5-8.5). Under both scenarios of climatic change, an expansion of suitable areas can be observed followed by a strong reduction for the medium–long future under the worst scenario. The multivariate environmental similarity surface analysis indicated future novel climates outside the current range. However, none of the species would occur in those areas. Even if many challenges remain in improving methods for forecasting species responses to global climate change and arbovirus transmission, ENM has strong potential to be applied to the geographic characterization of these systems. Our study can be used for the monitoring of Culex (Melanoconion) species populations and their associated arboviruses, contributing to develop region-specific public health surveillance programs.

Graphical Abstract

Open in new tabDownload slide

Introduction

Many diseases caused by arboviruses (arthropod-borne viruses) are emerging or reemerging worldwide, largely due to human impact on ecosystems and climate change (Contigiani et al. 2016). Ongoing anthropogenic disturbances in the Neotropics, such as deforestation, conversion of natural landscapes to crops and pasture, and resulting forest fragmentation are major environmental concerns in the phytogeographic provinces (PP) of South America (e.g., Cerrado, Imerí, Madeira, Pantepui, Pará, Rôndonia, Roraima) (Brown et al. 2006, Brando et al. 2014). Fire-related tree mortality in the Amazon and Chaco basins is increasing, while urban spread in the Pampas has led to major reductions in natural areas (Brown et al. 2006, Brando et al. 2014). Most of the anthropogenic activities mentioned alter ecosystem stability and resilience, triggering processes leading to the desertification, and savannization of biomes (Lapola et al. 2009).

Change in the landscape and climatic conditions can affect vector distributions and disease transmission dynamics (Rabitsch et al. 2017, Ryan et al. 2019). It is estimated that there are hundreds of species of mosquito-borne arboviruses in the Neotropics, at least 100 of which can infect humans (Forattini 2002). Between 1954 and 1988, 187 species of arboviruses were isolated by the Evandro Chagas Institute (Vasconcelos et al. 2001) indicating the magnitude of potential health problems. For example, dengue, Saint Louis encephalitis, and West Nile viruses, family Flaviviridae, are perhaps the best-studied arboviruses. Most arboviruses circulate in endemic/enzootic cycles in forest habitats and rarely, if ever, cause human disease (Contigiani et al. 2016). The genus Culex L., particularly the subgenus Melanoconion Theobald, contains important vectors of several arboviruses in the Neotropics (Mitchell et al. 1985, Sabattini et al. 1998, Ferro et al. 2003, Weaver et al. 2004, Turell et al. 2006). Culex ocossa Dyar and Knab, Cx. pedroi Sirivanakarn and Belkin, and Cx. taeniopus Dyar and Knab, are important vectors of Venezuelan equine encephalitis virus (VEEV), family Togaviridae, in Central and South America (Weaver et al. 2004, Aguilar et al. 2011). Several other species in the subgenus Melanoconion, including Cx. bastagarius Dyar and Knab, Cx. educator Dyar and Knab, Cx. intrincatus Brèthes, and Cx. delpontei Duret, have been found to be naturally infected with Bunyamwera virus (BUNV) (Peribunyaviridae) (Mitchell et al. 1985, Mitchell et al. 1987, Gallardo et al. 2019) and Cx. delpontei has been found naturally infected with VEEV (Mitchell 1985). In addition, VEEV and Western equine encephalitis virus (WEEV) have been isolated from specimens of Cx. ocossa and Cx. delpontei (Mitchell et al. 1985, Mitchell et al. 1987) and VEEV from Culex vomerifer Komp and Cx. gnomatos Sallum, Hutchings, Leila and Ferreira (Turell et al. 2006).

Because most arboviruses are RNA viruses, they have a very high mutation frequency (Holland and Domingo 1998), which allows them to alter host and vector ranges (Weaver and Barrett 2004). Several studies have shown that arboviruses, which are considered sylvatic are able to enter urban ecosystems and cause epidemics due to environmental and climate change (Contigiani et al. 2016). The future spread of arboviruses and parasitic diseases is projected to depend on the spread of effective vectors and their responses to current and future climate change. Recent environmental changes have created favorable habitat conditions for vector species in the Brazilian Amazon, for example, increasing the risk of malaria (Castro and Singer 2011). To this end, ecological niche modeling (ENM) and the estimation of vector species potential distributions in the Americas have received increasing attention (Larson et al. 2010, Altamiranda-Saavedra et al. 2017, Baak-Baak et al. 2017, Wiese et al. 2019, Figueroa et al. 2020, Sloyer et al. 2022).

Understanding how these processes and dynamics change is important because of public health implications. Therefore, our main objective was to estimate the current potential distribution of 9 species of Culex (Melanoconion) using an ENM approach and to assess their spatiotemporal variation under future climate change with health implications for the Neotropics.

Materials and Methods

Occurrence Data and Calibration Area

Presence records of 9 species of Culex (Melanoconion) were obtained from the literature, the Global Biodiversity Information Facility database (GBIF, http://www.gbif.org), previous field collections by ML, and unpublished records by colleagues (Supplementary Table S1). Given that morphological identification of females of Culex (Melanoconion) is problematic (Mureb-Sallum and Forattini 1996, Torres-Gutierrez and Mureb-Sallum 2015), records from the GBIF database were carefully reviewed to discard dubious identifications, considering the person who made the specific identification as well as the life stage reported. The study area (M) was based on the user-defined model calibration area (Barve et al. 2011), following Soberon and Peterson (2005) and Peterson et al. (2011) (i.e., the accessibility and mobility area based on the BAM diagram framework). This area was delineated using the Neotropical region proposed by Löwenberg-Neto (2014), which follows the biogeographical regionalization of Morrone (2014). This M was used to generate current and future species distribution models based on the BAM diagram approach (Soberón and Peterson 2005, Soberón et al. 2017).

Environmental Variables

We obtained 19 bioclimatic variables from the Chelsa v2.1 database (Karger et al. 2017) (https://chelsa-climate.org/), with a spatial resolution of 30 arc-s (~1 km) for both current (1979–2013) and future scenarios (2041–2070 and 2071–2100, see current and future species potential distribution estimations section for more detail). Before analyzing the species distributions and their future projections, we first selected the bioclimatic variables following an exploratory analysis using the “jackknife test” in Maxent v.3.4.4 (Phillips et al. 2006), with a random 70:30 split (70% for training and 30% for validation) via a bootstrap method using the total presence records for each species. Then we analyzed potential correlations between variables across the M area using a Spearman correlation test in R (R Development Core Team 2011). The final variable selection (Table 1) was made considering the species biology, variable importance, and contribution value according to the jackknife test, and correlation scores under 80% (i.e., r < 0.8) (Elith et al. 2006). In the end, we used a single set of 7 variables that best-represented climatic variation in the ecological space for each of the 9 species.

Current and Future Species Potential Distribution Estimations

We generated a species distribution model (SDM) based on an ENM approach (Elith and Leathwick 2009, Cobos et al. 2019) using the Maxent algorithm available in the KUENM package in R (Cobos et al. 2019). The model calibration process to estimate the currently suitable areas for each species was performed as follows. Before using KUENM, the species records were divided into calibration and testing data using the checkerboard 1 method via the ENMeval package implemented in the software Wallace 1.0.6.2 (Kass et al. 2018). Then, in KUENM, we created candidate models following Cobos et al. (2019) using the calibration occurrences datasets, 6 regularization multiplier values (RM) (0.5, 1.0, 1.5, 2.0, 2.5, 3.0), 7 combinations of feature classes (FC) (“l,” “q,” “p,” “qp,” “lp,” “lq,” “lqp”; linear = l, quadratic = q, and product = p), and the single set of environmental variables for each species (Table 1) with 10,000 background points. Candidate models were evaluated using the testing dataset and based on statistical significance (partial ROC, with 500 iterations and 50% data for bootstrapping), omission rate (OR), and model complexity according to Akaike information criterion corrected for small sample sizes (AICc) (Warren et al. 2010, Cobos et al. 2019). Best and final models were selected according to 3 criteria: (i) significant models, (ii) OR ≤ 15%, and (iii) models with delta AICc ≤ 2. Candidate models were created using the kuenm_cal function and evaluation, and best model selection was done using the kuenm_ceval function (Cobos et al. 2019). Finally, to estimate the future suitable areas (as a proxy of the future species potential distribution), we transferred the best current model of each species to the future general circulation model IPSL-CM6A-LR included in the CMIP6 from Chelsa v2.1 database (https://chelsa-climate.org/). This model was selected according to its capacity to represent the effects of ENSO (El Niño-Southern Oscillation) variability following Bellenger et al. (2014). Specifically, we used 2 shared socioeconomic pathway scenarios (ssp1-2.6, ssp5-8.5), with 2 time periods: 2041–2070 and 2071–2100. The least drastic future scenario, ssp1-2.6, shows an increase in global mean surface air temperature (GMST) relative to 1850–1899 ranging from 1.8 to 2.6 °C between 2041 and 2070 and 2.1–2.7 °C between 2071 and 2100; and atmospheric CO₂ concentration ranging from 417 to 459 ppm between 2041 and 2070 and 437–459 ppm between 2071 and 2100 (Boucher at al. 2020). The most drastic future scenario, ssp5-8.5, shows a higher increase in GMST relative to 1850–1899 ranging from 2.5 to 4.5 °C between 2041 and 2070 and 4.4–6.6 °C between 2071 and 2100, and atmospheric CO₂ concentration ranging from 500 to 722 ppm between 2041 and 2070 and 722–1,100 ppm between 2071 and 2100 (Boucher at al. 2020).

Models were produced not allowing extrapolation or clamping, so as not to predict beyond the range of environmental values available in the training range (Elith et al. 2010, Merow et al. 2013). Future models were created using the kuenm_mod function (Cobos et al. 2019). To identify extrapolation risks in model transferability due to no-analog climates (i.e., current and future climates), we performed a multivariate environmental similarity surface analysis (MESS, Elith et al. 2010, Phillips et al. 2017). The final model evaluations were based on statistical significance and OR using independent data (unpublished records from colleagues not used during model calibration).

Results

A total of 409 occurrence records (Supplementary Fig. 1) for the 9 Culex (Melanoconion) species were analyzed (Table 1). The calibration process generated a total of 84 models, with 42 candidate models for each species, most of which were statistically significant (p < 0.01) (Supplementary Table S2) (except 4 of Cx. pedroi). Applying the 3 evaluation criteria, few models met the criteria for the best model selection for each species (Table 1). The output of the final model selected for each species is shown in Fig. 1. Mean temperature of the warmest season around 27 °C, precipitation during the driest month around 50 mm, and precipitation during the warmest season over 200 mm contributed greatest to model performance of all species except Cx. bastagarius, Cx. delpontei and Cx. intrincatus (Supplementary Fig. 2). The environmental predictors contributing the greatest to model performance under current environmental conditions for Cx. bastagarius are the annual mean temperature between 26 and 30 °C (4.12% of contribution), the mean temperature of the warmest season around 30 °C (12.17% of contribution), and precipitation of the driest month over 60 mm (21.95% of contribution); for Cx. delpontei are temperature seasonality around 54.2% of contribution, mean temperature of the wettest season of 23 °C (21.3% of contribution), and precipitation of the warmest season over 200 mm (21.3% of contribution); for Cx. intrincatus are mean temperature of the warmest season of 26 °C (8.92% of contribution) and precipitation of the wettest and the warmest seasons, around 250 and 450 mm, respectively (0.06 and 0.02% of contribution, respectively). The areas predicted most suitable under current environmental conditions are shared by almost all species (Fig. 1C and E–I). These areas are located geographically in the southeastern Chaco Phytogeographic Province (PP), the northeastern Pampean PP, the southwestern area of the Parana Forest PP, and the southwestern Cerrado PP (Fig. 1C and E–I), except for Cx. bastagarius, Cx. delpontei, and Cx. intrincatus (Fig. 1A, B, and D). The suitability area of Cx. bastagarius is divided in 2 regions: the northern Pampean, southeastern Chaco, and southern Parana Forest PPs; and Napo, Imerí, and Chocó-Darién PPs (Fig. 1A). The suitability area of Cx. delpontei is located in the south-central Chaco PP, northeastern Pampean PP, and a small remnant in the central Monte PP (Fig. 1B). The suitability area of Cx. intrincatus covers an extensive area that almost reaches the extreme north and south of the Chaco PP, north of Pampean and Araucaria Forest PPs, and south of Parana Forest and Atlantic PPs (Fig. 1D). Culex ocossa and Cx. pedroi share suitability areas and are potentially distributed in the Imerí, Madeira, Napo, and Ucayali PPs, the eastern part of the Atlantic and Caatinga PPs, and in the Veracruzan and Yucatán Península PPs (Fig. 1E and F). The suitability areas of Cx. taeniopus are presented in a disjunct manner along the Amazon River and its tributaries in the Roraima, Madeira, Imerí, Napo, and Ucayali PPs, as well as in the coastal areas of the Guiana Lowlands, Guajira, Veracruzan, Yucatán, Pacific Lowlands, Cuban, and Hispaniola PPs (Fig. 1H).

In the near future (2041–2070) and under both scenarios of climatic change (ssp1-2.6 and ssp5-8.5), a slight and more marked enlargement of suitability areas can be noticed for Cx. bastagarius, Cx. delpontei, Cx. educator, Cx. ocossa, and Cx. taeniopus (Figs. 2–4, 6, and 9). For Cx. bastagarius, Cx. delpontei, Cx. educator, Cx. ocossa, and Cx. taeniopus, the suitability areas may be slightly reduced for the periods 2071–2100, given the conditions of the least dramatic scenario (ssp1-2.6) (Figs. 2–4, 6, and 9). Under the worst scenario (ssp5-8.5), a strong reduction of the suitability area is expected for 2071–2100, with a southward (and eastward in some cases) shift of these species distributions. The Chaco PP, apparently, would no longer be suitable for Cx. ocossa and Cx. taeniopus (Figs. 6 and 9) and only would have a very small suitable remnant toward the south suitable for Cx. bastagarius, Cx. delpontei, and Cx. educator (Figs. 2–4).

On the other hand, the Monte PP would present favorable environmental conditions for Cx. bastagarius and Cx. educator (Figs. 2 and 4). A first contraction followed by an increase in the suitability area of Cx. intrincatus is expected under the least drastic scenario (ssp1-2.6) (Fig. 5). The most extreme scenario (ssp5-8.5) shows a marked decreasing trend in the suitability area for Cx. intrincatus over time, reaching a maximum for the 2071–2100 period, with small remnants available at northeastern and southwestern parts of Pampean and Atlantic PPs, respectively (Fig. 5). A slight increase in the suitability area is expected for Cx. pedroi under both scenarios of climatic change over the time (Fig. 7). The suitability areas of Cx. pilosus and Cx. vaxus are not modified over time under the environmental conditions of the least drastic climatic change scenario (ssp1-2.6) and the period 2041–2070 of the worst conditions (spp5-8.5, Figs. 8 and 10). However, a marked contraction of the suitable areas of both species would be restricted for the extreme scenario south of the Monte PP, the entire Pampean, Araucaria and Atlantic Forest PPs, and in most of the Parana Forest and southeastern area of the Caatinga (Figs. 8 and 10). In relation to Cx. ocossa, the suitable area increases, in general, with a notable retraction in the Chaco PP in both scenarios for the 2041–2070 period (Fig. 6). For the most temporally distant period, a reduction in suitability is expected for both scenarios, with a stronger reduction for the worst conditions (spp5-8.5) and a marked shift southward into the Pampean and Monte PP (Fig. 6).

Some risk of nonanalogous climates (negative values shown in red) in much of the calibration area can be noted for all the species for the 2071–2100 period and worst conditions (spp5-8.5, Figs. 2–10). This MESS analysis indicates apparently novel climates, in the future outside the current range. These conditions may be found in the eastern slope of the Andes to the Atlantic coast under the most drastic and temporally distant scenario, except for the Monte, Pampean, Araucaria Forest, Atlantic, Parana, and the South of the Cerrado PPs (Figs. 2–10). These novel climate conditions would correspond to the mean temperature of the warmest season, the mean temperature of the wettest season, and to a lesser extent the maximum temperature of the warmest month. However, apparently, none of the species would occur in those areas.

Discussion

It is expected that the Culex (Melanoconion) species here assessed are potentially distributed based on the availability of suitable areas. These are consistent with what has been proposed for the subgenus Melanoconion, which includes 177 species (Harbach 2023) widely distributed from southern Canada southward through most of Brazil, Paraguay, Peru, on the western side of the Andes to central Argentina (Hunter et al. 2015, Torres-Gutierrez and Mureb-Sallum 2015). However, most species within the group appear to be adapted to tropical and subtropical environments such as the Amazon region and other forested environments in the Caribbean, central America, and northern and western South America (Forattini and Mureb-Sallum 1992, Hutchings et al. 2011, Rodrigues de Sá et al. 2022). Therefore, based on our modeling results, it will be informative to attempt to collect these species in those geographic areas where they have not been recorded, especially in areas of human development, such as cities and surrounding rural areas, where other species of the subgenus have been previously documented (Forattini and Mureb-Sallum 1987, Rossi et al. 2002, Pires et al. 2009, Stein et al. 2013, Burkett-Cadena 2022, Rodrigues de Sá et al. 2022, Sloyer et al. 2022). As pointed before, several arboviruses have expanded from a sylvatic cycle to include an urban cycle (Contigiani et al. 2016). However, it is important to note that although our models were calibrated using adult mosquito presence data using an ENM approximation rather than a mechanistic model, multiple conditions can contribute to their occurrence, and adaptation or toleration to anthropogenic environments and even blood feeding on humans. Therefore, it is expected that future arbovirus distributions may depend on the future distributions of their vectors, as shown in our results. For example, current environmental changes have increased the risk of malaria by creating suitable conditions for Anopheles species in the Brazilian Amazon (Castro and Singer 2011).

Despite this scenario, most of the studies about Melanoconion have focused on taxonomy (Forattini and Mureb-Sallum 1985, 1992, 1993, Mureb-Sallum and Forattini 1996), phylogenetic relationships (Navarro and Weaver 2004, Torres-Gutierrez et al. 2018, Bangher 2020), and natural infection (Mitchell et al. 1985, 1987, Turell et al. 2006). Except for the work of Sloyer et al. (2022), none have focused on the potential geographical distribution and conditions that could be suitable for these species’ biology. Our results show that suitable conditions for these species are associated with environmental predictors, such as mean temperature of the warmest season between 26 and 30 °C (except for Cx. delpontei), mean temperature of the wettest season, precipitation of the warmest season over 200 mm, and precipitation of the driest month. In this sense, suitable areas for most species under both scenarios of climatic change would be larger, but with a marked contraction of suitable areas and a shift to higher latitudes and coastal zones under the worst scenario for the medium–long future (2071–2100). Although our models were calibrated based on adult mosquito occurrence, a plausible explanation may be that increases in temperature and atmospheric CO₂ concentration would allow species to retain distributions within optimal ranges for their development and survival (immature and adults). This is consistent with the relationship between mosquito development and temperature proposed by Paaijmans et al. (2010) being one of the key interactions to interpret the current and future distribution of vector-borne diseases. This highlights the importance of the upper tolerance thermal limits for a species in the context of global warming (Mordecai et al. 2019). The main environmental predictor for Melanoconion species distributions was related to temperature, in agreement with Sloyer et al. (2022) who found that mean annual temperature contributed the greatest to model performance to predict the potential distribution of Culex (Mel.) cedecei (Stone and Hair) in Florida. Other authors also have reported temperature bioclimatic variables as important factors to predict the geographic distribution of other Culex species (Masuoka et al. 2010, Gardner et al. 2012, Figueroa et al. 2016, 2020) as well as Aedes species in suburban/rural areas (Wiese et al. 2019).

The potential suitability areas for the species seem to be stable in the near future (2030), but a contraction is anticipated following year 2070, a pattern that is also predicted for other mosquito species like Aedes aegypti (L.) and Nyssorhynchus darlingi (Root) (Khormi and Kumar 2014, Laporta et al. 2015). These contractions could be related to expected increases in the concentration of atmospheric CO₂ that would reduce cloud cover and therefore precipitation (Vila-Guerau de Arellano et al. 2012) and that would ultimately increase the heat stress at high temperatures due to the loss of humid areas. Rising temperatures (>4.4 °C) and atmospheric CO₂ concentration (>722 ppm) were herein predicted to have a strong impact on the distribution of most Melanoconion species, with a contraction of their distribution and marked displacement toward coastal areas. This trend is also consistent with that proposed for Anophelinae distributions, which seem to be driven by precipitation, particularly by the length of the dry season (Fu et al. 2013). Likewise, for Ny. darling, the projected climate change for 2070 would significantly reduce suitable habitat and favor the expansion of suitable areas for the albitarsis Complex in South America (Laporta et al. 2015).

Arbovirus transmission systems often represent complex interactions among multiple species (e.g., vectors, hosts, pathogens) (Peterson 2006). Therefore, it is important to treat the results of ENMs with caution when trying to predict species ranges and their changes in distribution within the context of climate change. We need to understand that these climate change scenarios are not facts, but rather plausible scenarios that depend on the socioeconomic paths that the planet takes. Thus, if we keep this in mind, ENMs have considerable potential for information that may be applied to the problem of geographic characterization of vector-borne disease systems, especially in climate change scenarios designed to predict possibly high-risk areas (Peterson 2006). In this sense, many challenges remain for improving the forecasting each species’ response to global climate change, as suggested by Peterson (2006). Prudent choice of environmental datasets and the assumptions, algorithms, and parameterizations of different methods are some of the most critical challenges to utilizing ENMs for these purposes (Peterson et al. 2001, Peterson et al. 20011, Beaumont et al. 2005). For example, the risk of nonanalogous climates and extrapolating (when it is made) to those areas/times representing environmental values outside the environmental range of the study area is an important issue. Therefore, it is important to interpret prediction results with caution, as ENMs do not necessarily identify directly physiological mechanistic processes between species and the novel environments in the future (Araújo and Rahbek 2006, Peterson et al. 2011).

It is crucial to understand the assumptions made in achieving model transferability to other landscapes or to other time periods to make appropriate interpretations (Peterson et al. 2011). Nonetheless, these approaches offer attractive possibilities for understanding the ecological distributions of species in a vector-borne disease context (Peterson and Shaw 2003). In this sense, despite the assumptions and limitations of the models developed here, the results could be used for future monitoring of Culex (Melanoconion) species populations and BUNV, VEEV, and WEEV, which would greatly contribute to the synthesis of region-specific surveillance programs that include other vectors and pathogen species of public health interest.

Acknowledgments

We are grateful to GC Rossi (Centro de Estudios Parasitológicos y de Vectores (CEPAVE), CCT CONICET La Plata) and DN Bangher (Universidad Nacional del Nordeste, CONICET) for providing records of Culex (Melanoconion) specimens from their unpublished databases. We thank ND Burkett-Cadena, PhD (University of Florida) for their useful comments on the article. M.L. is member of the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET).

Author Contributions

Magdalena Laurito (Conceptualization [Equal], Data curation [Lead], Formal analysis [Equal], Investigation [Equal], Methodology [Equal], Supervision [Equal], Validation [Equal], Visualization [Equal], Writing—original draft [Lead], Writing—review & editing [Equal]), and Andrés Arias-Alzate (Conceptualization [Equal], Formal analysis [Equal], Investigation [Equal], Methodology [Equal], Supervision [Equal], Validation [Equal], Visualization [Equal], Writing—review & editing [Equal])

References