How important is it to consider lineage diversification heterogeneity in macroevolutionary studies? Lessons from the lizard family Liolaemidae

Macroevolutionary and biogeographical studies commonly apply multiple models to test state‐dependent diversification. These models track the association between states of interest along a phylogeny, although many of them do not consider whether different clades might be evolving under different evolutionary drivers. Yet, they are still commonly applied to empirical studies without careful consideration of possible lineage diversification heterogeneity along the phylogenetic tree. A recent biogeographic study has suggested that orogenic uplift of the southern Andes has acted as a species pump, driving diversification of the lizard family Liolaemidae (307 described species), native to temperate southern South America. Here, we argue against the Andean uplift as main driver of evolution in this group. We show that there is a clear pattern of heterogeneous diversification in the Liolaemidae, which biases state‐ and environment‐dependent analyses in, respectively, the GeoSSE and RPANDA programs. We show here that there are two shifts to accelerated speciation rates involving two clades that have both been classified as having “Andean” distributions. We incorporated the Geographic Hidden‐State Speciation and Extinction model (GeoHiSSE) to accommodate unrelated diversification shifts, and also re‐analyzed the data in RPANDA program after splitting biologically distinct clades for separate analyses, as well as including a more appropriate set of models. We demonstrate that the “Andean uplift” hypothesis is not supported when the heterogeneous diversification histories among these lizards is considered. We use the Liolaemidae as an ideal system to demonstrate potential risks of ignoring clade‐specific differences in diversification patterns in macroevolutionary studies. We also implemented simulations to show that, in agreement with previous findings, the HiSSE approach can effectively and substantially reduce the level of distribution‐dependent models receiving the highest AIC weights in such scenarios. However, we still find a relatively high rate (15%) of distribution‐dependent models receiving the highest AIC weights, and provide recommendations related to the set of models included in the analyses that reduce these rates by half. Finally, we demonstrate that trees including clades following different dependent‐drivers affect RPANDA analyses by producing different outcomes, ranging from partially correct models to completely misleading results. We provide recommendations for the implementation of both programs.

| 1287 PERSPECTIVE macroevolutionary rates (Rabosky & Glor, 2010). These models track associations between the states of interest and speciation and extinction rates along a phylogenetic tree, and they have been implemented in hundreds of empirical studies (Rabosky & Goldberg, 2015). In addition, interest has also increased in the application of correlative algorithms (Condamine et al., 2013;Condamine, Sperling, Wahlberg, Rasplus, & Kergoat, 2012;Morlon et al., 2016;Steeman et al., 2009;Winkler et al., 2010), as implemented in the software RPANDA (Morlon et al., 2016), to use environment-dependent models to test whether gradual changes in palaeoenvironments have significantly influenced speciation and extinction rates. These state-and environment-dependent models have become popular, but important concerns have been raised (at least) for the SSE family of models that do not consider whether unrelated traits are associated with shifts in diversification rates (Beaulieu & O'Meara, 2016;Maddison & FitzJohn, 2014;Rabosky & Goldberg, 2015. For example, strong correlations with diversification are sometimes inferred from rate shifts, leading to dramatically high rates of false positives even if the shift is unrelated to the targeted state (FitzJohn, 2010;Maddison & FitzJohn, 2014;Maddison, Midford, & Otto, 2007;Rabosky & Goldberg, 2015. These issues may also be affecting environment-dependent correlative models in such scenarios . Consequently, while larger trees are preferred due to the presumed increase in power, this also increases the risk of including clades that differ in factors that can affect diversification rates along a tree (ecological requirements, dispersal abilities, and life history traits [Li, Huang, Sukumaran, & Knowles, 2018]). Hence, at least for the case of SSE models, new algorithms that include "hidden states" (HiSSE) have been proposed for binary traits (BiHiSSE; Beaulieu & O'Meara, 2016), and more recently geographic-dependent diversification hypotheses (GeoHiSSE; Caetano, O'Meara, & Beaulieu, 2018) and multiple states (SecSSE; Herrera-Alsina, Els, & Etienne, 2019). Hidden states refer to unsampled traits that are related to diversification shifts in a phylogenetic tree, thus incorporating heterogeneous diversification into the original SSE.
The HiSSE has improved the SSE model family and potentially resolved the problem of unrelated shifts in the phylogeny when appropriate sets of models are tested (Beaulieu & O'Meara, 2016;Caetano et al., 2018;Herrera-Alsina et al., 2019). On the other hand, there are no currently available solutions for environment-dependent models (i.e. hidden states have not been implemented so far), besides splitting clades into separate analyses (Lewitus, Bittner, Malviya, Bowler, & Morlon, 2018). Even though models that do not include hidden states are not suitable for trees including groups with heterogeneous diversification patterns, they are still commonly applied in macroevolutionary studies without careful assessment. For instance, the Geographic State Speciation and Extinction model (GeoSSE; Goldberg et al., 2011) was recently used in biogeographic studies of lizards (Esquerré, Brennan, Catullo, Torres-Pérez, & Keogh, 2019), tanagers and tortoises (Román-Palacios & Wiens, 2018), plants (Canal et al., 2019), including palms (Bacon, Velásquez-Puentes, Hoorn, & Antonelli, 2018) and oaks (Hipp et al., 2018), and a recent mega-phylogeny of mushrooms (Varga et al., 2019). Given the above-noted limitations of this model, these empirical cases should be revisited (e.g. Harrington & Reeder, 2017). Here, we re-analyse the recently published phylogeny of the lizard family Liolaemidae (Esquerré et al., 2019) (Cei, 1986

| MATERIAL S AND ME THODS
Here, we briefly describe our methods, but for further details, see the Extended Materials and Methods in the Data S1.

| Identifying possible diversification shifts along the tree and exploration of the relationship between speciation rates versus altitude
We explored the association between speciation (and extinction) rates with the maximum altitude of occurrence known for any species (phylogenetic tree and altitude data taken from Esquerré et al., 2019). The phylogenetic tree includes the mono-

| Hypothesis testing: role of the
Andean orogeny in diversification of the Liolaemidae

| Geographic Hidden-State Speciation and Extinction (GeoHiSSE)
We then implemented the GeoHiSSE model (Caetano et al., 2018) using the hisse package (Beaulieu & O'Meara, 2016). Given that we do not fully agree with the original classification of Andean and non-Andean species by Esquerré et al. (2019), here we propose an alternative classification, and performed these analyses using both (Table S8). As an example, the "Patagonia" group is distributed across a huge area that was assumed to be "Andean" in Esquerré area-independent models (CID) indicating no relationship between diversification rates and geographic areas, and 10 area-dependent diversification models (all models listed in Table S7). Given the relatively high rates distribution-dependent models receiving the highest AIC weights found in our simulations (see results), and issues of overparametrization (Caetano et al., 2018)

| Correlative environment-dependent models in RPANDA
Esquerré et al. tested whether speciation and extinction rates are correlated with the Andean uplift using the R package RPANDA (Morlon et al., 2016). Similar to the description above, they implemented this model on the complete phylogeny of Liolaemidae and did not analyse the two different genera into separate analyses.
They fitted a total of 10 models, including eight for Andean uplift dependency versus two null models where rates are constant through time. We believe that this set of models is lacking in alternative scenarios, and that Andean uplift dependency could have been selected simply because null constant rate models are unrealistic.
Specifically, the study fails to include time-dependent models (see Morlon, Parsons, & Plotkin, 2011

| RE SULTS AND D ISCUSS I ON
See also extended results in Data S1.  Figure S4). Specifically, the genus Phymaturus has the highest speciation rate that is also associated with a high extinction rate.

| Disparate patterns of diversification in
This result is concordant with another study using a different phylogenetic tree (Olave et al., in press), and also supported by similar results recovered using ClaDS ( Figure S2). We note that both programmes find two shifts to accelerated speciation rates involving clades that were classified entirely as "Andean" by Esquerré et al.

| Exploration of the relationship between speciation and extinction rates versus maximum altitude
We constructed linear models between the maximum altitude (MA) of species occurrence records, and the species-specific speciation and extinction rates. Linear models reveal non-significant correlations of the speciation/extinction rates with the MA for all target clades (Table S5). Scatter plots show no relationship between rates and MA (Figure 2; Figures S5-S11).
Other CID models showed lower support (M3, M19, M9, M20, M5 and M11; Figure 3a; Table S7). On the other hand, results based on the original geographic classification by Esquerré et al. (Figure 3b;  Figures 1, 3 and S2) and that, when F I G U R E 1 Color-coded phylogenetic trees for the speciation (a) and extinction (b) rates through time for the Liolaemidae inferred by BAMM. Our preferred geographic classification of species is shown on the left tree and the original geographic classification on the right tree; corresponding to Andean, non-Andean and widespread, shown at each tip in white, black and grey respectively. See Figure S1 for diversification rates [Colour figure can be viewed at wileyonlinelibrary.com] incorporating the HiSSE into the models, results show no apparent signal of Andean uplift increasing speciation rates in the Liolaemidae ( Figure 3; Table S7). This is true in both cases, when considering either our preferred geographic classification or the original classification ( Figure 3; Table S7). These results show that Esquerré et al.
have confounded clade-specific rate accelerations in their distribution-dependent diversification results.
In addition, the GeoHiSSE analyses also return estimates of ancestral distributions (white-black colour gradient in phylogenetic trees in Figure 3). Given the fact that we have changed the geographic classification, these results show that our preferred clas- palluma: widespread; and P. patagonicus: non-Andean; Figure 3a).
These results also argue against the Andes acting as "species pump".

| Correlative environment-dependent models in RPANDA applied to study the lizard family Liolaemidae
We fitted 42 models to each genus within Liolaemidae using RPANDA. This algorithm does not incorporate "hidden states", so the only proper way to approach such heterogeneous scenarios is by running separate analyses for each genus (but see extended results in Data S1 and Table S9 for analyses of Liolaemidae as a whole). The best model to describe evolutionary rates in Liolaemus is a birth-death model in which speciation and extinction rates are linearly correlated with time (i.e. our model 18; Table S9; Figure 4a). This model received an AIC weight average of 0.26 ± 0.17 (lambda = 0.098, alpha = 0.14; Table S9) among 100 trees sampled from the posterior, followed by a model in which speciation and extinction rates are linearly correlated with global temperature changes (AIC weight mean = 0.08 ± 0.07; model 15; see Table S9 TA B L E 2 Species counts for the geographic classification from Esquerré et al. (2019); taken from their Data S1. Note that some species were classified as both Andean and non-Andean as "widespread" (classification details in Table S8)

Considered "Andean species"
Considered "Non-Andean species" is remarkably similar across species, suggesting that it may be evolutionarily or ecologically constrained (Cruz et al., 2009;Ibargüengoytía et al., 2008). Thus, it is biologically plausible that once the temperature conditions became more favourable for this group, then Phymaturus' speciation rates increased (Figure 4b,c). Different abiotic factors have been suggested to affect South American climates, including Andean uplift as one of the possible promotors (e.g. Blisniuk et al., 2006;Gregory-Wodzicki, 2000). Is it possible that the Andes uplift has indirectly driven the increase in speciation rates in Phymaturus by promoting local climate change? There has been considerable debate over the role of the Andes relative to South American climate change, including whether the Andes promoted climate change, or global climate change promoted the Andean uplift (Lamb & Davis, 2003). However, post-Cenozoic temperature decrease was a global phenomenon (Zachos et al., 2008;, and a combination of relevant factors also acted as locally cooling promotors, such as the Humboldt Current generated by the closure of the Central America Seaway (Garreaud, Molina, & Farias, 2010;Hartley, 2003). We emphasize that multiple factors have likely driven Phymaturus evolution, and assuming a temperature-only dependent diversification is unnecessarily reductionist. We also cannot discard other important aspects in the peculiar biology of this genus. For example, all >40 cold-adapted species only occur in isolated patches of rock outcrops; therefore, migration between populations is likely severely limited relative to Liolaemus (Vicenzi, Corbalán, Miles, Sinervo, & Ibargüengoytía, 2017). The strong fidelity of Phymaturus species to specific microhabitats, "islands" of big boulders with deep crevices in volcanic cliffs, peaks and plateaus (Cei, 1986), in efficient selection for adaptation to a stable and narrow niche, reducing the cost of trade-offs by abandoning traits needed to utilize a wider range of resources (Futuyma & Moreno, 1988).
The contrasting biology between the sister groups Liolaemus and Phymaturus (Table 1)  3.2 | The importance of considering lineage diversification heterogeneity in macroevolutionary studies: warnings in implementation of GeoHiSSE and RPANDA programs

| Assessing the power of GeoHiSSE models: are hidden states improving state-dependent diversification models?
Previously, it has been shown that including hidden states to the model is a plausible solution to account for diversification heterogeneity and improve the accuracy of the original SSE models (Caetano et al., 2018). We assessed the power of GeoHiSSE analyses by simulating random permutations of geographic areas, and show that most of the CID models received the highest weights over distribution-dependent models (~85%). This is undoubtedly an improvement over the high level of false positives reported for the original SSE models (Rabosky & Goldberg, 2015. However, we still found a relatively high rate of distribution-dependent models receiving the highest weights (=15%) when considering all 35 models and permutations based on the original geographic classifications ( Figure S14a; Table S7), and 16% with simulations based on our preferred classification (Table S7). Interestingly, all state-dependent models receiving the highest weights are equal or greater than the maximum number of free parameters (=15 -21) in the most complex CID model (M11 and M17). Furthermore, most of these models are supported by high AIC weights ( Figure S14b). We show that excluding models with >15 free parameters (i.e. selection among 32 models) reduced these rates to 7% ( Figure S14c) following the original classification, and to a ~ 10% rate given the permutations on our preferred classification (Table S7; see also extended results in Data S1 for further details).
Given the recent documentation that a proper set of state-independent models for HiSSE methods should, necessarily, have the same number of free diversification parameters than the state-dependent models for a "fair" comparison (Caetano et al., 2018), we think that the original set of 35 models should be taken with caution. This issue was paradoxically shown in Caetano et al. (2018) for BiHiSSE implementation, but not corrected among the 35 proposed models for GeoHiSSE in the same paper. Our study does not mean that the same power estimations will necessarily apply to another empirical phylogeny, but we have demonstrated the importance of assessing power before applying HiSSE models, as well as reinforced the notion of potential risks related to the set of models (and number of free parameters) implemented for HiSSE hypothesis tests.

| Impact of heterogeneous trees in RPANDA: trees including clades following different dependentdrivers affect results
Inspired in our empirical results, we performed a simulation study to address how RPANDA analyses could be affected when two different clades in a phylogenetic tree evolve following different drivers. An extensive simulation study to test the performance of RPANDA was performed by , including exploration of different proportions of missing tips, accurate model selection and also rate shifts along the phylogeny.
Although the authors show that RPANDA seem to be robust when rate shifts are present along the phylogeny under their simulated conditions, here we focused on a different type of heterogeneous scenario.
Our results show that, even though RPANDA can accurately recover the true scenario when there is a single driver for speciation rates in a phylogenetic tree ( Figure S15; Table S10)  ing clades with important biological differences is key to assess potential differences a priori. As we have shown here, the two contrasting groups of Liolaemus and Phymaturus not only differ in speciation rates but also in several aspects of their biology (Table 1). Such biological differences are expected to affect the macroevolutionary patterns (Li et al., 2018). Empirical studies must explore such possibilities when using models that do not allow heterogeneous trees, as is it the case of the current version of RPANDA. showing that: (a) there is no support for increased speciation rate in Andean species (Figure 3 and Figure S7), (b) a time-dependent model better explains diversification of the genus Liolaemus than Andean uplift (Figure 4; Table S9) and (c) a model considering the interaction between time and global temperature changes is a better fit for the genus Phymaturus than Andean uplift (Figure 4; Table S9). As a final remark, our preferred classification of species distributions has now changed inferences in ancestral-range reconstructions, and suggests that ancestors of both genera and all four subgenera were non-Andean or widespread (i.e. both Andean and non-Andean; Figure 3a), thus, challenging the Andean "spe- have demonstrated that trees including clades following different dependent drivers affect RPANDA results ( Figure 5), and the only way to circumvent such problem is to separate clades for independent analyses, given that this method does not account for hidden states. We also provide recommendations for the implementation of GeoHiSSE models to account for fair comparisons among models included in the analysis. We reinforce the previously known case of removing the distribution-dependent overparameterized models in GeoHiSSE ( Figure S14; Table S7), to prevent exceeding the number of free parameters in the most complex null model.

| CON CLUS IONS
We note that, although here we focused on the recently published Liolaemidae phylogeny as a model system, it is very likely that many earlier empirical studies may have limitations similar to those we report here, and these should be re-assessed.

ACK N OWLED G M ENTS
We thank Richard Ree and two anonymous reviewers for useful comments made on earlier versions of this article. We also thank D.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data to support the findings of this study will be openly available