Remote sensing and citizen science to characterize the ecological niche of an endemic and endangered Costa Rican poison frog

. Habitat encroachment can have devastating effects upon biodiversity, especially amphibians. Phyllobates vittatus is an endemic frog from Costa Rica, where land cover has seen signiﬁcant changes over recent decades. Here we use remote sensing to create a land cover map of the region and carry out ecological niche modelling to identify the main abiotic factors associated to the distribution of this species. We have informed our models based on our own ﬁeld observations, those from other researchers


Introduction *
Amphibians constitute the most threatened group of vertebrates, with ∼50% of described species facing extinction (González-del-Pliego et al., 2019). Habitat change is one of the main threats, alongside the disease chrytridiomicosis caused by pathogens such as Batrachochytrium dendrobatidis and B. salaman-ner et al., 2007). Many amphibians, particularly endemic species, have small ranges and are therefore more vulnerable to habitat loss (Staude et al., 2019). In many cases, little is known about rare species ecology and demography, making it difficult to develop well-planned and informed conservation plans (Purvis et al., 2000;Bishop et al., 2012, Conde et al., 2019. Here we focus on one such species, Phyllobates vittatus (Cope, 1893) known as the Golfo Dulce Poison Frog, an amphibian species of the family Dendrobatidae, endemic to the southern Pacific rainforest of Costa Rica.
Phyllobates vittatus has been recently catalogued as Vulnerable (VU) by the IUCN (International Union for Conservation of Nature) (IUCN SSC Amphibian Specialist Group, 2020), but it has been catalogued as Endangered (EN) for the previous three revisions. Moreover, the IUCN Red List lists the population trend as decreasing and distribution notes highlight the severe fragmentation of the populations (De Plecker and Aguilar Montenegro, 2020;Mebs et al., 2014;Ryan, 2002). Its main threats are identified as anthropogenic habitat encroachment in the form of deforestation followed by transformation into plantations or other agricultural lands, water pollution due to the contamination from gold mining activities, and the illegal pet trade (IUCN SSC Amphibian Specialist Group, 2020).
Costa Rica has suffered high rates of deforestation, leading to severely fragmented forests (Sánchez-Azofeifa et al., 2001), in many cases reaching right up to the boundaries of the national parks (Sánchez-Azofeifa et al., 2002a), thus isolating many of the protected areas. On the Osa Peninsula, the core of the distribution of P. vitattus, 44% of the forest cover is degraded or affected by fragmentation (Jiménez, 2000;Sánchez-Azofeifa et al., 2002b). The continuous changes in the land use of the Southern Pacific coast of Costa Rica, may have led to a reduction of suitable habitat for P. vittatus, leading to the severe fragmentation of the populations. Given these circumstances, it has become particularly critical to characterize the ecological niche of P. vittatus and assess potential suitable habitat for the species. Knowledge on the ecological niche of a species is key to inform suitable conservation measures (Thorn et al., 2009;Urbina-Cardona and Flores-Villela, 2010). It can provide information on population connectivity (Thornton and Peers, 2020), indicating which populations might be under threat of becoming isolated. It can also determine potential suitable areas for future surveys (Brito et al., 2011) or reintroductions (Kuemmerle et al., 2010, assessment of conservation areas (Peers et al., 2016), evaluation of the effect of climate change, and even lay out further research questions (Thornton and Peers, 2020). Furthermore, the combination of citizen science and other sources to expand the database (Robinson et al., 2020), and the use of remote sensing to obtain land cover, topographic and climatic variables can improve the predictions from these models (Thornton and Peers, 2020).
In this study (1) we create an up-to-date land cover map for the Southern Pacific coast of Costa Rica using remote sensing and (2), characterize the realized niche of P. vittatus using ecological niche modeling (Sillero, 2011). Since the habitat of P. vittatus has been described to be associated to humid lowland forests and streams (Savage, 2002;Silverstone, 1976), we hypothesize that forest cover and distance to rivers will be the main contributors to the predicted distribution. Moreover, with the continued changes in land cover in the South Pacific, we would expect to find most of the extent of its potential geographic distribution within protected areas. Thirdly, we aim to extend the general knowledge of the biology of this species and identify suitable habitats in regions currently not under protection to inform conservation efforts for this threatened and endemic species.

Species locality data
To analyse the ecological niche of P. vittatus, we collected occurrence records from field surveys (our own conducted for this study and those of colleagues), citizen science and online databases. Field surveys: we compiled encounters (visual and/or acoustic) recorded between 2018 and 2019 from personal surveys (n = 19 encounters). We also gathered the data from occurrences shared by colleagues from 2018 to 2020 (n = 80) and published works from 2001 to 2019 (n = 11).
Citizen science: to multiply the effort, we involved local communities and showed them how to identify the target species with the open-source smartphone application iNaturalist (https://www.inaturalist.org). Other studies have previously shown how citizen science can work as a tool to multiply researchers' efforts, contribute to expand biodiversity knowledge (Chandler et al., 2017;Kobori et al., 2016), and to improve distribution and ecological niche models (Miller et al., 2019;Reich et al., 2018;Robinson et al., 2020). Phone apps like iNaturalist are curated by peers, including experts that confirm the species identification (Nugent, 2008). Authenticity of records can further be verified, as these records are generally accompanied by their respective photographs. We carried out 13 workshops with 243 participants. In June 2020, we collected all the observations available for P. vittatus from iNaturalist (n = 155; from 2002 to 2020), ∼50% of them recorded during the project. Since iNaturalist does not have a protocol for researchers to obtain obscured coordinates of endangered species, we reached out to users to share the coordinates from their P. vittatus observations.
Online databases: we also acquired data from the open-source Global Biodiversity Information Facility database (www.GBIF.org, accessed 6 July 2020; GBIF Occurrence Download https://doi.org/10.15468/dl.e5ev5u). All records (n = 141; from 1961 to 2013) were assessed and vetted as follows. Observations of P. lugubris, mistaken for P. vittatus were discarded based on images provided alongside the observations and their location (Caribbean not Pacific coast). Replicated observations that were recorded on the exact same coordinates, date and time, were discarded; as well as those with uncertainty regarding the exact location over 1 km 2 or with invalid coordinates (na values). Six occurrences recorded by SMNS in 1989 were reinterpreted from latitude −38 to −83, since −38 was located in the Atlantic Ocean and −83 is the latitude where the Osa Peninsula is located. Moreover, when we re-interpreted the coordinates, these occurrences were located in an area where the presence of P. vittatus has been recorded by many researchers. Overall, we were able to collect 406 occurrences, of these, 254 were considered valid occurrences.
Records from all three sources were combined and then filtered to discard any duplicated occurrences (based on identical GPS coordinates). Finally, in order for the records to be comparable with remote sensing data, we kept only records collected between 2016 and 2020. This data processing resulted in a total of vetted 156 records, from a total of 406 crude records. To reduce issues due to spatial biases during niche modelling , we reduced this dataset further to one point per km 2 , resulting in 67 independent occurrences for the analysis. Due to the high level of vulnerability of the species and because illegal pet trade is a major threat, the specific coordinates of the field observations collected for this study will not be made publicly available. Moreover, some stakeholders agreed to collaborate on this study on the condition that the exact locations would not be made public, given the risk of possible poaching. However, a map with obscured locations of the 67 independent occurrences used in this study can be seen in supplementary fig S1 for a better understanding of this study.

Land cover map
Since the habitat of P. vittatus is related to forested areas (Savage, 2002), the forest percentage is likely a strong feature associated to its ecological niche. To get an up-to-date forest cover of the study area, we mapped the land cover using data from the Sentinel-2 satellite, processed using QGIS (QGIS.org, 2020. QGIS Geographic Information System. QGIS Association). For a detailed description of the land cover classification, see supplementary text S1 and supplementary fig. S2. Our study area covered four tiles (remote sensing imagery stored as large images): L16PKK, L17PKK, L16PKL, L17PKL, which were stacked as mosaics and then clipped to our study area. Tiles are 100 × 100 km 2 ortho-images in Universal Transverse Mercator projection (UTM). We used the data of these four tiles for two dates (2 March 2019 and 10 February 2019) to reduce the cloud cover interference. To differentiate vegetation from non-dense vegetated areas, we calculated the Normalized Difference Vegetation Index (NDVI) from the images (Tucker, 1979). Then, we carried out an unsupervised 'k-means' classification (Likas et al., 2003) on each image, clustering the pixels into 35 groups. Each cluster was inspected to determine which land cover type was the defining feature. Clusters that were characterized by the same land cover type were reclassified and grouped together. Clusters corresponding to forested areas were processed further, since palm, white teak and pineapple plantations are all classified as dense vegetated areas. These plantations were distinguished from natural forests by manually inspecting these areas with high-resolution Google Satellite orthophotos. The final clustering resulted in nine categories: Forest, tall/dense grassland, pasture, mangrove, palm plantation, pineapple plantation, white teak plantation, wetland and non-dense vegetated areas. We evaluated the land cover classification estimating accuracy and area from the sample data in an error matrix. We applied the proportional approaches sample allocations methodologies suggested by Olofsson et al. (2014), using the Cochran´s (1977) sample size formula for stratified random sampling. Using this land cover map, we calculated the corresponding percentage of forest for every 1 km 2 pixel.

Topographic and climatic variables
Several other abiotic variables of interest, in addition to the land cover, were prepared for our ecological niche modelling. Elevation was estimated based on the digital elevation model (DEM) acquired from Aster (imagery courtesy of "NASA/METI/AIS/Japan Spacesystems, and U.S./Japan ASTER Science Team, Aster"; 30 meters resolution). Aspect and slope were also considered as predictors as they can influence the ground and water temperatures at a given site (Brown and Krygier, 1970). Aspect and slope were obtained through spatial analysis of the DEM using GDAL 3.1.9 (GDAL/OGR contributors 2021). Nineteen bioclimatic variables were obtained from WorldClim 2.1 (Fick and Hijmans, 2017). Considering the little information available about the ecology of the Golfo Dulce Poison Frog (Savage, 2002;Silverstone, 1976), choosing a priori which of the nineteen bioclimatic variables would be most important to its biology was not contemplated for this study as it would be subjective to the authors' opinion. WorldClim 2.1 climatic data is at a spatial resolution of 1 km 2 and corresponds to the times series for 1970-2000 (data available at the time of analysis, January 2020). This represented a limitation due to this being a different time frame as our occurrence's records.
Phyllobates vittatus is often found near streams. We therefore generated a spatial raster for: distance to rivers and water bodies (lakes and community-managed aqueducts known as ASADAS [Administrative Associations for Aqueducts and Sewers]), extracted from the Costa Rican digital atlas from 2014 (Ortiz-Malavasi, 2016). Unfortunately, many ephemeral and narrow creeks and streams were not incorporated in our model, as they had not been mapped at the time of our study. To standardize the spatial resolution, all variables were transformed to a resolution of 1 km 2 .

Ecological niche model fitting
Model fitting to determine the ecological niche of P. vittatus was done using the machinelearning algorithm of Maximum Entropy (Max-Ent) (Philips et al., 2006) in the programming environment RStudio (RStudio Team, 2020), through the dismo 1.1-4 (Phillips et al., 2006) and ENMeval 2.0.4 packages (Kass et al., 2021), run with rJava 0.9.13 (Urbanek, 2020). Occurrence data was previsualized in R, and the land cover, topographic and climatic variables layers were stacked and used as predictors. To reduce collinearity in the predictors, we carried out an exploratory analysis to visualize correlations among predictors (supplementary table S1). We found no significant correlation between the land cover forest percentage and topographic variables. We did not find any correlation between the climatic variables and most of the other variables, either. However, we did find a strong correlation between some bioclimatic variables and the elevation (>0.7, supplementary table S1). To reduce multicollinearity in our models, we excluded annual precipitation, precipitation of wettest month and precipitation of driest month, as they are already within precipitation of wettest quarter, driest quarter, warmest quarter and driest quarter. We kept these last four variables and precipitation seasonality, as they had no correlation among them. Regarding temperature-related variables, we only included annual range and elevation, since most of the other variables were explained by these two variables. Ultimately, we reduced the variables with which we ran the models to thirteen (supplementary table S1).
To test the effect of background point selection for model validation, we ran models with 300 background points sampled from differentsized background areas (Elith et al., 2011;Merow et al., 2013;VanDerWal et al., 2009). To exclude areas that had not been sampled in our study (Elith et al., 2011), we created two models where the background points were drawn from buffer areas at increasing distances around the occurrences (VanDerWal et al., 2009;Feng et al., 2017). The distances chosen were 5, and 10 km in order to take background points from other areas over the uncertainty of the occurrences, and still not cover the whole extent of our study area ( fig. 1). To test the effect that the total extent could have (Elith et al., 2011;Merow et al., 2013;VanDerWal et al., 2009), we also added a model where the background points were drawn from our total study area (from the coastline to 1,500 altitude) (Mebs et al., 2014;Ryan, 2002;Savage, 2002).
The occurrence data was divided into training (80%) and testing (20%) data using 5-fold cross-validation (Merow et al., 2013). For each background extent, we tested different combinations of the parameters "feature class" and "regularization multiplier" (herein β; Elith et al., 2011;Merow et al., 2013). We tested all the combinations of the "feature class" that Max-Ent supports: linear (L), quadratic (Q), product (P), threshold (T), and hinge (H) (Morales et al., 2017;Perkins-Taylor and Frey, 2020). We also tested these models at different levels of β: 1, 2, 3, 4, 5, 10, 15, and 20 (Morales et al., 2017;Perkins-Taylor and Frey, 2020). For each of the three background extents, we chose the best model based on the Akaike information criterion or AICc (Warren and Seifert, 2011). The best model for each background extent was the one scoring the smallest AICc value and which delta AICc was 0 (Morales et al., 2017). Furthermore, we compared the average AUCtest and AUCtrain between the background extents. Finally, we performed a prediction of each model for the whole extent of the study area in a map showing the relative suitability of the sites on a scale from 0 to 1.

Land cover map
The study area encompasses a total of 10,226.85 km 2 ; 8,966.64 km 2 from Costa Rica and 1,260.21 km 2 from Panama. The land cover map of our study area was developed at 10 m spatial resolution, in UTM zone 17 north in World Geodesic System 1984 datum ( fig. 2). Our land cover classification has a 92.1% overall accuracy with a confidence level of 97.5%. The classes of the land cover map corresponded to forest (4001.94 ± 106.84 km 2 ); urban, bare soil and others (3855.4 ± 89.6 km 2 ); dense grassland (963.4 ± 121.31 km 2 ); palm plantation (864.25 ± 53.65 km 2 ); pasture (382.82 ± 56.73 km 2 %); mangrove (159386 ± 0 km 2 ); pineapple plantation (83.17 ± 0 km 2 ); teak plantation (45.47 ± 0 km 2 ); and the Corcovado lagoon (11.32 ± 0 km 2 ). The main forest cover remains within Corcovado National Park, the Golfo Dulce Forest Reserve and Piedras Blancas National Park. Outside the Osa Peninsula, forest cover is highly segregated, resulting in a complex matrix of different classes. Palm plantations are abundant in the south, fragmenting the forest in this area and, subsequently, isolating the forest patch in the area of Punta Burica. The forest corridor corresponding to the biological corridor Paso de la Danta runs through Uvita and Dominical. Some forest remnants that appear as fragmented intrusions into urban and pasture areas in the Terraba river valley were difficult to identify in the land cover classification, due to their small area. However, these forest remnants are still key for species like P. vittatus, as the observations records from this valley were located in these patches. Of the total 156 occurrences, 133 were in areas with a percentage of forest cover greater than 47%. The remaining 23 occurrences in sites below this value were the ones located in fragmented forest remnants.

Ecological models
The model built from the total extent with β = 3 and linear and quadratic autofeatures was the model chosen as it had the lowest AICc value and the highest AUC value of the three (Table 1). According to the permutation assessment of the importance of each variable, we observed that the main variables for this model were elevation (57.5%, supplementary  fig. S3A), distance to lakes (13.7%, supplementary fig. S3B), forest percentage (13.4%, supplementary fig. S3C) and annual temperature range (bio 7, 6.9%, supplementary fig.  S3D). The other variables, which contributed less than 5%, can be seen in supplementary table S3). Relative suitability decreased logarithmically as we move to higher elevations (supplementary fig. S4A). It increased linearly with distance from lakes (supplementary fig. S4B), whereas it decreased with distance to rivers (supplementary fig. S4C). Increased forest percentage increased linearly the relative suitability (supplementary fig. S4D). The response to annual temperature range was a negative linear response, with relative niche suitability going from over 0.6 to under 0.35 within a temperature increase of 4°C (supplementary fig. S4E).
Based on this best-performing model, Osa Peninsula, Punta Burica, and the biological corridor Paso de la Danta ( fig. 3) are predicted to be key hotspots for potentially suitable P. vittatus habitat ( fig. 3). Nonetheless, the prediction also shows how the relative suitable sites are segregated in patches. The valley of Terraba river is also predicted to be of high suitability ( fig. 3), clearly highlighted in these models, but also highly fragmented.

Discussion
We used remote sensing to create a land cover map of the region and carry out ecological niche modelling to identify the main abiotic factors associated to the distribution of P. vittatus. We informed our models based on our own field observations, those from other researchers, and from citizen science participants to obtain a comprehensive database of P. vittatus occurrences. Elevation, forest percentage, distance to lakes and community-managed aqueducts were found to shape the ecological niche of P. vittatus -mostly located within protected areas. We identify populations that might be isolated, and areas where presence has not yet been verified or that have not been occupied by the species and calculated the area of occupancy. Finally, we suggest reassessing P. vittatus' as "Endangered". We discuss all of this, along with model performance and caveats of the study below.

Evaluating model performance and variable contributions
MaxEnt has been shown to be sensitive to different extents/background areas when constructing Figure 3. Relative niche suitability for Phyllobates vittatus across its distribution in the South Pacific of Costa Rica. The suitability map was generated during 2020 through the combination of eleven different environmental predictors, with elevation, forest percentage, distance to lakes and distance to ASADAS explaining the greatest proportion of the variance. The legend shows niche suitability ranging from low suitability (0; white) to high suitability (1; dark green). We represent in white high-altitude areas (>1500 m) and large plantations identified during the classification of the land cover, which were not included in the model. This prediction was made with a model built using a 5 km buffer area around the known occurrence points. Large protected areas are represented by their management category: 1: Corcovado National Park; 2: Piedras Blancas National Park; 3: Golfo Dulce Forest Reserve; 4: Paso de la Danta Biological Corridor. Sources: Costa Rica Atlas 2014, Costa Rica Conservation Areas National System (SINAC) 2016 and 2020. This map is at a 1:1,150,000 scale and CRS WGS84. models (Barbet-Massin et al., 2012;Elith et al., 2011;VanDerWal et al., 2009). Here too, we find that the extent from which the background points were sampled affected the performance and predictions of the models. Following Lobo et al. (2008) and Hijmans (2012), we opted for the model that performs well statistically based on the AICc. The model built from the total extent was deemed to be the most appropriate based on this evaluation.
Elevation, forest percentage, distance to lakes and rivers, annual temperature range and precipitation variables contributed to characterize the ecological niche of P. vittatus according to the most suitable model. Contrary to what we expected, forest cover and distance to rivers are not the two main drivers. Forest percentage is still one of the major ones and, although distance to rivers is not one of the major ones, it still contributes to characterize its niche. Elevation turned out to be one of the major abiotic factors influencing habitat suitability. Both climate and biota in Costa Rica change rapidly with altitude. In the south Pacific, the lowland rainforest gives way to the pre-montane and cloud forest, which typically have lower temperatures and in the latter, a constant mist that increases humidity (Kappelle, 2016a). It is also notable that highly suitable areas do not correspond to low-lying flat areas (see fig. 3 and  supplementary fig. S3A). This relates to other abiotic drivers, like distance to lakes and distance to rivers. These two variables seem to indicate that P. vittatus avoids large streams or rivers that flow into lakes or directly into the sea. The three abiotic factors together indicate that P. vittatus prefers areas of minor elevation where narrow creeks and streams occur. This is congruent with the habitat in which it has been typically found, as well as with its breeding biology (Leenders, 2016;Savage, 2002). Moreover, this would also fit the positive relation that the slope (o.7% contribution) might be having on characterizing the niche of the species, as narrow creeks and streams occur in areas with greater slopes.
Another main driver of habitat suitability for P. vittatus is increased forest density. We believe that it was important to conduct an additional identification of palm, pineapple and teak plantations with orthophotos when classifying land cover as these densely vegetated areas would have otherwise been erroneously classified as forest, which might have undermined the relative importance of dense forest vegetation. Forest percentage is one of typical metrics of landscape composition (Arroyo-Rodríguez et al., 2020), which are used to establish areas of conservation interest and in conservation planning (de Mello et al., 2016). This metric can help us to evaluate forest changes, which are important to monitor in areas under strong anthropogenic pressures. However, we do not recommend ignoring other relevant aspects of forest habitats nor use it as the only metric for habitat conservation.
From the climatic variables used, annual temperature range (bio 7), precipitation of warmest and coldest quarters (bio 18 and bio 19, with a contribution of 0.9% and 2.4%, respectively) and precipitation seasonality (bio 15, contribution of 3.9%) were significantly contributing variables to the modelled ecological niche of P. vittatus. We expected precipitation during the rainy season (coldest quarter) to play a significant role since this corresponds to the peak of the reproductive season of P. vittatus is during the rainy season (Silverstone, 1976). Our model also suggests that precipitation of the warmest quarter (dry season; Kappelle, 2016b) and high temperatures, might also be limiting suitability for P. vittatus. Low levels of precipitation during this dry period result in the drying out of the small streams this species prefers, forcing them to migrate and/or increase intra-and inter-specific competitive interactions. During the dry season of 2019, we observed P. vittatus individuals using the shelter that two leaking artificial water tanks (for human activities purposes) provided, as the nearby streams were completely dry, and sharing it with Dendrobates auratus. In areas of lower precipitation, P. vitattus may therefore also experience an increase in intra-and inter-specific competitive interactions. These variables already pinpoint a potential concern for the conservation of P. vittatus in the future, especially since climatic models forecast increased dryness in Central America, with longer and dryer summers (Aguilar et al., 2005;AlMutairi et al., 2019;Rauscher et al., 2008). Unfortunately, the effects on Costa Rica's biodiversity are already noticeable, with amphibians as one of the main groups beginning to be affected by this threat as global warming favours the spread of the pathogenic chytrid fungus (Batrachochytrium dendrobatidis) (Karmalkar et al., 2008;Pounds et al., 2006).

Model prediction and conservation implications
Conservation action requires reliable ecological data. Yet for many tropical amphibian species, such information is scarce. Here we constructed an updated land cover map of South Pacific Costa Rica and constructed ecological niche models to weigh the importance of different habitat components for the occurrence of a threatened and endemic species of poison dart frog of Costa Rica, Phyllobates vittatus. Our model prediction projected over the study area ( fig. 3) indicates multiple areas of high relative niche suitability for P. vittatus. The largest continuous area of high suitability is within the eastern side of the Osa Peninsula. Although there are few pixels highly suitable (over 0.8) within Corcovado National Park, the largest area with high suitability seems to be within the Golfo Dulce Forest Reserve. The Osa Peninsula is a key hotspot for biodiversity (Cortés and Jiménez, 1996;Kappelle, 2016b;Quesadaalpízar et al., 2006). It is listed as 'essential life support area' (ELSAs, UNDP 2020), and is clearly vital for the conservation of P. vittatus. High suitability areas can also be seen in Piedras Blancas National Park and its surroundings. However, not many occurrences have been recorded in this area. We recommend future surveys within Piedras Blancas National Park and its neighbouring areas to confirm putative populations of P. vittatus.
Outside the Osa Peninsula, highly suitable areas can be found within the biological corridor Paso de la Danta. This corridor seems to be of high relevance for P. vittatus in the northern part of our study area, where it is the only area with high suitability for the species. This corridor is also considered an ESLA region 'essential life support areas' (UNDP, 2020). In the Terraba's river valley, where recent occurrences have been recorded, we find a small area with moderate predicted suitable niche (0.4-0.6). Despite not being part of any large protected area (Sistema Nacional de Áreas de Conservación de Costa Rica (SINAC), 2020), this valley is considered an ESLA area to restore (UNDP, 2020). Moreover, it could be connected with another area with high suitability for P. vittatus, the biological corridor Paso de la Danta. In the area around Punta Burica there is an area with high suitability for P. vittatus. Populations have been confirmed in this area (Mebs et al., 2014) and we suggest further surveying efforts to search for additional populations, especially since this area could have become isolated by the growth of surrounding palm plantations. These plantations could mean another potential threat, as Costa Rica is a major user of pesticides in agriculture (Claydon, 2017;Alvarado, 2018), and monitoring programs in the country have shown how these pesticides run off into the nearby streams and leave the aquatic-dependant organisms exposed to a high toxicity (Carazo-Rojas et al., 2018;Echevarría-Sáenz et al., 2018;Pérez-Villanueva et al., 2022;Weiss et al., 2023). This potential threat has not been identified within the assessment by the IUCN (IUCN SSC Amphibian Specialist Group, 2020), but it could be a major threat for P. vittatus and other amphibians which niche is surrounded by plantations. Furthermore, there are currently no government protected areas within this region in either Costa Rica or Panama (Sistema Nacional de Áreas de Conservación de Costa Rica (SINAC), 2020; Sistema Nacional de Áreas Protegidas de Panamá (SINAP), 2005), and apart from a very small area considered restored (UNDP, 2020), there seems to be no planned conservation actions in the area. The fact that this area spans across country borders may also complicate the declaration of protected areas for conservation purposes. Considering the species was unknown in Térraba and Punta Burica until recently (2020), and that they are outside of protected areas, these populations might be the most vulnerable. If isolated, the genetic diversity of these populations could be affected by inbreeding and drift, increasing their probability of local extirpation. Further studies regarding population connectivity and gene flow between populations is needed to get a better insight about these populations and the best measures to protect them.
Based on the number of 2 km 2 grids with at least one occurrence record, the area of occupancy (AOO) of P. vittatus is 106 km 2 . Species with areas below 500 km 2 are classified "Endangered" by the IUCN Red List (IUCN criteria B2). In this study we showed how the populations outside of the Osa Peninsula seems to be extremely fragmented. Furthermore, there are extreme fluctuations in the subpopulations and there seems to be a decrease in the habitat quality (IUCN SSC Amphibian Specialist Group, 2020). Hence, these conditions of the criteria B of the IUCN Red List are sufficient to reassess the status of P. vittatus and catalogued it as "Endangered" and not "Vulnerable".

Conclusions
The use of citizen science in biodiversity monitoring has grown in the last decades, especially via the use of modern technology applications that help facilitate participation and ease of use. This serves not only as a tool to expand biodiversity data but to also engage people with nature conservation efforts (Pocock et al., 2018). However, it is necessary to remark that previous training on the use of 'apps' is vital to ensure the correct recording of the encounters, as well as post-processing to ensure data quality. Despite the caveats of this study that could influence the accuracy of our models, such as having to use unsupervised land cover classifications or relying on the climatic data just for 1970-2000, the citizen science data permitted us to expand our occurrences' database and discover locations that were not previously confirmed.
This work is a much needed first step towards a better understanding of the habitat requirements of P. vittatus. It has given us a deeper insight into the main abiotic factors that render a site suitable for P. vittatus. Additionally, we find that a few known populations could be potentially isolated under habitat encroachment in small forest fragments, and uncover some potentially suitable areas where its presence has not yet been verified. These isolated populations may be particularly vulnerable to other potential threats like climate change and/or chytridiomycosis.
Furthermore, we also analyzed the AOO of P. vittatus, showcasing how P. vittatus complies three conditions of the criteria B in the IUCN Red List, indicating that this species should be listed as "Endangered". Future surveys on these areas to confirm the species presence, periodic monitoring of the populations and genetic studies to analyse populations' connectivity are recommended to have a proper and updated baseline to create suitable conservation measures to protect an endangered species that symbolizes the endemism of an iconic region.

Acknowledgements.
We thank to all the citizen scientists that collaborated in this study and all the stakeholders that joined it. Special thanks to M. López Morales for the logistics and fieldwork support across the Osa Peninsula. Special thanks to M. Sánchez, M. López Morales, R. Núñez E., J. González, P. Rojas, F. Protti-Sánchez, C. Barrio-Amorós, R. De Plecker, J. Vargas, and G. Köhler for sharing their occurrences of Phyllobates vittatus and to all the iNaturalist users for sharing the coordinates of their observations. Thanks to M. Monge for helping to raise funds to carry out this study. SierpeFrogs, Tamandua Biological Station, Tropenstation La Gamba, the communities of Dos Brazos de Río Tigre and Rancho Quemado, the local guides associations of Aguilosa and Aguinadra, Lapa Rios Rainforest Lodge, Nueva Tierra de Osa, Danta Corcovado Lodge, Ecoturístico La Tarde, Golfo Dulce Retreat, and Saladero Lodge logistically supported this study. We thank the financial and technical support of the Rufford Foundation, Alongside Wildlife, Stiftung Artenschutz and the Amphibian Conservation Fund by Zoos and private participants in the German-speaking region and Osa Conservation. We also thank the reviewers for their insightful and thorough comments, which have helped to greatly improve this manuscript.