The shadow of the past: Convergence of young and old South American desert lizards as measured by head shape traits

Abstract Convergence is a pervasive phenomenon in the Tree of Life, and evolution of similar phenotypes sharing the same environmental conditions is expected in phylogenetically closely related species. In contrast, contingent factors are probably more influential in shaping phenotypic diversity for distantly related taxa. Here, we test putative convergent evolution of lizard head morphologies among relatively closely related desert dwelling Liolaemus species, and the very distantly related Ctenoblepharys adspersa. We estimated a multilocus time‐calibrated phylogeny of 57 species of South American liolaemus lizards, based on seven molecular markers. We collected head shape data for 468 specimens, and used three phylogenetic comparative methods (SURFACE, CONVEVOL, and WHEATSHEAF index) to test for and estimate the strength of convergence. We found strong evidence for convergence among Pacific desert lizard C. adspersa, Liolaemus audivetulatus, Liolaemus insolitus, Liolaemus poconchilensis, Liolaemus stolzmanni, and a candidate species (Liolaemus “Moquegua”). Our results suggest that, despite the long divergence and phylogenetic distance of C. adspersa with respect to convergent Liolaemus species, natural selection was probably more important than historical contingency in shaping phenotypic evolution in these desert lizards.


| INTRODUC TI ON
Evolutionary convergence is a pervasive phenomenon in the Tree of Life and can be defined as the repeated, independent evolution of the same trait (or complex of traits) in two or more clades (McGhee, 2011). Two possible goals in the study of evolutionary convergence are its identification (whether convergence is present) and quantification (estimating its frequency and strength; Arbuckle & Speed, 2017). The frequency of convergence is quantified by enumerating the cases in a group of taxa, while the strength of convergence estimates how similar is/are the trait(s) of the convergent taxa while taking phylogenetic distance into account.
Another important issue in studying convergence is whether natural selection and constraints can erase the contingent nature of evolution. Natural selection and constraints should be more prominent in shaping similar adaptive phenotypes of related species (e.g. species of the same genus), but historical contingencies may lead to a more significant imprint when taxa are not closely related (Ord & Summers, 2015). However, because the number of possible forms is finite, even phylogenetically distant taxa will evolve the same adaptations when exposed to the same selective pressures. So, under similar environmental conditions, closely related as well as distantly related taxa may either (a) evolve similar phenotypes independently of past events; or (b) not show repeated evolution because of historical contingencies (Ord & Summers, 2015).
Here, we test for convergence in South American desert lizards of the Liolaemus montanus group (Figure 1). Some species in this group have toad-like ("phrynosauroid") head shapes and pronounced serrated combs formed by the projecting outer ciliary scales, which are lacking in other closely related members of this species group ( Figure 2). Other species and populations have a similar toad-like head, but lack the pronounced serrated combs. These "phrynosauroid" lizards and similar forms inhabit the extremely arid desert environments (mean annual rainfall ~1-15 mm) of the South American Pacific coast (Rundel, Villagra, Dillon, Roig-Juñent, & Debandi, 2007), in contrast to the remaining, mostly Andean, species of the L. montanus group. Moreover, these species resemble distantly related taxa present in the same arid deserts, the monotypic Ctenoblepharys adspersa (Tschudi, 1845) (Liolaemidae).
The taxonomic history of these desert species is a good example of how putative convergence has confused taxonomists.
In more recent studies, these "new" genera were later rejected, and all species except C. adspersa were returned to Liolaemus (Etheridge, 1995;Lobo, Espinoza, & Quinteros, 2010). Moreover, there is not a comprehensive phylogeny of the L. montanus group which includes "Phrynosauroid" lizards and related forms, thus preventing a formal test of convergence and its strength.
The arid conditions and sandy substrates of the Peruvian and Atacama deserts likely have exerted strong selective pressures for the evolution of convergent phenotypes among these taxa. Although toad-like head shapes with enlarged ciliaries in "Phrynosauroid" lizards have been mentioned together as convergent traits for living in desert conditions (Valladares, 2004), they may reflect different aspects of their natural histories. Toad-like head shapes may be related to a diet composed largely of small preys, as documented for North American Phrynosoma lizards (Meyers, Herrel, & Nishikawa, 2006;Meyers, Nishikawa, & Herrel, 2018), in comparison with other desert living and hard-preyed specialist lizards of the genera Gambelia and Crotaphytus (Modlin, 2018). However, dietary data for C. adspersa and putatively convergent Liolaemus are limited. A high percentage of Hymenoptera and small Coleoptera was found in C. adspersa and L. insolitus stomach contents, respectively (Cei & Pefaur, 1982;Pérez & Balta, 2007). Given that small preys comprise most of the diets in C. adspersa and the putatively convergent Liolaemus, the repeated evolution of their distinct head morphology may reflect a loss of morphological traits related to the acquisition and processing of a more generalist diet.
In contrast, strongly enlarged ciliaries may protect the eyes from sand (Etheridge, 2000) while burrowing, or inside a burrow.
Burrowing in loose sand is known in this convergent group in C. adspersa (CA, personal observation), and phrynosauroid Liolaemus may use the projecting ciliaries to deflect sand from their eyes when they use burrows (J. Troncoso-Palacios, personal observation; Díaz-Vega, 2015).
The aims of this paper are to: (a) test the monophyly of "phrynosauroid" lizards and related forms of the L. montanus group and estimate their divergence times using a fossil-calibrated multilocus dataset; and (b) test for phenotypic convergence in Southern Pacific desert lizards using quantitative head traits. If the phylogenetic signal of natural selection overrides historical contingency, we should expect C. adspersa and "phrynosauroid" forms of the L. montanus group be very similar in head traits. Alternatively, if stochastic events dominated the evolutionary history of this clade, we should expect "phrynosauroid" lizards and related forms of the L. montanus group to be more similar to each other than to C. adspersa.

| Phylogenetic analyses
Details of field sampling, laboratory procedures, specimen assignment and locality data are summarized in Supporting Information Appendix S1. Sequences were aligned in the MUSCLE (Edgar, 2004) plugin in GENEIOUS®PRO v5.6.6 (Kearse et al., 2012), and protein coding sequences were translated to check for premature stop codons. Bayesian information criteria in JMODELTEST v2.1.3 (Darriba, Taboada, Doallo, & Posada, 2012) was used to identify the best-fit models of molecular evolution. Maximum likelihood (ML) analyses were performed using RAXML v.8 (Stamatakis, 2014) partitioned by gene, and 1,000 standard bootstrap replications were estimated using the rapid hill-climbing algorithm (Stamatakis, Hoover, & Rougemont, 2008) in the Cyber Infrastructure for Phylogenetic Research (CIPRES; Miller, Pfeiffer, & Schwartz, 2010). The ML tree is shown in Supporting Information Figure S1.
To estimate divergence times, a concatenated Bayesian tree (BT) was generated and calibrated as in Aguilar, Wood, Belk, Duff, and Sites (2017). This analysis was implemented in BEAST v2.4.6 (Bouckaert et al., 2014) using the bModelTest (Bouckaert  & Drummond, 2017) to simultaneously generate and explore the model substitutions across space to estimate model parameters, and to generate a fossil-calibrated phylogeny. We implemented 10 independent runs for 100 million generations, sampling every ten thousand generations, and we checked for convergence of the runs using TRACER v1.6 (Drummond, Suchard, Xie, & Rambaut, 2012) to ensure effective samples sizes (ESS) were >200. A Yule speciation prior and a log-normal relaxed clock model were applied. We discarded 10% of the trees and log files as burn-in, and the remaining trees were combined using LOGCOMBINER v1.8.0 and sampled at a lower frequency, resulting in 10,000 trees. A maximum clade credibility tree (MCC) was then constructed using TREEANNOTATOR v1.8 , and keeping mean and 95% confident intervals for node ages.

| Morphological data and analyses
Shape analyses were performed using principal component (PC) analysis after a Generalized Procrustes approach. Procustes and PCA analyses were performed using MORPHOJ v1.03d (Klingenberg, 2011), and PC scores were extracted using the GEOMORPH package (Adams & Otarola-Castillo, 2013) in R (R Development Core Team, 2014). The first three principal components explained 74.8% of the variance and were retained for further analyses (Supporting Information Appendix S3). We also estimated average of PC1-PC3 scores for specimens representing each species (Supporting Information Appendix S3). The use of PC scores in phylogenetic comparative studies has been criticized when only a few PC axes have been selected (Adams & Collyer, 2017;Uyeda, Caetano, & Pennell, 2015), but it is the best available approach for our purposes. Geometric morphometric data are necessarily reflecting one or more multivariate traits (Collyer, Sekora, & Adams, 2015); these are usually reduced using PC analysis, and the first PC axes that explain most of the variance are employed instead of the original data (e.g. Muschick, Indermaur, & Salzburger, 2012;Esquerre & Keogh, 2016).
To examine convergence traits in the L. montanus group, the BT was combined with head shape data. Ten landmarks on the dorsal head view (Aguilar et al., 2017; Supporting Information Appendix S3) of 468 lizards representing 57 species (Supporting Information Appendix S3) were set on digital pictures using TPSdig v1.4 (Rohlf, 2004). Landmarks (number in parentheses) were set on the tip of the snout (1), nostrils in each side (2, 3), beginning of the first superciliary scale in each side (4, 6), end of the last superciliary scale in each side (5, 7), interparietal scale (8), dorsal widest part of the head at the level of anterior margin of the external ear in each side (9, 10). These landmarks were selected to represent shape differences in different, but homologous areas of the head. Additional homologous landmarks were difficult to set because they were not repeatable across the very distantly related taxa (e.g. C. adspersa and other Liolaemus species).
We tested for deviation of head symmetry using Procrustes ANOVA in MORPHOJ v1.03d (see Supporting Information Appendix S3 for details). For each species, only adults of each sex (as assessed by larger snout vent length) were used to avoid ontogenetic allometric bias. When sample sizes within species allowed comparisons between sexes, we tested for sexual dimorphism for each species using the PCA scores of the first two PCs. Because males and females did not cluster separately in each species tested, data of both sexes were pooled together for further analyses. We also performed a discriminant function analysis using Procrustes coordinates to test for shape differences between sexes and all species, but we did not find strong significant differences between males and females (Supporting Information Appendix S3). Because specimens of Liolaemus lentus and Phymaturus sitesi were not available, we used another species of the same group (Liolaemus pseudoanomalus) or same genus (Phymaturus patagonicus), respectively.

| Convergence analyses
The BT was pruned using R package APE v5.1 (Paradis, Claude, & Strimmer, 2004) to include only taxa for which morphological data were available (see above). Evolutionary changes in head shape were reconstructed by using the PC scores onto the BT phylogeny, and then employing a squared-change parsimony analysis (Maddison, 1991) method in MORPHOJ. This method was selected to visualize the evolutionary changes between the reconstructed head shape of the focal taxa and their most recent common ancestor by deformation grids and displacements of landmarks.
Three convergence analyses were performed using PC1-PC3 axes. The first analysis was performed using the R package SURFACE (Ingram & Mahler, 2013); this algorithm employs an Ornstein-Uhlenbeck (OU) process to identify cases of convergence without the a priori designation of convergent taxa. The method has a forward phase in which selective regimes are inferred using a phylogenetic tree and quantitative traits, and a reverse phase in which taxa having the same (convergent) regime are identified. In the forward phase, selective regimes are added to a Hansen model (Hansen, 1997) and then further regime shifts (models) are added in a stepwise process. Model performance is evaluated using corrected Akaike information criterion (AICc). In the reverse phase, all selective regimes obtained in the first phase are combined in a pairwise manner and collapsed into a shared regime. This procedure is repeated until no more stepwise combinations improve the models, and convergent (collapsed) regimes are estimated (again using AIC).
The SURFACE model with convergent regimes (OUc) was compared with a model with multiple non-convergent regimes (OUnc), and simpler stochastic models such a single regime (OU1) and Brownian motion (BM) models.
Focal taxa (C. adspersa, L. audituvelatus, L. insolitus, L. stolzmanni, L. poconchilensis and Liolaemus "Moquegua"), and two subsets (one excluding L. "Moquegua," and another excluding L. insolitus and L. "Moquegua") were used to estimate convergent metrics with CONVEVOL (Stayton, 2015a). This method estimates four distances (C1, C2, C3, C4) and one frequency-based (C5) measure of convergence. C1 is based on the idea that the more dissimilar the ancestors, and the more similar the descendants, the stronger is the convergence. C1 represents the proportion of the maximum distance between two lineages that has been brought together by subsequent evolution, and ranges from 0 to 1 as convergence increases.
A value of 1 indicates that lineages are fully convergent, and a value of 0 means that lineages are phenotypically different, and convergence is absent. C2 is another measure representing the absolute amount of evolution that has occurred during convergence, with larger values indicating greater convergence. C3 and C4 are based on C2 and allow comparison between datasets (in contrast to within datasets). C3 is the proportion between C2 and the total amount of evolutionary change along the lineages leading from the common ancestor of the convergent taxa to those taxa. C4 is the proportion between C2 and the total amount of evolution in the entire clade defined by the common ancestor of the convergent taxa (Stayton, 2015a).
C5 is a frequency-based measure and is defined as the number of focal taxa that reside within a limited but convergent region of a phylomorpho-space (the phylogenetic connections between taxa represented graphically in a plot of morphological space).
Statistical tests of convergence as measured by C1, C2, C3, and C4 were evaluated using 500 evolutionary simulations via a BM model. Specifically, we tested whether the simulated measures are significantly different from the observed values. In the same way, F I G U R E 3 Divergence time tree including species of the Liolaemus montanus group and Ctenoblepharys adspersa. Putative convergent taxa are in red. Red and black dots on nodes indicate <0.95 and ≥0.95 posterior probabilities respectively. Grey bars are 95% confidence intervals for node ages the statistical significance of convergence as measured by C5 was tested using 500 simulations. Results of all tests were considered significant at a p-value ≤0.05.
We implemented a WHEATSHEAF analysis to measure the strength of convergent evolution in focal taxa and subsets (see above), as implemented in the R package WINDEX (Arbuckle & Minter, 2015;Arbuckle, Bennett, & Speed, 2014). This index calculates the similarity of focal (convergent) species to each other and the separation in phenotypic space of the focal group from non-convergent species, all corrected for phylogenetic relatedness.
Convergence is stronger when focal species are more phenotypically similar to each other but phylogenetically more distant, relative to the non-focal species. Convergent focal taxa in SURFACE and CONVEVOL were used to estimate the WHEATSHEAF Index and 95% confidence intervals. The null hypothesis that the observed WHEATSHEAF index is no higher than expected by chance is rejected when p ≤ 0.05 (indicating exceptionally strong convergence).
Expected WHEATSHEAF indexes were derived from 1,000 bootstrap replications.

| Phylogenetic relationships and divergence times
Main results of the Bayesian divergence time tree (BT; Figure 3) with posterior probabilities (PP) are shown, and differences with the maximum likelihood (ML) tree (Supporting Information Figure S2) are summarized below. Detailed results are given together with Supporting Information Figure S2    than the other subsets ( Table 2). The C5 metric shows that C. adspersa, L. audituvelatus, L. insolitus, L. poconchilensis, L. stolzmanni and L. "Moquegua" and the two subsets (Table 2) significantly cluster in a distinct region of the phylomorphospace (Table 2, Figure 6).

| Convergence analyses
The WHEATSHEAF index was higher for the subset C. adspersa, L. audituvelatus, L. poconchilensis and L. stolzmanni (Table 2) than for the other subset and all focal taxa. However, results were not significant for all of them. F I G U R E 6 CONVEVOL phylomorphospace of 57 species in principal components (PC) 1 and PC2. The black lines connect only non-focal species and red arrows connect nodes to six focal convergent species. Ctenoblepharys adspersa (Ca), Liolaemus audituvelatus (La), Liolaemus insolitus (Li), Liolaemus poconchilensis (Lp), Liolaemus stolzmanni (Ls) and Liolaemus "Moquegua" (LM). The violet curve defines a distinct region of the phylomorphospace where these six taxa are present convergence implies that natural selection has produced the same phenotype in similar environments in unrelated taxa (Losos, 2011).

| D ISCUSS I ON
Constraints (developmental and phylogenetic) imposed by ancestors are enforced by natural selection and thus by the circumstances under which organisms develop and evolve (Vermeij, 2015).
Natural selection and constraint are probably more important than historical contingency in explaining the repeated adaptations to the same habitat in phylogenetically close taxa (Hagman & Ord, 2016). In contrast, historical contingency may be more important in producing unique adaptations in phylogenetically distant taxa (Hagman & Ord, 2016). Natural selection is thought to erase the stamp of historical contingency when comparing species of the same genus or populations within a species, but not so in highly divergent taxa (Ord & Summers, 2015). In addition, evolutionary paths that converge on a local fitness optimum over a short time period will not necessarily achieve the highest adaptive peak over a longer period because of stochastic events (Orgogozo, 2015).
At least one of the three methods applied in our study, SURFACE, is based on the assumption that convergence is adaptive (Arbuckle et al., 2014;Ingram & Mahler, 2013). SURFACE employs the metaphor of Wright and Simpson that evolution is a local search by fitter geno/ phenotypes climbing higher peaks in an adaptive landscape (Arnold, Pfrender, & Jones, 2001;Niklas, 1995). This method identifies convergent taxa via evolutionary models and results are interpreted as occupation of the same or very similar adaptive peaks, which is then taken to characterize a selective regime (Ingram & Mahler, 2013;Mahler et al., 2017). In contrast, the C1-C5 metrics test whether convergence in putative taxa is different from chance without invoking any evolutionary process a priori (Stayton, 2015a(Stayton, , 2015b. However, there is still controversy about whether phylogenetic comparative models used to identify and measure convergence, are or are not process free (Mahler et al., 2017;Stayton, 2015b).
Our SURFACE analysis identified 13 convergent peaks and one of them reached by our six species focal species: C. adspersa, L. audituvelatus, L. insolitus, L. poconchilensis, L. stolzmanni and L. "Moquegua" ( Figure 5; Table 1). Convergence of these six taxa is strongly supported by higher C1-C4 metrics, in comparison with a subset including only taxa with highly modified head shape (L. insolitus excluded) and another including only taxa with serrate combs in the eyes (L. insolitus and L. "Moquegua" excluded; Table 2). Both results suggest a significant case convergence for these Pacific desert species reaching the same adaptive peak.
Convergent Liolaemus and C. adspersa show a lateral widening, and posteromedial stretching of the head (Figure 4c-h). Assuming evolutionary modifications in head shape in these convergent taxa are related to a loss of bite force and reduction in the morphology of prey-processing features, this reduction in traits may be a response of disuse driven by natural selection (Fong, Kane, & Culver, 1995;Porter & Crandall, 2003). The benefit of an unused character may be outweighed by its fitness cost in the current environment, but reduction of an unused trait may be associated with increased fitness (Hall & Colegrave, 2008). In environments with low resources, the fitness cost imposed by an unused trait is greatest (Hall & Colegrave, 2008 Liolaemus are morphologically similar in head shape, but not identical (Figure 4c-h). In addition, Pacific desert Liolaemus are not phylogenetically distant between each other (Figure 3).
The origin of the hyper-aridity in the Peruvian and Atacama deserts (~25 My) is consistent with the ages of focal Liolaemus clades in our time-calibrated tree (Figure 3; Rundel et al., 2007). In addition, our BT resolves the origin of this head shape first in C. adspersa and three times in the L. montanus group ( Figure 5). This suggests that, in addition to selection, phylogenetic constraint may have also been involved in producing the same phenotype in these desert lizards. However, information is lacking on whether similar or identical mechanisms (e.g. same developmental pathway, same genes, or same specific-site mutation in the same genes) have been responsible for the independent origin of these same phenotypic traits. Results of such research would clarify the role of constraint in the evolution of adaptation (Agrawal, 2017;Arendt & Reznick, 2007).
Our results do suggest, however, that despite its deep divergence and thus potential exposure to different contingent factors, natural selection has been stronger than historical contingency in producing the same convergent traits in the phylogenetically distant C. adspersa and younger species of the L. montanus group (for traits examined in this study). However, historical contingencies probably were also important when considering these traits at higher taxonomic levels. For instance, the convergent taxa in our study have been hypothesized to be phenotypically similar to the phylogenetically very distantly related lizard Phrynocephalus (family Agamidae; Valladares, 2004), suggesting than a toad-like head may have been present in an older lizard ancestor (~168 My; Zheng & Wiens, 2016).
Although a formal test of convergence is needed, both lizard families have probably been exposed to different contingent events and may also have unique adaptations for living in different desert environments (Arnold, 1995;Melville, Harmon, & Losos, 2006).

ACK N OWLED G M ENTS
We thank Herman Núñez (

DATA ACCE SS I B I LIT Y
New sequences were deposited in GenBank and their numbers are given in Supporting Information Appendix S2.