Getting Its Feet on the Ground: Elucidating Paralouatta’s Semi-Terrestriality Using the Virtual Morpho-Functional Toolbox

Püschel, Thomas A.; Marcé-Nogué, Jordi; Gladman, Justin; Patel, Biren A.; Almécija, Sergio; Sellers, William I.

doi:10.3389/feart.2020.00079

ORIGINAL RESEARCH article

Front. Earth Sci., 25 March 2020

Sec. Paleontology

Volume 8 - 2020 | https://doi.org/10.3389/feart.2020.00079

This article is part of the Research Topic Evolving Virtual and Computational Paleontology View all 16 articles

Getting Its Feet on the Ground: Elucidating Paralouatta’s Semi-Terrestriality Using the Virtual Morpho-Functional Toolbox

$\r\nThomas A. Püschel*$ Thomas A. Püschel^1*

Jordi Marcé-Nogué^2,3,4

Justin Gladman⁵

Biren A. Patel^6,7

Sergio Almécija^8,9,4

William I. Sellers¹⁰

¹Primate Models for Behavioural Evolution Lab, Institute of Cognitive and Evolutionary Anthropology, School of Anthropology, University of Oxford, Oxford, United Kingdom
²Department of Pathology and Anatomical Sciences, Jacobs School of Medicine and Biomedical Sciences, University of Buffalo, State University of New York, NY, United States
³Department of Mechanical Engineering, Universitat Rovira i Virgili, Tarragona, Spain
⁴Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Barcelona, Spain
⁵Shared Materials Instrumentation Facility (SMIF), School of Engineering, Duke University, Durham, NC, United States
⁶Department of Integrative Anatomical Sciences, Keck School of Medicine, University of Southern California, Los Angeles, CA, United States
⁷Human and Evolutionary Biology Section, Department of Biological Sciences, University of Southern California, Los Angeles, CA, United States
⁸Division of Anthropology, American Museum of Natural History, New York, NY, United States
⁹New York Consortium in Evolutionary Primatology, New York, NY, United States
¹⁰Department of Earth and Environmental Sciences, University of Manchester, Manchester, United Kingdom

Currently, there are no living platyrrhine primates inhabiting the main Caribbean islands. Nevertheless, the fossil record of this area has provided outstanding findings of different New World monkeys that were part of a diverse radiation exhibiting remarkably unusual morphologies. Among these, the Cuban genus Paralouatta corresponds to one of the most enigmatic primates ever found in the Greater Antilles. Some researchers have argued that Paralouatta’s post-cranium shows evidence of semi-terrestriality, a locomotor adaptation that is unusual, if not unique, in platyrrhine evolutionary history. Whether or not Paralouatta was truly semi-terrestrial remains uncertain, however, due to a lack of more sophisticated functional analyses on its morphology. Using novel virtual morpho-functional techniques on a comparative sample of 3D talar models belonging to diverse primate species representing three substrate preferences, this study aims to further evaluate whether Paralouatta was a semi-terrestrial genus or not. Geometric morphometrics and finite element analysis were used to empirically assess shape and biomechanical performance, respectively, and then several machine-learning (ML) classification algorithms were trained using both morphometric and biomechanical data to elucidate the substrate preference of the fossils. The ML algorithms categorized the Paralouatta specimens as either arboreal or as species commonly active on both ground and in trees. These mixed results are suggestive of some level of semi-terrestriality, thus representing the only known example of this locomotor behavior in platyrrhine evolutionary history.

Introduction

No extant platyrrhine primates are known in any of the Greater Antilles (i.e., Jamaica, Hispaniola, Cuba; Fleagle, 2013). However, during the Quaternary these Caribbean islands were home of some the most mysterious neotropical primates (Cooke et al., 2011), different from all extant species currently found in the New World (Cachel, 2015). At least four different genera belonging to this endemic radiation of unique platyrrhines have been described, including the following species: Antillothrix bernensis from the Dominican Republic (Rímoli, 1977; MacPhee and Iturralde-Vinent, 1995; Kay et al., 2011; Rosenberger et al., 2011); Xenothrix mcgregori from Jamaica (Williams and Koopman, 1952; Rosenberger, 1977; MacPhee and Fleagle, 1991; MacPhee and Meldrum, 2006); Insulacebus toussaintiana from Haiti (Cooke et al., 2011); and Paralouatta varonai and Paralouatta marianae from Cuba (Rivero and Arredondo, 1991; Horovitz and MacPhee, 1999; MacPhee and Meldrum, 2006). Most of these Antillean fossil primates are from the Quaternary, with some even surviving until the human settlement of these islands or even the time of the European conquest of the Caribbean (MacPhee and Fleagle, 1991; Turvey, 2009; Cooke et al., 2017a, b), at which time they went extinct along with most Quaternary non-volant mammal species, most likely due to anthropic pressures (Morgan and Woods, 1986; Cooke et al., 2017a, b).

These Caribbean primates represent a diverse radiation that displayed particularly unusual morphologies comprising several autapomorphies and seemingly derived conditions which exist across different platyrrhine groups. Even though clearly platyrrhines, they also display traits and feature combinations that are atypical or absent in any extant species (e.g., Xenothrix lacks third molars, and its first molars are much larger than the second ones). This morphological distinctiveness coupled with their noticeable differences when compared to extant continental taxa and the paucity of their remains has fueled an intense debate regarding the phylobiogeographic models that could explain the origin and adaptive radiation of this group. In fact, it has been challenging to determine the phylogenetic relationships of the Caribbean platyrrhine genera, and until recently they were not particularly clear. At different times, they have been associated with nearly all known platyrrhine sub-families (e.g., Williams and Koopman, 1952; Rosenberger, 1977; Rímoli, 1977; Ford and Morgan, 1986; Rivero and Arredondo, 1991; Woods and Ottenwalder, 1992; Horovitz, 1999; Rosenberger et al., 2011), classified as a sister group of living platyrrhines (Kay et al., 2011), or as stem platyrrhines that had established in the Greater Antilles in the Early Miocene (i.e., prior to the origin of modern New World monkey families; Kay, 2015). Several authors consider that the Caribbean endemic primates correspond to a monophyletic group, sister to Callicebus sensu lato (MacPhee and Horovitz, 2004), or to all crown platyrrhines (Kay, 2015), whereas others regard the striking morphological variation between the different species as evidence of multiple colonization events by different mainland lineages, or alternatively one multilineage colonization episode (Rosenberger, 2002). Fortunately, Woods et al. (2018) were able to retrieve DNA from a Xenothrix sample, providing several insights concerning this issue. They found that instead of being remotely related to living platyrrhines, Xenothrix corresponds to an extremely derived titi monkey (i.e., a member of the sub-family Callicebinae) that experienced significant morphological changes after arriving to Jamaica. Furthermore, they also found that based on the estimated splitting date between titi monkeys and Xenothrix (∼11 Ma), platyrrhines colonized the Greater Antilles more than once (Woods et al., 2018), since the oldest known Caribbean primate, P. marianae, was found in Miocene sediments dated to 17.5–18.5 Ma (MacPhee et al., 2003). This indicates that the Caribbean primates cannot be monophyletic and that there were at least two platyrrhine colonization events of the Greater Antilles at different times during the Neogene.

Among their many unique array of traits, Caribbean primates exhibit some of their most unusual in their post-cranial skeletons. For example, MacPhee and Fleagle (1991) considered that some aspects of the post-cranial morphology of Xenothrix (e.g., the adductor process of the femur) were so peculiarly dissimilar when compared to the post-crania of other platyrrhines, that they proposed that this species could have been a slow arboreal quadruped, thus representing a unique locomotion among platyrrhines. However, even more strikingly different are the set of post-cranial features of P. varonai which are thought to be related to a semi-terrestrial lifestyle (MacPhee and Meldrum, 2006). If correct, this interpretation would not only imply a locomotor convergence between some Old and a New World monkey species, but could also represent the first case of a locomotor behavior that does not exist in New World monkeys (among which all modern species are habitually arboreal) and that has no known analog in the platyrrhine fossil record. The amount of time that P. varonai may have spent on the ground, as well as to what degree it might have been similar to cercopithecids in postural and locomotor behaviors is yet to be resolved (MacPhee and Meldrum, 2006). Additionally, it remains unclear whether the older P. marianae exhibited any semi-terrestrial adaptations. Recently, Püschel et al. (2017, 2018) inferred the main locomotor behavior of P. marianae as one comparable to alouattines (i.e., exhibiting different levels of arboreal quadrupedalism, clambering and climbing), but it is important to bear in mind that they were not able to rule out possible semi-terrestrial adaptations, since they had not consider this category in their analyses.

Paralouatta marianae was originally described based on a single talus (MacPhee et al., 2003). It was described as only slightly different in morphology from that of P. varonai, in spite of the 17–18 Ma that allegedly separates them; the primary difference noted was that P. marianae’s talus is noticeably smaller (MacPhee and Meldrum, 2006). It has been argued that there is no appropriate morphological analog for the talus of Paralouatta amongst living platyrrhines (MacPhee and Iturralde-Vinent, 1995; MacPhee and Meldrum, 2006). MacPhee and Iturralde-Vinent (1995) described that atelid tali are different from Paralouatta in showing a “wedged” trochlea with a low trochlear relief (i.e., a trait associated with an increased mobility at the talocrural articulation), whereas Paralouatta displays a talus associated with increased stability. Additionally, the talus of Paralouatta has an evident cotylar fossa, hence providing a stable articulation for the medial malleolus (MacPhee and Iturralde-Vinent, 1995). This feature, that usually does not exist in large New World monkeys, is commonly seen in highly terrestrial cercopithecines e.g., Theropithecus (MacPhee and Meldrum, 2006). MacPhee and Meldrum (2006) used twelve linear measurements to compute a principal component analysis (PCA) of different tali. They found that Paralouatta’s talus is particularly distinct from the tali of other New World monkeys, especially because the absence of trochlear “wedging” distinguishes Paralouatta from all other large New World monkeys. In contrast, Püschel et al. (2017) performed a PCA of landmark data using geometric morphometrics (GM) and found that Paralouatta occupied a position in the resulting morphospace near Alouatta, as well as to some of the oldest platyrrhines from Patagonia (i.e., Dolichocebus, Soriacebus, and Carlocebus). Similarly, they applied a hierarchical clustering analysis that placed this fossil near Cebus/Sapajus and Dolichocebus, Soriacebus, and Carlocebus, hence showing that, at least from a morphological perspective, the talus of Paralouatta is not as unusual as initially thought (Püschel et al., 2017).

Paralouatta varonai was discovered in association with Late Quaternary fauna in Western Cuba, whilst P. marianae was discovered in Early Miocene (17.5–18.5 Ma) deposits, hence representing the oldest Caribbean platyrrhine known and establishing an early date for the arrival of platyrrhines to the Greater Antilles. P. varonai was originally considered to be similar to Alouatta (hence the genus name), since its estimated size was close to that of an extant alouattine, and the cranium was somewhat similar to that of a howler monkey (Rivero and Arredondo, 1991; Rosenberger et al., 2011). Further phylogenetic studies have either confirmed a possible alouattine connection (Rosenberger, 2002; Rosenberger et al., 2015), or have classified Paralouatta as related to the Callicebinae from South America (Horovitz and MacPhee, 1999; MacPhee and Horovitz, 2002, 2004). Meanwhile and as mentioned above, classical morpho-functional analyses (i.e., simple observations, linear measurements, ratio computations, etc.) of the post-cranial remains of Paralouatta suggested a semi-terrestrial locomotor mode for the genus, based on traits such as its short digits for its size, combination of deep olecranon, narrow trochlea, fossa and retroflexed medial epicondyle, among other traits. However, it is important to stress that MacPhee and Meldrum (2006) admit that Paralouatta exhibits an unforeseen mix of traits that when taken together differentiate this genus from all other known extant and extinct New World monkeys. Among these distinguishing features there are some that can be interpreted as evidence of some degree of terrestriality (i.e., based on the morphology of ground dwelling cercopithecoids), whereas other traits do not show this behavioral signal. Whether the traits that seem adaptive to terrestriality are actually indicative of this locomotor behavior, or rather represent another form of locomotor adaption not observed in extant platyrrhines or even anthropoids is currently not clear. Certainly, a more complete post-cranial fossil record of Paralouatta would contribute to our understanding of this issue. In the absence of more fossils, it is also possible to utilize virtual morpho-functional toolkits (i.e., “engineering toolbox”), which are increasingly being used to analyze both living and fossil functional morphology (Rayfield, 2019). Accordingly, the primary goal of the present study is to perform a battery of different computational analyses in order to test the hypothesis that Paralouatta was habitually semi-terrestrial. This is of importance because if semi-terrestriality is confirmed, this would represent the first case of such behavior among both extant and extinct New World monkeys.

This study specifically focused on the talus as it is the only post-cranial element representing both P. marianae and P. varonai. Additionally, the talus plays an informative arthrodial connection between the foot and the leg (i.e., ankle joint), because of its important role in joint mobility and stability, as well as weight support (Püschel et al., 2018). In fact, it has been demonstrated that there is a strong and significant covariation between locomotor behaviors and talar shape (Püschel et al., 2017). As a result, we decided to apply a set of different techniques to elucidate whether or not Paralouatta represents the first known semi-terrestrial platyrrhine species. First, GM was applied to characterize the talar morphology of Paralouatta among a large and diverse anthropoid sample (Püschel et al., 2017). Second, we analyzed the mechanical strength of the extant anthropoid talar sample by applying finite-element analysis (FEA), as well as carrying out the simulation in the two Paraloutta tali available, since it has been shown that talar biomechanical performance can also be used as locomotor proxy (Püschel et al., 2018). Finally, given that the main aim of the present work is to classify the fossil tali within a substrate preference, several machine-learning (ML) algorithms were trained using the morphometric and biomechanical data from the extant species.

Materials and Methods

Sample

The extant anthropoid sample included 3D surface renderings of tali from 109 individuals of 85 species representing all anthropoid sub-families; a large portion of the 3D data are available at https://www.morphosource.org/ (Boyer et al., 2016; Copes et al., 2016). Further details about the sample, including the Morphosource projects where the specimens can be found are provided in Supplementary Table S1. The fossil sample includes tali belonging to both Paralouatta species (Figure 1 and Supplementary Table S2). The extant sample was classified according to their main substrate preference based on the database of Galán-Acedo et al. (2019). This database provides some important ecological traits, including substrate preference (i.e., locomotion type in the database) for 497 primate species. This information was compiled through a meticulous review of 1,216 studies published between 1941 and 2018 (Galán-Acedo et al., 2019). The substrate preference categorization scheme is provided in Supplementary Table S1 and classifies each taxon as arboreal, terrestrial or as both substrate preferences (i.e., semi-terrestrial). This classification scheme refers to primary substrate preferences which reflects preferred or habitual environmental niche. The arboreal substrate preference comprises primate species that are strictly arboreal, which in undisturbed environments would seldom go to the ground; the terrestrial category considers primates that are mostly terrestrial (i.e., they carry out most of their daily activity on the ground), whilst the category “both” involves species that are regularly active on both substrates (i.e., ground and trees). Although rather simplistic, this classification system can be considered as a first approximation to locomotor behavior that can be easily applied to the fossil record without making further assumptions. In addition, this classification scheme primarily deals with our main goal which is to elucidate the possible nature of Paralouatta’s terrestriality.

FIGURE 1

Figure 1. Fossil sample under study.

3D Model Rendering

The refinement and smoothing tools available in Geomagic Studio^® (3D Systems, v. 2014.3.0), Rock Hill, SC, United States) were used to repair the irregularities observed in some of the 3D models. All the tali were aligned with respect to the same reference position (further details about the alignment are available in Supplementary Material S3). The talus of Paralouatta varonai is relatively incomplete (i.e., the talar head is missing, and other areas exhibit minor damage), whereas the one for Paralouatta marianae is fairly complete (only the talar head exhibits some minimal erosion). Therefore, different reconstruction procedures were applied to generate models suitable for both GM and FEA analyses (these reconstruction procedures are described in Supplementary Material S4).

Geometric Morphometrics

Thirty 3D landmarks were collected using Landmark editor v. 3.6 (Wiley et al., 2005) on the surface of the virtual tali in order to represent their morphological variation (Figure 2; Harcourt-Smith, 2002; Turley and Frost, 2013). These raw coordinates are available in Supplementary Material S5. These landmarks were analyzed using the “geomorph” R package (Adams et al., 2018), where a generalized Procrustes analysis (GPA) was performed (Rohlf and Slice, 1990). The GPA translates the landmark configurations to the origin, scales them to a standard size (i.e., unit-centroid size), and rotates them (by employing a least-squares criterion) against each other until a best-fit of corresponding landmarks to each other is achieved. The aim of the GPA is to remove non-shape information from anatomical objects. Hence, the resulting aligned Procrustes coordinates correspond to the shape variables of each specimen under analysis. These shape variables were used to carry out a PCA to decompose shape variation into orthogonal axes of maximum variation. In addition, a multivariate phylogenetic generalized least square regression (PGLS) was used to assess the association between shape and the logarithm of centroid size for the extant dataset (centroid size corresponds to the square root of the sum of squared distances of a set of landmarks from their centroid) (Bookstein, 1997). The PGLS was performed using the procD.pgls() function of the “geomorph” R package, which performs ANOVA and regression models in a phylogenetic context assuming Brownian motion, in such a way that can cope with high-dimensional data (Adams, 2014). This allowed us to evaluate the influence of size on shape when taking into account phylogenetic non-independence. Since the PGLS requires a phylogeny, we downloaded a consensus phylogeny from https://10ktrees.nunn-lab.org/Primates/ that was computed from 10,000 Bayesian trees that included most of the species present in our dataset (Arnold et al., 2010). This phylogeny was slightly adjusted to incorporate some species that were originally not present. The resulting phylogeny is available in Supplementary Material S6. All statistical analyses in this work were performed in R v. 3.6.0 (R Core Team, 2019).

FIGURE 2

Figure 2. A Chiropotes satanas talus (AMNH 95760) showing the 30 landmarks used in this study.

Finite Element Analysis

The virtual talar models were imported into ANSYS^® (Ansys, Inc., v. 17.1, Canonsburg, PA, United States)¹ to perform FEA. This engineering method reconstructs deformation, strain and stress in a digital structure after simulating a loading scenario (Richmond et al., 2005), and it is routinely used in vertebrate paleontology to tackle questions of organismal form, function, and evolution (Rayfield, 2007). In the present work we applied FEA to analyze talar mechanical strength, as it has been previously shown that talar biomechanical behavior can serve as locomotor proxy (Püschel et al., 2018). Each talus was modeled as a surface exclusively made of cortical bone; hence this required the use of shell elements in the finite element mesh. Homogeneous, linear and elastic material properties were used. Cortical bone values from a Homo sapiens talus were applied (Young’s modulus: 20.7 GPa; Poisson’s ratio: 0.3; Parr et al., 2013). The models were meshed using an adaptive mesh where the thickness of the shell elements representing the cortical bone was assumed to be constant. These values and further information about the FE models can be found in Supplementary Table S3. Cortical thickness values were obtained by measuring CT-data in some specimens and then, via a linear extrapolation of this data (see details of this procedure in Supplementary Material S8).

Loading Scenario and Boundary Conditions

Following the approach taken in Püschel et al. (2018), we applied a load on the trochlear surface of each talus in order to simulate a basic quadrupedal scenario. We decided to only model a generalized standing posture for all taxa since talar arthrokinematics are unknown in nearly all primates. The load was directed in the z-axis on the aligned talar models to simulate gravity and distributed on the trochlea to simulate a compressive force. The talus was constrained as indicated in Figure 3.

FIGURE 3

Figure 3. Loading scenario tested in the FEA illustrated using a talus of Alouatta caraya (AMNH 211513). The large arrows represent the applied load on the trochlear surface.

In this study, the values of the compressive force of each model were calculated using a quasi-homothetic transformation for planar models (Marcé Nogué et al., 2013). This methodology is also appropriate for shell models because the scaling of the forces is a function of the difference between the surface and overall thickness of the model. We used the talus of Alouatta caraya as a reference and applied a force of F = 15.8186 N. These values were obtained following Püschel et al. (2018) which computed the applied force as the 30% of the average body mass of the species multiplied by standard acceleration due to gravity g = 9.81 ms^–2. The scaled values of the force for the other species can be found in Supplementary Table S3.

Stress Values and Intervals’ Method

We computed the von Mises stress because it combines all Cartesian components of the stress tensor into a single value (equivalent stress; Zienkiewicz et al., 2005). This enables easier and more understandable comparisons when assessing different models, and it has been used successfully as a proxy to compare the strength of bony structures (Marcé-Nogué et al., 2017b), where species with lower values of stress represented the stronger structures. Furthermore, it has been demonstrated that if the bone is modeled using isotropic material properties, then the von Mises criterion is the most suitable one when comparing of stress in bones (Doblaré et al., 2004).

New variables corresponding to different intervals of stress values were computed following the Intervals’ Method proposed by Marcé-Nogué et al. (2017a). These interval variables were then used to analyze the FEA results. The Intervals’ Method generates a set of variables, each one defined by an interval of stress values. Once all the stress values of a single specimen were obtained, all elements of the model are sorted according to their von Mises stress values. These are then grouped into a finite number of intervals, defined by their range of stress values. Each one of the intervals represents the amount of the volume of the original model (as a percentage value) exhibiting a specific range of stress values. This method allowed us to analyze the data from finite element models in a comparative multivariate framework. The method of Marcé-Nogué et al. (2017a) needs the definition of a Fixed Upper Threshold (i.e., FTupper = 12 MPa) and a number of intervals (i.e., N = 100). This value is obtained based on a convergence procedure based on a PCA that was performed to define the minimum number of intervals (Marcé-Nogué et al., 2017a). We considered that convergence has been achieved when the correlation values of the PCs of the intervals were higher than 0.99. The values of each vector for stress interval when N = 100 can be found in Supplementary Table S4. These newly generated variables were analyzed using a PCA performed on the correlation matrix.

Fossil Substrate Preference Classification

Previous studies have indicated that it was possible to distinguish between different locomotor modes using talar shape or stress information (Püschel et al., 2017, 2018). Consequently, we applied the same approach here but using different substrate preferences (i.e., arboreal, terrestrial, or both substrate preferences). Two different datasets were analyzed and used to elucidate the main substrate preference of Paralouatta: (a) morphometric and (b) biomechanical data. Each one of these datasets corresponded to the PCs that accounted for 90% of the variance of the sample using the shape and interval variables, respectively. Two pairwise PERMANOVA tests with a Holm correction for multiple comparisons were performed to assess for differences between the three substrate preferences using both the morphometric and biomechanical datasets. In both cases, Euclidean distances were selected as similarity index.

Six well-known supervised algorithms were chosen as they correspond to a diverse range of different classification techniques: (a) linear discriminant analysis (LDA); (b) classification and regression tree (CART); (c) k-nearest neighbors (KNN); (d) naïve Bayes (NB); (e) support vector machine (SVM); and (f) random forest (RF) (further details about this models can be found in Püschel et al., 2018). All the models were run using the “caret” package for R (Kuhn, 2008, 2015). Performance was calculated using the confusion matrix from which the overall classification accuracy was computed, as well as Cohen’s Kappa (Kuhn and Johnson, 2013a). The complete dataset was resampled using a “leave-group-out” cross-validation (Kuhn and Johnson, 2013b), as a way to asses classification performance. This cross-validation strategy generates multiple splits of the data into modeling and prediction sets. This process was carried out 200 times and the data were split into a modeling sub-set comprising 75% of randomly assigned observations, whereas the testing sub-set considered the remaining 25%. The number of repeats was chosen to get a consistent classification performance and to minimize uncertainty. The models with the best classification performance for the morphometric and biomechanical data were then applied to deduce the principal substrate preference of both Paralouatta specimens by calculating their category probabilities.

Finally, the Euclidean distances obtained for the morphometric and biomechanical datasets were used to visualize morphological and biomechanical affinities between the extant species and the Paralouatta fossils. An unweighted pair-group average (UPGMA) algorithm for agglomerative hierarchical cluster analysis was used to generate two dendrograms (i.e., biomechanical and morphometric affinity dendrograms) that allowed us to assess general affinities (Sokal and Rohlf, 1962). Cophenetic correlation coefficients (CPCC) were computed as a way of measuring how closely the obtained dendrograms preserved the pairwise distances between the specimens (Sokal and Rohlf, 1962).

Results

Geometric Morphometrics Results

The PCA performed using the shape variables obtained through a GPA returned 90 PCs. The first 24 PCs accounted for 90% of the total variance of the sample, hence offering a reasonable estimate of the total amount of talar shape variation, which were then used in the classification analyses (i.e., morphometric dataset). The PCA shows the main regions of occupied shape space (Figure 4A). Platyrrhines (which are almost exclusively arboreal; only Cebus albifrons is considered into the “both” category) are located on the extreme positive side of PC1 (i.e., lower right and extreme right of the upper right quadrants), whilst most cercopithecid monkeys occupy the upper left quadrant showing mixed substrate preferences. Apes are located on the lower left quadrant, with gorilla and chimpanzees displaying the lowest scores on PC2. PC1 seems to mostly distinguish apes from platyrrhines, with cercopithecids occupying an intermediate position. Knuckle-walking terrestrial African apes (i.e., chimpanzees and gorillas) show the minimum values along PC2, followed by habitually suspensory genera including Pongo, Ateles, Lagothrix, and Hoolock. The most positive PC2 values are shown by cercopithecines and colobines. Interestingly, Paralouatta fossils are located in an intermediate position between cercopithecids and platyrrhines, with P. varonai being closer to the former rather than to the latter. The variation on the negative side of PC1 can be related to a shorter anterior and a longer posterior calcaneal facet, and a broader and shorter talar head. The morphological variation on the positive side of PC1 can be associated with an increased anterior calcaneal facet and an antero-posteriorly shorter trochlea.

FIGURE 4

Figure 4. (A) Principal component analysis (PCA) of the talar shape variables (only the two first PCs are shown). The models closest to the mean shape (i.e., Cercopithecus mitis) was warped to match the multivariate mean using the thin plate spline method. The obtained average model was then warped to represent the variation along the two plotted PC axes in both analyses (mag = 2); (B) Principal component analysis (PCA) of the stress data (i.e., intervals) obtained from the finite element simulations. The models located at the extremes of the PC axes correspond to the individuals showing the maximum or minimum values along the two first PCs.

The PGLS regression showed that there is no association between talar shape and the logarithm of talar centroid size (R²: 0.015; F: 1.325; Z: 1.027; p-value: 0.151; 9,999 permutations). This means that talar shape variation cannot be attributed to evolutionary allometric effects.

Finite-Element Analysis Results

Figure 5 displays stress maps for all the tali under analysis. By representing stress in a visual manner, it is possible to make qualitative inferences about biomechanical behavior, since the observed stress patterns can be read in terms of relative strength (i.e., individuals showing higher stress levels are weaker under the applied load). The correlation based PCA carried out using the 100 variables generated using the Intervals’ method returned 100 PCs. The first nine PCs that accounted for 90% of the total variance of the sample were used in the classification analyses (i.e., biomechanical dataset) as a way to avoid collinearity and to reduce the dimensionality of our data. Figure 4B shows the first two PCs of this correlation based PCA.

FIGURE 5

Figure 5. von Mises stress distribution for all the analyzed specimens.

When comparing locomotor behaviors in extant species, most specimens display moderate values, which makes it difficult to stablish a clear pattern with this proxy. However, the weakest tali -with high average values of stress- correspond to some of the arboreal species. Figure 4B shows the PCA of the stress data (i.e., intervals) obtained from the FE simulations. Mostly terrestrial and arboreal/terrestrial (i.e., “both” category) occupy the lower left quadrant of the PCA, whereas exclusively arboreal species are spread over the three remaining quadrants. This is probably because the “arboreal” category comprises several different locomotor styles (e.g., climbing/clambering, leaping, arboreal quadrupedalism, etc.). On the positive side of PC1 are located those individuals with intervals corresponding to higher stress values, whilst those exhibiting intervals with lower stress values are located at the opposite extreme of the axis. In general, suspensory species tend to show higher stress values than terrestrial and arboreal/terrestrial (i.e., “both” category) specimens. The two Paralouatta specimens present low average values of stress in coincidence with both arboreal and terrestrial species and are located in the PCA close to the limit between the upper and the lower left quadrant, exhibiting lower stress values in an area occupied by many species using both ground and arboreal substrates.

Fossil Substrate Classification Results

There were significant differences between all substrate preferences when using the morphometric dataset (i.e., 24 PCs from the PCA carried out using the shape variables) (Table 1a.). When the same categories are compared but using the biomechanical dataset (nine PCs from the PCA performed using the stress intervals), there were no significant differences between “both” and “terrestrial” categories. Nonetheless, there were significant differences between the arboreal group and the other two categories (Table 1b).

TABLE 1

Table 1. Pairwise PERMANOVA results.

Figure 6 displays the accuracy and Cohen’s Kappa results for all the applied models, for both the morphometric and biomechanical data after carrying out the “leave-group-out” cross-validation procedure and after applying an automatic grid search. The morphometric data slightly outperformed the biomechanical data when classifying substrate preferences in both accuracy and Cohen’s Kappa metrics. The best model for the morphometric data was a simple LDA model (Table 2a), whilst in the case of the biomechanical data, the best model was the KNN (Table 2b), although other models (i.e., LDA and NB) showed similar performance levels.

TABLE 2

Table 2. Summary statistics of the “leave-group-out” cross-validation procedure for the accuracy and Cohen’s kappa values for all the tested models (best models in bold).

FIGURE 6

Figure 6. Dot-plot comparing the accuracy and Cohen’s Kappa values of the different classification models applied to morphometric and biomechanical data. The dots represent the average accuracy and Cohen’s Kappa values after performing the “leave-group-out” cross-validation (200 repeats), while the whiskers display their respective 0.95 confidence level. Model acronyms: linear discriminant analysis (LDA); classification and regression tree (CART); k-nearest neighbors (KNN); naïve Bayes (NB); support vector machine (SVM) and random forest (RF).

There are no extra parameters in the morphometric LDA, so no further tuning was required (average accuracy: 0.824; average Cohen’s Kappa: 0.623; Figure 7A). The variables that contributed the most to the separation between categories were PC1 and PC5, followed by PC10 and PC13 (Figure 7B). Then, this model was used to classify the Paralouatta sample into the analyzed substrate preferences (i.e., arboreal, terrestrial and both). Using morphometric data, P. marianae was classified as an arboreal individual, whereas P. varonai was classified as “both” (i.e., arboreal/terrestrial species; Table 3). The obtained KNN model for the biomechanical data achieved its best performance with k = 7 (number of neighbors ranging from 5 to 23 were tested; average accuracy: 0.74; average Cohen’s Kappa: 0.42; Figure 7C). The variables that contributed the most to the separation between categories were PC1 and PC2, followed by PC4 and PC7 (Figure 7D). By applying the final KNN model, both Paralouatta specimens were classified into the “both” category (Table 3). This likely indicates that they engaged in activities both in the ground as in trees.

FIGURE 7

Figure 7. Decision boundary plots for (A) morphometric and (C) biomechanical data. The two variables that contributed the most to each one of the models are displayed. The space is colored depending on what substrate preference the (A) LDA or the (B) KNN algorithm predict that region belongs to, whereas the lines between colored areas represent the decision boundaries. Color intensity indicates the certainty of the prediction in a particular graph area (i.e., darker colors mean a higher probability of belonging to a certain category). Symbols surrounded by a white rim correspond to misclassified specimens. In addition, variable importance scores for the predictors used are provided for the (B) morphometric (only the top 20 variables are shown) and (D) biomechanical models.

TABLE 3

Table 3. Prediction results for the fossil sample.

The agglomerative-hierarchical cluster analysis of the PCs using the UPGMA algorithm displayed the affinities between living species and Paralouatta fossils for both the morphometric and biomechanical data (Figures 8A,B, respectively). The CPCCs of both datasets indicate a reasonable agreement between the cophenetic distances obtained from the trees, and the original Euclidean distances (morphometric CPCC: 0.78; biomechanical CPCC: 0.79) (Farris, 1969; Saraçli et al., 2013). In the case of the morphometric data, the clusters seem to predominantly reflect phylogenetic relatedness (i.e., platyrrhines, cercopithecids, and apes). In fact, the two Paralouatta fossils are located within the platyrrhine cluster next to Chiropotes satanas. The biomechanical data seems to mostly reflect a combination of substrate preference categories (i.e., most terrestrial and “both” species are located closer to each other), as well as phylogenetic relatedness. Interestingly, the Paralouatta fossils clustered next terrestrial (e.g., Papio and Theropithecus) or arboreal/terrestrial (e.g., Macaca) individuals.

FIGURE 8

Figure 8. UPGMA hierarchical clustering analyses of (A) shape and (B) biomechanical PCs. Fossils are in bold.

Discussion

The diverse post-cranial adaptations exhibited by extinct platyrrhine primates demonstrates that they occupied a wide variety of habitats and environments in the Americas (MacPhee and Horovitz, 2002). Among the different fossil New World monkeys, the origin and relationships of the endemic Caribbean primates are possibly one of the least understood aspects of platyrrhine evolution. These Antillean monkeys are of special interest for different reasons, including their intricate phylogenetic relationships to mainland forms, as well as due to their mysterious biogeographic history (Halenar et al., 2017). Nonetheless, the fact that they display a variety of features that are rare or absent in platyrrhines from the mainland is perhaps one of the most intriguing and enigmatic ones. It is likely that these traits appeared as a response to selective pressures that are particular to island environments, since insular species tend to significantly differentiate from mainland forms (e.g., non-volant small mammal species tend to evolve larger bodies whereas non-volant large mammal species tend to evolve smaller bodies as has been summarized by Foster’s rule) (Foster, 1964; Case, 1978; Lomolino, 2005; Whittaker and Fernández-Palacios, 2007). Among the features distinguishing the extinct Caribbean species are several post-cranial traits that seem to indicate locomotor adaptations not seen in extant mainland taxa. The original morphological description provided by MacPhee and Meldrum (2006) suggested that Paralouatta was a primate adapted to semi-terrestrial locomotion, perhaps similar to some Old World monkeys. This is particularly striking since no extant platyrrhine species within this adaptive radiation is known to exhibit frequent terrestrial behaviors as part of their locomotor repertoire. Nonetheless, MacPhee and Meldrum (2006) also mentioned that some of these terrestrial features observed in Paralouatta could be merely plesiomorphic traits shared with other anthropoid ancestors. In order to elucidate this issue, we applied a combined approach from the virtual morpho-functional toolkit (i.e., “engineering toolbox”), analyzing both morphometric and biomechanical data. We have arranged what is, to our knowledge, the largest number of 3D FE models to carry out meaningful comparisons that can contribute to a better understanding of the problem. Nevertheless, to properly address the polarity of changes in Paralouatta (e.g., to establish whether the talar traits are plesiomorphic or not) phylogenetically informed methods need to be applied (Almécija et al., 2019).

From a morphometric perspective our results indicate mixed locomotor behaviors. The morphometric PCA located the two Paralouatta specimens between the cercopithecids and the platyrrhines. This indicates that at least for the main axes of morphological variation (i.e., PC1 and PC2) Paralouatta’s talus seems quite distinctive, occupying an area of the morphospace almost vacated by other species. The classification algorithm (i.e., LDA) using the morphometric data classified P. marianae as an arboreal species, while P. varonai was categorized into the “both” category signifying a preference for a semi-terrestrial lifestyle. This is in agreement with Püschel et al. (2018), where it was found that P. marianae was most likely a clamber/suspensory species. Nevertheless, from the observation of Figure 7A, which displays the two variables that contribute the most to the classification models, it is evident that both Paralouatta specimens are located quite closely to the boundary between the arboreal and terrestrial categories. In general, it can be concluded that for P. varonai the traits that indicate semi-terrestriality are more pronounced when compared to P. marianae. If we take into account the millions of years separating these two species, one can speculate that the terrestrial behaviors which were more incipient (less frequent) in P. marianae became more ubiquitous in the later species P. varonai. The clustering analysis of the morphometric data mostly showed broad phylogenetic relatedness. The two Paralouatta specimens where located close to Chiropotes satanas a species known for its usual above branch quadrupedal locomotion, as well as some suspensory postures (Fleagle and Meldrum, 1988). However, it is important to notice that when carrying out the same analysis but using fewer morphometric PCs (e.g., 14 PCs that account for 80% of the variance of the sample), the fossils are located within the predominantly cercopithecid cluster, close to terrestrial or arboreal/terrestrial (i.e., “both”) specimens (results not shown), thus indicating that additional PCs provide additional phylogenetic (rather than “functional”) information.

The biomechanical results also indicate some mixed locomotor behaviors. The biomechanical PCA using the intervals showed that the two Paralouatta species are close to the “both” and terrestrial categories. However, it is important to stress that the “arboreal” species occupy most of the biomechanical-space and that many “arboreal” species are also close to the fossil specimens, probably because this category encompasses several other locomotor categories (e.g., leaping, climbing/clambering, arboreal quadrupedalism, etc.). The KNN algorithm classified both fossil species into the “both” category, hence indicating some levels of terrestriality. Interestingly, P. marianae also shows a high posterior probability for the “arboreal” category which suggests that the terrestrial traits observed in Paralouatta are more incipient in the early evolution of this genus (i.e., semi-terrestriality is certainly a derived character). The clustering analysis shows that the Paralouatta fossils grouped close to terrestrial or arboreal/terrestrial (i.e., “both”) individuals, which again shows some level of terrestriality.

Although the obtained results shown here are highly informative, there is certainly room for improvement. One limitation of our approach is the classification scheme applied. In fact, it is evident from both morphometric and biomechanical analyses that the arboreal category encompasses several locomotor styles that can vary greatly. Further studies could refine this classification scheme to provide a finer perspective when carrying out locomotor classifications. An alternative option, which would not require to force any species into discrete categories, would be to establish major patterns of covariation between a given shape data matrix (e.g., talar shape) and a locomotor behavior data matrix (see for e.g., Table 3 in Hunt, 2016). This approach has been successfully applied using two-block partial least squares analysis in a few studies (e.g., Almécija et al., 2015; Püschel et al., 2017) and it is certainly worth further exploration. Another limitation in our study arises from the fact that we only simulated one simple loading scenario (i.e., quadrupedal standing), which may not reflect the most realistic loading scenarios of the talus in order to distinguish between habitual substrate preferences adopted by primates. Hence, further works could simulate different loading scenarios and test their different discriminatory capabilities when elucidating locomotor behaviors.

Paralouatta was probably an island-adapted large-bodied genus that most likely diverged from other platyrrhines during the early Miocene. This would certainly explain the similarities of Paralouatta to the other platyrrhines, as well as many traits that are evidently unique to this genus and that seem to be exaggerated in the later species P. varonai. The talar morphology of Paralouatta combines some more primitive morphological aspects (both anthropoid and platyrrhine) with derived features associated to some terrestriality levels, as initially though by MacPhee and Meldrum (2006). Given that selection pressures and ecologies can vary significantly between islands and the mainland, different adaptations associated with species that are endemic to islands are to be expected. Paralouatta adapted to a different environment and probably employed a different locomotion, which based on our results it is highly likely to have involved a significant level of terrestrial activity, while still retaining arboreal behaviors as shown by Püschel et al. (2017, 2018). By the island rule (i.e., Foster’s rule) small sized species are expected to become larger, whilst large species tend to become smaller (Foster, 1964). It is surely suggestive to think that when Paralouatta arrived in Cuba it might have become larger due to a more relaxed predation pressure, and therefore able to shift from an arboreal lifestyle to a more terrestrial one. However, this interpretation should be treated with caution, since more refined phylogenetically informed analyses are required to better establish the polarity of changes. Further trait-evolution phylogenetic comparative analyses could certainly contribute toward settling this issue, but firstly the phylogenetic position of Paralouatta needs to be fully resolved.

Finally, this study has shown that a combined virtual morpho-functional approach can help to the understanding of locomotor behaviors in other fossil taxa. By combining morphometrics, biomechanics and ML methods it is possible to provide a broader perspective regarding the locomotor behaviors of fossils species by analyzing different aspects of their functional morphology. The proposed methodological approach can certainly be beneficial when figuring out not only locomotion in fossil species, but also when assessing any other past behaviors that can be inferred from their morphology.

Data Availability Statement

All datasets generated for this study are included in the article/Supplementary Material.

Ethics Statement

Ethical review and approval was not required for the animal study because we exclusively analyzed bones and fossils belonging to different museum collections. We did not perform any experiments on living animals.

Author Contributions

TP and JM-N designed the study and carried out the different analyses. TP, SA, BP, and JG collected and analyzed data. All authors interpreted the data and wrote the manuscript.

Funding

This project was funded by the NERC (NE/R011168/1). TP was funded by the Leverhulme Trust Early Career Fellowship, ECF-2018-264, SA by the NSF-BCS 1316947, AEI/FEDER EU (CGL2017-82654-P), and BP by the National Science Foundation (BCS-1317047) and the Leakey Foundation. JM acknowledges Dr. Juan Liu and the Laboratory for the Evolution and Anatomy of Fish of the University at Buffalo. Both SA and JM-N acknowledge the CERCA Programe (Generalitat de Catalunya).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with one of the authors SA.

Acknowledgments

We are grateful to the editors who invited us to take part of the “Evolving Virtual and Computational Palaeontology Research Topic” and to the two reviewers whose helpful comments improved this manuscript. A version of this work has been submitted as an abstract to the Program of the 89th Annual Meeting of the American Association of Physical Anthropologists (Püschel et al., 2020).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.00079/full#supplementary-material

References

Adams, D. C. (2014). A method for assessing phylogenetic least squares models for shape and other high-dimensional multivariate data. Evolution 68, 2675–2688. doi: 10.1111/evo.12463

PubMed Abstract | CrossRef Full Text | Google Scholar

Adams, D. C., Collyer, M. L., and Kaliontzopoulou, A. (2018). Geomorph: Software for Geometric Morphometric Analysis. R Package Version 3.0.6. Available online at: http://CRAN.R-project.org/package=geomorph (accessed November 1, 2019).