Predicting climate effects on aquatic true bugs in a tropical biodiversity hotspot

Climate change is a matter of worldwide concern with severe predicted impacts on biodiversity. Here, we analysed the potential impacts of current and future climates on aquatic true bugs (Heteroptera) in relation to their distribution patterns and ecological preferences (based on a database generated from existing literature references and field collections). We considered the traits as ‘species thermal range’ and ‘emergence period’ to evaluate the future climate change impacts on the distributions of aquatic true bugs in the riverine regions of a tropical biodiversity hotspot, the Western Ghats of India. We used Species Distribution Models (SDMs) to evaluate the potential impacts of climate change on the distributions of aquatic true bugs. We modelled the distributions of twenty-six species of aquatic true bugs using different modelling tools through a carefully examined set of occurrence records to generate potential present distributions and to project these distributions into future scenarios of climate change. We observed increasing/decreasing range sizes of the species in the current and future scenarios. We found losses and increases of species' ranges in some regions, but not much variation in species richness. Similarly, no significant effect was observed in the distribution ranges for species with different duration of emergence period and thermal range in current and future climatic scenarios. Losses and gains in species richness would be concentrated in the mountainous area of the Western Ghats, whereas loss of species and the bigger difference between current and future richness will occur in the adjacent lowlands and towards central regions, including the network of protected areas of the Western Ghats. These areas are critical to buffer regional species loss in the future. Given the importance of aquatic true bugs as bioindicators and biological control agents, monitoring their range shifts should be routinely addressed in the conservation contexts of the Western Ghats.


Introduction
Climate change is projected to become one of the most predominant drivers of biodiversity change by the end of this century (IPCC 2013). In particular, freshwater biodiversity is increasingly threatened by human activities, resulting in rapid and major declines in the diversity of aquatic ecosystems, habitats, and species over a wide range of scales from global to local (e.g., Dudgeon et al. 2006;Reid et al. 2019;Albert et al. 2020). Climate warming will profoundly impact the distribution, abundance, and composition of aquatic insect species and communities because it will affect species through almost all life-history traits, including emergence, growth rate, and voltinism (Bale et al. 2002;Villalpando et al. 2009). Nevertheless, the knowledge of climate change impacts on aquatic insects is poorly developed, particularly based on findings from biodiversity hotspots (Hering et al. 2009;Sundar et al. 2020).
It is expected that global warming will affect aquatic insects in many ways, ranging from changes in their distributions to local, regional, and global extinctions . However, our capability to predict the effect of climate change on different aquatic insect groups, based on field data, is variable. Considering the significant threat posed by mosquitoes to human health, the relationships between climate change and these aquatic insects are perhaps the best studied, in both experimental and modelling work, and current knowledge suggests that the ranges of mosquito-borne diseases could expand dramatically in response to climate change (Reiter 2001). Despite the increasing availability of digitized information on biodiversity data and species occurrence records (Graham et al. 2004;Newbold 2010;Pyke and Ehrlich 2010), our knowledge about most aquatic insect species' geographic distribution is woefully incomplete to allow accurate analysis of climate change effects. This is particularly true for aquatic insects in tropical biodiversity hotspots around the world. Up to date, only one study has addressed the potential impact of climate change on stoneflies in a tropical biodiversity hotspot, i.e., the Atlantic Forest of Brazil (Silva et al. 2019). Lack of studies causes concerns because some groups of aquatic insects, such as the cold stenothermic species, are expected to be highly threatened by climate change (Hering et al. 2009). Moreover, there is recent evidence that climate change could affect each species to different degrees and that the protected areas network would not be effective to conserve species under climate change mainly because it would not cover future species distribution of aquatic insects (Sánchez-Fernández et al. 2013;Rosso et al. 2018).
In this study, we attempted to predict the distribution patterns of aquatic true bugs (Heteroptera) based on Species Distribution Modelling (SDM) tools (Elith and Leathwick 2009;Elith et al. 2010). We compared distribution changes between current and future climatic conditions in species with different emergence times and thermal ranges. Specifically, we expected that species with short emergence periods would be particularly sensitive to alterations in temperature (Kotiaho et al. 2015). Concerning the thermal range, species with narrow thermal range are expected to be more strongly threatened by climate warming compared to the species with wide thermal range (Heino et al. 2009;Silva et al. 2019). The generated database of aquatic true bugs we used in this study is from a poorly-studied tropical biodiversity hotspot, the Western Ghats of peninsular India.

Study area
The Western Ghats (latitude: 8-21 °N), a UNESCO World Heritage Site, among the 36 global biodiversity hotspots in the world, occupies 5th position in the economic potential of its biological resources (Shameer et al. 2016). Western Ghats' major vegetation types are moist deciduous, dry deciduous, dry shrub vegetation, semi-evergreen, evergreen, the shola and grasslands forests (Singh et al. 2002). It is a chain of mountain ranges with approximately 140,000 square kilometers (a length of 1600 km) and extends from the west coast of peninsular India from the River Tapti in the north to Kanyakumari in the south. Some important rivers originate from the mountains of Western Ghats, such as Godavari, Cauvery, Krishna, Thamiraparani, and Tungabhadra (Anbalagan et al. 2012;Shameer et al. 2016).
The Western Ghats' climate is primarily influenced by the South-West monsoon between June and September when it receives the maximum rain (> 80%). The annual rainfall varies from 2350 mm in the north to 7450 mm in the south, even though the states Karnataka, Tamil Nadu, and Kerala that comprise the Western Ghats, receive significant rainfall from the North-East monsoon from October to December (Anbalagan et al. 2012). The Western Ghats experience a tropical climate that is warm and humid during most of the year. The average annual temperature for the whole Western Ghats is around 15 °C. During winter, frost is common in some parts, and the temperature drops to freezing point. Mean temperature ranges from 20 °C in the south to 24 °C in the north. The coldest periods in the South Western Ghats are also the wettest and receive the most precipitation (Dahanukar et al. 2004).
Large-scale land use changes leading to deforestation have contributed to an increase in mean temperature by 0.5 °C and a decline in the number of rainy days (Ramachandra and Setturu 2019), which necessitates prudent landscape management strategies involving all stakeholders for conservation of the ecologically fragile Western Ghats. The primary forests have been transformed into other land-uses in the past four decades due to commercial establishments, hydroelectric projects, industries, and monoculture plantations (Ramachandra and Setturu 2019). The main threats affecting freshwater biodiversity in the Western Ghats include pollution, biological resource use, residential and commercial development, dams and other natural system modifications, alien invasive species, agriculture, aquaculture, energy production, and mining. Stream vegetation and riparian forests have also been subjected to similar threats with direct impacts on biodiversity and indirect effects through increased sediment loads, erosion, flash floods, loss of habitats such as stagnant pools, inconsistent flow, and disappearance of primary and secondary streams (Molur et al. 2011).

Species occurrence data
First, we gathered occurrence data on aquatic true bugs in the Western Ghats from all available published literature sources (Appendix), and the first author's (S.S.) field experience data on the distributions and ecological preferences of the aquatic bugs in the Western Ghats river basins were pooled to generate baseline data. Although this database is subjected to the prejudice of lack of uniformly designed spatial and temporal surveys, a genuine attempt was made to arrive at a tentative to standardize the information on the distribution of aquatic true bugs in the Western Ghats. In our study, six species are endemic from Western Ghats, and then our dataset and modelling strategy do not cause any bias in projecting their potential distribution across India. In the case of species that are not restricted to Western Ghats, we recognize that our approach of using only data from this region could increase the risk of artificially assigning low habitat suitability to environmental conditions in which the species could be actually living (Raes 2012). However, as Western Ghats comprises a huge environmental variation of multiple ecoregions, we believe that our dataset is adequate to project area for most parts of India.
A potential source of criticism to the amassed occurrence dataset concerns the low availability of individual species records of occurrences (Table 1). A second criticism may concern the number of modelled species we are focusing on in this study. Mainly, this occurs because, when compared to other biological groups, insect species are generally less charismatic and the species charisma is an important driver in conservation (Mammola et al. 2020). Consequently, they are more strongly affected by the several data-related shortfalls that usually affect data for biological groups (Cardoso et al. 2011;Hortal et al. 2015). Despite potential criticisms regarding how comprehensive and extensive this dataset is, it is currently the best occurrence dataset for aquatic true bugs available in the Western Ghats. Even though a lack of large numbers of occurrences may affect modelling results, several past studies were able to produce a preliminary assessment of a given species' distribution ranges with a low number of records. We used only reliably georeferenced records of aquatic true bug species in this study. Records with strikingly incorrect georeferences (e.g., points on the sea, inverted latitude and longitude, inverted signs) were excluded from all analyses.
Model per for mance depends on both sample size and species' prevalence, being the fraction of the study area occupied by the species. Previous studies involving the SDM approach were able to produce reliable species distribution prediction with a minimum of one occurrence (De Siqueira et al. 2009). Yet, a higher number of records (e.g., 10 to 50 occurrences per species) would allow us to obtain more reliable range predictions for the modelled species. Therefore, we focused only on species with at least 10 spatially unique occurrences to allow a minimally reliable distribution range to be modelled at the grid cell-size resolution we used (see below). Following the criteria mentioned above, we modelled the final set of 26 aquatic true bug species depicted in Table 1.

Predictor variable processing and modelling procedures
We obtained the 19 variables for both current and future climate conditions from Worldclim 1.4 database (Hijmans et al. 2005), considering a 2.5 arcmin (~ 0.041º or 4 km at the equator. We standardized all variables from all scenarios to have their means equal to zero and unit variances. For the current scenario, we used Principal Component Analysis (PCA) to generate orthogonal/independent variables to be used as new variables to predict the species' distribution. Using this method, we avoided model overfitting (Jiménez-Valverde and Lobo 2011) and decreased both the number of variables and variable collinearity in our models (Dormann et al. 2013). The first seven principal components (PCs) accounted for ~ 97% of the original climate variables' variation. The contribution of each climatic variable to each one of the seven PCs is shown in Table 2.
The variables for the future scenarios were available from 17 atmosphere-ocean global circulation models (AOGCMs) available from Worldclim, considering the  0, BCC-CSM1-1, CCSM4, CNRMCM5,  GFDL-CM3, GISS-E2-R, HadGEM2-AO, HadGEM2-CC,  HadGEM2-ES, INMCM4, IPSL-CM5A-LR, MIROC-ESM-CHEM, MIROC-ESM, MIROC5, MPI-ESM-LR,  MRICGCM3, NorESM1-M). Under these RCP conditions, the global human population is expected to grow up to 9-12 billion people continuously. No significant socio-politiceconomical changes are expected to reduce the emission of greenhouse gases significantly. According to the reports from the Intergovernmental Panel on Climate Change (IPCC), the average temperatures are expected to rise 3.7 °C, with variation from 2.6 to 4.8 °C). All AOGCMs had the same resolution than the current scenario and were also standardized to have the average equal to zero and unit variances. Then, we also applied a PCA to each of our 17 future scenarios by projecting the linear coefficients obtained for the current scenario on each one of them to create a dependency of each future scenario with the current one. Later, we applied a PCA to each one of the 17 AOGCMs, and also selected the first seven PCs to be used as predictor variables of these species in each future.
To predict the distribution patterns of aquatic true bugs, we used Species Distribution Modelling (SDM) methods (Elith and Leathwick 2009;Elith et al. 2010). These methods are based on correlations from climate change models with the known current occurrence records to project onto the geographic space climatically similar areas with no available distributional information of the species. These methods are being systematically used to detect potential new areas of occurrence for insect species (Almeida et al. 2010), aiming to delimit and support active conservation actions (Nóbrega and De Marco 2011), predict suitable areas for invasive species (Silva et al. 2014;, and to predict climate effects upon biodiversity for both the past (Peres et al. 2015) and future (Martins et al. 2015) climatic scenarios. These methods are based on the Biotic-Abiotic-Movement (BAM) diagram proposed by Soberón and Peterson (2005) and Soberón (2007). According to the theory involving the BAM diagram, the realized niche of species is related to the locations where the abiotic and biotic conditions would allow its occurrence and are accessible for the species.
Since many of the species we modelled had very few occurrence records, we used the Maximum Entropy method (abbreviated as "MaxEnt" hereafter; Phillips et al. 2006;Phillips and Dudík 2008) as our modelling method. The MaxEnt is a machine-learning presence/pseudo-absence method that is reliable when the number of occurrences available for modelling a species distribution is limited (Hernandez et al. 2006;Pearson et al. 2007). In our modelling procedures with MaxEnt, we only used both linear and quadratic features to simplify the biological interpretation of distribution models (Elith et al. 2011;Anderson and Raza 2010).
We bootstrapped our species occurrence dataset into two subsets in our modelling procedures, containing 70% of the occurrences for model training and the remaining 30% as testing datasets. We also used pseudo-absences along the geographic space by first establishing a multivariate bioclimatic environmental space. Then, we allocated the pseudoabsences in geographic areas that were climatically outside of such environmental space as previously done in other studies (VanDerWal et al. 2009;Lobo and Tognelli 2011). This allows MaxEnt to have higher discriminatory and explanatory capacities and a lower probability of climatic conditions impossible to be accessed by the species to affect the modelled species' response curves (Barve et al. 2011). Finally, we also restricted the available climatic information to train the models with the ecoregions shape file obtained from World Wildlife Fund website (https ://www.world wildl ife.org/biome s) to the ecoregions within the Western Ghats' known occurrences of the modelled species. This approach avoids model over prediction and unreliable predictions for the species distribution range (VanDerWal et al. 2009;Giovanelli et al. 2010;Acevedo et al. 2017).
In the modelling, we considered the "lowest presence training" (LPT hereon; Pearson et al. 2007) as the threshold value to cut the suitability matrices obtained for all modelled species in all scenarios into presence/absence maps. In ecological terms, the LPT threshold identifies pixels that are predicted as minimally suitable for the species' occurrences as its actual observed records. Since the species' realized ecological niches may be underestimated, given the small number of occurrences, we intended to estimate a distribution range for these species that could be near their true distribution range by using this threshold.
We considered the Jaccard similarity index, as proposed by Leroy et al. (2018), to be used as the evaluation metric in our study. This metric varies from 0 to 1, where values equal to or higher than 0.7 are considered acceptable. One of the main advantages of using this metric is that it does not need pseudo-absences to evaluate the produced models, being able to produce more reliable predictions than those that consider pseudo-absences to generate contingency tables for model evaluation [e.g. area under the curve (Fielding and Bell 1997) or true skill statistics (Allouche et al. 2006)]. To generate the final ensemble distribution for each one of the species we modelled, we used a weighted mean based on the Jaccard similarity index to ensemble the final distribution range for all species in all scenarios. We ran all analyses with the ENMTML R package by Andrade et al. (2020).
We follow general classification proposed by Hering et al. (2009) to classify ecological and distributional traits, i.e. species thermal range (narrow < 20 °C) and duration of emergence period (short < 50 days) for the occurrences from the Western Ghats. Unfortunately, there is little information available to figure out if aquatic true bug species that occur within a narrow thermal range are necessarily thermal specialists in terms of their fundamental thermal niche. So, to avoid this problem, we classified them as "species with narrow thermal range" and "species with wide thermal range". The species thermal range was defined according to the temperature range in which aquatic true bugs occur in nature because no physiological information was available. We checked whether the species with differences in their (1) duration of emergence period (Short < 50 days) and (2) thermal range (< 20 °C) (Table 1) showed any differences in distribution ranges in current and future scenarios. Therefore, we used Factorial Repeated Measures ANOVAs using the species' analyzed traits and climatic scenario as predictor variables and their distribution range size in each scenario as the resulting variable. Since we are using the very same occurrences to generate both current and future predictions of our insect species, the predictions are dependent one another, what requires the use of a repeated measures ANOVA analysis. We also used a Pearson correlation to check if the species' ranges in both current and future scenarios were correlated. Finally, we also computed the species richness and evaluated the species richness change in both scenarios, by summing the final modelled ranges for all species in specific each scenario.

Results
In general, our models reached highly acceptable Jaccard values (0.947 ± 0.051; average ± standard deviation), a result that indicates the predicted ranges for all 26 species we modelled here were adequate. Several species reached Jaccard values equal to 1, whereas the one which reached the lowest value was D. rusticus (Table 1). The overall contributions of each variable to the PCs we used as new environmental variables in each one of the climatic scenarios we modelled may be checked in Table 2.
After comparing the distributions of the species in both current and future scenarios (Figs. 1 and 2), we observed that there were practically no differences in the species distributing range size in each scenario. Yet, some patterns were noticeable. For instance, most species (n = 16) showed some range increase, while very few of them showed decreased ranges when current and future ranges were compared (Fig. 3a). Among these species, A. hyalinipennis (~ + 86%), H. indicus (~ + 70%), L. grossus (~ + 41%), and D. rusticus (~ + 33%) were the species that showed ranges increasing the most from the current to the future climatic scenario. However, some species (n = 9) showed distribution ranges decreasing from the current to the future scenarios (e.g. A. barbatus, A. s. sardeus, C. asiaticus, C. productus, D. molestus, H. rotundatus, L. maculatus, M. communis, and M. s. scutellaris). Among these last species, H. rotundatus (~ -80%), C. asiaticus (~ 31%), A. exiguus (~ − 30%), and L. maculatus (~ − 20%) were the species that were predicted to exhibit the highest range loss. In general, all species had their current range without much difference in relation to their future range, although many species had a narrow distribution range in both scenarios, while very few species had a broad distribution range in both scenarios (Fig. 3b). We did not observe any effect of the species duration of the emergence period (F1,23 = 1.119; p = 0.301) and thermal range (F1,23 = 1.119; p = 0.301) on the species ranges in both current and future scenarios.
Finally, considering the modelled species richness in both current (Fig. 4a) and future (Fig. 4b) scenarios, we did not observe a strong difference between the scenarios. Yet, when we subtracted the current species richness from the future species richness, we observed that there was loss of species richness in some regions, whereas there were increases in the modelled species richness in others (Fig. 4c).

Discussion
Our study adds evidence from the Western Ghats that some aquatic insects will face significant threats of declining habitat suitability for species due to climate change in the coming years. However, our findings did not support the idea that their response to climate change can be easily predicted based on ecological traits, such as species thermal range and emergence period. Our results also indicated that the altitudinal and latitudinal ranges of the Western Ghats are key characteristics to retain the balance between gains and loss in the distribution of aquatic true bugs.
The emergence period of aquatic true bug species largely depends on temperature and can vary from less than 50 days in some species at high temperatures to more than one year under colder temperature conditions (Saulich and Musolin 2007). It is expected that species with short emergence periods are particularly sensitive to temperature rise, leading to accelerated development, more generations per year, and, consequently, changes in their distribution patterns. Contrary to this expectation, we did not find an effect of duration emergence time in both current and future scenarios. A possible explanation could be related to the intensity of climate change in the region that may not be strong enough to elicit clear responses of species with different emergence periods to warming. However, considering that, in the 2030s, the mean annual temperatures of the Western Ghats region are projected to rise from 26. 8 ± 0.4 °C to 27.5 ± 0.4 °C (MoE-FCC 2010), we believe that the lack of relationship can also be attributed to the high variability in the ecological requirements and behavioral responses of various taxa to climate warming, which are not represented by simple categories of emergence periods we used in this study.
The thermal range of aquatic true bugs did not affect our findings. It is a surprising finding given that there are many species with narrow thermal range in the Western Ghats and some species with narrow thermal range could have wider thermal tolerances and the lack of differences in distribution patterns observed in this study could be affected by this limitation. The topographic, altitudinal, and latitudinal variation may promote opportunities to some species to change their distribution, resulting in no net effect on their distributional range size. Moreover, aquatic true bug species in the Western Ghats may have experienced strong temperature variation along the latitudinal and altitudinal gradient (increasing and decreasing cycles) in the Quaternary, which may generate resilient species with wide thermal range. Variable spatio-temporal patterns of precipitation (and the resultant expansion and contraction of suitable habitats, droughts, floods and so on) are also essential to promote opportunities for species resilience and co-existence through niche partitioning (Brown et al. 2020). In this context, new studies about the effect of climate change on aquatic insects' traits should address the potential impacts of both temperature and precipitation. Our results showed that climate warming would decrease the distribution ranges of nine species of aquatic true bugs. All species predicted to experience range loss in the coming years would be restricted to the mountainous part of the Western Ghats, and some species, such as H. rotundatus, C. asiaticus, would be limited to the extreme south of the region (close to the Sri Lanka Bridge). Previous studies evaluating the effects of climate change on the distribution of aquatic insects (e.g. Silva et al. 2019;Tierno de Figueroa et al. 2010;Domisch et al. 2013) have already shown that some species will be impacted. This trend is even more concerning if we think about the predicted hydrological changes expected to occur worldwide (Rodell et al. 2018), where many catchments around the planet will exhibit drier conditions (Rolls et al. 2018). Although abundance changes at the community level may not necessarily occur (Jourdan et al. 2018), sensitive species may suffer range size decreases in some areas, especially those surrounded by intensified agriculture. The changes in the distribution ranges for our modelled species were not as evident as observed for other species in previous studies involving SDMs (e.g. Tierno de Figueroa et al. 2010;Domisch et al. 2013;Shah et al. 2014). Such differences may be related to the fact that the species we modelled here are more adapted to the climate conditions near the Tropics than the species modelled in previous studies conducted in subtropical and temperate areas. Therefore, these species may be considered more climate generalists, and such a feature may indicate a better ability to cope with climate change. Yet, considering the ecological and economic importance of aquatic true bugs as bioindicators and biological control agents (Schaefer and Panizzi 2000), careful monitoring of their range shifts can pinpoint the species requiring particular conservation attention.
Losses and gains in species richness would be concentrated in the mountainous area of Western Ghats, whereas loss of species and the bigger difference between current and future richness will occur in the adjacent lowlands and towards central regions. This finding suggests that the Western Ghats' topography and latitudinal range can buffer regional species loss in the future. Furthermore, in India, anthropogenic pressures on biodiversity are concentrated in the lowlands, including near big cities (e.g., Mumbai, Pune), which will exacerbate the loss of biodiversity in the adjacent lowlands of the Western Ghats. Besides, the network of protected areas and complementary conservation strategies should include a range of lowlands and uplands in the Western Ghats, particularly the transitional areas between uplands and lowlands where considerable changes of species distributions are expected. Most protected areas are located in the mountainous regions of the Western Ghats, and they cover a sizeable latitudinal variation in the region. However, little information about the effective conservation of aquatic biodiversity is available in this region.
Despite of the results we obtained in this study, it is of utmost important to stress that SDMs are not silver bullets. These methods are being significantly used to determine distribution range hypotheses concerning several different objectives, such as: (1) delimitation of conservation areas (e.g., Nóbrega and De Marco 2011); (2) determine the impacts of invasive species and whether their niches change during the invasion process (e.g., Faleiro et al. 2015;Montalva et al. 2017;Silva et al. 2016Silva et al. , 2019; (3) evaluate the potential impacts of past and future climatic changes on biodiversity (e.g., Peres et al. 2015;Carvalho et al. 2017); and (4) pinpoint interesting areas for future sampling and to indicate species under-sampled regions. Still, considering the biotic-Abiotic-Movement (aka BAM) diagram (Soberón and Peterson 2005;Soberón 2007), these methods are mere simplification of the reality, since important ecological processes, such as migration and biotic interactions are not fully accessed, while the importance of abiotic processes may be overestimated. Considering the technology and theoretical background we know by now, we are not able to solve and include and all of these ecological processes into our models. Although our models are wrong, they are still useful (Box 1979), and with our results we provided a first glimpse on how the distribution range of these modelled species may be now and in near future. We also hope that our research serves to incentivize future research to be produced in this area of the globe, with these and other insect species, in order to fight the Wallacean shortfall that severely affect our knowledge on insect species (Hortal et al. 2015;Cardoso et al. 2011).
To conclude, our current observations suggested that climate change is a potentially significant threat to aquatic true bugs, although habitat fragmentation and degradation may have larger effects than climate change on these species' distributions in the tropics. Also, a considerable portion of aquatic true bugs may be driven to local or regional extinction in India because of rapid changes in temperatures or extreme climate events, even if the possibility of altitudinal and latitudinal range shifts cannot be fully neglected. We hope that our findings spur additional studies on the status, diversity, and distributions of aquatic true bugs in the tropics.