Macroevolutionary diversification of glands for chemical communication in squamate reptiles

Chemical communication plays a central role in social, sexual and ecological interactions among animals. However, the macroevolutionary diversification of traits responsible for chemical signaling remains fundamentally unknown. Most research investigating evolutionary diversification of glands responsible for the production of chemical signals has focused on arthropods, while its study among vertebrates remains neglected. Using a global-scale dataset covering > 80% (7,904 species) of the living diversity of lizards and snakes (squamates), we investigate rates, trajectories and phylogenetic patterns of diversification of their follicular glands for chemical communication. We observed these glands in 13.66% of species, that their expression has varying phylogenetic signal among lineages, and that the crown squamate ancestor lacked follicular glands, which therefore originated and diversified subsequently during their evolutionary history. Additionally, our findings challenge the longstanding view that within squamates the Iguania are visually oriented while Scleroglossa are chemically-oriented, given that Iguania doubles Scleroglossa in the frequency of glands. Our phylogenetic analyses identified stabilizing selection as the best model describing follicular gland diversification, and revealed high rates of disparity. We provide the first global-scale analysis investigating the diversification of one of the main forms of communication among reptiles, presenting a macroevolutionary angle to questions traditionally explored at microevolutionary scale.

has influenced our understanding of animal interactions. For example, although the role of female mate choice in sexual competition has been shown to be widespread among animals 15,16 , this form of sexual selection has been difficult to detect in reptiles, leading to the idea that competition over mates if fundamentally mediated by male-male contests 17 . However, emerging evidence suggest that female lizards can in fact choose males based on the chemical signals they produce [18][19][20] . Similar studies have shown the role for chemical signals in other forms of social interactions in these vertebrates, including territoriality 21,22 , social recognition 23 and conspecific assessment 24,25 .
Even though reptiles deliver scents in multiple ways (e.g., via the skin, or feces) 13 , chemical communication among these organisms seems to be predominantly based on signals produced by follicular epidermal glands (FG) [26][27][28] . These glands are specialized tubular structures embedded in the dermis that discharge waxy secretions that are delivered into the external environment through epidermal pores 26,29 . The numbers and body locations of FG differ widely across species and lineages 26 , varying from zero to ~130, while they can be situated around the cloaca ('precloacal'), on the ventral surface of the thighs ('femoral'), or in both regions 26 . In the majority of cases, these FG are restricted to, or show more complexity in males. Surprisingly, however, only a few studies restricted to specific lineages have investigated the evolution of these glands. Some studies have hypothesized that the wide interspecific diversity of FG could be the result of selective mechanisms mediating the enhancement of efficacy of the chemical signals as a function of the environments that species reside in [30][31][32] . Given the energetically costly demands associated with the development and functioning of these structures, selective pressures (e.g., allocation of energy, environmental factors) could shape adaptive variation of FG. On the other hand, a comparative study conducted in the hyper-diverse lizard genus Liolaemus revealed that, despite a previous conclusion of adaptive evolution 32 , the primary factor explaining patterns of FG numbers across lineages was shared ancestry 30 . Likewise, another study on lacertid lizards showed a marginal environmental effect on FG number and strong phylogenetic inertia 31 . However, and despite these phylogenetic effects within lizard families, FG have extensively evolved along the squamate phylogenetic history. Therefore, although these studies have shed some light on the factors influencing variation in FG, no research investigating the evolutionary trajectories and rates of these organs across a truly diverse range of lineages and areas of the world exists.
In this paper, we present the first global-scale study investigating the evolution of FG across the squamate tree of life. Based on an originally created dataset for 7,904 species and a molecular phylogeny containing 3,533 of the species for which FG data are available, we implement a model-selection approach aimed to quantitatively establish the evolutionary trajectories, tempo and mode of diversification of these organs.

Results
Diversity of follicular epidermal glands. We found that FG are present in 13.66% of the 7,904 species in our dataset (Supplementary Table S1). After excluding snakes (which entirely lack FG), a total of 24.8% of the squamates have FG. These proportions vary across clades, FG being present in 35.2% of the Gekkota species, 26.82% of Iguania, and in 96.8% of Lacertoidea. In contrast, lower proportions of species with FG were found in Dibamidae (14.28%) and Scincoidae (1.11%) ( Fig. 1B; Table 1). FG are entirely lacking from the clade Anguimorpha (173 species). Importantly, the proportion of species with FG is higher in Iguania (26.82%) than in Scleroglossa (which contains all other squamate lineages, with an 11.11%). The gecko Mniarogekko chaoua was identified as the species with the highest mean number of FG (95), while multiple species of different subclades presented only 1 or 2 of them. Regarding the location of FG, we found that 47.82% of species only have precloacal FG, 36.49% only femoral FG, and 15.69% have both types (Fig. 1B).
Analyses performed among species grouped according to whether they have precloacal, femoral or both forms of FG revealed significant differences in their mean number (phylANOVA; F 2,605 = 172.58, P = 0.016). Species with precloacal FG have on average a lower number of these glands (mean = 8.51) compared to species with femoral FG (mean = 29.98; P = 0.024), or with glands in both locations (mean = 28.81; P = 0.034).

Phylogenetic signal and ancestral state.
Overall, the number of FG showed a moderate phylogenetic signal in Squamata, with a high λ and intermediate K ( Table 2). A qualitatively similar degree of phylogenetic signal was observed in the subclades Gekkota and Lacertoidea. In contrast, K-values were very high in Iguania and especially in Scincoidea, meaning that species from these clades resemble each other more in their mean number of FG than expected under Brownian motion of evolution. The body location of FG (cloacal, femoral or both) was found to be highly phylogenetically conserved in squamates, and in lizards in particular ( Table 2). Macro-evolutionary patterns and models. Our model-selection analyses investigating macroevolutionary diversification dynamics of FG across squamates identified the Ornstein Uhlenbeck model (OU) as the best approximation describing the tempo and mode of evolution of these structures in squamates in general, as well as in Gekkota, Iguania and Lacertoidea when analyzed separately. The Delta, Brownian Motion (BM) and Early-Burst (EB) models were ranked in decreasing order based on their Akaike Information Criterion (AIC) values after the OU model (Table 3).
Disparity-through-time analyses (DTT) performed on the squamate tree returned a positive morphological disparity index (MDI = 0.07), indicating high levels of subclade disparity: i.e. the subclades have diversified considerably and the ranges of their FG numbers overlap extensively (Figs 2B and 4). We also obtained positive relative disparity indices for Gekkota (MDI = 0.23), Iguania (MDI = 0.16) and Lacertoidea (MDI = 0.10), when analyzed separately (Fig. 4), which, again, indicates that most of the variation in FG numbers occurs within, not among subclades (Figs 2 and 3).

Discussion
Our study provides the first global-scale analysis investigating the phenotypic diversity, phylogenetic distribution and evolutionary diversification dynamics of the FG used by squamate reptiles in chemical communication. Our  results reveal that the degree of phylogenetic signal underlying variation in the number of FG is moderate across the squamate phylogeny. However, the strength of this phylogenetic signal varies among subclades, with some squamate lineages showing higher evolutionary lability than other clades. The diversification dynamics of FG numbers across lineages is described by a stabilizing selection model of evolution (OU). The evolution of relative disparity in FG numbers lies above the values expected from a Brownian motion model, which implies that subclades overlap with one another in FG morphospace, indicating a substantial degree of evolution of similar 'gland strategies' across different squamate clades. Our ancestral state reconstructions coupled with diversification analyses suggest that despite the observed phylogenetic signal, FG have extensively diversified among lizards over time and across lineages, after stemming from a basal ancestor that likely lacked these glands, and from which these structures emerged and disappeared in repeated evolutionary episodes during the squamate evolutionary history. Squamate follicular glands play a paramount role in mediating and shaping conspecific and heterospecific interactions via the delivery of chemical signals for social and sexual communication [see reviews 13,28,33 ]. However, most information about the evolution of these glands known so far only comes from two lizard groups, lacertids and Liolaemus 30,34,35 . Indeed, patterns of morphological (e.g. tongue shape, number of sensory cells in their vomeronasal organ) and behavioral (number of average tongue-flicks and foraging mode) variation suggests that the evolution of chemo-sensation in squamates experienced a drastic episode of divergence early in the evolutionary history of this reptile group, followed by phylogenetic stability in the diversification of modes of chemical communication in the sublineages that subsequently originated from these crown ancestors [36][37][38][39] . Thus, Iguania is known to mostly be a 'visually-oriented' lineage, while the Scleroglossa are traditionally known to be 'chemically-oriented' organisms 36,38 . Remarkably, however, our findings posit a challenge to this longstanding hypothesis, given that the proportion of species with FG is in fact much higher in Iguania than in Scleroglossa.

FG number
Anatomical FG location  Therefore, our global-scale analysis reveals the need to hypothesize that chemical communication is likely to have played a central role during the evolution of social and interspecific interactions within the cosmopolitan clade Iguania and not only Scleroglossa 30,35,39,40 . Furthermore, some exceptionally-diverse lineages within Scleroglossa (e.g. skinks and snakes) have few or no species with FG (1.11% and 0, respectively), reinforcing the need to re-assess the currently accepted conclusions about the role of chemical communication in the global radiation of squamate reptiles as a whole (which accounts for >96% of living reptiles 41 ). However, it is important to point-out that delivery of information via chemical signals is rarely (if ever) the result of one single source of structures that produce and deliver chemical signals into the environment. In fact, vertebrates in general have evolved multiple approaches to produce signals through the skin, feces, body fluids (e.g. saliva, urine), and vaginal secretions, among others. However, despite the wide diversity of systems for chemical signaling that have independently proliferated among different lineages of vertebrates, it is interesting to note that glands such as FG and the equivalent structures in a range of mammals show some important level of phenotypic conservatism within lineages. For example, while a range of studies have reported extensive variation in the chemical composition of FG signals across closely related species of lizards 32, 42 , our study confirms that there is a considerably degree of phylogenetic   signal that underlies the expression of the glands' structures across species. The same remains true for other distantly related lineages such as rodents. Therefore, although our findings reveal extreme diversification of FG across the squamate tree of life, it seems clear that selection operating on these structures is independent from selective pressures driving diversification of the chemical structure of the signals themselves. Scincoidea and Iguania, in particular, were found to show phylogenetic conservatism in the variation of FG numbers across subclades, a finding that is consistence with a previous study 30 . In Gekkota and Lacertoidea, the phylogenetic signal is only moderate, which highlights the asymmetries in the strength of phylogenetic inertia and adaptive lability of chemical glands across the squamate tree of life. In fact, a recent study conducted on lacertid lizards revealed a moderate effect of phylogenetic ancestry shaping variation in FG, but also, a moderate predictability of variation in FG as a function of differential occupation of substrates among species 31 . In addition, our results identified stabilizing selection as the best approximation explaining the diversification of FG numbers, a finding that is further supported by the strong subclade overlap in morphospace also revealed by the DTT analysis. Therefore, and although shared ancestry remains an important factor underlying FG numbers, our results suggest that the number of these glands gravitates towards optimal values in both squamates as a whole, and in Gekkota, Iguania and Lacertoidea in particular.
Our analyses reveal that the ancestral state in squamates (and in lizards) is the lack of FG (see Table 1, for multiple ancestral states with and without FG across squamates). In view that the absence of FG seems to be the ancestral state for this clade, the subsequent appearance of these glands might have arisen in response to the need to engage in communication via an alternative signaling channel. Thus, given the high extent of diversity of signaling methods known among squamates (sounds, colours, behavioural displays, chemicals, among others), the differential degree of investment of each species into chemical signaling could condition the final expression (i.e. number and location) of the FG. For example, some lineages (e.g. snakes), in which the lack of FG is well known 13,27 , are likely to have satisfied their communication demands via the specialization of alternative sources of chemical signals, such as the skin or feces. Indeed, multiple studies have shown the high physiological costs associated with the production and maintenance of systems for chemical communication in squamates 27,33,34 . Therefore, species with multimodal signaling systems (acoustic, visual or/and chemical), might be more constrained in the production and maintenance of these signaling traits than those that base their communication on one mode of communication only. In addition, to be dependent on other communication channels, FG are also closely interrelated with their chemical secretions. Accumulating evidence reveals the extraordinary qualitative and quantitative diversity of compounds present in these secretions 13,43,44 . Therefore, the production and maintenance of chemical signaling structures that improve the efficiency of the chemosensory system do not only depend on the FG, but also on the specific level of energetic allocation into the secretions 42 . Consequently, the evolution of chemical communication may be primarily targeted by selection operating on the efficiency of the signal itself rather than on the FG.
While ancestry is revealed as an important factor underlying the global patterns of variation in FG (including both number and location), the factors (e.g. environmental demands) underlying the observed tempo and mode of diversification presented by our study remain largely an open question. The strength of our global-scale study provides the first compelling overview that reinforces the view that different lineages have evolved different degrees of adaptive lability underlying the diversification of these glands. Some authors have suggested the hypothesis of between-channel compensation 31 in which, given harsh environmental conditions, species might increase investment in additional or alternative signaling channels that are likely to promote changes in the evolutionary direction of the existing sensory channel (e.g. FG), leading to shifts in numbers, origins or losses of them. The field of chemical communication in reptiles offers a plethora of open questions and our study provides the fundamental empirical basis to guide the directions of further studies investigating the evolutionary dynamics of these interactions among one of the most diverse groups of tetrapods on Earth.

Material and Methods
Data Collection. We assembled a global dataset on the presence, number and location of FG for 7,904 species of squamates from the literature. These data cover 94% of all squamate families and over 80% of all species (see Table S1). To guarantee a comprehensive account of the phylogenetic distribution of the variation of these glands, and to inform the phylogenetic models about where the traits exist and where they have been lost, our data include species with and without them. Also, given that females of many species lack FG, we focused on males only. For each species, we obtained the mean pore number (i.e. pores on both left and right thigh) calculated from the average of multiple samples or as the midpoint between the minimum and maximum number of pores depending on the kind of information made available in the literature. When multiple independent sources provided information for the same species, we averaged data provided by all published sources for a species. FG were classified based on their location as 'precloacal' (located on the edge of the cloacae), 'femoral' (on the ventral surface of the thighs), both (when a continuous row of glands expands from one hind limb to the other through the cloacae area), or neither (glands are absent).

Comparative analyses.
In all our analyses, we used Pyron et al. 's 45 molecular phylogenetic supertree for 4,161 squamate species. Of the 7,094 species for which we had FG data, we were able to include 3,533 species into this tree.
To assess evolutionary patterns in the number and location of FG, we first tested if the average number of FG differed depending on their anatomical location, using phylogenetic analyses of variance ('phylANOVA') 46 . We then examined potential differences in FG number among the four main taxa where FG are known to be present (i.e. Gekkota, Iguania, Lacertoidea and Scincoidea) based on non-phylogenetic GLMs. For these analyses, we deliberately excluded the Dibamidae family given that this small family of squamates has only one species with FG. Pairwise comparisons were based on Tukey's HSD tests in all cases 47 . All statistical tests were performed using R 3.2.2 and SPSS 20.0.0 software. Phylogenetic signal and ancestral state. We estimated the phylogenetic signal for FG number and location in all squamates. Generally, phylogenetic signal is recognized to be the tendency for related species to resemble one another for a specific trait, and Pagel's λ and Blomberg K are two quantitative measure of this pattern 48,49 . A λ-value close to 0 indicates no phylogenetic structure in the trait, whereas a λ-value close to 1 indicates an increasingly stronger effect of shared ancestry on the expression of the trait across the species in the tree 48 . A K-value lower than 1 implies that relatives resemble each other less than expected under Brownian motion evolution along the hypothesized tree, whereas K > 1 implies that closely related species are more similar than expected under Brownian motion evolution. The signal of the continuous variables was assessed using Pagel's λ and Blomberg's K (with nsim = 1,000) and the 'phylosig' function in the 'phytools' package 46 . For the discrete variable 'anatomical location of the pores' , we calculated Pagel's λ only, using the 'fitDiscrete function' in the 'geiger' package 50 .
Additionally, we used ancestral character state reconstructions ('ace' function in the 'ape' package 51 ) to quantitatively estimate ancestral states of FG location at the root of the Squamata tree, and at specific internal nodes of the different lineages within Squamata (i.e. Sauria, Gekkota, Iguania, Lacertoidea and Scincoidea). Finally, we employed the approach implemented by Revell and Freckleton 52 to visualize the dynamic trajectories of ancestral states of FG numbers along branches of the phylogenic tree. These analyses use a Maximum-Likelihood approach to reconstruct ancestral states.

Modelling FG evolution.
To examine the diversification of FG number across Squamata, we used two approaches aimed to elucidate both global and lineage-specific macro-evolutionary diversification dynamics. We firstly analyzed the whole global dataset and then, we focused on the three lizard lineages with the largest number of species with FG in our dataset (i.e. Gekkota, Iguania and Lacertoidea). We fitted four alternative evolutionary models that describe different regimes of phenotypic evolution: i) the Brownian-motion model (BM) according to which evolution proceeds as a random walk through trait space; the expected phenotypic difference between two taxa grows proportional to the time since their divergence 53 ; ii) the Ornstein-Uhlenbeck model (OU) which describes the evolution of traits under stabilizing selection -it modifies the BM model by considering one or several adaptive optima to which trait evolution is attracted as diversification proceeds during the phylogenetic history of the lineage 54 ; iii) the Early-Burst model (EB) or 'accelerating-decelerating model' , which assumes that SCiENtiFiC RepoRts | 7: 9288 | DOI:10.1038/s41598-017-09083-7 species evolve under a density-dependent availability of niche space, predicting rapid early diversification followed by decreases in evolutionary rates as a result of saturation of niche space over time [55][56][57] ; iv) the Delta model which allows the evolutionary rate to either decrease (δ < 1) or increase (δ > 1) through time, from root to tips. In a model with δ < 1, the trait changes rapidly early in the history of a clade and then slows down through time. If δ > 1, trait change accelerates through time 48,53 . We performed the comparisons of goodness of fit for these four models based on Akaike Information Criterion (AIC) 58,59 . We use the bias-corrected version of AIC (AICc) to determine the ΔAICc values, which result from the difference between the lowest AICc and the AICc of each alternative model. Therefore, the best-fit model has ΔAICc = 0 53, 60 . All these analyses were conducted using the 'geiger' package in R 50 .
We further investigated the disparity of the FG number over time. We performed disparity-through-time analyses (DTT) by plotting fluctuations in average relative disparity over time [61][62][63] . In such analyses, the mean disparity in the trait (here: FG number) at any node along the clades' history is compared to the disparity under the null model of Brownian motion (estimated by the mean of 10,000 simulations). The average of these two values (i.e. FG number disparity from dataset and simulated data from BM) are plotted against node age to obtain the morphological disparity index (MDI) 53,62,64 . Thus, values below zero (i.e. those values lower than expected under BM model) indicate that most of disparity is among subclades, which are distributed in smaller and isolated morphospace regions. Instead, MDI values above zero mean that disparity among subclades is highly overlapped in the morphospace 53 . We carried out DTT analyses using the R package 'geiger' 50 . Subsequently, we finally built a "traitgram" plotting the Squamata phylogeny onto the FG number morphospace over time since the root origin. The resulting projection is based on ancestral node estimations using maximum likelihood approaches 65 . These analyses were performed using the R package 'phytools' 46 .