Evolution of developmental sequences in lepidosaurs

Background Lepidosaurs, a group including rhynchocephalians and squamates, are one of the major clades of extant vertebrates. Although there has been extensive phylogenetic work on this clade, its interrelationships are a matter of debate. Morphological and molecular data suggest very different relationships within squamates. Despite this, relatively few studies have assessed the utility of other types of data for inferring squamate phylogeny. Methods We used developmental sequences of 20 events in 29 species of lepidosaurs. These sequences were analysed using event-pairing and continuous analysis. They were transformed into cladistic characters and analysed in TNT. Ancestral state reconstructions were performed on two main phylogenetic hypotheses of squamates (morphological and molecular). Results Cladistic analyses conducted using characters generated by these methods do not resemble any previously published phylogeny. Ancestral state reconstructions are equally consistent with both morphological and molecular hypotheses of squamate phylogeny. Only several inferred heterochronic events are common to all methods and phylogenies. Discussion Results of the cladistic analyses, and the fact that reconstructions of heterochronic events show more similarities between certain methods rather than phylogenetic hypotheses, suggest that phylogenetic signal is at best weak in the studied developmental events. Possibly the developmental sequences analysed here evolve too quickly to recover deep divergences within Squamata.


INTRODUCTION
With over 10,000 species, Squamata (lizards, snakes and amphisbaenians) are one of the most species-rich extant tetrapod lineages (Uetz, Freed & Hošek, 2016). However, our understanding of their evolutionary history is confounded by the conflict between phylogenetic hypotheses based on morphology and molecular data (e.g., Losos, Hillis & Greene, 2012). Morphological analyses suggest that the first divergence within Squamata was between Iguania (iguanas, agamas, chameleons and kin) and Scleroglossa (all other lizards and snakes) (e.g., Estes, de Queiroz & Gauthier, 1988;Conrad, 2008;Gauthier et al., 2012), while molecular studies indicate that iguanians are highly derived lizards, closely related to anguimorphs (e.g., monitor lizards) and snakes, and that limbless dibamids or gekkotans (geckos and kin, sometimes also including dibamids) are the first-diverging branch of squamates (e.g., Townsend et al., 2004;Vidal & Hedges, 2005;Wiens et al., 2010;Wiens et al., 2012;Pyron, Burbrink & Wiens, 2013). Increasing the number of taxa and characters in these analyses has not led to an improvement of our understanding of squamate phylogeny, but rather has only increased the discordance between the hypotheses based on those two lines of evidence. Combined morphological and molecular analyses (e.g., Wiens et al., 2010;Reeder et al., 2015) generally favour the molecular topology (but see Lee, 2005). However, some authors argue that molecular data may not be ideal for resolving the higher-level phylogeny of squamates because of the large genetic distance between squamates and their closest living relative-the tuatara (Sphenodon punctatus), the only extant rhynchocephalian-and thus the only reasonable proximal outgroup to Squamata in phylogenetic analyses (McMahan et al., 2015). Despite numerous publications on this subject (Gauthier et al., 2012;Losos, Hillis & Greene, 2012;Reeder et al., 2015), the debate continues and still new approaches to the problem are being taken (McMahan et al., 2015;Harrington, Leavitt & Reeder, 2016;Pyron, 2017).
Developmental data may be useful for phylogenetic inference (e.g., Laurin & Germain, 2011) but they rarely have been used in squamate phylogenetics. Notable exceptions are the studies of Maisano (2002) and Werneburg & Sánchez-Villagra (2015), using ossification sequences. The former found that these sequences are useful for determininig relatively shallow divergences but failed to recover deeper nodes, possibly because of their high rate of evolution (Maisano, 2002). Werneburg & Sánchez-Villagra (2015) found that developmental data were most congruent with the close relationship between snakes and varanids, as postulated by some morphological studies (e.g., Lee, 1997) but also some combined morphological and molecular analyses (e.g., Lee, 2005). Sequences of other developmental traits were studied by Andrews, Brandley & Greene (2013) but the authors regarded relationships of squamates as ''well defined'' and reconstructed the ancestral states only on the molecular topologies. Moreover, their study did not consider the tuatara, a taxon critical in studying lepidosaur evolution. We attempt to supplement their data with the developmental sequence of the tuatara and reconstruct ancestral states using both molecular and morphological topologies. We also assess phylogenetic utility of timing of organogenesis using several different methods.

Character construction and cladistic analyses
Developmental sequences of 20 characters in 21 species representing most major squamate lineages (Tables 1 and 2) were obtained from Andrews, Brandley & Greene (2013). Developmental sequences of seven other squamate species were taken from the literature ( Table 1). The developmental sequence of the tuatara was compiled from Dendy (1899) and Sanger, Gredler & Cohn (2015) (see also Moffat, 1985). These sequences were transformed into continuous characters, where the first event has a value of 0, and the last one-a value of 1 (Germain & Laurin, 2009;Laurin & Germain, 2011). These values constituted the basis for cladistic characters, which were created following Werneburg & Sánchez-Villagra (2015)-values between 0 and 0.09 were coded as 0, between 0.1 and 0.19 were coded as 1, and so on. The missing data were coded as unknown (?), while limb characters in snakes were coded as inapplicable (-). Cladistic analyses employing these characters were conducted in TNT v. 1.1 (Goloboff, Farris & Nixon, 2003;Goloboff, Farris & Nixon, 2008), using the traditional search option, with 10 replications of Wagner trees. These trees were held in RAM and subjected to tree bisection reconnection, holding 10 trees per replicate.
In the first analysis, all characters were unordered (non-additive), and in the second one, all were ordered (additive) (see Werneburg & Sánchez-Villagra, 2009;Laurin & Germain, 2011). In both analyses, Sphenodon was used as the outgroup. Another set of cladistic characters was created using the event-pairing method (Smith, 1997;Velhagen Jr, 1997;Jeffery et al., 2002a;Jeffery et al., 2005). Comparing 20 developmental events in 29 species resulted in 190 event pairs. These characters were analysed in the same way as continuous characters.
With these cladistic characters and files with both molecular and morphological topology in memory, Templeton test (Templeton, 1983) was performed in TNT (using a script written by Alexander Schmidt-Lebuhn: https://www.anbg.gov.au/cpbr/tools/templetontest.tnt). Four replications were conducted: using either ordered or unordered characters; and employing continuous or event-paired characters.

Ancestral state reconstruction and heterochronic events
Reconstruction of ancestral states was performed in Mesquite v. 3.2 . Developmental sequences were mapped on two competing phylogenetic hypotheses of lepidosaurs-first one, from Pyron, Burbrink & Wiens (2013), using seven nuclear and five mitochondrial genes, and the second one, from Gauthier et al. (2012), the largest morphological analysis to date. Ancestral states were reconstructed using both maximum parsimony and maximum likelihood for event-paired data and square-changed parsimony for continuous data. The branch length may have a significant effect on reconstruction of ancestral states (e.g., Andrews, Brandley & Greene, 2013;Boyd, 2015), so analyses using maximum likelihood and square-changed parsimony were performed on both molecular and morphological topologies. In the first analysis, all branches were given an equal  (Nydam & Caldwell, 2015) but teiids are universally accepted as gymnophthalmid sister group, so the oldest known teiid is used to provide a calibration point for gymnophthalmids in the analyses.
length (=1), while in the second, the branch lengths were calibrated to reflect the fossil record of a given group. The oldest-known fossil of a total group was used to calibrate the tree rather than that of a crown group (Table 3). Only fossils unquestionably placed within a given group were included. When the fossil record of a group was unknown (mostly in relatively recently diverged species), the branch length was set, arbitrarily, as 3. Square-changed parsimony reconstruction using continuous data was performed using root node reconstruction in PDAP:PDTREE module of Mesquite (Midford, Garland Jr & Maddison, 2011). This module calculates 95% confidence intervals (Garland Jr & Ives, 2000) for each character of a hypothetical ancestor of all taxa included in a tree (in this case, ancestral lepidosaur). A statistically significant heterochronic event occurs when a value of character state of a given taxon is beyond the confidence interval. In the second analysis, Sphenodon was pruned from the tree, and reconstruction was made for the ancestral squamate and compared to the values of terminal taxa. Event-pair synapomorphies were mapped on both topologies using synapomorphy mapping in TNT. These synapomorphies were subjected to event-pair cracking, following the procedure described in detail by Jeffery et al. (2002a). Only deviations from their methods are described below. Clades supported by only one event-pair synapomorphy, two synapomorphies involving four different events and so on were excluded because the number of developmental changes was insufficient for determining the background pattern and heterochronies. In the ordered dataset, when degree of change was ambiguous (e.g., from 0 to 1 or 2), a mean was taken (in this example, 1.5). Characters in which the direction of change could not be unambiguously reconstructed (i.e., from 1 to 0 or 2) were excluded from further analysis. This should not have significant effect on the analysis, as there was only a few such characters (Tables S1-S8). Only events with total relative change (TRC) beyond the 95% confidence interval calculated for the mean TRCs at a given node were regarded as heterochronic. This is more a conservative approach than the one taken by Jeffery et al. (2002a) but will make the analysis more comparable to the continuous analysis described above.

Cladistic analyses
Cladistic analyses conducted using the transformed continuous data generated trees that are not similar to trees obtained in either morphological or molecular analyses. Analysis using unordered characters yielded 214 most parsimonious trees (MPT; tree length = 109, consistency index = 0.560, retention index = 0.628), the strict consensus tree of which is almost completely unresolved. This analysis failed to recover clades of very closely related species such as Liolaemus (Fig. 1A). When all characters were ordered, it resulted in 174 most parsimonious trees (TL = 133, CI = 0.459, RI = 0.625). The strict consensus tree is mostly unresolved-the only groups that were monophyletic in all MPTs are Liolaemus, Tropidurus + Strophurus, Calyptommatus + Anolis and a clade including Uta, Agama, Furcifer, Mabuya, Gehyra, Chamaeleo and Zootoca. A 50% majority rule tree does not resemble published morphological or molecular phylogenies (Fig. 1B).
Similar to the continuous dataset, the event-paired data did not result in a topology matching any previously published phylogeny. Analysis of unordered characters generated 10 MPTs (TL = 185, CI = 0.530, RI = 0.552). In the strict consensus tree Furcifer and Varanus indicus are in trichotomy with the clade including all other squamates. This clade is divided into a group containing seven species of iguanians, gekkotan Strophurus, snake Thamnophis, scincoid Mabuya and lacertiform Zootoca, and the second group to which all other squamates belong ( Fig. 2A). Analysis using ordered characters yielded 16 MPTs (TL = 220, CI = 0.464, RI = 0.599). The strict consensus tree is poorly resolved but excluding Varanus indicus from it significantly improves resolution. After this, squamates are divided into two clades-the first one includes eight species of iguanians, Thamnophis and Mabuya, while the second group includes all other squamates (Fig. 2B).
Mapping of continuous characters indicates slight differences in tree length between morphological and molecular topologies. With all branches being assigned equal length (=1), the former is 1.49630768 steps long and the latter-1.51610078. With the fossilcalibrated tree, the morphological topology is 0.19257679 steps long and molecular-0.17638729. Mapping of unordered event-paired characters gives the molecular topology a length of 250 steps and the morphological-252 steps. With ordered characters, the molecular topology is 322 steps long, while the morphological is 327 steps long.
Neither replication of the Templeton test detected any statistically significant differences between morphological and molecular phylogenies under both present continuous and event-paired character datasets (p > 0.05 in all cases).

Developmental diagnoses
There are several event-pair synapomorphies diagnosing some higher-level taxa (i.e., family-level clades or higher). However, at least some of these groups are represented by only a few members (e.g., Anguimorpha, Scincoidea), so these apomorphies may in fact diagnose less inclusive clades (Table 4).

Heterochronic events
Inferred heterochronic events show more consistency between given methods than between phylogenies (e.g., event-paired data for morphological phylogeny are more similar to event-paired data for molecular topology than to continuous data for morphological tree). Only a few of these events are common to all methods and phylogenies (Figs. 3-14).

DISCUSSION
Developmental cladistic characters failed to recover topology similar to those based on other data (i.e., molecular or morphological). This was also found in similar studies (Maisano, 2002;Werneburg & Sánchez-Villagra, 2009;Werneburg & Sánchez-Villagra, 2015). This may be a consequence of uneven sampling of different squamate clades in the present Table 4 Event-paired developmental synapomorphies of higher-level squamate clades. Asterisk denotes synapomorphies present only in analysis using ordered characters, while plus denotes synapomorphies present only in analysis employing unordered characters.

Serpentes
(1) pharyngeal slits closed earlier than eyelid forms as a thin ribbon-like sheet of tissue*

Gymnophthalmidae
(1) jaw initiated simultaneous with maximum pharyngeal slits, (2) jaw initiated earlier than hemipenal buds form on cloacal lips, (3) jaw initiated earlier than three-segmented limb*, (4) pharyngeal slits closed simultaneous with hemipenal buds form on cloacal lips*, (5) pharyngeal slits closed simultaneous with threesegmented limb*, (6) jaw completion simultaneous with digits differentiated in limb paddle Figure 3 Heterochronic events in lepidosaur evolution. Mapped onto molecular phylogeny, using continuous data, in relation to the ancestral lepidosaur. Length of all branches equals 1. Numbers within boxes refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 4 Heterochronic events in lepidosaur evolution.
Mapped onto molecular phylogeny, using continuous data, in relation to the ancestral squamate. Length of all branches equals 1. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 5 Heterochronic events in lepidosaur evolution.
Mapped onto molecular, stratigraphically calibrated phylogeny, using continuous data, in relation to the ancestral lepidosaur. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 7
Heterochronic events in lepidosaur evolution. Mapped onto morphological phylogeny, using continuous data, in relation to the ancestral lepidosaur. Length of all branches equals 1. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development. ( (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 9
Heterochronic events in lepidosaur evolution. Mapped onto morphological, stratigraphically calibrated phylogeny, using continuous data, in relation to the ancestral lepidosaur. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.  (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 11
Heterochronic events in lepidosaur evolution. Mapped onto molecular phylogeny, using unordered event-paired characters. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 12
Heterochronic events in lepidosaur evolution. Mapped onto molecular phylogeny, using ordered event-paired characters. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Figure 13
Heterochronic events in lepidosaur evolution. Mapped onto morphological phylogeny, using unordered event-paired characters. Numbers refer to developmental events (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.  (Table 2). Down arrow denotes earlier development of a given structure, while up arrow represents later development.

Skawiński and
analysis-out of 28 included species, 11 are iguanians and six are gekkotans, while there are only three anguimorphs (and all of them belong to a single clade, Varanus) and one scincoid. Members of other important clades, like Amphisbaenia and Dibamidae, were not included. Some of these groups only recently were studied in terms of development (e.g., Gregorovicova et al., 2012). Moreover, development of lepidosaurs included in this analysis is incompletely known. Thorough study of developmental sequences of these and other members of these diverse clades will be beneficial to future analyses. However, it may be that homoplasies are very common in developmental sequences of squamates. Moreover, the phylogenetic signal in organogenetic events (at least those used in this study) may be weak or detectable only in deeper nodes of the phylogenetic tree (cf. Jeffery et al., 2002b;Maisano, 2002). This may be indicated by higher congruence between methods in reconstructing heterochronic events than between given phylogenies. The only cladistic analyses that slightly resembled published phylogenies employed event-paired characters, especially ordered ones (Fig. 1B). In this analysis, eight of eleven included iguanian species formed a monophyletic group with Thamnophis and Mabuya that was sister to all other squamates. This resembles the morphological topology, where iguanians are sister group to all other squamates (e.g., Estes, de Queiroz & Gauthier, 1988;Conrad, 2008;Gauthier et al., 2012). This may suggest that developmental sequences of most iguanians and the tuatara are relatively similar. Under morphological topology, these similarities would represent symplesiomorphies but under molecular one, would be considered homoplasies. Reeder et al. (2015) suggested that support for basal placement of Iguania comes from the cranial characters. This is not the case in the present analysis. Character mapping and ancestral states reconstructions of event-paired data suggest that potential symplesiomorphies between the tuatara and iguanians (as a whole or one of their major subgroups-Acrodonta and Pleurodonta) are connected with the relatively later torsion completion, rather than of some events concerning head development. Other groups recognized by morphological analyses also receive some support. For example, Scleroglossa are supported by earlier occurrence of torsion completion (simultaneous with occurrence of hyomandibular slit and allantois bud), unlike in tuatara and Iguania. Scincomorpha are supported by simultaneous development of otic placode, allantois bud and secondary optic vesicle.
Gekkotans differ from other squamates in later development of the allantois (Andrews, Brandley & Greene, 2013) but in that trait they resemble the tuatara. Under molecular topology, earlier development of the allantois bud supports the Unidentata (Table 4). This may represent a genuine signal of monophyly of that group, however, caution is warranted. Gekkotans display many paedomorphic features, including their morphology (e.g., Daza, Bauer & Snively, 2014) and development (Jonasson, Russell & Vickaryous, 2012). Thus, the condition in gekkotans may represent reversal to the primitive condition (presumably, as displayed by the tuatara) rather than plesiomorphy. This situation is similar to the development of a single egg tooth, which purportedly supports the monophyly of Unidentata (see discussion in Assis & Rieppel, 2011). To gain more insight into that matter, it would be crucial to sample development of dibamids, the only other non-unidentate squamates.
In the fossil time-calibrated continuous analysis, only one event in two species is inferred to show heterochrony in relation to the ancestral lepidosaur. This may seem surprising, as some squamates show heterochrony to the ancestral squamate (much closer phylogenetically). However, if all studied taxa are extant (as is the case in the present analysis), the long branches would result in wider confidence intervals and thus ancestral state reconstructions for deep nodes of the phylogenetic tree would be less certain (Germain & Laurin, 2009). Integration of data from fossils would be useful in that regard but it seems highly unlikely that information on organogenesis can be preserved in the fossil record, despite recent significant advances in developmental palaeobiology (e.g., Skawiński &Tałanda, 2015).
In the continuous analyses (both calibrated and uncalibrated and using either molecular or morphological topology), values of all developmental events of the tuatara are located within the confidence interval of the ancestral squamate. This suggests that present data are equally consistent with either hypothesis of squamate phylogeny (cf. Germain & Laurin, 2009).
In this study only two major phylogenetic hypotheses of squamates were used. It is not beyond imagination that neither of these phylogenies is fully correct. For example, in the analysis combining morphological and molecular data conducted by Lee (2005) the ''fossorial group'' is polyphyletic, as suggested by molecular analyses (e.g., Wiens et al., 2012;Pyron, Burbrink & Wiens, 2013), but division of squamates into Iguania and Scleroglossa is retained, as in morphological analyses (e.g., Conrad, 2008;Gauthier et al., 2012). This could, to some extent, explain the discrepancies in reconstructions of heterochronic events, as none of these would be done on the basis of the correct tree.

CONCLUSIONS
Cladistic analyses conducted using characters generated by event-pairing and continuous analysis do not resemble any previously published phylogeny. Ancestral state reconstructions are equally consistent with both morphological and molecular hypotheses of squamate phylogeny. Results of the cladistic analyses, and the fact that reconstructions of heterochronic events show more similarities between certain methods than phylogenetic hypotheses, suggest that phylogenetic signal is at best weak in the studied developmental events.