1For the past 30 years, virtual imaging has increasingly been applied to non-medical fields of research. Especially in palaeontology and palaeoanthropology (e.g., Balzeau et al., 2012). This methodological revolution has been accompanied by the development of geometric morphometric methods (GMM, i.e., quantitative analysis that allows to explain and visualise the differences in shape (i.e., form without size) that have been mathematically analysed (Zelditch et al., 2004; Adams et al., 2011). Both methodologies have been adopted by taxonomists and are now widely used to test hypotheses on how to classify specimens of extinct species (Mounier et al., 2011; Mounier, 2012; Profico et al., 2016).
2Nevertheless, the rise of these new tools largely ignored a previous methodological improvement in taxonomy: cladistics. Cladistics renewed systematics by introducing in the analysis the concept of ancestor-to-descendant relationships (Matile et al., 1987), while previous taxonomic analyses relied on a phenetic theoretical framework that focused solely on the morphological similarities (Sneath and Sokal, 1973). Morphological data obtained from virtual imaging and geometric morphometrics can be directly analysed within a phenetic theoretical framework but need to be adapted in order to work within a cladistic phylogenetic framework (MacLeod et al., 2002).
3Nevertheless, identifying phylogenetic signals from geometric morphometric data has been the focus of several scientific works (MacLeod et al., 2002; Klingenberg and Gidaszewski, 2010) and a few protocols for integrating this type of data into cladistic frameworks have been developed (González-José et al., 2008; Clouse et al., 2011). However, such protocols remain anecdotal, are heavily criticised (Catalano et al., 2010; Adams et al., 2011) and only occasionally used in palaeontology (Smith and Hendricks, 2013). Three main issues are raised by researchers as to why the use of cladistics with GMM should be avoided.
4First, should shape data be considered as one trait or as a set of traits, and, if so, how should the partitioning into separate characters or characters’ states be performed (Bookstein, 1994; MacLeod et al., 2002; Zelditch et al., 2004)? For analytical procedures operating in a phenetic framework, such as neighbour-joining or those based on likelihood (e.g., maximum likelihood), the question of partitioning does not arise, because the phenetic inferences are made on the distances between taxa (Swofford et al, 1996; Felsenstein, 2004), so that the choice of a particular set of shape variables does not affect the result as long as information on the relative positions of taxa within a multidimensional space is used (Klingenberg and Gidaszewski, 2010). Thus, the definition of a trait (i.e., its coding) is not central and the very notion can be abandoned. On the other hand, for cladistic methods based on Wagner’s parsimony (i.e., those which only take into account apomorphies (derived) and whose nodes are built on the degree of shared derived traits), the character need to be well-defined and are assumed to be independent from one another. Such is not the case for 3D morphometric data as the landmark coordinates which represent the extracted shape variables cannot be considered as independent of each other and their definition as discrete traits is, at best, debatable. Each landmark does contain morphological information, but the information itself is only revealed when landmarks are associated with each other.
5In order to address these questions, several protocols have been tested in the past. One consisted of using partial distortion scores for landmarks that were "identical" between all the specimens/taxa included in an analysis (e.g., Zelditch et al., 1995; 2000). This method was highly controversial (e.g., Bookstein, 1994; 2002; Monteiro, 2000; MacLeod et al., 2002) and was subsequently abandoned (Zelditch et al., 2004), as the choice of data was considered to be too arbitrary. Another approach subdivided the global shape of the analysed operational taxonomic units (OTU) into smaller parts and considered each of these parts as discrete traits (MacLeod et al., 2002; González-José et al., 2008). While this approach mathematically addressed the issue regarding independence of the variables used to perform the cladistic search – the variables in the cladistic analysis are imputed from the coordinates of the mathematically independent Principal Components (PCs) obtained after a Principal Component Analysis (PCA) –, the way each variable was calculated raised additional concerns. Indeed, the global shapes of the OTUs were divided into shape blocks which were aligned separately through multiple Generalised Procrustes Analyses (GPA). Unfortunately, this creates distortion in the database, as the PCs are derived from several GPAs and PCAs (Adams et al., 2011; Klingenberg and Gidaszewski, 2010). The use of uncorrelated shape variables, as can be obtained from PCA based on terminal taxa (González-José et al., 2011; MacLeod et al., 2002), has also been criticised, as the covariation between terminal taxa is generally not the same as the covariation between evolutionary changes in conformational variables (Felsenstein, 1985; 2004).
6While the issue of trait definition and independence is at the root of many controversies surrounding the use of morphometric data (Rohlf, 1998; Monteiro, 2000; Bookstein, 2002), it should be noted that it is also a problem when using "traditional" discrete morphological character because their definitions vary widely (e.g., Thiele, 1993) as it is often difficult to assess their independence (e.g., MacLeod et al., 2002; Mounier et al., 2009).
7Finally, according to Adams et al. (2011), methods consisting of forcing multivariate shape variables into a form compatible with cladistic software (Zelditch et al., 1995; Gonzalez-Jose et al., 2008) would distort the information present in these variables and this would severely limit the useful biological conclusions that can be made. This is in fact the case for all implementations of discontinuous coding from a continuous dataset and is therefore not specific to geometric morphometric data.
8In order to further contribute to this long-lasting debate, the present study tests two additional protocols that could help to merge the methodological strengths of both approaches: geometric morphometrics to describe and analyse morphological variations, with the aim to infer taxonomy and phenetic relationships from 3D morphological data, and cladistics to reconstruct phylogenies. These methodological procedures aim at analysing a large corpus of 3D data describing the calvarium (i.e., the skull without the mandible and upper face) of specimens belonging to the genera Pongo, Gorilla, Pan and Homo, while addressing the methodological concerns raised in the past regarding cladistics based on geometric morphometric data. More specifically, the first procedure uses PC coordinates as variables for the cladistic analyses obtained from a unique GPA alignment of the OTU 3D coordinates, hence addressing both criticisms regarding mathematical independence of the variables used as traits and the distortion of the database due to multiple GPA alignments. Another raised issue – distortion of the morphological information of the data due to the transformation of multivariate shape variables into a form compatible with cladistic software – is also addressed in the present study, as the TNT software used was developed to accommodate different types of continuous data, including 3D morphometric data (Goloboff et al., 2008).
9Finally, additional possible issues regarding the use of GMM in cladistics will be considered through a comparison of the first protocol with a second approach, using the 3D coordinates as variables, and with results from "traditional" numerical taxonomic methods, neighbour-joining and UPGMA, allowing the evaluation of the methodological limitations and strengths of the different protocols.
10The 3D models were obtained by one of the authors (AM) through three different procedures depending on the availability of the equipment/specimens: 1) medical computed tomographic scans processed in Amira (v5.5, FEI); 2) photogrammetry using the software Metashape (v1.8.1); and 3) 3D surface scans using an optical scanner (HDI Advance, 45µ accuracy, LMI) and the software Flexscan (V.3.3 LMI) (see Supplementary Information table S1). The shape of the calvarium was described using landmarks and semi-landmarks, the latter being allowed to slide (Bookstein, 1997), according to three protocols: A) the first contains 27 landmarks and 121 semi-landmarks (figure 1A; table 1); B) the second contains 27 landmarks and 320 semi-landmarks (figure 1B; table 1); and C) the third contains 27 landmarks and 608 semi-landmarks (figure 1C; table 1). Those datasets were collected with the software Landmark (IDAV, Wiley, 2005).
11The data used as variables in our study are composed, for the first set of analyses (tables S2, S4 and S6), of all the coordinates of the PCs generated from a PCA (Abdi and Williams, 2010) led on the aligned 3D landmarks (i.e., after GPA) and, for the second set of analyses, of aligned (i.e., after GPA) 3D coordinates of anatomical landmarks describing the morphology of the calvarium of cast and original specimens (tables S3, S5 and S7).
Figure 1
Position of the landmarks on the calvarium for each dataset (A. 148 landmarks; B. 347 landmarks; C. 636 landmarks). From left to right: norma facialis; 3/4 view; and norma basalis |
Position des points-repères sur le calvarium pour chaque base de données (A. 148 points-repères ; B. 347 points-repères ; C. 636 points-repères). De gauche à droite : norma facialis, vue de trois-quarts, et norma basalis
Table 1
Summary of the 27 landmarks used in the study and their abbreviations |
Présentation des 27 points-repères utilisés dans l’étude
12To implement and test phylogenetic reconstructions based on 3D geometric morphometrics, different taxa belonging to the family Hominidae Gray, 1825 (e.g., Lockwood and Tobias, 1999) were analysed (table S1). The outgroups are composed of three specimens for each of the genera Pongo, Gorilla and Pan. We used female individuals in order to avoid any distortion of the results due to the presence of strongly expressed secondary sexual characters (e.g., the sagittal crest) in some ape males. The 20 other OTU are formed by representatives of the Homo and Australopithecus genera.
13The specimens STS5 ("Mrs. Ples") and STS71 constitute the OTUs corresponding to the species A. africanus Dart, 1925. This species was first considered to be the common ancestor of two later hominid taxa, A. robustus/boisei and Homo according to Tobias (1980). The specimen STS5 has been the subject of much debate, including its taxonomic attribution, with some authors classifying it as A. africanus (Lockwood and Tobias, 1999; Villmoare et al., 2015) and others as A. afarensis (Clarke, 2008).
14Several taxa of the genus Homo representing a broad spectrum of morphological and taxonomical variations were selected (table S1). The two specimens KNM-ER 1470 and KNM-ER 1813 form the OTU representing the species H. habilis Leakey, Tobias and Napier, 1964 sensu lato. The attribution of these two specimens to the same species has been the subject of much debate (Prat, 2022). Indeed, KNM-ER 1470 has been attributed to H. habilis (Walker and Leakey, 1978), H. rudolfensis Alexeev, 1986 (Zeitoun, 2000; Argue, 2017) and to the genus Australopithecus (Walker and Leakey, 1978), while some others conclude it was a sister group of H. habilis and H. rudolfensis (Groves, 1989). The KNM-ER 1813 skull has sometimes been assigned to the species A. africanus (Walker and Leakey, 1978) and H. ergaster Groves and Mazák, 1975 (Groves, 1989; Zeitoun, 2000), but is generally classified as H. habilis (Lieberman et al., 1996).
15Two specimens from Dmanisi constitute the OTU representing the species H. georgicus Gabounia, de Lumley, Vekua and Lordkipanidze, 2002: D2282 and D4500. Analysis of the Dmanisi skulls, particularly that of D4500, has also been used to suggest that some previously recognised species or subspecies within the genus Homo correspond to a single evolutionary lineage: H. erectus (Lordkipanidze et al., 2013).
16Another OTU was formed by combining two specimens, attributed to the species H. ergaster Groves and Mazák, 1975: KNM-ER 3733 and 3883. Initially assigned to the species H. habilis (White et al., 1981), these specimens were assigned later to H. erectus (Walker et al., 1993; Grine et al., 1996) before being incorporated into H. ergaster (Groves, 1989; Tattersall, 2013) or even to be affiliated each to two new species (Zeitoun, 2000).
17Regarding the species H. erectus Dubois, 1893, we used three pairs of specimens belonging to the Sangiran, Sinanthropus and Ngandong series, in order to illustrate the chronological and geographic diversity of this polytypic species. The Ngandong 7 and 14 specimens, first described as belonging to a new species H. (Javanthropus) soloensis by Oppenoorth in 1932, were later assigned to different categories: "archaic" H. sapiens or evolved H. erectus (Zeitoun, 2000; Zeitoun et al., 2010).
18The specimens Sangiran 2 and 17 were initially assigned to the species Pithecanthropus erectus Dubois, 1894 by Koenigswald (1940) before being considered archaic H. erectus (Sartono, 1982).
19Finally, the specimens Sinanthropus III and XI were used. Initially described as belonging to Sinanthropus pekinensis Black, 1927, Sinanthropus XI was included later in an H. sapiens clade (Zeitoun, 2000) before being affiliated with H. erectus.
20Regarding the species H. neanderthalensis King, 1864, the choice was made to consider taxonomic units that reflect populations recognised as distinct based on genetic data and otherwise belonging to different periods and/or geographical areas (Fabre et al., 2009). Thus, the taxonomic operational unit "early Neanderthal" includes two specimens: Saccopastore 1 and Ehringsdorf H. A "Near East Neanderthal" taxonomic operational unit consists of the specimens Amud 1 and Shanidar 1. The Monte Circeo, La Ferrassie 1 and La Chapelle-aux-Saints specimens constitute the "Classic Neanderthal" taxonomic operational unit.
21Finally, to take into account the variability within the species H. sapiens, different OTUs have been constituted. The first, called "early sapiens" corresponds to fossil specimens of H. sapiens: Qafzeh 6, Qafzeh 9 and Skhül V. Then eight OTUs, composed of five recent individuals each, were formed in order to reflect the evolutionary history of human populations as indicated by the genetic diversification between different metapopulations (table S1) (Mounier and Mirazón Lahr, 2016).
22A total of 78 specimens were used in the study. We chose to run the cladistics analysis at the populational/species level, therefore the shape describing the specimens included in the study were averaged per populations/species (see "Modality of data processing") to form the 23 operational taxonomic units of the cladistic analyses.
23We estimated the position of missing landmarks, by mirroring the existing landmarks onto the other side and, when not possible, by estimating them by thin-plate-spline interpolation (i.e., TPS, Bookstein, 1989). Then, the 78 landmark configurations containing 636 landmarks describing our samples were aligned using GPA (Gower, 1975; Rohlf and Slice, 1990; Goodall, 1991) and bilateral asymmetry was removed from the data (Klingenberg et al., 2002). The 23 OTUs (including the three outgroups) were then computed as an average of the specimens forming each group described above and six different matrices were calculated. Matrices A.1 (table S2) and A.2 (table S3) were based on the analysis of 148 landmarks. Matrix A.1 used the PC coordinates for each OTU obtained after running a PCA on the Procrustes residuals resulting from the alignment of the 23 averaged configurations (table S2). Matrix A.2 was directly formed with the Procrustes-aligned 3D coordinates of the 148 landmarks averaged per OTU (table S3). Matrices B and C were based on 347 and 636 landmarks respectively. The same protocol used for the above configuration of 148 landmarks was applied to matrices B and C: B.1 (table S4) and C.1 (table S6) use the PC coordinates obtained after GPA and PCA, while B.2 (table S5) and C.2 (table S7) directly use the aligned 3D coordinates of the landmark describing each OTU. Variables presenting negative values cannot be processed by TNT software version 1.5 (Goloboff and Catalano, 2016). The protocol using the PC coordinates present negative values and to allow the software to process the data and perform the analysis, we added ten units to each variable of the three matrices. To verify that the order of the taxa did not impact the results several analyses were performed with different outgroups, changing the order of the taxa in the matrix.
24While the current debate on the use of geometric morphometric data in cladistics remains open, few technical means are available to test the approach implemented here. Among all the available cladistic software (e.g., PAUP*, Hennig68), only TNT of the Willi Hennig Society presents algorithms allowing the direct use of continuous (overlapping) quantitative data. These algorithms will attempt to reconstruct ancestral trait states by establishing the shape, estimated from landmarks, with the goal of parsimoniously minimising the morphological evolution of the tree based on optimisation logic (Catalano et al., 2010; Goloboff and Catalano, 2011; Ascarrunz et al., 2019; 2021).
25The approach used here is that of Goloboff and Catalano (2010; 2011). It operates directly on the coordinates of the PCs or of those of the landmarks, attempting to minimise the individual displacement of landmarks between ancestors and descendants. To do this, it will estimate the ancestral configurations at internal nodes in a manner analogous to the estimates of trait states in traditional parsimony, i.e., with discrete (non-overlapping) qualitative data. The ancestral patterns will be "spatial" and the coordinates of the landmarks will be optimised as a multivariate trait. The protocol aims to find the ancestral position of a landmark, by trying a limited set of positions determined by a grid encompassing the range of observed landmark positions in descendants. The ancestral position can be refined by increasing the number of grid subdivisions or by repeating the procedure using a smaller grid around the position of the first approximation.
26Therefore, six parsimony phylogenetic analyses were run, one for each of the protocols for both matrices. The use of new samples of landmarks for matrices B and C was driven by the need to best describe the morphology of the calvarium, while addressing possible technical problems, such as the limit of the number of analysable landmarks by the TNT software and a prohibitive computational time above a certain number of landmarks. Indeed, the parsimony analysis of the second protocol of matrix C could not be completed as the high number of variables (i.e., 636 landmarks) is beyond the computational capacity of the TNT software. In each case, 20 OTUs formed the ingroup (the genera Homo and Australopithecus) and three OTUs served as outgroups (the genera Gorilla, Pan and Pongo) to root the search and polarise the transformations of the variables. Phylogenetic analyses were performed with TNT version 1.5 parsimony software (Goloboff and Catalano, 2016), allowing the comparison of multiple hypotheses of relatedness and the quick selection of the most parsimonious. Parsimony has been chosen because it is the only phylogenetic reconstruction method that allows for the identification and discussion of the number and nature of synapomorphies observed at individual nodes (Matile et al., 1987).
27In the analyses presented here, characters were treated in an unordered manner and no transformation cost matrix was applied. The analyses of the matrices were performed with a heuristic algorithm, with random addition sequence (RAS), and tree bisection and reconnection (TBR) branch swapping, i.e., the tree is split into subtrees, which are then moved and reconnected to each other until the final tree, with the least number of steps, is found.
28The consistency index (CI) gives us an idea on the robustness of the trees with an indication of the degree of homoplasy present in each tree. However, if the number of autapomorphies is high in a tree, this index value will be high (tending towards 1), regardless of the homoplasy content of the other traits. To solve this issue, we also calculated the retention index (RI), which represents the ratio between the maximum number of homoplasies observed and the number of homoplasies that can be observed.
29In the case of the first protocol (analyses A.1, B.1 and C.1), a weighting of the variables, using the RIs of each PC, was applied before running a subsequent tree search, which used a non-exhaustive exact "Branch and Bound" algorithm, guarantying an optimal solution (Mounier and Caparros, 2015). The pre-variable RI gives information about the presence or absence of phylogenetic information, i.e., an RI=1 indicates the presence of a synapomorphy, on the opposite an RI=0 represents a homoplasic character (Farris, 1989).
30The robustness of each of the branches was finally evaluated via 100 bootstrap replications (following Felsenstein, 1985).
31The morphological variation present in the initial data sets (i.e., with configurations of 148, 347 and 636 landmarks) was first explored at the level of individual specimens using PCA analyses (table S8). In the resultant morphospaces, specimens occupy positions consistent with their usual taxonomic assignment in palaeoanthropology (tables S9-S11). When the average landmark configurations (i.e., used as variable for the OTUs) were used to compute the PCAs, the morphospaces described by PCs 1 to 3 indicate, as expected, an even stronger pattern of taxonomic groupings than what could be observed initially (tables S12-S14).
32The initial heuristic search on analyses A.1, B.1 and C.1 (i.e., composed of 22 PCs) resulted in several trees (figures 2-4; tables S15-S17). Out of the 22 PCs, three presented a RI superior to 0.35 for matrix A (figure 5), six PCs for matrix B (figure 6) and four PCs for matrix C (figure 7).
Figure 2
Parsimony phylogenetic trees based on a configuration of 148 landmarks: A) Matrix A.1 analysis of 22 PC coordinates resulting in a unique tree (2,630 steps, CI=0.569 and RI=0.806); B) Matrix A.2 analysis of the aligned 3D coordinates of the landmarks resulting in a unique tree (2,65227 steps and CI=0.468). Values shown above the branches represent bootstrap scores over 100 replications (node values below 50% are not shown). Red letters are used to name the nodes |
Arbres phylogénétiques les plus parcimonieux, basés sur l’analyse de la conformation de 148 points-repères : A) Arbre unique (2,630 pas, et IC=0,468) issu de l’analyse de la matrice A.1 (coordonnées des 22 CPs) ; B) Arbre unique (2,65227 pas, IC=0,468 et IR=0806) issu de l’analyse de la matrice A.2 (coordonnées alignées des points-repères) Les valeurs indiquées sur les branches représentent les scores de bootstrap après 100 réplications (les valeurs inférieures à 50 % ne sont pas montrées). Les lettres rouges correspondent aux noms des nœuds
Figure 3
Parsimony phylogenetic trees based on a configuration of 347 landmarks: A) Matrix B.1 analysis of 22 PC coordinates resulting in a unique tree (968,626 steps, CI=0.561 and RI=0.771); B) Matrix B.2 analysis of the aligned 3D coordinates of the landmarks resulting in a unique tree (2,70438 steps and CI=0.468). Values shown above the branches represent bootstrap scores over 100 replications (node values below 50% are not shown). Red letters are used to name the nodes |
Arbres phylogénétiques les plus parcimonieux, basés sur l’analyse de la conformation de 347 points-repères : A) Arbre unique (968,626 pas, IC=0,561 et IR=0,771) issu de l’analyse de la matrice B.1 (coordonnées des 22 CPs) ; B) Arbre unique (2,70438 pas, IC=0,468) issu de l’analyse de la matrice B.2 (coordonnées alignées des points-repères) Les valeurs indiquées sur les branches représentent les scores de bootstrap après 100 réplications (les valeurs inférieures à 50 % ne sont pas montrées). Les lettres rouges correspondent aux noms des nœuds
Figure 4
Parsimony phylogenetic tree (931,203 steps, CI=0.575 and RI=0.781) based on the analysis of a configuration of 636 landmarks. Only the analysis of matrix C.1 (22 PC coordinates) yielded results. The values shown above the branches represent bootstrap scores over 100 replications (node values below 50% are not shown). Red letters are used to name the nodes |
Arbre phylogénétique le plus parcimonieux (931,203 pas, IC=0,575 et IR=0,781), basé sur l’analyse de la conformation de 636 points-repères. Seule l’analyse de la matrice C.1 (coordonnées de 22 CP) a permis d’obtenir des résultats. Les valeurs indiquées sur les branches représentent les scores de bootstrap après 100 réplications (les valeurs inférieures à 50% ne sont pas montrées). Les lettres rouges correspondent aux noms des nœuds
Figure 5
Cranial deformations represented by PCs obtained from the analysis of 148 landmarks. Only PCs with RI≥0.35 are shown. Values shown in percent represent variance |
Déformations crâniennes représentées par les composantes principales (CP) obtenues après l’analyse de 148 points-repères. Seules les CP présentant une valeur de IR≥0.35 sont représentées. La variance de chaque CP est representée en pourcentage
Figure 6
Cranial deformations represented by PCs obtained from the analysis of 347 landmarks. Only PCs with RI≥0.35 are shown. Values shown in percent represent variance |
Déformations crâniennes représentées par les composantes principales (CP) obtenues après l’analyse de 347 points-repères. Seules les CP présentant une valeur de IR≥0.35 sont représentées. La variance de chaque CP est representée en pourcentage
Figure 7
Cranial deformations represented by PCs obtained from the analysis of 636 landmarks. Only PCs with RI≥0.35 are shown. Values shown in percent represent variance / Déformations crâniennes représentées par les composantes principales (CP) obtenues après l’analyse de 636 points-repères. Seules les CP présentant une valeur de IR≥0.35 sont représentées. La variance de chaque CP est representée en pourcentage
33For analysis A.1, heuristic analysis at 200 replications – i.e., the construction of the Wagner tree (the most parsimonious tree) and the application of TBR branch-swapping are replicated 200 times independently – produces two equally parsimonious trees that differ in the relationships between the modern human populations and in the positions of H. habilis, H. ergaster and Sinanthropus (table S15, node U of figure 5). Both trees are 2,630 steps long with a CI of 0.422 and a RI of 0.575 (table S15). After reweighting the 22 PCs using the RI values of each trait, a subsequent Branch-and-Bound analysis produces a single tree of 931,403 steps (CI=0.569 and RI=0.806) (figure 2A). The heuristic analysis of matrix B.1 results in three equally parsimonious trees of 2.456 steps, with a CI of 0.437 and a RI of 0.576 (table S16). After reweighting the 22 characters using the RI values and submitting them to a Branch-and-Bound analysis, we obtain a unique tree of 968.626 steps (CI=0.561 and RI=0.771) (figure 3A). Heuristic analysis of matrix C.1 results in two equally parsimonious trees of 2,384 steps with a CI of 0.443 and a RI of 0.57 (table S17). After reweighting the 22 traits using the RI values, a Branch-and-Bound analysis generates a unique tree of 915,513 steps (CI=0.575 and RI=0.781, figure 4).
34Analysis A.2, corresponding to 148 landmarks, results in a unique tree of 2.65227 steps (figure 2B), with a CI of 0.468. Nevertheless, the RI could not be computed in TNT and is, therefore, unknown. Analysis B.2, i.e., 347 landmarks, results in a unique tree of 2.70438 steps (figure 3B) with a CI of 0.468. Finally, the heuristic analysis C.2, corresponding to the 636 landmarks, could not be processed by TNT (see Methods). While matrices with several million characters, such as genomic sequences, have been processed successfully before (Torres et al., 2022), the use of landmarks has been found to be currently limited to 541.
35The unique most parsimonious tree resulting from the two-step analysis of matrix A.1 tree (figure 2A) shows that, at node D, the shift towards a value of -0.1 in PC1 explains the formation of the node. Node H supports a monophyletic group composed of the series: Sinanthropus, Sangiran and H. georgicus. Sangiran and H. georgicus constitute a terminal clade at node I. At node H, it is the values greater than -0.042 in PC1 and less than -0.012 in PC2 that explain the formation of the node. The Neanderthals appear as a paraphyletic group, but the Near East Neanderthal and the earliest early Neanderthal forms are grouped into one clade. At node M, the shift to a value of 0.02 in PC1 explains the formation of the node. At node N – i.e., the OTUs belonging to H. sapiens ("early H. sapiens" and the various modern human populations today) –, it is the shift to values greater than 0.025 in PC1 and -0.013 in PC2 that explain the formation of the node.
36In the tree from the A.2-matrix analysis (figure 2B), the monophyletic group including Sinanthropus and Sangiran/ H. georgicus at nodes H and I is supported by bootstrap values of 95% and 100%, respectively. The monophyletic H. sapiens group appears at node N, with a bootstrap score of 100%. In contrast, Neanderthals also appear as a paraphyletic group.
37The unique most parsimonious tree resulting from the two-step analysis of matrix B.1 tree (figure 3A) shows that, at node D, the shift towards a value of -0.1 in PC1 explains the formation of the node. Contrary to the tree obtained from matrix A.1, the one obtained with matrix B.1 shows that, at node E, we observe that the taxon A. africanus forms a clade with the taxon H. habilis sensu lato, supported by a bootstrap of 53%. It is a shift towards a value of 0.25 in PC2 that explains the formation of this node. Node H supports a monophyletic group composed of the series: Sinanthropus, Sangiran and H. georgicus. The latter two taxa constitute a terminal clade at node I. At node H and I, it is a value inferior to -0.04 in PC3 and inferior to -0.005 in PC4, respectively, that explains the formation of the nodes. The Neanderthals appear as a paraphyletic group. At node N, i.e., the OTUs belonging to H. sapiens ("early H. sapiens" and the various modern human populations today), it is the shift to values greater than 0.04 in PC1 and -0.025 in PC2 that explains the formation of the node.
38In the tree from the B.2-matrix analysis (figure 3B), the monophyletic group including Sinanthropus and Sangiran/ H. georgicus at nodes H and I is supported by bootstrap values of 100% in each case. The monophyletic H. sapiens group appears at node N, with a bootstrap score of 100%. In contrast, Neanderthals also appear as a paraphyletic group. During the bootstrap analyse, the formation of a clade A. africanus/H. habilis, with a bootstrap of 55%, can be observed.
39As with the tree obtained from matrix B.1, the one obtained with matrix C.1 (figure 4) also shows that, at node E, the clade composed with the taxon A. africanus and H. habilis sensu lato, is supported by a bootstrap of 51%. It is the shift to a value of 0.029 in PC2 which explains the formation of the A. africanus and H. habilis clade, while the other taxa keep a value less than 0.029. At node D, the shift towards a value of -0.05 in PC1 explains the formation of the node. On the other hand, we find the same Sinanthropus and Sangiran/H. georgicus clade at nodes H and I, supported respectively by bootstraps of 56% and 76%. At node H, it is a value smaller than -0.03 on PC3 that explains the formation of the node. The Neanderthals appear as a paraphyletic group. At node N, i.e., the OTUs belonging to H. sapiens, it is the shift to a value of 0.03 in PC1 which explains the formation of the node.
40Regarding the morphological changes associated with the nodes of the tree obtained from matrix A.1 (figures 2A and 5), at node D, it is the shift towards a value of -0.1 on PC1 that explains the formation of the node. This translates into various anatomical changes: the cranial outline in norma occipitalis, the definition of which depends on the position of the greatest cranial width (bi-euryon), is of triangular shape ("tent-like") (Grimaud, 1982); the post-orbital constriction is less marked compared to the one observed in apes; in norma verticalis, the profile of the supraorbital area goes from concave in its medial part at the level of the glabella to rectilinear; in norma lateralis, the profile of the antero-posterior frontal bone is more convex than the one of apes; and in norma lateralis, the transition between the planum occipital and the planum nucale is rounded, compared to the angular one of apes. At node H and I, it is values greater than -0.042 on PC1 and less than -0.012 on PC2 which explain the formation of the node regrouping H. erectus (Sinanthropus and Sangiran) and H. georgicus. This shows that the cranial outline in norma occipitalis is of triangular shape ("tent-like") (Gri-maud, 1982) compared to that of Neanderthals, for whom it is circular ("bomb-like" shaped) (Boule, 1911, 1912, 1913; Vandermeersch, 1981; Condemi, 1992); in norma verticalis, the profile of the supraorbital area goes from rectilinear for Sinanthropus to concave for Sangiran and H. georgicus; and in norma lateralis, the angle between the planum occipital and the planum nucale is more salient, compared to the other OTUs; we can also observe the fusion of the supraciliary arches in the medial portion of the supraorbital area. In this tree, the Neanderthals appear as a paraphyletic group, with the particularity that the Near East Neanderthal and the early Neanderthal are grouped into one clade. At node M, the shift to a value of 0.02 on PC1 explains the formation of the node. This shows that the cranial outline in norma occipitalis is more similar to that of H. sapiens, for whom it is pentagonal ("house-like") (Broca, 1868; Olivier, 1960), than to that of the "classic Neanderthal" OTU which presents a circular shape ("bomb-like") (Boule, 1911; 1912; 1913; Vandermeersch, 1981; Condemi, 1992); in norma verticalis, the profile of the supra-orbital region is convex and does not present the concavity in its medial part at the level of the glabella; in norma lateralis, the convexity of the antero-posterior frontal bone is more convex than for the "classic Neanderthal" OTU, but not as much as that of H. sapiens; and in norma lateralis, the angle formed by the planum occipital and the planum nucale is almost absent. The grouping of early and Near East Neanderthal could correspond to morphological similarity between the oldest Neanderthals whose morphology may not yet show the complete Neanderthal morphological pattern and others that show a slightly different morphotype, closer to the complete Neanderthal morphological pattern (Hublin, 1998). The N-node supports the grouping of all H. sapiens: "early H. sapiens" and the various modern human populations today. At this node, it is the shift to values greater than 0.025 on PC1 and -0.013 on PC2 that explain the formation of the node. This means that the cranial outline in norma occipitalis is pentagonal shape ("house-like") (Broca, 1868; Olivier, 1960); in norma verticalis, the profile of supra-orbital region is convex; in norma lateralis, the convexity of the antero-posterior frontal bone is very pronounced; and in norma lateralis, the angle between the planum occipital and the planum nucale is totally absent.
41As for the morphological traits of the tree of matrix B.1 (figures 3A and 6), at node D, it is the shift towards a value of -0.1 in PC1 that explains the formation of this node. The morphological changes indicated by this value are, for the most, similar to those observed on the matrix A.1 tree: the cranial outline will present the same shape; the post-orbital constriction is also less marked compared to the one observed in apes; in norma verticalis, the profile of the supraorbital area goes from concave to rectilinear in its medial part at the level of the glabella; in norma lateralis, the convexity of the antero-posterior frontal bone is more convex than the one of the great apes; and in norma lateralis, the angle between the planum occipital and the planum nucale is similar to the one observed in the matrix A.1 tree. At node E, it is the shift to a value of 0.25 on PC2, which explains the formation of the A. africanus and H. habilis clade. This indicates that the cranial outline in norma occipitalis presents a more accentuated triangular shape ("tent-like") (Grimaud, 1982) than the rest of the OTUs; the post-orbital constriction is important compared to the other specimens; in norma verticalis, the profile of the supraorbital area goes from rectilinear for H. habilis to convex for A. africanus; and in norma lateralis, the angle between the planum occipital and the planum nucale is more rounded, compared to the strong angle observed in H. ergaster. At node H and I, it is values greater than -0.04 on PC3 and less than -0.005 on PC4 respectively, which explain the formation of the node regrouping H. erectus (Sinanthropus and Sangiran) and H. georgicus. In this tree, the Neanderthals appear as a paraphyletic group, with the disappearance of the Near East and early Neanderthal clade. As for node N, supporting the grouping of all H. sapiens, i.e., "early H. sapiens" and the various modern human populations, it is the shift to values greater than 0.04 on PC1 and -0.025 on PC2 that explain the formation of the node. The morphological characteristics of this node are very similar to those from the preceding trees. We can observe the same pentagonal shape for the cranial outline in norma occipitalis; in norma verticalis, the profile of the supraorbital area will also be convex; in norma lateralis, the convexity of the antero-posterior frontal and the angulation between the planum occipital and the planum nucale are identical to those formerly observed in tree A.1.
42As mentioned above, the tree obtained from matrix C.1 is very similar to those of matrices A.1 and A.2. In the tree of matrix C.1 (figures 4 and 7), except for some minor changes in the H. sapiens clade, there is only the presence of the node A. africanus/H. habilis (as in matrix B.1, figure 3A) and the absence of the clade Near East/early Neanderthal (in contrast with matrix A.1, figure 2A). At node D, it is the shift to a value of about -0.05 on PC1 that explains the formation of the node. The distribution of morphological characteristics is similar to those of the other trees. At node E, the shift to a value of 0.29 on PC2 explains the formation of the A. africanus and H. habilis clade. This references the same characteristics as those of node E in the matrix B.1 tree: the cranial outline in norma occipitalis is in a more accentuated triangular shape ("tent-like") (Grimaud, 1982) than in the rest of the OTUs; the post-orbital constriction is important compared to the other specimens; in norma verticalis, the profile of the supraorbital area goes from rectilinear for H. habilis to convex for A. africanus; and in norma lateralis, the transition between the planum occipital and the planum nucale is more rounded, compared to the angle in H. ergaster. At nodes H and I, it is a value inferior to -0.03 on PC3, which explains the formation of the node regrouping H. erectus (Sinanthropus and Sangiran) and H. georgicus. The morphological traits for those two nodes are similar to those described for the former trees including the same clade. In this tree, the Neanderthals still appears as a paraphyletic group, but, as in the B.1-matrix tree, the Near East Neanderthal and the early Neanderthal are not grouped into one clade. In node N, i.e., supporting the grouping of all H. sapiens, it is the shift to a value of 0.03 on PC1 which explains the formation of the node, with morphological characteristics similar to those allocated to H. sapiens (e.g., a globular braincase; a short and high cranium; a long and high parietal arch; a high frontal arch; a rounder frontal bone; etc.).
43In this study, the CI and RI indices of the parsimony trees in matrices A.1, B.1 and C.1 are relatively high and reflect a robust phylogenetic signal. However, the RI for the matrix A.2 tree could not be obtained despite several attempts, and the literature does not explain this either.
44The main objective of this study was to test different protocols to combine 3D geometric morphometrics and cladistics. We ran several analyses using different types of data: coordinates of shape PCs and 3D coordinates of aligned landmark configurations.
45Such methodological approaches raise questions regarding inherent technical constraints in carrying out these analyses. First, these issues will be addressed, before focusing on the palaeoanthropological implications raised by the results of the present study, including the impact that the initial step of grouping the specimens into morphometric combinations designed to form OTUs may have had on the interpretation of the phylogenetic hypotheses that have been formulated.
46The methodology also presents biases related to several constraints that appeared in the processing of the information and the exploitation of morphometric units taken as OTUs of different taxonomic ranks. For analysis 1 (A.1, B.1 and C.1), the reduction of dimensionality through a PCA transformed the original variables into independent, orthogonal variables, at least from a mathematical point of view, thus, addressing possible problems linked to trait dependency generally present in the case of 3D geometric morphometric data. This issue, nevertheless, remains for the second set of analyses (A.2, B.2 and C.2) which use the correlated untransformed morphological contents as variables (i.e., 3D landmark coordinates). Comparisons of the results obtained from both approaches bring insight into the importance of the definition and the independence of the characters, especially when they can be compared to results obtained through "classical" numerical taxonomy. Indeed, the phylogenies obtained from both cladistic approaches are remarkably congruent and, except for a few differences such as the grouping of A. africanus and H. habilis, the obtained phylogenies are very similar. Moreover, the phenetic approaches which were tested on the PC coordinate database show a much higher degree of instability, producing a phenetic tree similar to the obtained phylogenies when 148 landmarks are used to describe the calvaria, while the phenetic trees from 347 and 636 landmarks produce very different typologies from the phylogenies (tables S19, S20). Therefore, and despite the possible theoretical issues related to using 3D geometric morphometric data in a cladistic framework, the obtained results in the present study, clearly indicate that the cladistic analysis better handles large 3D databases than classic numerical taxonomic approaches.
47Regarding the phylogenetic and palaeoanthropological content of the results obtained, firstly, the relative position of the OTUs in the different phylogenetic trees computed from five different matrices (A1, A2, B1, B2 and C1) are both congruent with the established phylogenies/chronologies and with the grouping observed in the morphospaces resulting from the PCA analyses (tables S9-S14), including the position of the outgroups. We can also notice the monophyly of H. sapiens, with slight differences in the relationships among modern human populations.
48Hence, our results contribute to the discussion regarding long standing debates in palaeoanthropology. For instance, trees B.1 and C.1, computed respectively from the PC coordinates derived from the analyses of 347 and 636 landmark datasets, show the formation of a clade A. africanus/H. habilis sensu lato, supported by bootstrap values of 55% and 51%. This result refers to the debate regarding the validity of the species H. habilis and its belonging to the genus Homo (Stringer, 1986; Prat, 1997). As we explained before, the two specimens forming the present study’s OTU of early Homo, KNM-ER 1470 and 1813, have seen their taxonomic attribution challenged over the years. While a number of studies are discussing the possibility that these specimens belong to different species of the genus Homo (i.e., H. rudolfensis and H. habilis (see Lieberman et al. 1996; Zeitoun, 2000; Argue, 2017), others question the attribution of these two fossils to the genus Homo to begin with and would advocate for their inclusion into the genus Australopithecus (Walker and Leakey, 1978; Wood and Collard, 1999). Therefore, the formation of a monophyletic clade A. africanus/H. habilis could be interpreted as support for the inclusion of the specimens attributed to H. habilis into the genus Australopithecus (Walker and Leakey, 1978). Alternatively, our results are based on the analysis of an OTU associating two specimens from the same site, Koobi Fora, but which may belong to different species (e.g., Alexeev, 1978; 1986; Walker and Leakey, 1978; Stringer, 1986; Chamberlain and Wood, 1987; Lieberman et al., 1988; 1996; Wood, 1991; 1992; Rightmire, 1993; Ferguson, 1995; Donnelly, 1996; Prat, 1997; 2004; Zeitoun, 2000). This could explain the apparition of this node in the trees of our analyses and raise the question: would this clade appear if we perform the same analyses exclusively using each specimen as representative of an OTU? The tree obtained with matrix A.1 (figure 4A) does not present the monophyletic clade associating H. habilis and A. africanus and would support the more traditional view of the inclusion of the H. habilis within the genus Homo (Prat, 2022). Another possible explanation for the formation of an A. africanus/H. habilis clade could be related to the anatomical region of the skull analysed here. Indeed, numerous morphological traits used to define A. africanus and H. habilis are based on facial and mandibular features (Rightmire, 1993; Prat, 1997; 2004; 2022), while the present study focuses exclusively on the calvarium, which is known to present less derived morphologies (e.g., Johanson et al., 1987; Rightmire, 1993; Prat, 1997; 2004; 2022), and could possibly explain the discrepancies in our results.
49The constitution of a monophyletic group composed by the OTUs Sinanthropus-Sangiran-H. georgicus, with Sangiran and H. georgicus forming a terminal clade, also contributes to the current debates in palaeoanthropology. The occurrence of these two clades refers to the still unresolved argument about the definition of H. erectus as one or more species. Our results appear to support a possible exclusion of the Ngandong specimens from a H. erectus clade which would include all the erectus-like Eurasian OTUs (specimens from Dmanisi, Zoukoudian and Sangiran) to the exclusion of the African erectus-like OTUs (specimens from Koobi Fora). The above results, unlike the clade A. africanus/H. habilis, can be observed in all the trees with relatively high bootstrap values (from 51% to 55%) and are therefore better supported. It is important to note that the question of the attribution of the Ngandong sample (possibly together with the specimens from Sambungmacan) to the species H. erectus remains open to debate, with some authors classifying it as a distinct species: i.e., H. soloensis (Zeitoun, 2000; Durband, 2007; Tattersall and Schwartz, 2009; Zeitoun et al., 2010). Regarding the Dmanisi sample, it is interesting to note that they cluster with populations recognised as H. erectus sensu stricto and are separate from African H. erectus or H. ergaster. These results are inconsistent with previous phenetic and phylogenetic interpretations of the Georgian hominin sample which underlined their morphological similarities with sub-contemporary African fossils designated as H. ergaster (Rightmire et al., 2006; Lordkipanidze et al., 2013). In this new scenario, the Dmanisi, Sangiran and Zhoukoudian fossils could be grouped into a unique H. erectus sensu stricto species, different from both the African H. ergaster and the late Asian H. erectus such as Ngandong.
50While the debate regarding the taxonomic position of H. neanderthalensis and H. sapiens has been somewhat fuelled by the recent genomic data indicating the existence of admixture between the two species (Green et al., 2010; Fu et al., 2015; Prüfer et al., 2021), most palaeoanthropologists now agree on the existence of two palaeoanthropologically valid species (e. g. Rak, 1993; Hublin, 2000; Harvati et al., 2004; White et al., 2014; Hajdinjak et al., 2018; Ni et al., 2021). However, the present study’s results regarding H. neanderthalensis could question, yet again, the taxonomic status of this species: in all of our results, the Neanderthal and H. sapiens OTUs are not distinguished as two separate clades, suggesting that they could belong to one group. This phylogenetic result is in disagreement with recent cladistic studies based on discrete morphological features (Mounier and Caparros, 2015; Mounier and Mirazón Lahr, 2016) but is supported by numerous studies where Neanderthals and modern humans are considered variants of a same single species (e.g., Smith and Trinkaus, 1991; Trinkaus, 1991; Relethford, 2001; Ahern et al., 2002; 2004; Curnoe and Thorne, 2003).
51There are a range of available interpretations that could explain such a result. First, the morphological differences between western and eastern populations of Neanderthals, as well as their relationship with H. sapiens, could be the consequences of a speciation by circular overlap or "ring species" (Voisin, 2006), i.e., "ring species provide dramatic evidence that normal genetic divergence within one species can build up to a sufficient level to generate two species" (Ridley, 2004:388). This would mean that "from Western Europe to the Near East, there was a succession of human populations that developed, over time, Neanderthal characters that became more and more marked from east to west; and in Western European Neanderthal populations, differentiation reached a level that did not allow interbreeding with modern humans. […] The farther those populations lived to the west, the more they displayed pronounced Neanderthal characters." (Voisin, 2006:311-318). However, it seems now unlikely that the classic Neanderthal populations would not have been able to interbreed with H. sapiens (Ahern et al., 2004; White et al., 2014; Prüfer et al., 2017; 2021).
52Finally, the topology of our tree is consistent in clustering first the classic Neanderthals, followed by the early and Near East Neanderthals. This could be due to the more derived Neanderthal morphology observed in the late Western European Neanderthals (Dean et al., 1998; Condemi, 2001; Profico et al., 2023), while a less derived morphology would explain both the formation of the Near East and early Neanderthal clade in the first phylogenetic tree (i.e., Analysis A.1) and a possible stronger phylogenetic relationship with early H. sapiens.
53An unequivocal interpretation of the results regarding H. neanderthalensis remains difficult at best, we cannot rule out that the use of different datasets could have had an impact on the building of the different phylogenetic hypotheses in this particular case. Additionally, most studies considering the phylogenetic and/or taxonomic status of H. neanderthalensis include other parts of the skull, such as the upper face and the mandible, where numerous specific Neanderthal features have been described (Hublin, 1998), and do not focus solely on the calvarium. Thus, we have to keep in mind that this region presents only a part of the total of morphological characteristics ("total morphological pattern") attributed to Neanderthals (e.g., Wood and Richmond, 2000; Smith et al., 2005; Tattersall and Schwartz, 2006).
54Our results, while questioning the status of Neanderthals as they are grouped within a paraphyletic sequence, also raise the question of the status of the Ngandong specimens (see node J of figures 2-4), whose position in the phylogenies are a reminder of the proposed H. sapiens soloensis trinomial by Dubois (1940) and later taken up by several authors (Dubois, 1940; Campbell, 1963; Jelinek, 1981; Tobias, 1985; Stringer, 1987; Bonde, 1989; Bräuer and Mbua, 1992; Widianto and Zeitoun, 2003).
55Regarding modern humans, in our results, they always form a monophyletic group, despite ambiguous relationships within the different populations of the species. It is important to note that H. sapiens are split in different OTUs representing relatively cohesive populations. This increases the morphological variation present in the database providing additional possibilities for the analysis to build phylogenetic hypotheses. In such a context, obtaining coherent and converging phylogenies, as is the case for the cladistic approaches developed in this study, advocates in favour of the reliability of the tested approach.
56Finally, cladistics is concerned with the structure of phylogenetic relationships, while morphometry is a tool designed to describe morphology and shape which can be used in a phylogenetic framework (e.g., the reconstruction of ancestral shapes at nodes, Mounier and Mirazón Lahr, 2016; 2019).
57To conclude on this point of reflection between interpretation and discussion on processes and taxonomy in the context of a phylogenetic analysis, let us recall that Tassy (1991) specifies that if the concept of phylogeny is conceived as the history of ruptures in the gene "pools" and the appearance of new "pools" isolated between populations, then we should indeed speak of tokogeny (Wheeler and Meier, 2000). This term describes infraspecific relationships according to Hennig (1966), who adds that relationships between terminal taxa are not necessarily hierarchical. Thus, if the aim of phylogeny is to put forward hypotheses of relationships between units seen as recognisable, isolated and closed sets, these units may be species but not necessarily identical to the evolutionary units. They can be populations that may be reduced to a few individuals and not necessarily isolated from other populations of the species (Zeitoun, 2000).
58In the present paper, we applied two protocols using geometric morphometric data in a cladistic framework on three different datasets, with the aim of testing the potential of such methodology. We obtained coherent results, describing possible phylogenetic and taxonomic hypotheses for the evolution of the genus Homo (i.e., a clade of A. africanus/H. habilis; a grouping of Sangiran/Sinanthropus/H. georgicus with the exclusion of Ngandong, paraphyly of Neanderthals and monophyly of H. sapiens). We are fully aware that caution is needed when it comes to interpreting such results. Many factors (i.e., the grouping of the OTUs, the anatomical region considered, the number and distribution of landmarks) could, indeed, contribute to shape the outcome of this kind of analysis and therefore the ability to understand the evolutionary mechanisms implied. Nevertheless, these results demonstrate the interest and potentiality of such an approach in palaeoanthropology, as they both confirm well-established ideas about the phylogeny of the genus Homo while bringing challenging ideas which can contribute to some of the most heated debates among paleoanthropologists today.
59The supplementary information is available in .xlsx from the BMSAP website (https://journals.openedition.org/bmsap/14080?file=1).
Acknowledgements: For permission to study specimens in their care, we thank directors and curators of the following institutions: Duckworth Laboratory, Cambridge; Distong National Museum of Natural History, Pretoria; Institut de Paléontologie Humaine, Paris; Musée de l’Homme, Paris; Museo di Antropologia, Sapienza Università di Roma; Museo Pigorini, Rome; National Museums of Kenya, Nairobi; ORSA database, Penn Museum; Peabody Museum, Cambridge; Senkenberg Forschungsinstiut und Naturmeseum, Frankfurt; and American Museum of Natural History, New York.
We thank K. Gerdau for the language proofing, Y. Heuzé, the Editor, and two anonymous reviewers for valuable comments and criticisms of earlier drafts that contributed to the improvement of this paper.
Finally, the program TNT is freely available, thanks to a subsidy from the Willi Hennig Society.
Flux RSS
