Promising abundance and composition of forest dwelling anurans in cashew plantations in a tropical semi-evergreen forest landscape

Plantation crops in tropical human-modified landscapes provide alternative habitat to biodiversity outside protected areas. The Western Ghats of India is home to a mosaic of closely spaced habitats, including forests and agroecosystems. Cashew is a widely grown plantation crop in the northern Western Ghats and is known to provide economic and societal benefits. However, its role as a supplementary habitat for anurans is not well understood. We assessed the factors that influence understorey anuran composition and abundance in cashew plantations, forest edges, and forest interiors in Tillari Conservation Reserve (TCR), Maharashtra. Species composition of cashew plantations differed significantly from forests and was positively influenced by understorey and canopy cover. Cashew plantations had near equal abundance of anurans as that of forest edges and interiors, which could be due to the preponderance of habitat generalists. Understorey positively influenced anuran abundance while ambient temperature had a negative influence. Reduced understorey and low canopy cover represent habitat modifications that occur in cashew plantations. Such structural changes could lead to reduced environmental refuges for anurans, thereby exposing them to large variation in temperature and moisture. Cashew plantations in TCR serve as supplementary habitat for anurans. That said, cashew cultivation practices and markets must be understood before biodiversity-friendly plantation practices are proposed.  

Methods
2.1 Study area

The study was carried out at Tillari Conservation Reserve (TCR) in Sindhudurg District, Maharashtra state, which falls under the Northern Western Ghats, India (Figure 1a). The region’s elevation ranges between 50 m and 1,030 m a.s.l. It receives an average annual rainfall of 3,000 mm. Almost all the rainfall is received during the monsoon, which spans from mid-June to October. It has tropical semi-evergreen and moist-deciduous forest types (Champion and Seth 1968). Besides these natural habitats, TCR is part of a human-modified landscape, that includes hydroelectric projects, human habitation and secondary evergreen forest. Cash crops include: cashew, banana, pineapple, and rubber. Cashew plantations cover roughly 537 km2 in Sindhudurg district, and that constituted 31 % of all area under cashew cultivation in the entire state of Maharashtra (Sengar, Mohod, and Khandetod 2012). In the landscape comprising of forests, plantations, reservoir and human habitations, the state owned primary and secondary forests that are part of TCR encompasses 29.53 km2 (Patil 2021). 

The sampling for anurans was conducted across forest interiors, forest edges, and cashew plantations. The cashew plantations are typically small land holdings, locally owned and predominantly characterized by 2–3 m-high cashew trees spaced 15 ft apart. They are interspersed with a few remnant native trees such as Terminalia paniculata, Mangifera indica, and* Artocarpus* heterophyllus, and native trees are occasionally logged. Cashew requires sunlight and is therefore grown primarily on gentle hill slopes in TCR. The age of the plantations in the area ranged between 7 and 10 years. The understorey of the plantations is cleared twice a year, in summer and after monsoon. No irrigation is provided to the cashew plantation during the dry period. The sampling was done in cashew plantations owned by several people of the Hevale village and prior permission for this work was obtained from land owners and village committee. The forest edge transects were laid at a distance of not more than 1.5 km from the bordering cashew plantations. Seven of the ten forest edge transects were present in a Sacred Grove. They are protected by the local communities and no modification to the understorey is made. The remaining three transects were laid in habitats that were not sacred groves but were bordering cashew plantations. The transects in forest interior habitats were laid about 2 km from the plantation border inside TCR The forest interior sites had no recent anthropogenic disturbance. However, locals practiced shifting cultivation and selective logging about four decades ago. 

2.2 Sampling anurans and environmental covariates

We laid thirty 50 × 2 m transect strips, with ten transects in each of the three habitats. Sampling was limited to a maximum of 3 m vertically. Each transect had a minimum of 50 m distance between them to avoid double counting. We employed a nocturnal Visual Encounter Survey (VES) for sampling anurans (Crump and Scott 1994), wherein, two people with headlamps searched for terrestrial and arboreal anurans. We sampled only on land and avoided aquatic habitats and canopy. Each transect was thoroughly scouted for anurans after sundown for 1 hour, typically between 1900–2100 h. Since we sampled anurans during the period when they were most active, we made the assumption that their passive encounters were representative of their relative abundances. Anurans were sampled from August 2020 to October 2020 during the monsoon. We conducted a pilot survey during the monsoon of 2019 and made an inventory of anuran species in the study area. Six out of nine species observed during this study were taxonomically identified based on our 2019 effort. The anurans with ambiguous taxonomic identity were carefully examined and identified using DNA barcoding by matching partial 16S rRNA sequences obtained from oral swabs (Kurabayashi et al. 2005). In addition, we recorded microhabitats used by anurans which included, bare ground, leaf litter, and bush. 

Environmental covariates such as elevation, canopy cover, percent understorey, leaf litter depth (LLD), and air temperature (Ta) were measured to determine their contribution in influencing the anuran composition and abundance across sampled habitat types. Each variable was measured at every 10 m from the beginning to the end of the transect and then averaged. The elevation was measured in meters above sea level using a Garmin-Etrex 10 handheld GPS with an accuracy of 5 m. A vertical densitometer was used to measure the percent canopy cover. The understorey was estimated by measuring the amount of 1–2 m-high vegetation in a 10 × 2 m plot inside the transect. Leaf litter depth was measured using a steel ruler inserted through the litter until it touched the soil or humus layer following Rowley and Alford, 2007, and the air temperature was measured using a digital thermometer, Extech™ RH101.

2.3 Measuring operative body temperature and prey availability

Operative temperature (Te) is the temperature attained by a non-thermo-regulating animal from radiation, conduction, and convection. It was measured by placing a physical model or operative temperature model (OTM) mimicking the size and shape in the microhabitat typically used by the species (Taylor et al. 2020). In this study, the OTM was made of 3% agar gel shaped into a frog of 6 cm in size weighing 10–13 g using a gypsum frog cast, mimicking the body size of an adult frog and covered by a plastic mesh (Figure 1b). The models had a DS1921G Thermochron® iButton temperature data logger embedded in them to record hourly temperature data in degrees Celsius. 

We installed pitfall traps coupled with drift fences guiding the pitfalls to measure the prey availability. Invertebrates found in the pitfall were identified to order level and prey availability was measured by counting the number of insect orders. The traps were made of 5 cm diameter cylindrical bottles with an inverted conical top to avoid escape of the trapped insects and small holes at the bottom to avoid waterlogging. The models and pitfalls were placed from 1800–0700 h in the three microhabitats: bare ground, leaf litter, and bush. We placed the models and pitfalls for two consecutive days in every transect for all habitats except the forest edge, where the Te and prey availability was not measured in one microhabitat i.e., leaf litter, as we did not find anurans in the leaf litter in this habitat. The model temperature and prey availability values were averaged at the transect level. We checked the pitfalls every day to avoid water logging or damage by animals. 

2.4 Statistical analysis

The abundance of anurans was the count of anurans of a species encountered in a transect in the three habitats. A Spearman correlation test was used to check for correlation among environmental covariates variables (Figure S1). We used Moran's I test to check for spatial autocorrelation (Kurz et al. 2016). To test for differences in the environmental covariates between habitats, we used a Kruskal-Wallis test coupled with effect sizes followed by a Dunn’s pairwise comparison test. We used non-parametric tests as some variables collected during the study (canopy cover, ambient temperature, and prey availability) did not follow normal distribution. These tests were performed using kruskal_test, kruskal_effsize, and dunn_test functions from the rstatix package (Kassambara 2020), and Mean ± SD was calculated. Generalized linear mixed models (GLMMs) with Poisson distribution and log link function were used to assess the influence of environmental covariates on anuran abundance (Bolker et al. 2009). We used the glmer function from the lme4 package (Bates et al. 2015) to perform this analysis. Habitat type, canopy cover, understorey, LLD, Ta, Te, and prey availability were used as fixed effects, and intercept varying among transects was used as a random effect. Elevation was removed for this analysis as it was positively correlated with canopy cover (ρ > 0.7, P < 0.001). The continuous fixed effect variables were scaled using the function scale in package base by dividing the values by standard deviation before performing the GLMM (R Core Team 2019). We used the dredge function from the MuMIn package (Barton 2020) to run all possible models. We performed model averaging instead of a single parsimonious model following Burnham and Anderson (2002), by selecting candidate models with ΔAICc<2 and using the function aictab from the AICcmodavgpackage for model averaging. A total of 128 candidate models were built by the dredge function, of which 6 models were selected for averaging. Following this, we extracted the model average predictions for all variables using the modavgfunction. We used Permutational multivariate analysis of variance (PERMANOVA) to assess the patterns in anuran species composition across habitat types using the adonis2 function (Anderson 2001). Bray-Curtis dissimilarity measure was computed from the species frequency data matrix using the vegdist function. Non-metric Multidimensional Scaling (NMDS) was performed using the metaMDS function to visualize species composition patterns across habitat types (Borcard, Gillet, and Legendre 2011). Following NMDS, the envfit function was used to fit environmental covariates as vectors to the ordination. Species composition analyses were performed using the vegan package in R (Oksanen et al. 2020).

Further, to characterize the value of the habitats, we conducted an indicator species analysis (ISA) to identify species that select a particular habitat type (De Cáceres, Legendre, and Moretti 2010). The analysis generates an indicator value (IV) for every species based on the specificity (A = mean abundance) and fidelity (B = relative frequency of occurrence). An IV of a species is a predictive value of a species as an indicator for a site or a combination of sites (De Cáceres, Legendre, and Moretti 2010). The ISA was carried out using a multipatt function with 999 permutations from the indicspecies package (De Cáceres, Jansen, and Dell 2020). All the analysis was conducted using R version 3.6.2 and plotted using the package ggplot2 (Wickham 2016).