In the neotropics wide areas are planted in grasses, being Urochloa spp. P. Beauv. (syn. Brachiaria spp.) the most extensive forage monoculture (Ghimire et al., 2015; Worthington et al., 2021). Its economic impact is estimated at USD12.4 million in Mexico, Central America, Colombia, and Brazil, the largest contribution comes from U. brizantha (Hochst. ex A. Rich.) R.D. Webster cv. Marandu in Brazil with USD 6.3 million (White et al., 2013). The major biotic stress affecting forage production and its quality in this region is caused by spittlebugs (Hemiptera: Cercopidae). A large group of species causes severe damage in susceptible grasslands (Cardona et al., 2004) with economic losses estimated at USD 840–2,100 million per year in all host crops (Thompson, 2004).

Although spittlebugs are found in most terrestrial ecosystems, the tropics are the most diverse ecozone harboring 70% of known species (Thompson, 2004; Dietrich, 2009). In the Neotropics, species are reported from the southeastern United States to northern Argentina (Peck and Thompson, 2008). Different spittlebug species coincide in each country. The main species that occur in Brazil are from the genus Mahanarva (Distant, 1909), Notozulia (Berg, 1879) and Deois (Fennah, 1949) (Resende et al., 2012). In Mexico, the species Aeneolamia albofasciata (Lallemand, 1939), A. contigua (Walker, 1851), and A. postica (Walker, 1858) are major pests of sugarcane and grasses (Cardona et al., 2004; Thompson and León González, 2005; Parada Domínguez et al., 2019). Whereas in Colombia the predominant species are A. varia (Fabricius, 1787), A. lepidior (Fowler, 1897), A. reducta (Lallemand, 1924), Prosapia simulans (Walker, 1858), Zulia carbonaria (Lallemand, 1924), and Z. pubescens (Fabricius, 1803) (Peck, 2001).

Climate change can modify the distribution of species by expanding their presence to new locations and disappearing from previously suitable areas (Hughes, 2000). Anthropic movement, land-use change, environmental degradation (e.g., habitat loss and fragmentation) and biotic interactions (e.g., competition, species introduction, and plant host distribution) produced by the on-going climate change are factors that influence this distribution (Wagner et al., 2021). Insects are well-known for being particularly susceptible to environmental changes of temperature, humidity, radiation, and resource availability driven by those factors (Larson et al., 2019). Processes that homogenize and simplify the landscapes as extensive agriculture, allow the growth of pests over native species (Cardoso et al., 2020). Several studies in recent years have warned about the decline of insect populations to extinction caused by changes in the seasonality and, consequently, in their life cycles. This reduction in the populations has great impact over the ecosystems as the loss of abundance and richness of species continue to occur (Hallmann et al., 2017; Goulson, 2019; Halsch et al., 2021).

Despite insect pest outbreaks are expected for the short term (Heeb et al., 2019; Liu and Shi, 2020), its severity may not be evenly increased due to the narrow environmental niche requirements, physiological tolerances of insects, and differential effects of climate variables on their life cycle (Lehmann et al., 2020). Previous models show an increase in suitable areas for pest species in Europe, e.g., Helicoverpa zea (Lepidoptera: Noctuidae), Aleurocanthus spiniferus (Hemiptera: Aleyrodidae), under climate change scenarios (Grünig et al., 2020). Thus, characterizing the effect of climate change in Colombian spittlebugs geographic distribution and identifying niches where these species would become key pests is important in the transition to more sustainable livestock systems.

In this context, ecological niche models (ENMs) provide an approximation to estimate potential geographical zones with environmental conditions that a species requires to maintain its populations (Peterson et al., 2011). This tool is widely used in insect pest management programs to anticipate unknown distributional areas, geographic potential of invasive species, and response to changing environmental conditions (Peterson and Soberón, 2012). ENMs can be built based on occurrence data (inductive or correlative niche models; Elith and Leathwick, 2009) or based on physiological data [deductive or mechanistic niche models; (Kearney and Porter, 2009)]. For spittlebugs associated with grasses, we identified only two studies focused on changes in suitability of geographical areas under climatic change scenarios. The first, based on physiological data of Mahanarva spectabilis (Distant) (Fonseca et al., 2016), and the second, based on occurrence data of four Mahanarva species (Schöbel and Carvalho, 2020).

This paper responds to the need to know whether A. varia, A. lepidior, A. reducta, P. simulans, Z. carbonaria, and Z. pubescens are potential key pests in new sites under climate change scenarios that consider the impact of human activities. Hence, spittlebug ENMs contribute to the development of adaptation strategies for tropical America climate-smart perennial grasslands, and sugarcane production, by addressing the need for shift toward more sustainable pest management practices (Macfadyen et al., 2018). For instance, adoption of cultivars with host plant resistance incorporated in high suitability predicted areas, or establishment of susceptible crops in low suitability sites, within intensive livestock and agriculture systems.

Our main objective was to determine the current distribution of these six species and estimate the potential distribution under two future climate scenarios via ecological niche methods based on presence-only data.

In total 590 occurrence records were obtained: 115 from human observations, 299 from CIATARC, 108 from GBIF, 24 from SpeciesLink, and 44 from literature review. After data cleansing, 48, 186, 19, 71, 55, and 120 points were used for A. lepidior, A. reducta, A. varia, P. simulans, Z. carbonaria, and Z. pubescens, respectively. Maps showing the occurrence records, estimation of current distribution and future potential distribution (2041–2060) under SSP126 - SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585) are presented in Figures 1–6. Suitable areas and suitability values as well as bioclimatic layers selected based on multicollinearity analysis differed according to the species in the study (Figure 7). Consequently, probability of occurrence (i.e., suitability) in the niches of each species as a function of two most representative biovariables (Figure 8) varied according to species. In general, the ensembled models reached acceptable values for metrics used to evaluate ENMs accuracy (see Supplementary Table S1). The most used, AUC and TSS metrics, showed high scores for all species under study indicating robust performance (Figure 9).

A. lepidior occurred in southern and central Costa Rica, central Panama, and northern Colombia. The ENM estimated a suitable area in central and north Colombia and some areas of Venezuela (AUC 0.97 ± 0.05, TSS 0.80 ± 0.1) (Figure 1). Bioclimatic layers with high contribution were isothermality (bio_3) and temperature seasonality (bio_4), showing high suitability with high values of bio_3 (>70 %) and low values of bio_4 (<77.45%) (Figure 7). Averages of AUC and TSS (± SD) were 0.97 ± 0.05 and 0.80 ± 0.1, respectively (Figure 9, Supplementary Table S1). A considerable increase in suitability is expected for large areas of Amazonas ecoregion of Peru, Venezuela, and the north of Brazil even in the optimistic scenario, with possible invasions in those sites and western Ecuador, northeastern Peru and northern Bolivia (Figure 1). Also in Panama, Costa Rica, and, in the pessimistic scenario, in Guatemala and Belize. Small areas in a few sites of the Pacific coast of Central America and tropical South America show a decrease in suitable areas for this species.

A. reducta occurred in Costa Rica, central Panama, and central and northern Caribbean Colombia. Fewer records were obtained in northwestern Venezuela and northern Brazil. The ENM estimated a suitable area in southern Costa Rica and Panama, as well as Eastern Ranges and Caribbean coast in Colombia, and Andean Venezuela (Figure 2). Bioclimatic layers with high contribution were minimum temperature of coldest month (bio_6) and isothermality (bio_3), showing high suitability with high values of both bio_6 and bio_3 (Figure 7). Average of AUC and TSS (± SD) was 0.94 ± 0.01 and 0.88 ± 0.05, respectively. An increase in suitable areas and possible invasions are expected for the future optimistic scenario in Colombian and Venezuelan Llanos and Colombian Caribbean region. In the pessimistic scenario, Amazonas ecoregion of Peru and Brazil, along with some sites in southern Costa Rica, Panama, Dominican Republic, and Mexico are predicted to be susceptible to new invasions (Figure 2).

A. varia occurred in central and southwestern Colombia and northwestern Venezuela. The ENM estimated a suitable area in Amazonas ecoregion of Colombia, Venezuela, and northern Brazil, and a smaller region in northern Peru (Figure 3). Bioclimatic layers with high contribution were precipitation of the coldest quarter (bio_19), temperature seasonality (bio_4), and precipitation seasonality (bio_15) (Figure 7). Average AUC and TSS (± SD) was 0.97 ± 0.01 and 0.89 ± 0.05, respectively. A decrease in suitable areas is expected for future scenarios compared to the same sites in current sites. Also, extinction is predicted in a few areas of Colombian and Venezuelan Llanos (Figure 3).

P. simulans was the most widespread species in this study. Occurrence records were obtained mostly from North America (Mexico) and Central America, with fewer records in western Colombia (Figure 4). Bioclimatic layers with high contribution were precipitation of the wettest month (bio_13) and precipitation of the coldest quarter (bio_18), showing high suitability with values <1,060 of bio_18 and values between 468 and 900 of bio_13 (Figure 7). Average of AUC and TSS (± SD) was 0.91 ± 0.06 and 0.73 ± 0.12, respectively. ENMs showed more habitats in South America and a small area in the Pacific Coast of Central America but with low suitability. An increase in suitability and possible invasions for small areas of Brazil Cerrado in both scenarios, along with Venezuelan Llanos in the optimistic scenario, and a noticeable decrease in Costa Rica is expected (Figure 4).

Z. carbonaria has been recorded only in western Colombia, across central Andes. The ENM estimated higher suitability in Colombian and Ecuadorian Andes (middle tropic) and the Amazonian Piedmont of Colombia, decreasing its values to zero in Colombian and Venezuelan Llanos (low tropic) (Figure 5). Bioclimatic layers with high contribution were isothermality (bio_3) and precipitation seasonality (bio_15), showing high suitability with values close to 40 of bio_15 and high values of bio_15 (Figure 7). Average of AUC and TSS (± SD) was 0.99 ± 0.02 and 0.93 ± 0.07, respectively. A decrease in suitability for the Amazonian Piedmont of Colombia and the Andes is expected (Figure 5).

Finally, Z. pubescens occurred widely in western and central Andes of Colombia, northern Ecuador and western Brazil, including Amazon and Cerrado biogeographic zones. Fewer records were obtained in southern Peru and northern Suriname (Figure 6). Bioclimatic layers with high contribution were temperature seasonality (bio_4) and precipitation seasonality (bio_15), showing high suitability with low values of bio_4 (<10) and values close to 40 of bio_15 (Figure 7). Average of AUC and TSS (± SD) was 0.89 ± 0.03 and 0.66 ± 0.09, respectively. An increase in suitability is expected for some areas of Ecuador, Peru, and Brazil in both climate change scenarios, being greater in the pessimistic scenario (Figure 6).

In our study ENMs of the occurrence data had a high grade of accuracy given the sample size of five species, except for A. varia, for modeling (>25 records) (van Proosdij et al., 2016; Schöbel and Carvalho, 2020). Despite small sample sizes methodologies based on calculation p-values through Jackknife are implemented in the SDM R package used in this study (Naimi and Araújo, 2016), more records may increase the model accuracy (van Proosdij et al., 2016). Low records for spittlebugs were previously reported for Mahanarva in Brazil (Schöbel and Carvalho, 2020) being underrepresented in occurrence databases. This phenomenon was also observed for the six species studied as most of the records were obtained from CIATARC collection and expert's reports through the years (human observation).

The ENMs also revealed differences in the distribution and ecological niche of the six spittlebug species in South America showing that these species ecological niche varies widely in the Neotropic, and has the potential to invade large areas, where livestock systems coincide. A. reducta y A. lepidior have great potential to impact grassland mainly in Colombian and Venezuelan Llanos where susceptible pastures (e.g., Urochloa decumbens) and sugarcane are planted in large areas. Another ecoregion where these two species have high suitability is the Amazonian ecoregion in Colombia and Brazil, where livestock extensive systems are increasing indiscriminately.

The evidence showed that Z. pubescens is distributed in a wide altitudinal range (8–3225 m.a.s.l) but with a local reduced temperature seasonality. Elevation has been reported as the most important variable with the highest contribution in the ENMs in other spittlebugs (Schöbel and Carvalho, 2020). Few species have such a wide altitudinal range, which allows us to propose two hypotheses: (1) Z. pubescens presents extreme thermal limits and (2) the species presents geographically separated populations. A case of biotypes is observed for the spittlebug Calitettix versicolor in China, which diverged in two lineages consistent with biogeographical regions separated by Hengduan Mountains (Yang et al., 2016). Similarly, this could be happening with Z. pubescens influenced by the Colombian Andes. Although the species is reported in Brazil (27 occurrence records; average of 400 m.a.s.l), the suitability values are lower than in Colombia and Ecuador (93 occurrence records; average of 1079 m.a.s.l.). The higher number of records in the highlands of Colombia and Ecuador could be causing an overestimation of the occurrence probability at these areas over the records of Cerrado places in Brazil, this would explain the current potential distribution estimated, and also could be reflecting the possible existence of, at least, two populations with different ecological niches.

The position of a species within an ecosystem is determined by the interactions with their biotic and abiotic environment (Polechová and Storch, 2019). Tropical spittlebugs have a seasonal dynamic strongly synchronized with rainfall patterns. For instance, Z. carbonaria and A. reducta in Colombia, P. simulans in Colombia and Venezuela, D. flavopicta in Brazil, as well as A. contigua and A. contigua in Mexico, reduce diapause rates and a higher abundance of nymphs is observed after rain season start (Peck et al., 2001, 2002; Sujii et al., 2002; Olán-Hernández et al., 2016). Hence, a strong effect of the biovariables 12 to 19 in the models, related with precipitation, in the models was expected but in our estimations, the distribution of habitat suitability of these six species also involved environmental variables related to temperature suggesting that variables derived from temperature has a strong effect on the biology of these species. For P. simulans, precipitation was more important than temperature to determine its distribution with a relative importance over 0.4 for precipitation of the wettest month, thus, greater probabilities of occurrence happen in precipitation between 500 and 940 mm. In general, the habitat suitability estimated for two-dimensional niches was low as the biovariables' relative importance varied among all the species with values below 0.4 (Figure 8). Similar results were obtained by Schöbel and Carvalho (2020) in ENM of four Mahanarva species showing that most of the WorldClim variables did not contribute to their analysis and that for M. fimbriolata and M. spectabilis the biovariables had contribution percentages from 15 to 27%.

Regarding the climate change scenarios proposed, we found that these have a significant influence on the potential distribution of the species in study, increasing the suitability value and suitable area for some (mainly for A. reducta and A. lepidior) or decreasing them for others (A. varia). Previous studies showed a declining tendency in suitability for Mahanarva across Central and South America (Fonseca et al., 2016; Schöbel and Carvalho, 2020) and Philaenus spumarius in North America (Karban and Huntzinger, 2018). Global warming and longer drought periods contribute to accelerate this phenomenon as spittlebug biology is highly dependent on plant water status. Being xylem feeders, they require excessive amounts of sap which flow is subject to transpiration (Novotny and Wilson, 1997). Under water stress conditions transpiration rates decrease as well as food availability for spittlebugs, particularly in the nymphal stages. Besides, these conditions may affect nymph thermoregulation by foam or “spittle” production, composed mainly of excreted semi-digested plant fluid, fatty acids, carbohydrates, mucopolysaccharides, and proteins produced by Malpighian tubules (Rakitov, 2002; Tonelli et al., 2018). Since the six species are Urochloa spp. key pests, a future limitation of ecological niche in future scenarios in livestock production zones should be taken into account as improved resistant grasses to spittlebug attack and increase the number of forage species are considered a sustainable strategy for the livestock systems under climate change (Rao et al., 2016; Schiek et al., 2018). Competition can influence species future distribution as well. Despite reaching the spittlebug habitat's food limits is unlikely (Schöbel and Carvalho, 2020), the variation among species' life cycles may determine the success of one species over others. A. reducta was reported for the first time in 2019 in Cauca River Valley, Colombia (Hernandez et al., 2021) where A. varia is a key pest of sugarcane and P. simulans of signalgrass [Urochloa decumbens cv. Basilisk; (Rodriguez Chalarca et al., 2003; Gómez, 2007)]. In Colombian Caribbean coast, A. reducta's entire life cycle is shorter (45.2 days) compared with A. varia (62 days) or P. simulans (71.9 days) in Cauca River Valley conditions (Peck et al., 2002; Rodriguez Chalarca et al., 2003; Castro Valderrama et al., 2011). Thus, A. reducta can coexist or even displace these two species in sugarcane and signalgrass for potentially having more generations per year in the region where ~208 thousand ha of sugarcane was harvested in 2018 (Asocaña., 2019).

The current study contributes to the ecological knowledge of spittlebugs, which will be useful in the development of prevention and control strategies for this pest in South America. Finally, we suggest carrying out studies of physiology and genetics of populations to determine the thermal limits of the species and to corroborate if there are genetic divergences between geographically separated populations.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

LH and PE contributed to the data collection, data curation, analysis and interpretation of maps, manuscript preparation, and supervision. DF contributed to data collection and manuscript preparation. VC contributed to manuscript preparation. JC contributed to interpretation and manuscript preparation. MG-J contributed to data collection, interpretation of maps, and manuscript preparation. All authors contributed to the article and approved the submitted version.

This work was funded by the CGIAR Research Program on Livestock. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. We also acknowledge the financial assistance of GROW Colombia from the UK Research and Innovation (UKRI) Global Challenges Research Fund (GCRF) (BB/P028098/1).

DF was employed by the company Semillas Papalotla S.A. de C.V. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.