Statistical methodology for the evaluation of leukocyte data in wild reptile populations: A case study with the common wall lizard (Podarcis muralis)

The leukocyte profile has the potential to be a reliable method to measure health conditions and stress in wild animals, but limitations occur because current knowledge on reference intervals is largely incomplete, especially because data come from studies on captive animals involving few individuals from single populations. Here we propose a general framework for achieving reliable leukocyte reference intervals, encompassing a set of internal and external factors, potentially affecting the leukogram. To do so, we present a systematic survey of the hematology of the common wall lizard, Podarcis muralis, involving 794 lizards from 54 populations over the whole geographic range of the species in Italy. Reference intervals for white blood cell (WBC) and leukocyte differential count were obtained by using linear mixed models in a Bayesian framework. The application of the procedure clearly showed that both internal (sex and size) and external (latitude and season) factors are a source of variation of leukocyte profile. Furthermore, the leukogram of common wall lizard has a strong variability among populations, which accounts for more than 50% of the whole variation. Consequently, some common assumptions used in studies on captive individuals are no longer supported in wild populations, namely, i) any group of individuals is a representative sample, ii) any population is representative of all others, iii) geographic clines do not occur over the species range, and iv) seasonal variation has limited effects. We encourage researchers aimed at the definition of leukocyte reference intervals for wild populations of reptiles to involve a large number of populations over a wide geographic range in ad hoc statistical models to disentangle local and geographic effects on leukocyte profile variation.


Introduction
Characterizing the physiological responses of wild animals to stressors, including humaninduced landscape changes, is an important question with deep implications for both animal conservation and ecology [1]. Indeed, stress is a key factor when dealing with animal welfare in both wild and captive contexts [1,2], as well in determining how resource investment in immune-function shapes species life-history traits [3]. Leukocyte profile has the potential to be a reliable method to measure stress in vertebrates, as an alternative to hormone assay [1]. For example, neutrophils (heterophils in birds and reptiles) and lymphocytes respond in the opposite way to stress, so that their ratio is positively related to the circulating glucocorticoids and the magnitude of the stressor [4]. The increase of the total number of leukocytes (leukocytosis) has been occasionally used as a measure of stress [reviewed in 1] as well as glucocorticoid-induced stress has been shown to reduce the eosinophil number, even if only in humans and other mammals [reviewed in 1]. Moreover, the increase in the number of eosinophils (eosinophilia) has also been used as a measure of stress, although with contrasting results [1]. Since the lower end of eosinophil reference intervals frequently include zero, and eosinophils typically represent a minor part of the leukogram in many species, changes to eosinophils are typically not used to assess for stress in human or veterinary medicine. Compared to hormone assay, the enumeration of white blood cells from blood smears is cheaper and makes the assessment of the baseline much easier to do since the initial response of leukocytes begins within hours (or in a few days) after the stress. The downsides are essentially two: the interpretation of leukocyte counts is difficult, and the leukocyte profile does not necessarily supply information about the ability of individuals to mount an immune response [1].
The potential of the leukocyte profile as a tool in ecological and conservational research can be fully realized only after robust information on reference values for "normal" individuals is achieved. Reference values (RV) describe the variation in blood analyte concentrations in well characterized groups of healthy individuals, while Reference intervals (RI) represent estimated distributions of RV from healthy populations of comparable individuals [5]. Population-based RI have become one of the most common tools used in clinical decision-making processes [5]. To be reliable for wild individuals, RI should be assessed by adopting general guidelines for the selection of reference individuals, which must embrace the variability among populations of the focal species [5]. Without this information, disclosing if the cell count of a given individual is high or low compared to the leukocyte profile of healthy individuals is, in fact, impossible. Unfortunately, many of the RI proposed so far have been obtained with captive animals or, when done in the wild, by a survey of only a few (even only one) populations (Table 1). Indeed, the approach used to date in ecological and conservational research has been borrowed from veterinary research and has involved samples generally around 50-70 units, from a single locality, without replication (Table 1). Such a kind of approach can provide reliable data when dealing with pets or farm animals, which live in strictly controlled conditions that sensibly reduce all those sources of variation in leukocyte profile proper to natural conditions, such as habitat heterogeneity, climatic variability, food distribution, and competition with other species, just to name a few. The effect of these external factors on leukocyte profile in wild populations, however, cannot be disregarded, as they can be one of the main causes of variation among populations, all over the range of the species [6,7].
Although differences between captive and free-ranging individuals have been documented in reptiles [14,16,24,33], the lack of RI for wild populations is particularly severe for this clade [34]. Many authors have described the morphology of reptiles' blood cells, particularly for species of commercial interest, such tortoises, crocodiles, and iguanas [10,33,35]. Several attempts to determine RI for leukocyte differential counts have also been performed, but mostly involved captive individuals or one or a few populations ( Table 1). The first attempt dates back to the work by Duguy [8,9], whereas the first review on RI for leukocytes in reptiles is by Frye [10], who reported the leukocyte differential count for 24 species of reptiles (including four tortoises, fourteen lizards, five snakes, and one crocodile). However, in all but one case, no more than three individuals for species were considered. Since then, several studies have investigated reptilian hematology in wild individuals, but always with the limitation of using a few individuals collected from one/few populations (Table 1). During the last decades a large amount of data has also been collected in parallel for veterinary purposes, but even in this case, the large majority of studies involved few individuals and did not account for possible seasonal variation neither did cover a significant portion of the geographic range of the species (Table 1). To date, only a very restricted number of studies have focused on reptiles hematology in wild populations over large geographic ranges [6,26,30]. Dickinson et al. [6] supplied RI (with seasonal and annual variation) for 36 radio-tagged Gopherus agassizii sampled in three localities of the Sonoran desert, Ozzetti, et al. [26] supplied leukocyte baseline for Oxyrhopus guibei and Xenodon neuwiedii in Brazil collecting individuals opportunistically, leading to sample de facto only one individual for each locality, whereas Gangloff et al. [30] supplied a baseline for the heterophils to lymphocytes ratio on a sample of 86 Thamnophis sirtalis collected in 12 populations widespread over the distribution range in North America. Given the lack of reliable RI, there is the need to develop standardized techniques [5] to collect data from wild populations, to be able to properly use leukocyte profiles in ecological and conservation studies of reptiles. A further issue regarding the reliability of the RI concerns the statistical treatment of data used to define them [5]. Most studies compare groups of individuals (mostly males and females) by simple one-way analyses (i.e., t-test or one-way ANOVA), and established intervals reporting ranges as well as standard deviations/errors of the observed data. However, such a kind of intervals only describes the probability of obtaining the same result if the measure is replicated but does not allow any inference on the possibility of classifying new observations. No previous studies have explicitly implemented the effect of the breeding season in the analyses, nor the effect of the geographic variability.
We believe that reliable RI for ecological and conservational purposes can only be achieved through the application of a general, well-defined protocol encompassing the patterns of variation of leukocytes, due to both internal (age and sex) and external factors (season and latitude) potentially affecting leukocyte parameters. Furthermore, we also think that RI for leukocytes should be estimated through an appropriate statistical model and not directly from the raw data. Only in this way it is possible to rule out the effects of the external and internal perturbations to leukocyte parameters and estimate the baseline of RV. In summary, to overcome the above issues, the best solution would be to use: 1) large sample within the population including individual from both sexes and of different age classes; 2) repeated sampling within the population to account for seasonal variation of leukocytes' relative abundances; 3) sampling data from different populations settled all over a large portion (possibly the whole range) of the distribution of the species; 4) appropriate statistical models that allow to decompose the variance among individuals according to the internal and external factors involved, and to assess the among-populations variation.
The common wall lizard has been extensively studied in light of hematology and immune function [e.g., 20, [36][37][38][39][40]. It is a small lacertid lizard (snout-vent length, SVL, 45-75 mm) widespread in Southern and Central Europe, which produces on average two clutches per year [41]. The species is sexually dimorphic [42] and the breeding season starts in late February and ends in July [41]. Body temperature during activity is near 33˚C, being slightly higher (33-36 C) in warmer regions (e.g., Central Italy) and lower (32˚C) in mountain areas [43][44][45]. Morphology of blood cells has been repeatedly studied [20,36,37,39] as well the immune response has been previously investigated concerning polymorphic ventral coloration [40,46], immuno-competence handicap hypothesis [38], or temperature [40,47]. Despite this, RI for the leukocyte profile in wild populations still lack, though some pivotal studies have been performed in single, isolated populations [39].
In this paper, we applied the suggested framework to the analysis of the common wall lizard (Podarcis muralis) leukocyte profile, with a double aim in mind: i) quantifying the bias we could introduce in the RI estimates by ignoring (i.e., not measuring) the among-population, seasonal, and geographical variability; ii) measure the performance of our procedure when dealing with those sources of variation, and its ability to eliminate or control for the effects of outliers.

Materials and methods
Adult common wall lizards were collected by noosing during March-September between 2008 and 2017 in 54 sites within the whole species distribution in Italy (Fig 1; Table 2). Only lizards in apparently healthy and in good condition (external examination) were selected for this study and were handheld only for the time needed for measuring and blood collecting.  [41], and 59 mm if males [48].
Blood samples (15-20 μl) were collected in the field, from the postorbital plexus, using heparinized capillary tubes, after holding the head of the lizard firmly in one hand [49]. Soon after blood sampling individuals were released healthy in the capture site. Air-dried smears were stained with May-Grünwald/Giemsa stain and scanned using a light microscope (OPTIKA) at 60× following standard routines. In each microscope field, leukocytes were counted and classified as lymphocytes, monocytes, eosinophils, heterophils, or basophils. In each smear, we explored no less than 50 fields for no less than 150 leukocytes (i.e., if after 50 fields the count of leukocytes did not exceed 150, new fields were added until this value was reached). This ensured us to be able to detect leukocytes until a frequency of 0.7%. Total white blood cell counts (WBC) were estimated from May-Grünwald/Giemsa stained blood films by multiplying the mean number of leukocytes per microscopic field by the objective power squared [50].

Statistical analyses
The WBC were log-transformed to achieve normality. Leukocyte differential count is inherently a composition (with dimension D = 5, [51]), and cannot be used in linear model as it stands. Therefore, we used the isometric log-ratio transformation (ilr) to map the leukocyte differential count of each individual in a D-1 Euclidean space, thus generating a set of four independent variables [51]. To account for seasonal variation of leukocyte variables, we used multiple-component cosinor models, which were first developed by Halberg et al. [52] to model physiological processes that have a circadian rhythm. Those models use linear combinations of cosine functions to assess the amplitude (A) and phase (φ) of physiological rhythms around the average value (i.e. the midline estimating statistic of rhythm, MESOR) over the period (τ). We used a single-component cosinor model, i.e., Y(time) = MESOR + Acos(2π(time)/τ + ϕ) + e(time) including time expressed as Julian date (1 = January 1 st ) and τ = 365 to account for effect of circannual rhythms around the baseline of each leukocyte variable. This model was incorporated into a random intercept Linear Mixed-Model (LMM) to assess the combined effects of circannual rhythm, biotic, and abiotic factors on lizards' leukocyte profile. Cosinor, sex, body size (standardized SVL), and latitude (standardized UTM coordinate) entered the model as fixed effects. We also added the cosinor × sex, cosinor × body size, and cosinor × latitude interactions to account for possible sex, size, and geographical specific patterns in circannual rhythms on leukocyte variables. The population entered the model as a random intercept to account for unexplained variation at the population level (σ 2 pop ) when controlling for the explanatory variables. An independent model was carried out for each hematological variable. Models were fit in a Bayesian analytical framework available in the package R2jags [53] in R 4.3.3 [54], which uses the samplers implemented in JAGS 4.3.0. Uninformative normal priors (μ = 0 and σ = 0.001) were used for model's coefficients, and gamma priors (a = 0.001 and b = 0.001 corresponding to μ = 1 and σ = 1,000) were used for both the error (σ) and the random intercept (σ pop ). Three chains were run using randomly selected initial values for each parameter within a reasonable interval, and conventional convergence criteria were checked. The number of iterations was selected for each run to obtain at least 10,000 valid values for each chain after convergence and thinning. The posterior distributions were back-transformed to obtain the posterior distributions of the variables in their original scales. According to ASVCP guidelines [5], baselines and RI were extracted from the posterior distributions of the values predicted by the models, as the mode and 95% prediction intervals.

WBC
The posterior distribution for WBC was skewed, with a longer upper tail (Fig 2a), suggesting that high values, even three-time greater than the mode of the distribution, of these two parameters could occur in wild individuals, albeit with low probability. Consequently, the RI for WBC as derived by the mode and the 95% prediction intervals of posterior distributions are sensibly asymmetric compared to the mean (Table 3). Accordingly, the probability of finding individuals with WBC exceeding the mode was 0.70, but it falls to 0.21 for individuals with WBC twice as high as the mode, until 0.05 for individuals with WBC three times the mode (despite two and three times the mode were both within the 95% prediction intervals).

Leukocyte differential count
As for WBC, the posterior distributions of leukocytes' relative abundances were sensibly skewed too, with a longer and shorter tail (Fig 3). Consequently, the 95% prediction intervals (i.e., the RI) were highly asymmetric (Table 3).

Among populations variation
We found a relevant among-populations variability (σ 2 pop ) on all hematological variables, which accounted for a consistent proportion of the unexplained variation after we controlled for the explanatory variables. For WBC, it accounted for 57% of the total variation, whereas σ 2 pop accounted for between 39% (lymphocytes) and 98% (eosinophils) of the total unexplained variation. The effect of population was negligible for monocytes (σ 2 pop = 3%) and basophils (σ 2 pop = 7%).

Discussion
In ecological and conservation research, leukocyte profile has the potential to be a reliable method to measure stress and health condition in animals [1], as long as the reference values for wild and healthy individuals are fully available. The methodological approach used to assess. RI is of primary importance, since the reliability of these values depends precisely on how they were obtained [5], and, as a cascade effect, also the validity of all decisions that are taken from these intervals. To date, there is no single and shared protocol to obtain reliable RI for wild animals (but see [5] for veterinary species), as evidenced by the outcome of our survey on previous studies on reptiles (Table 1). What is crucial, is that no one has yet made any attempt to test the reliability of those procedures, or to assess how important the effects of the variability among individuals and populations of the same species can be. In this study, we proposed a general method for defining RI for leukograms having a dual advantage: i) encompassing a set of internal (linked to individuals) and external (linked to the populations) factors, potentially affecting leukogram's RI, and ii) assessing how much noise each of these factors causes to the RIs. These two goals can be achieved thanks to the use of linear mixed models for the definition of RIs instead of simple means and standard deviations on raw data as in previous studies. By using linear mixed models, baselines for RVs are obtained through the model intercept, while the effects, if any, of the external factors can be estimated through the model coefficients. In addition, the effect of the variability among populations can be assessed by inserting the population as a random effect into the model. Concerning previous approaches, the novelty of computing baselines and RIs through linear modeling is that upper and lower limits of the RI are estimated after removing the effects of both internal/external factors affecting leukocyte concentration and among-population variation over the distribution of the species. Furthermore, if models are computed under a Bayesian analytical framework, the RIs (i.e., 95% prediction intervals) represent the actual probability distribution of leukocytes of any individual in any population within the range of the species. In other words, intervals derived from the 95% prediction intervals gave us the range within which a new individual was predicted to fall with an accuracy of 95%. These ranges can be regarded as the natural variability of leukocyte parameters of the species and any deviation from it can be reasonably regarded as an indicator of possible diseases including inflammatory responses, stress, toxins, and leukemia.
The application of the procedure to wild populations of the common wall lizard data confirmed our hypothesis that both internal and external factors are a source of variation of the leukocyte profile. The direct implication is that some of the assumptions implicit in the studies in captive animals are no longer supported in studies involving wild populations, namely, i) any group of individuals is a representative sample, ii) any population is representative of all others, iii) lack of geographic clines over the species range, and that iv) seasonal variation has limited effects.
The assumption common to all previous studies that "all individuals involved were healthy" (i.e., the values obtained could be regarded as baselines) and, therefore, representative of the species baseline, was violated as demonstrated by our analyses. Indeed, the probability of including non-healthy lizards into the sample was far from being negligible, as, for example, 3.4% of males and 1.6% of females included in the sample exceeded the upper limit of the RI (leukocytosis). These are not trivial deviations, but rather they can be a serious problem when RIs are estimated as means of raw data and intervals are derived using the minimum and maximum of observed values [6,12]. Indeed, by using minimum and maximum of observed values, all sampled individuals showing extreme values are regarded as healthy, when they might actually suffering from some disease (false negative). Our results show that such an approach may supply not accurate reference and overestimate the range of variability of leukocytes in healthy individuals, leading to type II errors (i.e., false negative) when screening individuals for disease occurrence. We further point out that also the statement that "all individuals were in healthy condition given that none was showing any pathological condition", which is often used to improve data reliability, is not fully supported, given that all lizards exceeding RI did not show any morphological or behavioral anomaly. Accordingly, the ASVCP guidelines for the determination of RI point out the need to include procedures for outlier detection [5], and our statistical methodology represents a reliable approach to deal with outliers and their effect on RI estimations.
Our analyses show that also the assumption that "any population can be regarded as representative of the whole species" is not reasonable for the common wall lizard. Indeed, we found a huge variability among populations in most RVs, and variation among populations explained much more than 50% of the whole variation in the leukocyte profile. Since the population was modeled as a random effect, variability at the population level was not related to sex, age, latitude, and season, and should, therefore, be regarded as a property intrinsic to each population. We should not be surprised by that. Rather, it is exactly what we should observe if the leukocyte profile records the stress and health condition of individuals. Indeed, much of the variation in leukocytes is expected to depend on the site-specific conditions that individuals are experiencing, such as parasites, diseases, food availability and quality, predators, and competitions, just to name a few. We also believe that small samples with a lack of replication in more than one population are likely to be the major cause for the large variability in RIs for the same species among different research groups [1]. From a methodological point of view, our results indicate that RIs obtained from a single population should be regarded with caution, even if computed over a large sample size, and reliable leukogram baselines should only be achieved by sampling a large number of populations from the whole species' range.
The assumption that RIs are constant all over the range of the species did not find support in our analyses, as we detected a latitudinal cline of variation in leukocyte profile. From Southern to Northern populations, the WBC increases, the proportion of heterophils decreases while the opposite does occur in lymphocytes, leading the H:L ratio to decrease northwards. This pattern of variation is consistent with the ecology of the common wall lizard, which is a continental species of mainland Europe [55], and in Italy, it becomes increasingly rare from North to South, where it occurs only in mountains [56]. More in general, species that are spread over a large geographical area may experience different climatic conditions and may have consequently evolved local adaptations, which are reflected in the geographic pattern of variation of RIs in leukocyte differential count. Again, from a methodological point of view, the only way to account for such a pattern is to collect data from several populations, spread all over the species' range.
Finally, we clearly showed that several components of the leukogram are characterized by a regular pattern of variation above and below the reference values during spring-summer, possibly reflecting seasonal variation in temperature or for other factors, namely hormone plasma levels [57,58]. Indeed, reference values of WBC and lymphocytes sensibly decreased during the peak of the breeding season (May), whereas heterophils and eosinophils increased. Patterns of variation in leukogram over the breeding season appear to be a component common to several species of reptiles [33,35], and our data indicate that seasonal oscillations around the baselines may exceed RIs. Therefore, if the seasonal pattern is ignored, some healthy individuals might be regarded as sick, leading to type I error (false positive). Baseline leukocyte profile also varies according to body size, which, in lizards, is a proxy for the age. In general, adults had more WBC, more heterophils, eosinophils, and monocytes, and fewer lymphocytes compared with juveniles and subadults. Thus, the leukocyte profile in this species seems to show senescence, as previously reported for other species [35]. Seasonal and age patterns we detected in common wall lizards are only two examples of how internal/external factors may affect leukocyte profile, but point out the need to account for such effects when trying to define leukogram's RIs.
In conclusion, the application of WBC counts and leukocyte differential counts in conservation, as well as ecological research on reptiles, is still limited mainly because of the lack of basic information, but also by the inconsistencies among studies performed until now. These troubles mainly arise from the heterogeneity among methodological approaches and low standardization in experimental design. We argue that it is possible to disentangle local and geographic effects on leukocyte profile only by involving a large number of populations over a wide geographic range in ad hoc statistical models. Otherwise, the inconsistency among studies will carry on without a solution, and we may never get to leukocyte profiles reliable enough to be used to assess the health conditions of wild animals.