Population ecology of the spectacled caiman ( Caiman crocodilus ) in the Apaporis River middle basin

. Population ecology studies on spectacled caimans ( Caiman crocodilus ) in Colombia have been few and far between with many covering short periods and de ﬁ ning population parameters based on relative indices (i.e., individuals/km). This re ﬂ ects a lack of information on the general effects that environmental variables have on annual cycles of population dynamics, as well as a bias in abundance estimations due to the uncertainty of detection error. Keeping this in mind, we assessed the abundance and demographic structure of the spectacled caiman population inhabiting the Apaporis River middle basin over a year, based on robust hierarchical model that accounts for imperfect detection. We recorded a total of 1156 caiman observations between December 2018 and November 2019, estimating an average predicted value for abundance across all surveys of 29.99 (cid:1) 13.17 individuals, slowly increasing as the transect length increases and increasing variation as months passed by. The average detection probability was 0.69 (cid:1) 0.25 across all surveys, with no apparent effect as water temperature and relative humidity change across space-time and slowly decreasing as months go through. The population size estimated based on the top-performing model was 1763 (cid:1) 786 caimans across ~ 7.1 km 2 assessed. We estimate the commonly used relative abundance (encounter rate) index as well as a generalized linear model and discuss how those relate with the values predicted by N-mixture models. We also discuss the relevance and cautions researchers should have when using N-mixture models to better understand spectacled caiman ecology.


INTRODUCTION
The spectacled caiman (Caiman crocodilus, Linnaeus 1758) is one of the most widespread crocodylian species in the Neotropics. It naturally occurs from Mexico to Peru on the Pacific coast and from Mexico to Brazil on the Atlantic/Caribbean between 0 and~1000 m above sea level (Medem 1981, 1983, Balaguera-Reina and Velasco 2019. This species has also been deliberately introduced in countries such as the USA, Cuba, and Puerto Rico (Ellis 1980, Rodríguez-Soberón et al. 1996, Rivero 1998. In Colombia, C. crocodilus occurs in a wide variety of freshwater habitats across the Caribbean, Pacific, Orinoco, Andes, and Amazon regions and has also been introduced in several insular areas such as San Andres and Gorgona islands (Medem 1981, Forero-Medina et al. 2006. The spectacled caiman is a highly plastic species and can tolerate anthropogenic intrusion and habitat transformation. This species can be commonly found in man-made water reservoirs, around dams, and fish farms (De La Ossa and De La Ossa-Lacayo 2013). Despite this plasticity, the unsustainable exploitation of C. crocodilus throughout the last four decades of the 20th century caused a dramatic reduction of wild populations in regions such as the Caribbean and Andes, raising concerns about the future of the species in the country (Barahona et al. 1996, Rodríguez 2000. Currently, C. crocodilus is catalogued regionally as Least Concern (Balaguera-Reina and Velasco 2019) and is listed in appendix II of CITES (2019) due to a reduced pressure of wild populations and a more restricted control of farming (Balaguera-Reina and Densmore 2014). However, the lack of a country-level assessments has avoided to clearly define the current conservation status of the species in Colombia.
Population ecology studies on spectacled caimans in Colombia have been highly concentrated on regions such as the Caribbean (Barahona et al. 1996, Ulloa and Cavanzo 2003, Cavanzo 2004 (Forero-Medina et al. 2006), and the Amazon (Naranjo 1996, Balaguera-Reina 2019. Most of these studies were performed across periods of time less than a year, which implies they were missing at least part of the story due to the effects of environmental annual cyclical variables (e.g., temperature, precipitation) on population parameters (e.g., number of individuals, demographic structure; Balaguera-Reina et al. 2018).
Historically, the most common method of estimating abundances of spectacled caimans' populations (and crocodylians in general), has been based on relative indices (i.e., individuals/km; Chabreck 1963). However, these indices do not account for imperfect detection (missed/unseen individuals across transects) and do not include covariates that could affect the number of individuals reported by sites (i.e., different habitats across transects) and surveys (e.g., air temperature at the time sampling was done), which have clear effects on the population parameters reported (Kéry and Royle 2016). The use of hierarchical models, specifically, N-mixture models (Royle 2004), has allowed some of those issues to be resolved providing a way to estimate site-specific abundances distributed according to some mixing distribution (e.g., Poisson). At the same time, they can account for the measurement error (detection probability) and thus are a reliable and robust set of methods to account for imperfect detection (Kéry and Royle 2016).
Regarding the Apaporis basin, only two studies have been done since Medem (1955Medem ( , 1981 hypothesized the presence of an endemic spectacled caiman subspecies (C. c. apaporiensis) in the area, both of these were one-time surveys (Naranjo 1996, Balaguera-Reina 2019), leaving uncertainty about the actual conservation status of the species in the area. Herein, we assessed the abundance and demographic structure of the spectacled caiman population inhabiting the Apaporis River middle basin across a year (2019) based on a robust hierarchical model method (N-mixture models) that account for imperfect detection. We also estimated and compared the commonly used relative abundance (encounter rate) index to make it equivalent to most of the studies done on the species, discussing how it relates with the values predicted by models. We use the specific (C. crocodilus) rather than the subspecific (C. c. apaporiensis) taxonomy reported for the area due to recent molecular studies that have shown extremely low genetic differentiation between C. c. apaporiensis and C. c. crocodilus, suggesting that the former lies within the genetic spectrum of C. crocodilus Amazon lineage (Balaguera-Reina et al. 2020).

Study area
The Apaporis River originates at the junction of the Tunia and Ajajú rivers between the Caquetá and Guaviare departments (Dos Santos locality), flowing 1060 km across four departments before emptying into the Caquetá (Japurá) River (Medem 1981, IDEAM et al. 2007). It is a whitewater river with high levels of suspended sediments that result in near-neutral pH, high conductivity, and a pale, muddy color, characterized by the presence of prominent waterfalls, which acts as zoogeographic barriers for aquatic organisms (Medem 1981). The Apaporis River middle basin geologically corresponds to the upper tertiary with quaternary formations exposed around rivers and creeks (Cárdenas-López et al. 2010

Fieldwork
We carried out monthly spotlight surveys (Chabreck 1963) Fig. 1). Differences in transect lengths were due to waterbody extension and relative navigability. Surveys were done between 1900 and 0400 h on new moon nights using a 5000-lumen headlamp (InnoGear 5000 LED, Morgan Hill, California, USA) and a handcrafted canoe powered with paddles. Animals were visually located by the eyeshine reflected when hit by light and slowly approached as close as possible to assure accurate species identification (two other species inhabit the area, Paleosuchus palpebrosus and P. trigonatus) as well as to estimate total length (TL AE 10 cm). Physical variables such as air temperature (AT), water temperature (WT), and relative humidity (RH) were taken with a hygro-thermometer (Extech Instruments, Fotronic Corporation, Melrose, Massachusetts, USA) at the beginning and the end of each transect to define the initial conditions and how those changed throughout the survey (Hutton and Woolhouse 1989). Observations were grouped based on the TL estimated by biologically structured size classes as defined by Ayarzagüena (1983): class I (<50 cm; majority juveniles within its first year after hatched), class II (50-120 cm; subadults), class III (121-180 cm; all females and~20% of males reproductively active [adults]), and class IV (>180 cm; all individuals are adults).

Data analysis
Two sample t tests were performed to assess any significant differences between environmental variables measured at the beginning and the end of each transect via R (version 3.5.2; R Core Team 2018). We estimated the mean and standard deviation (SD) of each pair of environmental variables (i.e., start and end AT, WT, and RH) and performed linear regression analyses to test for collinearity. N-mixture models (Royle 2004) were used to account for imperfect detection simultaneously estimating detection probability (p) and abundance (λ) across space-time (transects and months). Even though the original formulation of this model assumed closed populations, subsequent extensions relaxed this assumption by modeling the population growth rate defining the current population size as a function of the population size in a previous step (Kéry et al. 2009, Dail andMadsen 2011). In our case, we assumed that λ was constant per transect per month but allowed it to vary as a linear function of month across the year, estimating the general trend and relaxing the assumption of population closure (Mazzotti et al. 2019).
We organized and tested the field dataset in two different ways: (1) by transect by month (as collected directly in the field, N = 60) and (2) by segment by transect by month, dividing each transect into small equal segments (four segments between 0.9 and 2 km in length, N = 240) as recommended by Fujisaki et al. (2011) and Mazzotti et al. (2019). We did this to understand the best way to organize data boosting spatial replication without affecting model general assumptions (i.e., data independency). We tested the most complex model for each dataset (~p month × AT mean × WT mean × RH mean~λ month × ilength × hab for the former and~p month × AT mean × WT mean × RH mean~λ month × transect × hab for the latter) to formally compare model performing. From the covariates, ilength is the inverse of transect length used to account for sampled distance bias (Royle and Dorazio 2006) and hab is the type of habitat (river or oxbow) where observations were made. We scaled and centered all continuous variables to improve model fitting by subtracting the mean from each measured value and then dividing by the SD Royle 2016, Mazzotti et al. 2019).
We tested whether the summation limit in the likelihood evaluation (K) was high enough at its default when modeling by running both complex models under the Poisson (P) distribution and then refitting them with K = 500 (Couturier et al. 2013, Dennis et al. 2015. Once knowing that this was true (Akaike index criterion [AIC] K default = K 500), we evaluated models under three distribution structures (P, zero-inflated Poisson [ZIP], and negative binomial [NB]) using the function pcount() from the package unmarked (Fiske and Chandler 2011). We defined the distribution model that best fit both of our datasets estimating the AIC and removing variables with large P values (>0.05; Kéry and Royle 2016) as well as running 10,000 replicates of the parametric bootstrap Goodness-of-Fit (GoF) test on a chisquare (χ 2 ) discrepancy (Kéry and Royle 2016), via the Nmix.gof.test() function from the package AICcmodavg (Mazerolle 2020). We selected the best dataset and distribution structure to be used for estimating p and λ based on the AIC and GoF χ 2 (the lowest) as well as testing the hypothesis of no spatial autocorrelation using the Moran.I() function from the package ape (Paradis and Schliep 2019). To run this latter analysis, we generated a matrix of inverse weights based on the UTM projected centroids of each transect/segment estimated via QGIS (QGIS.org 2020). Once the best dataset and distribution model was selected (see Results), we fit a total of 54 models under the NB distribution structure based on all possible combinations of site and observation covariates starting from no covariance (null model~1~1) to the best performing full model (month + WT mean + month:RH mean + WT mean:RH mean~month + ilength + month: ilength + month:hab + month:ilength:hab).
v www.esajournals.org We estimated p and λ averaging the top models with a ΔAIC ≤ 2 via the model.avg() function from the package MuMIn (Barton 2009) with the main idea of including the large majority of the variation capture by models. We used this cutoff limit based on Burnham and Anderson (2002) recommendations, considering that models with ΔAIC < 2 were well supported by the data, whereas those models with ΔAIC greater than 10 were not supported. We predicted p and λ keeping covariates constant except the one of interest using the modavgPred() function from the package AICcmodavg. Finally, we estimated the population size predicted by the model averaging for the area assessed based on the sum of the values predicted.
We estimated the relative abundance (encounter rate) per transect and month based on the common information used (individuals/km) running also generalized linear models (GLM) using the gml() and the gml.nb() functions from the package MASS (Venables and Ripley 2002). To choose the best model, we used the automated model selection dredge() function to generate all possible models based on the covariates used. We then averaged models via model.avg() function, sub setting it on the basis of a ΔAIC < 2. We did this to formally compare the results of the model that includes detection probability (N-mixture) with a model that does not account for it. We reported all values derived from models followed by one standard error (SE).

RESULTS
We recorded a total of 1156 observations between December 2018 and November 2019, of which Arriba Lagoon had the largest number of observations (418), followed by Churuco Lagoon (210), Inana Lagoon (201), the river segment in between La Victoria and Inana Lagoon (186) and the river segment in between Inana and Churuco Lagoons (141; Table 1, Fig. 1). March was the month with the largest number of observations (242), followed by February (224) and January (199) whereas July (0) and June (15) had the lowest number of observations. Demographically, 90.2% of the caiman observations were from classes I (291 yearlings or below), II (560 subadults), and III (192 mixed group of subadults and adults), with only 0.4% represented by class IV sightings (five adults) and 9.3% which were unassigned observations (eyes only, 108; Fig. 2). Adults were seen in December, January, February, and September whereas the majority of class III individuals (~61%) were observed in January, February, and March. We saw caimans belonging to these two groups across all transects except one (in between La Victoria and Inana Lagoon, no class IV and only~4% of all class III). Arriba Lagoon had the largest number of caimans seen in all classes except one (class IV).
The continuous model covariates were in average AT = 25.5°AE 0.8°C, WT = 28.6°AE 1.2°C, and RH = 84.5% AE 2.7% at the start point and AT = 25.7°AE 1.4°C, WT = 28.5°AE 1.2°C, and RH = 87.2% AE 2.0% at the end point, varying all significantly between them across the whole study (AT [t test = −3.54, df = 1716, P < 0.001], WT [t test = 2.02, df = 1716, P = 0.04], and RH [t test = −24.099, df = 1716, P < 0.001]). After estimating the mean, we found very low (almost no) correlation between the predictors (AT mean-WT mean r 2 = 0.06, AT mean-RH mean r 2 = 0.06, WT mean-RH mean r 2 = 0.00), allowing us using all variables to model p and λ without overfitting them. Moran's I analysis rejected the null hypothesis of no spatial autocorrelation in the dataset organized by segment/transect/month (Moran's I observed = 0.02, expected = −0.004, P < 0.001) but cannot reject this hypothesis when data were organized by transect/month (Moran's I observed = −0.03, expected = −0.01, P = 0.136). Concomitantly, model-fitting analyses showed that the full NB model from the transect/month dataset outperformed both the ZIP and the P distribution structure models from this dataset as well as all full models from the segment/transect/month dataset (Table 2). GoF χ 2 test also corroborated the NB distribution from the transect/month dataset as the best method to model our data due to no evidence of lack of fit (P = 0.62) or overdispersion (ĉ = 0.89) compared with ZIP and P distributions (P < 0.05,ĉ = 2.52 and 9.22, respectively).
The top nine performing models (out of the 54 ran total) varied <2 AIC units and included all the covariates used in the present study, receiving a total of 56% of the model weight (Table 3). Once all models with an AIC < 4 were included (22 total) the model weight increases up to 89%. Model averaging analysis showed four component models with a ΔAIC < 2 and 70% of model weight (Table 4). Model-averaged coefficients  Relative abundance (encounter rate) index (individuals/km) values were around five times lower than those predicted by models per transect (Table 5). The GLM showed the best performance under the NB variation structure (AIC = 450.6,ĉ = 1.34) compared with the P distribution (AIC = 775,ĉ = 9.94). Model averaging of the NB variation structure models with an AIC < 2 showed a significant positive association with all response variables (conditional average, P < 0.05) except with month in which it has a negative relationship. The model that included month, transect, AT, WT, and RH as response variables received the highest model weight (54%), followed by two other models with all these variables but alternating the inclusion of RH and transect (25% when RH absent and transect present and 21% when RH present and transect absent; Table 6).

DISCUSSION
As expected, abundance values predicted by models that account for imperfect detection (estimating detection probability) were much higher than those estimated assuming 100% detection (relative abundance index), reinforcing the idea of a generalized underestimation of this population parameter in spectacled caiman studies (and likely in many other crocodylians; Fujisaki et al. 2011) that employed the latter method. A recent study also showed how indices such as King's and Messel's sighting fraction (commonly used in crocodylian population studies) underestimated population size estimates due mainly to Notes: Notice the differences between the top-performing model (NB from Transect/Month) and the other five (ΔAIC) as well as the number of parameters (nPars) used. Notice also the critical values obtained after performed 10,000 parametric bootstrapping goodness-of-fit χ 2 test (GoF) and theĉ per model. Model selection criteria used in the present study showed how subdividing transects in small segments to increase spatial and temporal replication not only increase the AIC of models (Table 2) but also lead to a spatial autocorrelation of the data (Moran's I P < 0.05), violating one of the main assumptions of N-mixture models (data independency). We recommend always testing for spatial data dependency as a criterion to define whether dataset accomplish this requirement for estimating detection probability and abundances under N-mixture modeling. The selection of a NB distribution structure over the ZIP and P distribution was highly influenced by the high discrepancy of mean (19.26) and variance (442.97) in data collected from the field. However, we found that under the P and ZIP models mean abundance was estimated twice as high (45.07 and 41.21 individuals, respectively) than under the NB distribution (21.13 individuals) being the former the most conservative method to be used with our data.
Our top-performing component models obtained after model averaging (ΔAIC < 2; month, WT, and RH as observation covariates and month and habitat as site covariates; Table 4) explained 70% of the variation in spectacled caiman's detection probability and abundance in the study area, meaning that factors other than the ones used also affect p and λ. In a similar way, Cartagena-Otálvaro et al. (2020) found that vegetation cover and water depth across channels explained 44% of the variation in p and λ of spectacled caimans in the Magdalena River middle basin. In contrast, Mazzotti et al. (2019) found that air temperature, water temperature, moon phase, habitat type, and salinity per routes and years explained 30% of the variation of the same parameters for American crocodiles (Crocodylus acutus) in Florida. Other authors such as Evert (1999) and Rosenblatt and Heithaus (2011) studying alligators (Alligator mississippiensis) in Florida, reported that the abundance of this species correlates to the levels of nutrients and prey availability. Thus, relating the occurrence of caimans to the distribution and abundance of prey items as well as additional habitat attributes ( Notes: We used negative binomial variance structure and the dataset set up by Transect/Month to run all models as it was the top-performing model that for these data. The significant covariates (P < 0.05) selected and used for modeling were as follows: month (mt), water temperature (wt), relative humidity (rh), the combination of mt:rh and wt:rh, the inverse of transect length (ilength), habitat (hab), and the combination of mt:ilength, mt:hab, and mt:ilength:hab. Notice that the top-performing model received 12% of the model weight (AICwt) and the first nine models received 56% of the cumulative weight (CWt). Table 4. N-mixture model analysis using a negative binomial variance structure showing the component models derived from the model averaging analysis based on the top nine models from the complete model set (see Table 3). (e.g., vegetation and water cover, water depth, current speed, wind speed), should increase the total variation explained by models describing p and λ. Basic variables such as precipitation and water cover are highly relevant when counting crocodiles due to the increase of habitat heterogeneity, which directly reduces encounter rates (Hutton and Woolhouse 1989, Thorbjarnarson 1989, Kofron 1993. However, the lack of meteorological stations and reliable imagery for the study area limited their inclusion and analysis in the present study. The population size predicted by our averaged model (1763 AE 786 caimans) is lower than the one estimated by Cartagena-Otálvaro et al. (2020) in the San Juan channel (2493 caimans) but higher to the one estimated in Caño Negro channel (500 caimans) across the Magdalena River middle basin. However, the lack of information regarding the area sampled cover, the fact that they only sampled eight days in a row rather than across a period of time that would then allow them to accurately estimate the actual variation of the population (substantial sample size), the general lack of fit of the different models examined, and the apparent overdispersion of the dataset, make these data incomparable. To date, no other studies of spectacled caiman populations have used N-mixture models to estimate either abundances or populations sizes, which limits broader analyses and comparisons.
Relative abundances/densities reported across C. crocodilus distribution have shown high variability ranging from 0 to 58 caimans/km in lakes and from 0 to 7 caimans/km in canals in Amazon rainforest habitats (Da Silveira et al. 1997), and from 4.3 caimans/km to 11.6 caimans/ha in lotic and lentic environments in Magdalena River middle basin habitats (Moreno-Arias et al. 2013). This variation reflects a broad spectrum of relative abundances and population sizes influenced by anthropogenic pressures (i.e., hunting, habitat loss) but also likely reflects different carrying capacities of the various ecosystems where the species inhabits. Unfortunately, the lack of robust and comparable estimates of abundance and population size across the range Fig. 3. Relationship between the spectacled caiman predicted abundance (λ) and detection probability (p) via model averaging and four of the covariates assessed. Notice that the predicted abundance slowly increases as transect length increases and abundance and detection probability decreases as months pass through. Notice also how these models found an apparent no effect of the water temperature and relative humidity in the predicted detection probability. Gray lines represent 1 standard error and black lines the mean. The other covariates in each figure were held at a constant value for model predictions.
of C. crocodilus limits our understanding of this matter.
We observed in the field a prevalence of spectacled caimans in lagoons (oxbows) compared with transects on the river, which relates to what has been reported for the species in areas such as Cispata Bay (Ulloa and Cavanzo 2003) and the Salamanca Island Road Park (Balaguera-Reina and González-Maya 2009) in the Colombia Caribbean region or Amazon rainforest habitats (Da Silveira et al. 1997). However, this prevalence was not supported by the expected abundance after model averaging, predicting the same values in both types of habitat across the analysis. This contradicting outcome could derive from the effect caused by differences in transect length, that although normalized by its inclusion as a covariate, could increase the expected number of caimans to be observed in transect 3 (largest transect), even though observations in the field were lower. Then, we recommend trying to keep as low as possible the variation among transects to avoid undesirable effects on modeling. This prevalence observed and reported in literature in other regions could be related to habitat conditions generally present in rivers (such as differences in current and wind speeds, water depth, and availability of floating vegetation) that do not necessarily facilitate high concentrations of individuals (Medem 1981, Ayarzagüena 1983, Cerrato 1991, Da Silveira et al. 1997, Balaguera-Reina and González-Maya 2009). Superimposed on this habitat variability is a larger anthropogenic impact due to human encroachment and movements (Ron et al. 1998).
The demographic structure found in the area with a large number of juveniles (class I) and subadults (class II) and a lower number of adults (especially class IV) might indicate that this population still is under some level of hunting pressure, resembling the demographic structure described by Ayarzagüena (1983) for such populations. Ethno-zoological analyses done reported the use of all three species of crocodylians present in the area (i.e., C. crocodilus, P. palpebrosus, P. trigonatus), with some indigenous communities having a higher preference for those crocodylians of the genus Paleosuchus (Balaguera-Reina 2019). Future studies focusing on the use of crocodylians in the area should be done to help to understand and quantify the impact(s) caused by local, unquantified use of caimans.
We found that the average detection probability decreased, and the SE increased as months passed through, which for us, represents an indirect reflection of the effect of a set of variables (precipitation, water level, water cover) not included in the study. The Amazon region's climate can be classified as equatorial superhumid  Note: Notice that the first model, which included all covariates (at, air temperature; mt, month; rh, relative humidity; tr, transect; wt, water temperature), received 54% of model weight (AICwt). area with high levels of precipitation throughout the year (Guzmán et al. 2014), especially flooding large areas of land around rivers and lagoons when rainfall is at its peak (from April to August), and dramatically influencing spectacled caiman habitat availability. This same pattern has been reported in areas such as the Brazilian Amazon and Venezuelan savannas to heavily influence in the number of crocodylians observed (Velasco et al. 2003, Da Silveira et al. 2008. Future studies using hierarchical models that include these variables could help to understand their overall effect on the detection probability and the relative density parameters in caimans, putatively improving model fit. It is important to highlight that inference derived from N-mixture models are very sensitive to model assumptions (closure assumption, false-positive errors, independence of detection, homogeneity of detection among individuals, and parametric modeling assumptions; Kéry and Royle 2016) so researchers must be careful about its use when data are not adequate (i.e., overdispersed, many extreme values, high probability of double counting; Link et al. 2018). In our case, we relaxed the assumption of population closure and reduced the likelihood of falsepositive errors assuming that λ was constant per transect per month but allowed it to vary as a linear function of month across the year, leaving us with monthly abundance estimates within which these two assumptions should be accomplished. This is because monthly spotlight surveys have a low probability of counting twice the same individual in a transect due mainly to (1) transects are unidirectional (normally upriver), (2) the time collecting information once observed is short (within minutes), and (3) the cryptic-awaiting hunting behavior of crocodiles. However, the impossibility to know for real whether we double counted or not any animal brings some unmeasured uncertainty to the model (Link et al. 2018).
Overall, the spectacled caiman population ecology parameters assessed in the present study in the Apaporis River middle basin showed abundances and demographic structures similar to populations under some level of hunting pressure (Ayarzagüena 1983, Da Silveira et al. 1997, Velasco et al. 2003. This finding will require indepth studies to define management measurements that guarantee a sustainable use of the species. However, at this time we can say that the species is not under threat and only requires regulation that helps to use this biological resource in a sustainable manner.

ACKNOWLEDGMENTS
This project was sponsored and approved by the University of Ibague (18-489-ESP), Texas Tech University, IUCN/SSC/Crocodile Specialist Group, CrocFest, and CrocDocs. We thank Lucrecia Castillo, Jarvi Vargas, Johnatan Vargas, Frank J. Mazzotti, Brian J. Smith, Venetia Briggs-Gonzalez, and the CrocDocs team for their help in the field and in the data analyses. We finally want to thanks to the editor and reviewers for their valuable input to improve the manuscript.