Ecology and morphology of the dwarf bromeliad boa Ungaliophis panamensis (Squamata, Boidae, Ungaliophiinae) in Costa Rica and Panama

Ecological and morphological data on Ungaliophis panamensis is extremely limited as this species is rarely encountered. These knowledge gaps have been advanced in this study where data was analysed from a small sample of snakes collected in two tropical forested environments in Costa Rica and Panama. Standardised major axis testing and a Bayesian latent variable ordination revealed that the species is sexually dimorphic, closely associated with tree trunks in natural forested areas, and occasionally discovered in rural buildings. Although further investigation into its natural history is warranted, this study shows that even with just a few individuals it is possible to elucidate ecological information that is relevant to the conservation of snake species.


Introduction
The natural history of tropical snakes is often little understood. This is due in part to their cryptic behaviour and remarkable camouflage, which results in infrequent detection and is a frequent issue for snake biologists (Durso and Seigel 2015;Shelton et al. 2018). The genus Ungaliophis consists of two species that are restricted to tropical Mesoamerica; they are considered highly cryptic and very little data exists regarding their natural history (Bogert 1968;Savage 2002;Köhler 2008). The northern species, U. continentalis, is found from southern Mexico to northern Nicaragua, whereas the southern species, U. panamensis, occurs from southern Nicaragua to northwestern Colombia (Villa and Wilson 1990). The genus was first described in 1880 by Müller who designated an individual from Guatemala as Ungaliophis continentalis (Savage 2002). It was not until 1933 that Schmidt described the species U. panamensis. Furthermore, U. panamensis was not confirmed from Costa Rica until 1974 (Corn 1974).
Since then, very little has been reported about the ecology and natural history of either species. The little data that exists regarding the natural history of either species of Ungaliophis strongly suggests that they are arboreal species found from lowland rainforests to high elevation cloud forests (Köhler 2008). Corn (1974) found several individuals of U. panamensis in banana plants and bromeliads. As such, many authorities believe that Ungaliophis naturally occurs in dense epiphytic growth in the tree canopy. Both species of Ungaliophis are small snakes capable of reaching between 470 mm and 482 mm in total length (Köhler 2008). Historically, the diet of wild Ungaliophis was unknown; however, several recent publications regarding the diet of U. panamensis in Costa Rica have included observations of bats (Solórzano and Carillo 2017), birds (Dwyer 2017), and geckos (Espinoza and Barrio-Amorós 2018) being taken by the species. Ungaliophis continentalis have been recorded to give birth to between 5 to 6 neonates per clutch in captivity (Ross and Marzec 1990;Burger and Ford 2007). While there is no information regarding reproduction of U. panamensis, it is likely to be similar to its sister species.
Identifying the habitat preferences of a given species is critical to understanding the extent of functional habitat available to that species in the wider landscape (Ocampo-Peñuela et al. 2020). Effective conservation strategies can be built once the available habitat for a species, or assemblage, has been quantified and understood (Riva and Nielsen 2020). Herein, we describe the microhabitat preferences and morphometry of 13 individuals of U. panamensis, ten of which were encountered in Barra del Colorado Wildlife Refuge in northeastern Costa Rica, and three from the Coclé and Panamá Oeste provinces of Panama. We hypothesise that U. panamensis has an association with natural forest habitats and that its morphology does not influence its habitat preference.

Study area
Individuals of U. panamensis were recorded across multiple locations within the Barra del Colorado Wildlife Refuge (BCWR) in the NE region of Costa Rica (Fig. 1), and from Coclé and Panamá Oeste provinces, Panama (Fig. 2). The BCWR region is dominated by lowland tropical wet forest (Holdridge 1967) and Manicaria (Arecaceae) swamp forests (Myers 1990;Lewis et al. 2010). The regional habitat in both Coclé and Panamá Oeste, Panama, is mid-elevation mature secondary cloud forest and primary premontane cloud forest (Ray et al. 2012;Ray 2015).

Field methods
Data collected for this study involved a pool of surveys that were performed from multiple longer term studies spanning 1997-2012 (Ray 2009;Lewis et al. 2011;Lewis et al. 2013;Zipkin et al. 2020). Most individuals were detected using visual encounter surveys in forested environments (Heyer et al. 1994;Lovich et al. 2012). However, some were observed by chance encounter in ceiling rafters of buildings and others reported to us by local people. Surveys at both sites were conducted from November 2005 to September 2008, between 20:00 and 02:00 in a variety of habitats: canal edge areas, deep forest and riparian gallery forest. Surveys averaged three surveyors and were conducted four nights a week for up to six hours per survey.
Biometric data comprised sex determined by careful probing and examination of the anal spurs, snout-vent length (SVL, mm), tail length (TL, mm) and mass (g). Microhabitat data comprised observations where the snake was located. If the individual was in Forest, specific substrate classes were recorded (Tree, Palm, Shrub), along with relevant structural features (Trunk, Leaf, Branch, Twig). Although infrequent, when snakes were found in rural environments, often buildings, they were assigned their own identity (Buildings).

Data analysis
Encounter rates of U. panamensis at both sites were calculated by dividing the number of individuals by the number of surveys and multiplying by 100, a technique adapted from Kiszka et al. (2007). For the purposes of exploring patterns of association between the morphology of U. panamensis and its microhabitat preferences we considered individuals to be the sampling unit.
The morphometric relations of U. panamensis were analysed using Standardised Major Axis (SMA) estimation created by the function SMA within the package smatr within the program R (Warton et al. 2012a; R Core Team 2020). An advantage  of using SMA to test for common slope allometric relations is its inclusion of fitting factor groups for y against x (Warton and Weber 2002). We also utilised Huber's M estimation in place of a least squares method using the robust function in smatr in order to make the fit more inclusive for marginal outliers (Taskinen and Warton 2011). Likelihood ratio tests evaluated the slope fit.
We chose to use a multivariate Bayesian, instead of univariate or distance based, approach for ecological analysis. This was because transforming data to meet assumptions for univariate and distance based approaches is problematic for smaller data sets (Warton et al. 2012b;Warton 2017). Microhabitat data were therefore analysed using a Bayesian unconstrained ordination based on a latent variable model (LVM) (Hui et al. 2015) using the package Boral (Hui 2016). Bayesian LVM's are useful at explaining multivariate composition while accounting for residual correlation. They are superior to non-metric multidimensional scaling (NMDS) because they make provision for possible mean-variance relationships in data without confounding location with dispersion (Warton et al. 2012b;Hui et al. 2015). Boral fits three types of model; 1) Covariates with no latent variables (fitting independent response GLMs where columns of y are assumed independent; 2) With no covariates (a pure LVM that constructs unconstrained ordination); and 3) Combined covariates with latent variables (fitting correlated response GLMs, with latent variables). We explored U. panamensis microhabitat preference using modelling option 2. The model comprised the following; Option 2 (pure LVM), where µ ij is considered the mean response at microhabitat level i for an individual snake j, Ɵ 0j is the individual-specific intercept, ȥ i = (ȥ il , ȥ i2 ) T is a vector of two latent variables, and Ɵ j = (Ɵ j1 , Ɵ j2 ) T are the corresponding individual-specific coefficients. This modelling approach enabled biplots to visualise the data in a similar way to a two-dimensional NMDS. From the model we extracted posterior median values of latent variables and used these as coordinates on ordination axes to plot microhabitat association (Hui et al. 2015). For Boral models, estimation is performed using Bayesian Markov Chain Monte Carlo (MCMC) methods via software JAGS (Plummer 2003). Our model comprised two latent variables, used a negative binomial family, ran 40000 iterations with 10000 discarded for burn-in, and was thinned by 30. Priors were set using Boral's modest automated uniform normal distribution detected and set through JAGS (priors = ~ dnorm (0,0.1)). Convergence was assessed using MCMC trace plots retrieved from package Coda and inspection of Geweke convergence diagnostic (Geweke 1992), a test which is similar to the Gelman statistic potential scale reduction factor (PSRF) (Gelman et al. 2013), but applicable given Boral operates with only a single MCMC chain (Hui 2016). Model assumptions of mean-variance and log-linearity were examined using Dunn-Smyth residuals vs. fit plots and normal quantile plots (Dunn and Smyth 1996).
Correlation and residual correlation were checked by plotting a residual covariance matrix of latent variable regression coefficients using function get.residual.cor in Boral and package Corrplot. A strong residual covariance/correlation between factor variables can be interpreted as evidence of autocorrelation in a model; however, acceptable levels have been recognised as indicative of an interaction/association (Pollock et al. 2014). The residual precision matrix in Boral can be used to directly identify association between factors (in our case microhabitats and individual snakes) (Ovaskainen et al. 2016). Two factors exhibiting a zero result in such plots may remain correlated, indicating they do not directly interact, but also can remain correlated through other factors. Residual precision matrix results should not exhibit elements equal to exactly 1 or -1 (suggesting strong autocorrelation). Nevertheless, relatively large values between these limits of precision imply a useful indication of a correlated relationship between two factors.
All analyses were carried out in the program R version 4.0.0. (R Core Team 2020) which we attached to JAGS version 4.3.0 for performing Bayesian routines.

Results
Our dataset comprised microhabitat and morphological data for a small number of individuals (N = 13) ( Table 1). In BCWR 258 surveys over an 8-year period yielded 10 individuals (an encounter rate of: 3.88) (Fig. 3). In Coclé and Panamá Oeste 1107 surveys over a 15-year period yielded 3 individuals (an encounter rate of: 0.27) (Fig. 4). This highlights just how poorly detected the U. panamensis is in the wild. Eight were found during surveys, three by random encounter, and two were brought to us by locals. A small specimen encountered, from Costa Rica (possibly juvenile) was removed from the analysis to avoid bias in the models. This snake measured 162 mm SVL, 29 mm TL, and weighed 2.5 g, and was encountered on a palm frond at a height of 76 cm. Inspection of morphometric data via boxplots showed a difference in the means of SVL and TL between males and females. However, boxplot distribution overlapped, indicating that such mean difference may not be necessarily statistically significant. Adult male snakes were longer in both SVL and TL, but adult females had the greatest mass despite mean differences in mass being equal (Fig. 5, Table 1). SMA showed that variables grouped by sex were uncorrelated; Females (R 2 = 0.67, P < 0.05), Males (R 2 = 0.73, P < 0.05). SMA slopes were unequal (LRT 0.48, P = 0.48) indicating a sexual difference in body size. The SVL/Mass [log scale] plot revealed females were subtly larger than males due to having greater mass (Fig. 6). Slope fit was checked via residuals and q-plots and deemed acceptable given the small data size (Fig. 6).
The LVM in Boral successfully computed with useful convergence (see: Suppl. material 1: Convergence trace plots). The Option 2 model (correlated response) resulted in a fairly low DIC (122.6) and a sufficient residual trace correlation test result (104.79). All Geweke Z-score p-values for both models were >0.05 (results of <0.05 is approximately equivalent to a PSRF of 1.96 which would exceed MCMC Geweke upper values and indicate poor convergence). Residual plots of model fit showed good distribution of linear predictors indicating minimal over-dispersion (Fig. 7).
The primary latent variable ordination plot for the model showed preference by U. panamensis for natural microhabitat features and clustered latent variables proximate to corresponding snake microhabitats (Fig. 8). Whilst buildings were utilised by five U. panamensis individuals as refuge and foraging sites, they were distally associated in the ordination with those snakes detected in natural microhabitats. It is important to note that the buildings these five individuals were found in were situated in very remote areas, surrounded by natural habitat and not within higher density urban areas. The plot also showed discrimination of Tree, Trunk and Leaf       variables, and Branch, Twig and Palm variables to the right-hand side of the ordination amplifying their relations as a group away from Buildings. Forest and Shrub were categorically assigned in the model and, although not tightly associated with the clustered microhabitats, are spatially situated away from variable Building, thus adhering them to the right-hand side of the plot. This indicates their grouping with natural microhabitats. The residual correlation plot showed most variables with a weak-minimal residual correlation (Fig. 9). Therefore, we accepted model correlation had not compromised the model and interpreted correlated signals as underlying relationships. Tree and Trunk were the strongest positively correlated variables. Building and Forest were strongly negatively correlated which also aligns with the contrast reflected in Fig. 8.
Sex, SVL, Mass, and TL were plotted against the latent variables for both models but exhibited no discernable relation to microhabitats with individuals distributed liberally between variables and associated to a spread of sex and size among the plots.

Discussion
This study provides a unique insight into the ecology of a little known and understudied snake and confirms general assumptions that have been reported in the literature (Corn 1974;Villa and Wilson 1990;Savage 2002). We have shown that like other boid species, U. panamensis exhibits sexual dimorphism, and that females are, in general, shorter than males. Our data confirm the hypothesis that U. panamensis has a significant preference for natural forest habitats, although it will use rurally located buildings in proximity to natural forest cover as refuges, or possible hunting grounds, a finding that is consistent with recent observations (Dwyer 2017;Solórzano and Carillo 2017). In the current study, given the forest habitat in the immediate area surrounding rural buildings where we encountered U. panamensis, it is likely these individuals may not have normally chosen to utilise anthropogenic structures.
The strong positive correlation between Tree and Trunk is an ecologically useful descriptor as this species is strongly associated with being found on tree trunks when detected in forested environments. The high correlation between Tree and Trunk variables identifies a hidden relationship between these two variables in the data that is not initially clear within the primary ordination plot. When buildings were present, they signalled a negative correlation with most natural microhabitat variables confirming that natural microhabitat components feature as the habitat components of choice by U. panamensis. These correlations between habitat variables in relation to the presence of U. panamensis are also observed for other arboreal boid species in the region such as Corallus annulatus (Lewis et al. 2011). Whilst such correlated habitat variables might not be unexpected, they have, prior to the work herein, not been adequately described for U. panamensis.
Eco-morphology and microhabitats are interesting concepts to describe a reptile's habitat preferences (Bickel and Losos 2002;Pyron and Burbrink 2009;Lewis et al. 2013). These preferences could be used to guide conservation management when considering options for land management (Berriozabal-Islas et al. 2017;Todd et al. 2017). Our data was collected from only 13 individuals, over regular surveys during a ten-year period, which attests to the species' cryptic behavior, remarkable camouflage, and low detection rate. Prior to techniques such as latent variable modelling, deciphering interrelated associations between habitats for a species from small data sets often amounted to anecdotal information and univariate associations. This kind of information often rendered a species' preferred microhabitat structure inadequately addressed and at the mercy of descriptive natural history notation and surveyor opinion.
It is well documented that detection rates in snakes hamper efforts to better understand their ecology and conservation needs; in some cases decades are needed merely to understand the extent of the snake assemblage at a given location (Doan and Arizábal 2002;Duellman 2005). This difficulty in detection historically has led to a paucity in ecological data for many snake species, particularly in the tropics, due to the low numbers of encounters (MacKenzie et al. 2003;Bailey et al. 2004;Durso et al. 2011;Stroud and Thompson 2019). For example, Corallus annulatus, Calliophis bibroni, and Trimeresurus medoensis are three tropical snake species so infrequently detected that their detection is of conservation and distributional interest whenever they are discovered (Lewis et al. 2011;Griffin et al. 2012;Arockianathan et al. 2014). We believe this study shows that even with just a few individuals it is now possible to elucidate ecological information that is relevant to the conservation of snake species for which that information is currently lacking. Our analysis provides evidence that although the species can be found in rural environments such as buildings situated close to forested environments, it ultimately is a forest dwelling species that is detected well on tree trunks and among lower foliage on trees. With such low detection rates, we expect future research and observations of U. panamensis to further the knowledge of the ecology of the genus substantially. This novel approach to small datasets also provides the opportunity to better inform conservation and habitat management decisions in locations inhabited by rarely detected species.