Geographic Distribution of Colombian Spittlebugs (Hemiptera: Cercopidae) via Ecological Niche Modeling: A Prediction for the Main Tropical Forages' Pest in the Neotropics

Spittlebugs (Hemiptera: Cercopidae) are the main tropical pests in Central and South America of cultivated pastures. We aimed to estimate the potential distribution of Aeneolamia varia, A. lepidior, A. reducta, Prosapia simulans, Zulia carbonaria, and Z. pubescens throughout the Neotropics using ecological niche modeling. These six insect species are common in Colombia and cause large economic losses. Records of these species, prior to the year 2000, were compiled from human observations, specimens from CIAT Arthropod Reference Collection (CIATARC), Global Biodiversity Information Facility (GBIF), speciesLink (splink), and an extensive literature review. Different ecological niche models (ENMs) were generated for each species: Maximum Entropy (MaxEnt), generalized linear (GLM), multivariate adaptive regression spline (MARS), and random forest model (RF). Bioclimatic datasets were obtained from WorldClim and the 19 available variables were used as predictors. Future changes in the potential geographical distribution were simulated in ENMs generated based on climate change projections for 2050 in two scenarios: optimistic and pessimistic. The results suggest that (i) Colombian spittlebugs impose an important threat to Urochloa production in different South American countries, (ii) each spittlebug species has a unique geographic distribution pattern, (iii) in the future the six species are likely to invade new geographic areas even in an optimistic scenario, (iv) A. lepidior and A. reducta showed a higher number of suitable habitats across Colombia, Venezuela, Brazil, Peru, and Ecuador, where predicted risk is more severe. Our data will allow to (i) monitor the dispersion of these spittlebug species, (ii) design strategies for integrated spittlebug management that include resistant cultivars adoption to mitigate potential economic damage, and (iii) implement regulatory actions to prevent their introduction and spread in geographic areas where the species are not yet found.


INTRODUCTION
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).
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.

Climatic Data
Elevation layer and 19 bioclimatic layers (bio_1 to bio_19) were obtained from Worldclim from 1970 to 2000 using raster::getData function. For the current climate data, the Version 2 Bioclimatic variables with a spatial resolution of 2.5 min were selected (Fick and Hijmans, 2017) with the aim of maintaining the same spatial resolution of the species georeferenced (Sillero and Barbosa, 2021). To extract values from the bioclimatic layers, the extract function was used. Finally, the species names were combined with coordinates (latitude, longitude), bio_1 to bio_19, and elevation values into a single data.frame.

Collinear Variables Removal
To prevent any multicollinearity-related bias in the models, a collinearity test among bioclimatic variables was performed using the vifstep function. Collinearity describes the situation where two or more predictor variables in a statistical model are linearly correlated (Alin, 2010). Therefore, it could inflate both the standard error and the confidence intervals, and prevent the determination of the significance of each variable on the dependent variable (Quinn and Keough, 2002). Variables with VIF (Variance Inflation Factors; Chatterjee and Hadi, 2006) values < 0.7 were selected for the subsequent analyzes. We created a sdmData object including species and previously selected variables, which means low collinearity, as predictors. Approximately 1,000 'pseudo-absences' points were randomly selected over the study geographical area for each species using argument method='gRandom'. Pseudo-absence refers to cells in which the species has not yet been recorded, not to cells in which the species is necessarily absent (Phillips et al., 2009).

Model Fitting
We used four species distribution models to predict the distribution of each spittlebug species under study. All models were based on presence and pseudo-absence data: Maximum Entropy (MaxEnt), Generalized Linear Model (GLM), Multivariate Adaptive Regression Spline (MARS), and Random Forest (RF) models. MaxEnt was used as default settings since it has shown the ability to achieve good performance as a default (Phillips and Dudík, 2008). Models are fitted with sdm function using two replication techniques (subsampling and bootstrapping) establishing 70% of the occurrence data as training data and 30% as test data. This process was repeated 3 times. As a result of our methodological procedure, a total of 24 different projections (4 models * 2 replication techniques * 3 repetitions) were generated for each species.

Model Prediction and Ensemble
We consider the accessible area of species under study as the entire neotropical ecoregion and that the species do not have restrictions since in this ecoregion there is a large pasture monoculture for livestock and it has a wide sugarcane planted area where cercopids can be established (Jank et al., 2014;Schöbel and Carvalho, 2021). The hypothesis was that climate change will impact or lead to an increase of future potential distributions of the species under study. Models obtained were used to estimate the current distribution in South America using the predict function from the sdm package. This function allows making a raster object with predictions from several fitted models (Naimi and Araújo, 2016). All 24 predictions were ensemble in one using the ensemble function which provides a consensus of multiple models. By combining projections from different models, errors tend to be canceled out thus aiding predictive accuracy (Diniz-Filho et al., 2010).

Model Evaluation
To evaluate model outputs, we used the receiver operated characteristics, analyzing the area under curve (AUC) (Fielding and Bell, 1997) and the true skill statistic (TSS) (Allouche et al., 2006). The AUC value is a standard method to assess the accuracy of predictive distribution models, AUC values below 0.7 were considered poor, 0.7-0.9 moderate, and >0.9 good (Araújo et al., 2005). TSS compares the number of correct forecasts, minus those attributable to random guessing, to that of a hypothetical set of perfect forecasts. TSS values close to one denote an ideal prediction; values of zero or less denote a prediction that is not better than random (Allouche et al., 2006). For each species, the relative importance of bioclimatic variables selected based on multicollinearity analysis and AUC metric were plotted.

Future Distribution Model
To build future potential distribution, we used the BCC-CSM2-MR global climate model from the Coupled Model Intercomparison Project 6 [CMIP6; available for use in the WorldClim (https://www.worldclim.org/data/cmip6/ cmip6climate.html); (O'Neill et al., 2016)] and two shared socio-economic pathways [(SSP); (1) SSP126: an optimistic scenario increasing shift toward sustainable practices with low greenhouse gas concentration levels and (2) SSP585: a pessimistic scenario that assumes an energy intensive, fossil-based economy with increasing greenhouse gas emissions over time Riahi et al., 2017)] in a 2.5-min resolution. Habitat suitability was modeled using selected previously bioclimatic layers under each SSP scenario. In this study, only one time period was used for near future prediction: 2050 (average for 2041-2060). To quantify the change between current and future distribution, maps were converted from probability of occurrence to presence and absence. For this, the mean threshold (occurrence probability values) was used in the ifelse function which allows reviewing the probability values. If the probability values are greater than or equal to the average threshold, the new value assigned is 1 (presence) and if the probability value is less on the average threshold, the new value assigned is 0 (absence). Later, the current distribution raster was subtracted from the future distribution raster, as a result, possible extinction and invasion were plotted.

RESULTS
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 Frontiers in Sustainable Food Systems | www.frontiersin.org FIGURE 1 | Ecological niche models of Aeneolamia lepidior. Distribution records, current potential distribution, future potential distribution (2041-2060) under SSP126 and SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585). The scale shows the habitat suitability being 1 = higher suitability. Scale in change maps −1 = possible extinction and 1 = possible invasion.
Frontiers in Sustainable Food Systems | www.frontiersin.org FIGURE 2 | Ecological niche models of Aeneolamia reducta. Distribution records (red point indicates the most recent report in a new niche), current potential distribution, future potential distribution (2041-2060) under SSP126 and SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585). The scale shows the habitat suitability being 1 = higher suitability. Scale in change maps −1 = possible extinction and 1 = possible invasion.
Frontiers in Sustainable Food Systems | www.frontiersin.org FIGURE 3 | Ecological niche models of Aeneolamia varia. Distribution records, current potential distribution, future potential distribution (2041-2060) under SSP126 and SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585). The scale shows the habitat suitability being 1 = higher suitability. Scale in change maps −1 = possible extinction and 1 = possible invasion.
Frontiers in Sustainable Food Systems | www.frontiersin.org FIGURE 4 | Ecological niche models of Prosapia simulans. Distribution records, current potential distribution, future potential distribution (2041-2060) under SSP126 and SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585). The scale shows the habitat suitability being 1 = higher suitability. Scale in change maps −1 = possible extinction and 1 = possible invasion.
Frontiers in Sustainable Food Systems | www.frontiersin.org FIGURE 5 | Ecological niche models of Zulia carbonaria. Distribution records, current potential distribution, future potential distribution (2041-2060) under SSP126 and SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585). The scale shows the habitat suitability being 1 = higher suitability. Scale in change maps −1 = possible extinction and 1 = possible invasion.
Frontiers in Sustainable Food Systems | www.frontiersin.org FIGURE 6 | Ecological niche models of Zulia pubescens. Distribution records, current potential distribution, future potential distribution (2041-2060) under SSP126 and SSP585 scenarios, and comparison between current and future scenarios (change SSP126 and SSP585). The scale shows the habitat suitability being 1 = higher suitability. Scale in change maps −1 = possible extinction and 1 = possible invasion.
Frontiers in Sustainable Food Systems | www.frontiersin.org   ( 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).

DISCUSSION
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 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 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.

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

AUTHOR CONTRIBUTIONS
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.

FUNDING
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).

ACKNOWLEDGMENTS
We would like to thank the person who enriched the spittlebugs collection deposited in CIATARC: Daniel Peck. Also, to the organizations that allowed the use of their data for our research: GBIF and SpeciesLink. This work was carried out as part of the CGIAR Research Program on Livestock. We thank all donors who globally support our work through their contributions to the CGIAR System. CGIAR is a global research partnership for a food-secure future. Its science is carried out by 15 Research Centers in close collaboration with hundreds of partners across the globe. We also thank the reviewers for their constructive comments that helped to improve the manuscript.