Phylogenetic relationships of Neogene hamsters (Mammalia, Rodentia, Cricetinae) revealed under Bayesian inference and maximum parsimony

Taxon set

Included taxa depend on data availability and general completeness of the material. Within the seven extinct genera, Apocricetus, Collimys, Cricetulodon, Hattomys, Neocricetodon (including species assigned to Kowalskia), Pseudocricetus and Rotundomys, a total number of 41 species could be coded, which makes up around 77% of the total number of species within these genera (53), that were found in the literature. Additionally, two extant taxa were included as well, Cricetus cricetus (Linnaeus, 1758) and Nothocricetulus migratorius (Pallas, 1773). As outgroup, Eucricetodon wangae Li, Meng & Wang, 2016 was added, as coded in López-Antoñanzas & Peláez-Campomanes (2022). For additional information about the included taxa, e.g., age interval, references, observed material, see Supplemental Material S1.

Morphological characters and matrix construction

The matrix was constructed in Mesquite v. 3.81 (Maddison & Maddison, 2023). It is based on the morphological matrix from López-Antoñanzas & Peláez-Campomanes (2022), and expanded here from 82 characters to 116 characters, introducing additional characters corresponding to the structures allowing to differentiate cricetine genera and species. Four characters concern the whole molar row, six refer to morphometrics, and the remaining 106 are related to the morphology of each dental element (M1: 37; M2: 11; M3: 16; m1: 23; m2: 9; m3: 10). In cases of intraspecific variability between different locations, only the condition found in the type locality has been considered. In case of variability in the type location, only the character state present in most of the specimens was taken into account, except for specimens, for which no clear majorities were seen. The morphological matrix is provided in Supplemental Material S2.

Phylogenetic reconstructions

All final trees were annotated and visualized in R with the packages treeio, ggtree and deeptime (Wang et al., 2020; Yu, 2022; Gearty, 2023). For input and output files of the Bayesian inference and maximum parsimony analyses, see Supplemental Material S3.

Maximum parsimony analyses. All maximum parsimony analysis were run with TNT v. 1.6 (Goloboff & Morales, 2023) with all characters treated as unordered.

Equal weights analysis. The analysis under equal weights (MP-EW) was conducted with new technology algorithms using initial trees from 1,000 rounds of random addition sequence, with 100 iterations or rounds for sectorial search, ratchet, and tree fusing. The resulting 106 most parsimonious trees (256 steps, consistency index (CI): 0.414, retention index (RI): 0.715) were used to calculate a 50% majority consensus tree (258 steps, CI: 0.411, RI: 0.711). Clade support (given in %) was calculated based on 1,000 bootstrap (BS) replicates under the same parameters (Felsenstein, 1985).

Implied weights analysis. An additional analysis was run under the same options, as before but including implied weighting (MP-IW) (Goloboff et al., 2008). Following recent suggestions, a larger concavity index (k) of 12 was used (see Goloboff, Torres & Arias, 2018). The resulting two most parsimonious trees (257 steps, CI: 0.412, RI: 0.713) were used to calculate a consensus tree, (258 steps, CI: 0.411, RI: 0.711) with clade support based on 1,000 BS replicates (Felsenstein, 1985).

Bayesian inference analyses. Bayesian analyses were run with the parallel version of MrBayes v. 3.2.7a (Altekar et al., 2004; Ronquist et al., 2012b) using the Cyber-Infrastructure for Phylogenetic Research (CIPRES) Science Gateway version 3.3 (Miller, Pfeiffer & Schwartz, 2010).

Non-clock analysis. The Mkv model (Lewis, 2001) was used with among character rate heterogeneity modelled under a gamma distribution (Yang, 1993). All characters were treated as unordered. The analyses were run with four independent Metropolis-Coupled Markov chain Monte Carlo (MCMCMC) runs with six chains and 30,000,000 generations, sampling every 1,000 steps and a burn-in of 30%. Convergence and sufficient length of the runs were checked, using the R package Convenience v. 1.0.0 (Fabreti & Höhna, 2022). Based on the posterior tree sample a maximum clade compatibility (MCC) tree was calculated, as a consensus tree.

Time-calibrated relaxed-clock analysis. All settings of the non-clock analysis were adopted, except the number of MCMCMC generations, which was increased to 50,000,000. Time-calibrated relaxed-clock analyses were performed under a fossilized birth-death (FBD) tree prior (Stadler, 2010; Zhang et al., 2016). To model the way in which extant and extinct taxa are sampled in the construction of the tree, different strategies can be used (Simões, Caldwell & Pierce, 2020). To avoid problems when inferring speciation or extinction rates (Höhna et al., 2011), we have tested two of the three strategies, that are compatible with the FBD tree prior. The option ‘diversity’, that assumes a sampling strategy to maximize the diversity of extant taxa, was excluded, as our database only includes two extant species. Consequently, we have tested the two models that assume randomly sampled extant species. The first one, with sampled ancestors, SA-FBD (‘random’), allows the fossil taxa to be tips or ancestors of other taxa, while in the second one, the so-called noSA-FBD (‘fossiltip’), the fossil taxa have to be tips. The use of one or another can have an impact in the estimations of divergence times (Gavryushkina et al., 2014; Simões, Caldwell & Pierce, 2020). For the extant sampling probability, the number of included extant taxa (2) is divided by the total number of extant Cricetinae species (18 after Musser & Carleton, 2005).

In order to time-calibrate the tree, the tip-dating approach was used (Ronquist et al., 2012a; Ronquist, Lartillot & Phillips, 2016). Age ranges of the fossil taxa, resulting from age uncertainties of one or multiple locations of one taxon, were addressed by assigning uniform prior distributions to the tip calibrations, which can help to avoid erroneous divergence time estimations (O’Reilly, Dos Reis & Donoghue, 2015; Barido-Sottani et al., 2019). For the root age, an offset exponential distribution was set as a prior, with a minimum of 33 Ma (= minimal age of the oldest included fossil Eucricetodon wangae) and a mean of 41.2 Ma (following López-Antoñanzas & Peláez-Campomanes, 2022).

To give an informative prior to the base rate of the clock, the median tree length, calculated by a preceding non-clock analysis, was divided by the median of the root age prior (3.189768/37.1 = 0.085978) (following Simões et al., 2018, 2020). This estimated rate in natural log scale (= −2.45367) was used as the mean of a log-normal distribution with the exponent of the mean (e^0.085978 = 1.08978) as the standard deviation (following Pyron, 2017). To enforce proper rooting of the tree and facilitate reaching convergence, the ingroup was constrained to be monophyletic.

For relaxing the clock, two different models, compatible with the FBD prior are implemented in MrBayes v. 3.2.7a. The IGR (Independent gamma rate) model draws substitution rates from a gamma distribution, uncorrelated between branches, which allows more dramatic rate changes (more punctuated mode of evolution) (Drummond et al., 2006). The second model, TK02, samples from a lognormal distribution and is autocorrelated between branches, which represents a more gradual mode of evolution (Thorne & Kishino, 2002). Both models were used here, resulting in a total of four different models with all combinations of ‘fossiltip’ vs. ‘random’ and IGR vs. TK02. To choose the best fit model, stepping-stone sampling was done to estimate the marginal likelihoods (Xie et al., 2011). These can be used to calculate Bayes factors to compare the fit of two models to the data. For the stepping-stone sampling, the number of MCMCMC generations was increased by a factor of 10 to 500,000,000 as suggested by Ronquist et al. (2020).

All analyses were checked for convergence by the R package Convenience v. 1.0.0 (Fabreti & Höhna, 2022), as mentioned for the non-clock analysis.

Rogue taxon identification and tree set pruning. To improve posterior probabilities of the resulting trees, so-called rogue (‘wildcard’) taxa, were identified. These taxa are characterised by an unstable position in the tree, as they are resolved in varying clades in the trees of a tree set, e.g., the posterior tree sample of a Bayesian inference analysis. This leads to decreased posterior probabilities or even less resolved consensus trees. Equally, in the case of maximum parsimony analyses, they can affect the consensus tree calculated by several most parsimonious trees or the support values, given by a set of bootstrap trees. Deletion of the rogue taxa from the tree sets before calculating the consensus trees, can therefore lead to better resolved and supported trees, while they are still incorporated in the actual reconstruction. The deletion of the rogue taxa from the taxon set followed by a re-run of the analysis (e.g., in Aberer & Stamatakis, 2011; Simões, Caldwell & Pierce, 2020) is seen critically by some authors as it means disregarding available and potentially important information (as discussed in Goloboff & Szumik, 2015).

In this study, the posterior tree samples of both clock trees were used to identify rogue taxa utilising the R package Rogue v. 2.1.6 (Smith, 2022, 2023). An additional examination of the 106 most parsimonious trees of the equal weighting maximum parsimony analysis, using the web interface of RogueNaRok (Aberer, Krompass & Stamatakis, 2013), did not result in any identified rogue taxa. For the R-code used to identify the rogue taxa, see Supplemental Material S4.

This study does not include the Pleistocene cricetine taxa, for which a revision beyond the scope of this paper is needed. The lack of the youngest cricetine fossil taxa from our analysis, makes the inferred position of the two included extant species uncertain. For this reason, they were removed from the resultant trees (Figs. 1, 2, S5.1–3) but are shown in Figs. S5.4–8 together with a discussion on the identified rogue taxa in Supplemental Material S5.

In total, the tree sets of all analyses were pruned by the same rogue taxa and the extant taxa, to obtain comparable trees. These trees, that are based on the pruned taxon set were then used for the subsequent analyses.

Assessing stratigraphic congruence

To compare the fit of the topologies resulting from the different reconstruction methods, the R package strap v. 1.6 (Bell & Lloyd, 2015) was used to first time calibrate the non-clock trees, and then to calculate the following stratigraphic fit indices: (i) the relative completeness index (RCI) assesses the amount of gaps in the fossil record in relation to the observed fossil ranges in the tree (Benton & Storrs, 1994); (ii) the gap excess ratio (GER) indicates the sum of ghost ranges in the tree scaled in relation to the possible minimum and maximum sum of ghost ranges in theoretical topologies (Wills, 1999); (iii) the modified Manhattan stratigraphic measure (MSM*) indicates the sum of ghost ranges in the tree in relation to the sum of ghost ranges of the theoretical tree of best fit to the stratigraphy (Siddall, 1998; Pol & Norell, 2001). For all indices, significance tests were carried out, that resulted in very small p-values. Resulting indices of all analyses are listed in Supplemental Material S6.

Character mapping

To ensure reliable results, two methods and two trees were used to map ancestral character states. Following the results of the stratigraphic congruence analyses, the implied weighting maximum parsimony tree was used to map characters in TNT under a parsimonious approach (see Supplemental Material S7) and the time-calibrated tree under the IGR clock model was used for stochastic character mapping with the R package phytools v. 2.1.1 (Revell, 2024). For the stochastic character mapping, three different models were fit to each character, with rates between states being either equal, symmetrical, or all different. The respective best fitting model was chosen using the Akaike information criterion and then used to map the character on the tree. Finally, for each character, the posterior probabilities of the different states were plotted on the nodes of the tree. For the reconstructed maps of characters, that are mentioned in the results, see Supplemental Material S8.

Maximum parsimony analysis

Maximum parsimony analyses of morphological data sets including highly homoplastic characters can be improved by weighting characters according to their homoplasy (Goloboff et al., 2008). Consequently, implied weighting parsimony analyses can produce more resolved and accurate trees than standard equal weights (Smith, 2019). To assess which tree, from equal or implied weighting, fits better with the chronostratigraphy, stratigraphic congruence indices were calculated. While the topologies of both parsimony trees are quite similar (see Figs. 1 and S5.3), our results evidence a better stratigraphic congruence of the implied weighting tree than of the equal weighting tree (see Supplemental Material S5). Therefore, we discuss below the topology retrieved by applying implied weighting (Fig. 2).

The topology of the consensus tree of the two most parsimonious trees shows a major basal split into two main clades. The first one (stemming from node 28) consists of Collimys, Rotundomys, Cricetulodon bugesiensis and Cricetulodon sabadellensis whereas the other one (stemming from node 3) includes the remaining cricetines. The most basal taxon of this latter clade is Cricetulodon hartenbergeri, followed by Pseudocricetus (node 27), which is in turn sister to a large clade (node 5) that splits into two lineages. The first of them (stemming from node 6) includes Cricetulodon complicidens and most of the species belonging to the genus Neocricetodon. The second one (stemming from node 17) has as most basal taxa the sister species N. browni and N. nestori, as sister clade to the group stemming from node 18: Apocricetus plinii, followed by Cricetulodon meini plus Cricetulodon lucentensis, and the sister clades Hattomys and Apocricetus sensu stricto (s.s.).

Tip-dated Bayesian analysis

The results of the stepping stone sampling strongly evidence a better fit for the ‘fossiltip’ sampling strategy under both clock models (2 log_e(B10) > 10, see Kass & Raftery, 1995). Comparing Tk02 with IGR under ‘fossiltip’, the better fit model is the Tk02 (2 log_e(B₁₀) = 4.6) but without strong evidence. However, the analysis under the Tk02 model showed some problems in reaching convergence (ESS below 200 for two parameters, ‘convergence failed’ according to Convenience). The IGR model showed a considerably better performance (lowest ESS > 3,300, ‘convergence reached’ according to Convenience). Moreover, the results of the strap analyses show a better stratigraphic congruence for the IGR model than for the TK02 model (see Supplemental Material S6). As both calibrated Bayesian inference analyses resulted practically in the same topology of the trees (see Figs. 2 and S5.2) with only slight differences in the posterior probabilities (generally a bit higher in Tk02) and in the relationships among the species belonging to Neocricetodon (particularly in that concerning the clades with low posterior probabilities), we describe below only the results obtained by applying the ‘fossiltip’ strategy under the IGR model (Fig. 3).

The topology of our tree evidences four major clades. The most basal one (stemming from node 35) consists of all species belonging to Collimys. Its sister clade (stemming from node 3) includes all remaining cricetines, with its most basal clade (stemming from node 31) including all Rotundomys species on the one side and the remaining two major clades (stemming from nodes 5 and 18) on the other. The first of the latter (stemming from node 5) consists of all species of Neocricetodon, as well as Cricetulodon complicidens. The second one (stemming from node 18) has Cricetulodon hartenbergeri as sister of two clades: a small one (node 30) including Cricetulodon bugesiensis and Cricetulodon sabadellensis and a larger one (stemming from node 19) constituted by a succession of clades with Pseudocricetus (node 29) at the base, followed by the sister taxa Cricetulodon meini and Cricetulodon lucentensis (node 28). One node up (node 22) inserts A. plinii, followed by the sister lineages Apocricetus s.s. and Hattomys. The topology of the tree supports the monophyly of Collimys (node 35), Rotundomys (node 31), Neocricetodon (node 5), Pseudocricetus (node 29), Apocricetus s.s. (node 26), and Hattomys (node 24). In contrast, the genus Cricetulodon is paraphyletic and only the type species Cricetulodon sabadellensis and Cricetulodon bugesiensis should be included in this genus.

Divergence times. The divergence times estimated using the two different clock models IGR and TK02, vary only slightly, with differences of mostly less than 500,000 years. In the same way, uncertainties on divergence times (measured by 95% highest probability densities (HPD) ranges) are quite similar independently of the clock model applied. In the following, the results of the analysis under the IGR model are reported. They reveal a late Early Miocene age for the first split within the ingroup (16.54 Ma, 95% HPD: 14.53–19.12 Ma). Collimys is recovered as the oldest genus, diverging during the Middle Miocene (14.81 Ma, 95% HPD: 14.07–15.85 Ma). The remaining five monophyletic genera diverged later, during the Late Miocene: Neocricetodon (11.28 Ma, 95% HPD: 10.26–12.64 Ma), Rotundomys (10.91 Ma, 95% HPD: 10.01–12.23 Ma), Pseudocricetus (7.9 Ma, 95% HPD: 6.69–9.29 Ma), Hattomys (6.7 Ma, 95% HPD: 5.17–8.47 Ma), and Apocricetus s.s. (5.7 Ma, 95% HPD: 4.59–6.81 Ma).

Ancestral character state reconstructions

The slight differences in the topology of both trees (particularly regarding the phylogenetic position of Rotundomys and Collimys) result in different reconstructed synapomorphies for some clades. However, concerning the genera, most of the synapomorphies found when analysing the results of the stochastic character mapping are also retrieved by the parsimonious mapping of synapomorphies (see Supplemental Materials S7 and S8). Only a few synapomorphies are not found with the latter approach, due to the coding of polymorphic states as ‘ambiguous’ in the mapping process in TNT. Therefore, below we list mainly the synapomorphies obtained by applying stochastic character mapping.

All species belonging to Collimys (stemming from node 35) share the exclusive synapomorphy of having an ectomesolophid on the m1 (68: 0→1). Additional ambiguous synapomorphies are: e.g., a labial spur of the anterolophule reaching the labial border of the M1 (103: 0→2), a long mesoloph on the M3 (49: 2→0) and the presence of a labial spur of the anterolophulid on the m1 (64: 0→1).

The clade Rotundomys (stemming from node 31) is defined by the following exclusive synapomorphies: a very weak to absent mesoloph on M1 and M2 (20: 1→2, 37: 1→2), a posteroloph that is merged with the posterior metalophule on the M1 (26: 1→3) and the presence of a short or hanging labial anterolophid on the m3 (88: 0→1), as well as many non-exclusive synapomorphies.

Representatives of the clade stemming from node 5 including Neocricetodon and Cricetulodon complicidens share non-exclusive synapomorphies such as the presence of the anterior metalophule on the M2 (38: 2→1/0). Additional non-exclusive synapomorphies, only shared with Collimys, are the presence of a long mesolophid on the m1 (66: 0/1→2) and having a medium to long labial spur of the anterolophule on the M1 that reaches the molar border in most of the taxa belonging to Neocricetodon (103: 0→1/2).

Representatives of Pseudocricetus (stemming from node 29) are characterised by having a short but distinct posteroloph on the M3 (56: 0→1) and the tendency to form a very small mesolophid on the m2 (73: 2→1).

Apocricetus s.s. (stemming from node 26) is characterised by the two non-exclusive synapomorphies of having a multi-lobed, crestiform anteroconid (57: 2→4) with a poorly developed labial anterolophid (60: 0→1) on the m1.

Representatives of Hattomys (stemming from node 24) are clearly distinct from those belonging to its sister group Apocricetus s.s. They are characterised by sharing the two exclusive synapomorphies of having large M1, longer than 3.2 mm (1: 1→3) and an hypolophulid connected to a medium sized mesolophid on the m1 (109: 2→5).

Maximum parsimony vs. Bayesian trees

All maximum Parsimony and Bayesian tip-dated and undated trees (Figs. 1, 2, S5.1–3) support the monophyly of Collimys, Rotundomys, Pseudocricetus Apocricetus s.s. and Hattomys. Neocricetodon includes Cricetulodon complicidens in all trees. It is monophyletic in the Bayesian trees but not in the parsimonious ones, as two species branch outside of the clade (N. browni and N. nestori in the implied weighting analysis and N. moldavicus and N. fahlbuschi in the equal weighting one).

Cricetulodon splits up in the same groups in both, the IGR Bayesian clock analysis and the implied weighting parsimony analysis (Figs. 1 and 2). However, some differences are observed concerning the phylogenetic position of the clades. In this way, the topology of the parsimony tree shows Cricetulodon bugesiensis and Cricetulodon sabadellensis as the most basal taxa of the clade including Collimys and Rotundomys (stemming from node 28, Fig. 2), whereas in the Bayesian topology this lineage belongs to the clade that comprises the remaining species of Cricetulodon, Pseudocricetus, Apocricetus and Hattomys (stemming from node 18, Fig. 3). Moreover, A. plinii is basal to Cricetulodon meini and Cricetulodon lucentensis in the parsimony tree, but sister to Apocricetus s.s. and Hattomys in the Bayesian one (inserting at node 18, Fig. 2 vs. node 22, Fig. 3).

The most striking difference when comparing the Bayesian clock tree to the parsimony one concerns the relationship between Collimys and Rotundomys, for which the parsimony analysis found a sister relationship (see node 29, Fig. 2, BS = 76) whereas in the Bayesian tree they insert sequentially (at node 2 and 3, Fig. 3). Therefore, depending on the topology, contrasting ancestral character state reconstructions are found. The characters that are different in Collimys and Rotundomys compared to the remaining cricetines are the following: M1 anterocone not clearly divided, but more crestiform vs. divided in two (5: 6 vs. 2); M1 lingual anteroloph clearly present vs. weak or absent (11: 0 vs. 1); upper molars protolophule posterior vs. double (7: 0 vs. 1, 34: 2 vs. 1, 47: 3 vs. 0); m3 lingual anterolophid well-developed vs. weak or absent (76: 1 vs. 0, reversed in Apocricetus s.s. and Hattomys (76: 0 vs. 1)). According to the stochastic character mapping, all these characters serve as synapomorphies for the clade Neocricetodon + Cricetulodon + Pseudocricetus + Apocricetus + Hattomys (stemming from node 4, Fig. 3). In the parsimony analysis, on the other side, the results are reversed. Only having a double protolophule on the M1 is reconstructed as a synapomorphy for this clade (stemming from node 3, Fig. 2), while the other above-mentioned characters serve here as synapomorphies for the clade Collimys + Rotundomys (node 29, Fig. 2).

The genera Rotundomys and Collimys were informally grouped by Kälin (1999) on the basis of emerging hypsodonty, whereas Heissig (1995) excluded Collimys as potential ancestor of Rotundomys, stating that the tendency of acquiring hypsodonty evolved independently.

In the parsimony analysis, Cricetulodon sabadellensis and Cricetulodon bugesiensis are additionally positioned as sister taxa to Collimys and Rotundomys. Considering the much younger age of the Cricetulodon taxa compared to Collimys, the arrangement in the Bayesian tree seems more plausible.