Phylogenetic analysis of stemmed points from Patagonia: Shape change and morphospace evolution

This work is focused in the study of Patagonian lithic projectile points shape variation from a phylogenetic perspective pursuing three main aims: first, generate a model of projectile point shape diversification and morphospace evolution; second, estimate shape variation through time, and finally, assess the robustness of previous results using the same methods but in a larger sample with better spatial coverage. A previous work using geometric morphometric and cladistic methods suggested a pattern of general morphological diversification across Patagonia related, at least in part, to the spatial distance between cases, distinguishing two main clades in northern (43-45° S) and southern (50-52° S) Patagonia. In the present work to study this pattern in a more detailed level, a sample of ca. 1200 projectile points was used to obtain statistically different morphological classes performing unsupervised K-means searching. Shape characters were used to describe the different taxonomic units and to perform the phylogenetic analysis (through the Neighbor Joining and Maximun Parsimony methods) using as an ancestor the earliest point type known to the region (Fishtail point). The new results suggest that projectile points with longer and narrow blades and smaller stems evolved later in Patagonia and occupy a different sector of morphospace that could be related to the emergence of different technical systems, like the bow and arrow. However, these results do not support the previous ones of a projectile point diversification pattern mediated by spatial distance, maybe due to the reduction of contrast between the extreme north and south of Patagonia by the larger spatial coverage used in the present analysis.


Phylogenetic perspective of projectile point shape evolution
It was shown that phylogenetic reconstruction is a useful tool to generate hypotheses regarding tempo and mode of technological change (O'Brien et al. 2001;O'Brien & Lyman 2003;Lipo et al. 2005;O'Brien et al. 2005, Cardillo 2009García Rivero & O'Brien 2014;among others), under the assumption that culture conforms and evolutionary system with a hierarchy of genealogical units analogous to the genealogical hierarchy of organic evolution (Boyd et al. 1997). Due to cultural transmission processes, artefacts are able to evolve in lineages that can be documented by different phylogenetic methods. For example, approaches like maximum parsimony, distance-based, maximum likelihood and Bayesian statistics have been applied to explore hypothesis about the evolution of basketry (Jordan & Shennan 2009), tapestry motifs (Tehrani & Collard 2002), ceramics (Harmon et al. 2006;Pardo Gordó et al. 2018), lithics (Buchanan & Collard 2007;Cardillo 2009;Darwent & O'Brien 2006;Lycett 2009;Mesoudi & O'Brien 2008a;2008b) and linguistics (Atkinson et al. 2008;Gray & Atkinson 2003). While techniques of phylogenetic analysis are different -since some of them are related to specific hypotheses about rates of change-the basic principle for application is based on the observation that culture constitutes an independent system of inheritance (but in many cases related to genetic one) (Cavalli-Sforza & Feldman 1981;Durham 1992).
The main strength of the cladistic method is its dependence for phylogenetic reconstruction on homologous traits, separating them from those analogous (Kitching et al. 1998). By definition, homologous traits are those shared traits inherited from a common ancestor. Instead, analogous traits are those forms that evolved independently in unrelated lineages. In cladistic terms, the former constitute the so-called synapomorphies, or inherited similarities, whereas the latter are homoplasies (Kitching et al. 1998; see also Collard 2006;O'Brien & Lyman 2003). The purpose of cladistics is to build hypothetical evolutionary relationships among taxa by documenting the branching structure between these taxa, distributed in related groups called clades and forming a tree-like pattern. So, cladistics assumes that the evolutionary process takes place by the binary divergence between classes increasingly derived from a hypothetical ancestor. Thus the quantity of homoplasy (analogous variation) in a dataset informs the degree to which the evolution of certain taxa cannot be explained by the tree-like model of branching divergence. These homoplasies are the result of convergence and other evolutionary processes such as reversion and parallelism but not of inheritance. As homoplasy increases, the phylogenetic signal of a dataset declines. The strongest of the phylogenetic signal of a given dataset is assessed by different goodness of fit measures to the branching pattern.

Shape spaces and phylogenies
The form (size plus shape) of projectile points is usually used as a classification tool to characterize the variation of these artefacts in time and space (Beck 1998;Bettinger & Eerkens 1999;O'Brien et al. 2001;O'Brien & Lyman 2003;Okumura & Araujo 2014). Thus, morphological variations over time are often explained by changes in the strategies for obtaining resources (Hughes 1998;Ratto 2003;Restifo 2013, among others), or by the existence of particular traditions in design selection, without this necessarily having any functional implications.
In artefacts so directly related to subsistence through energy capture, like projectile points, it is expected that the morphological variability reflects (at least to a limited extent) functional restrictions (such as cutting capacity, penetration, etc.) as well as the interdependence between structural aspects (such as weight, raw materials, symmetry, hafting requirements), which will not be the same for the different technical systems (Ratto 1990;Hughes 1998;Shott 2011). In this scenario of correlation between the different structural and morphological factors, it is expected that the evolutionary trajectories preferentially follow some ones over others (Cardillo 2009). On the contrary, if these restrictions do not exist, it would be expected that all the potential morphological variation was made over time. Within a phylogenetic perspective like the one proposed in this paper, it is expected that in the absence of restrictions, the different clades have the same rate of diversification over evolutionary time when they are represented within the morphological space. Kindemberg (2010) refers to the first scenario as one of total or partial morphological restrictions, while the second would have no restrictions.
On the other hand, if there were restrictions, diversification would not occur in an equivalent way throughout the morphological space but would be greater around certain design combinations. This would be more in line with the existing evidence for a large number of technologies, as demonstrated by Basalla (1998). Neither these design spaces would be fixed, but it is expected that these types conform lineages (evolutionary trajectories of types connected by descent with modification) displaced over time in relation to changing functional requirements favoured by selection. Although in general it is expected, as mentioned above, that these changes are made in preferential directions or regions of the design space (Cardillo 2009).
One way to approach this phenomenon is from the study of morphological spaces (or morphospaces) as an approximation to the general design of the lithic projectile points. The morphological spaces are by their nature continuous and multidimensional, so it is common to generate them from multivariate methods (see for example, McGhee 1999). In the case they are estimated by a specific number of real cases, these spaces will be of an empirical nature (McGhee 1999) and its amplitude will be, at least in part, a function of the morphological variation present in the dataset. Another possibility is to generate theoretical spaces from morphologies defined by geometric functions. In this case the spaces are not determined by empirical variation and are especially useful for complex morphological variables that involve numerous dimensions (McGhee 1999;. In our case study, we will use empirical spaces generated from geometric morphometrics. Although empirical spaces have limitations, they will be more robust as the sample size increases. This allows representing the phylogeny within a morphological space and, in this way, generating a visual representation of the diversification path of the lithic projectile points over time. In this context, evolutionary change is represented, as Klinberg (2010) refers to as paths of ancestors to descendants within morphospace. The patterns of occupation and displacement of the ancestors-descendants of projectile points within the morphospace thus provide information about the evolutionary dynamics of shape. Some sectors of the total design space may present more restrictions than others, so various patterns are expected in relation to this (McGhee 1999;Sidlauskas 2008; see also Gould 1989 about constraints in evolution). Therefore, it is possible that some tree branches have more potential to generate new classes (greater diversity) than others, but with less global morphological variation. Within this scenario, certain clades will occupy more restricted spaces where less potential variation is feasible (for example, by functional constraints). In this case morphological channelling is expected and, therefore, an imbalance in the way in which this morphospace is occupied. Alternatively is possible that the morphological space is occupied homogeneously and all the sectors present the same probability of diversification. In this second scenario, the restrictions on the potential for morphological diversification are minor, irrespective of their diversity of classes ( Figure 1). This is in accordance with what was proposed by Sidlauskas (2008) about the evolution of the morphological space, in which there is a scenario where the lineages within the clades with high morphological diversity experienced a higher diversification rate per branch and a second scenario where the exchange rate is equivalent for all the clades but the greatest morphological diversity is linked to the exploration (mode) of new morphospace regions. This means that in one case the high diversity of classes is accompanied by a high morphological variation (disparity) and in the other both are decoupled.

Previous research on projectile point shape variation since a phylogenetic perspective
Patagonia is the southern tip of South America, covering a spatial scale of ca. 1500 km between 39º (Colorado river) and 52º (northern coast of strait of Magallanes) of South latitude. It is characterized by the presence of Andean Cordillera in the west and plateau and low plains in the east. The regional climate is conditioned by temperature gradient which decreases southward (Clapperton 1993). The earliest evidences for human occupations at different places across Patagonia is dated to ca. 10-12000 BP (Bird 1938;1946;Borrero 1994Borrero -1995Borrero & Franco 1997;Massone 1987;Miotti 1995;1996;Nami 1985Nami -1986Prieto 1991). These hunter-gatherer populations had a diet mainly centered in guanaco (Lama guanicoe) hunting and a lithic technology well-known by the presence of Fishtail points (Bird 1938;1946;Hermo & Terranova 2012;Hermo et al. 2015;Miotti 1995;1996;Massone 1987;Mengoni Goñalons 1987;Nami 1985Nami -1986Politis 1991, among others). From the Middle Holocene and especially during the Late Holocene stone points from Patagonia show a wide range of metric and morphological variation. Many of these changes are related with functional diversity, use-life, hafting techniques, and spatial and temporal variations, among others (Álvarez 2011;Banegas et al. 2014;Cardillo & Alberdi 2015;Charlin et al. 2013;Franco et al. 2005;Gómez Otero et al. 2009;Nami 1986;2003;Ratto 1990;.
From a phylogenetic framework, in a previous work we analyzed late Middle and Late Holocene stemmed points shape variation with the aim to explore how spatial dimension mediates on the process of point shape diversification . We studied a sample of 301 complete stone points recovered from continental Patagonia between 40º and 52º of South latitude, which was separated in six groups according to latitudinal strips. Through geometric morphometrics, mean shapes by strip were obtained, which were then used in cladistic analysis to model diversification trends.
These analyses showed a pattern of general morphological diversification related to the spatial distance between groups, showing a geographical gradient from north to south. Two large groups of morphologies with similar deformation patterns were distinguished in Northern (43º-45º S) and Southern Patagonia (50º-52º S). While point shape of higher latitudes showed a more uniform pattern, at middle and low latitudes a greater variability was observed.
These results suggested that variability in late Middle and Late Holocene point morphology could be explained by the occurrence of geographical (spatial) and historical macroscale-related mechanisms. The divergence into two large groups appeared as a phenomenon channelled by spatial distance and related to mobility and information flow among human populations since spatial model explained 79% of phylogenetic variability . This process was related with the Santa Cruz River (50º S) functioning as a biogeographic barrier, like the distribution of other lines of evidence had also suggested (Borrero 2001;Cardillo 2011;Charlin & Borrero 2012;Franco 2002;Orquera 1987).
Given the considerable environmental variability in Patagonia, the pattern we had observed might be also linked to ecological mechanisms since the Patagonian environment is highly conditioned by latitude (Clapperton 1993). Thus, it may be expected that point design was influenced by performance requirements in different environments. Pursuing this aim,  assessed the correlation between spatial and environmental variables (precipitation and temperature) and point morphological change, enlarging the sample up to 1445 stemmed points, including insular Patagonia (samples from the Isla Grande of Tierra del Fuego in southernmost Patagonia). A global trend for the distribution of shapes according to environment was not observed at this largest scale. Contrarily, the results showed a pattern of high morphological variation in lithic points in a local or micro-regional scale across overall Patagonia. This phenomenon is similar to that recorded by other lines of evidence, such as the distribution of raw materials Alberti & Fernández 2015;Borrazzo 2012;Charlin 2009;Cirigliano et al. 2018;Franco 2002), flake vs. blade technologies (Charlin et al. 2011;Franco 2008;Franco et al. 2016;, coastal technologies (Cardillo 2011), rock art (Charlin & Borrero 2012 and references therein) and diet breadth (Barberena 2002;Barberena et al. 2009;Borrero & Barberena 2006;Borrero et al. 2001), which suggest an increase in the regionalization of human populations, especially in the Late Holocene.

Aims
In order to study the Patagonian projectile points morphological diversification in a more detailed level, our aims here are: first, generate a quantitative model of Patagonian stemmed point evolution through phylogenetic reconstruction based on shape variation; second, estimate shape variation through time, and finally, assess the existence of correlation between the diversification pattern and the spatial distance between the classes, as previously research showed ).

Sample composition
The studied sample is composed by 1197 complete stemmed points from overall continental Patagonia (Figures 2 & 3). According to sample density, southern continental Patagonia, especially the Pali Aike volcanic field region (Santa Cruz Province, Argentina and Magallanes, Chile), is better represented than other northern areas since several previous works were focused there (Charlin &  The whole sample is composed by non-Fishtail stemmed points belonging to late Middle and Late Holocene and come from our research projects, an extensive survey of museum collections and published images taken from local literature.

Geometric morphometrics
Geometric morphometrics (GM hereafter) is the statistical analysis of form based on Cartesian landmark coordinates (Adams et al. 2013;Bookstein 1991;Mitteroecker & Gunz 2009;Slice 2007;Webster & Sheets 2010). In the last years this method has been increasingly applied to the study of stone tools form (  The application of GM methods on stone tools allow representing their physical configuration (their size and shape) as a mathematical object by means of Cartesian coordinates (Mitteroecker & Hutteger 2009). These methods also allow quantification of variation in size and shape as separate variables in the absence of allometry (Zelditch et al. 2004), which is a great advantage over traditional methods, since they usually studied the shape of artefacts through linear measurements that are in most of the cases highly correlated among them and with size (Bookstein 1991). Thus, these measures actually describe form (size+shape) rather than shape and provide largely redundant information (Iovită 2010;Shott & Trail 2010). In morphometrics, the term shape denotes the geometric properties of an object invariant to scale, position and orientation, whereas form comprises both its shape and size (Bookstein 1991;Mitteroecker & Gunz 2009;Mitteroecker & Huttegger 2009). Hence, through GM it is possible to perform multiple statistical analyses to assess projectile point size and shape change.
In GM the form of an object is captured by discrete points called landmarks and semilandmarks. The most important property of the former is their homology, either in light of a biological or geometrical principle (Bookstein 1991). When homologous points are difficult to identify, like in curved outlines (e.g., projectile point blades or end-scraper edges), semilandmarks are used (Bookstein 1997). They are arbitrary points defined in terms of its position on curves and surfaces used to capture homologous structures where isolated anatomical or geometric loci are not evident. The arbitrary spacing of semilandmarks can be controlled by sliding them following different procedures and algorithms (see Bookstein 1997;Gunz & Mitteroecker 2013;Gunz et al. 2005;Perez et al. 2006).
In landmark-based methods shape parameters are estimated by a Procrustes superimposition procedure (the Generalized Procrustes Analysis, GPA), that translates the original forms to a common origin, scales them to the same centroid size, and rotates them using a least-squares criterion (Rohlf 1990;Rohlf & Slice 1990). In this way, the GPA removes the effects of translation, rotation, and scaling, which results in shape coordinates free of variations in position, orientation, and size.
In this analysis we used 24 morphometric points located on the contour of the projectile points in order to achieve a good representation of their shape: they comprise seven landmarks located in homologous loci according to stemmed point design and 17 semilandmarks in projectile point curved sections (especially on blade outline) and in the middle-point between landmarks in shoulders and stem portions without distinguishable morpho-technical traits (Figure 3 A).
The arbitrary spacing of semilandmarks was removed by sliding them following the minimun bending energy criterion (Bookstein 1997).
Landmark configurations were superimposed performing a GPA using the tpsRelw (ver. 1.69) software (Rohlf 2017). After superimposition, pure shape information (named in general Procrustes or shape aligned coordinates) was obtained to be used in cladistic analyses.

Definition of morphological types
As many experimental and allometric studies have shown, artefact shape variation is a continuous phenomenon (Andrefsky 2006;Bettinger & Eerkens 1999;Bettinger et al. 1991;Bradbury & Carr 2003;Buchanan 2006;Buchanan & Collard 2010;Dibble 1984Dibble , 1987Flenniken & Raymond 1986;Flenniken & Wilke 1989;Hiscock 1994Hiscock , 2006Hiscock , 2007Hiscock & Attenbrow 2002, 2003Hiscock & Veth 1991;Hunzicker 2007;Morrow & Morrow 2002;Shott & Ballenger 2007;Towner & Warburton 1991). However, phylogenetic analysis require a set of types or classes described by characters (discrete or continuous) to build a tree. Therefore, to define different discrete entities (classes) in a continuous distribution could be a difficult task and, in some cases, a very subjective one. Moreover, in stone tools like arrowheads, discrete shape classes, if exist, are affected by replicative errors, mechanical differences in the knappable materials used, as well as by the life history. In fact, in a previous work  showed evidence that the morphological variation in stemmed projectile points from southern Patagonia was affected by rejuvenation processes (see also , which result -in evolutionary terms-in shape convergence, since they generate morphological equifinality among projectile points from different technical systems and chronological contexts . For all these reasons we considered classes as hypothetical operational taxonomic units (OTUs) in the sense of Dunnell (1989), which would carry information about life history. These kinds of units are useful in the cases in which obvious discontinuities are not directly observable, as often happens in continuous features. To derive OTUs, we perform unsupervised K-means group searching in the morphological space defined by all Procrustes aligned specimens (shape coordinates of projectile points). K-means is a simple automatic learning algorithm that is used to solve clustering problems. The goal of K-means clustering is to find groups in the dataset, where K represents the number of groups to be defined. The algorithm works iteratively to assign each data point to one of K groups based the observed similarities in the data matrix (in this case the shape matrix). Hence, we first perform a gap statistic to estimate the optimal number of groups (Charrad et al. 2014). Gap statistics is a method to find the gap in the continuous multivariate distribution which defines the minimum number of possible groups. To achieve a stable solution 500 bootstrap replications of the searching processes were made. A minimum of 1 and a maximum of 20 groups were set for K-means algorithm searches. Once the optimal number of clusters was found, 30 interactions of the K-means were allowed in order to accurately define group centroids and boundaries. After this, the mean difference between groups was tested by permutational MANOVA using 10000 bootstrap replications (Anderson 2001) at a significance level of α=0.05. In order to decrease the chance of committing type 1 error in-between-group comparison, p-values were adjusted with Bonferroni method at α/n, where n was the number of groups to being compared. More information about the analyses can be found in the attached R script file 1 and 2 of supplementary material.

Phylogenetic analyses
For phylogenetic reconstruction, the mean shape of each cluster was used as an operational unit to build a phylogenetic tree through two phylogenetic reconstruction methods: one based on distances and other on maximum parsimony. The Neighbor Joining method (NJ), which is a distance-based phylogenetic tool (Saitou & Nei 1987), was used in previous research by the authors  with overall good results. Although the method is similar to cluster analysis (in the sense that it uses total similarity as input of the clustering processes), NJ is considered a minimum evolution algorithm, because minimizes the total amount of change as well as other phylogenetic methods. The tree can also be polarized to indicate the direction of change. Moreover, experimental studies show that NJ is either effective in recovering the true phylogeny or in many cases is significantly closer to the actual tree (Atteson 1997;Gascuel & Steel 2006;Mihaescu et al. 2009). In this context homoplasies are related to additivity. For distances to fit into an NJ tree, they must achieve this condition (also called four point condition, Saitou & Nei 1987). When estimated distances generate lineages that go backwards, the method fails to produce a correct evolutionary tree.
The result of NJ reconstruction is only one fully resolved tree, while other methods such as Maximum Parsimony, can generate an indeterminate number of more parsimonious trees (see below). Bootstrap resampling is a common method to measure the uncertainly in tree reconstruction. In this case, the support of each clade was evaluated by bootstrap resampling 10000 times and subsequently the majority consensus was estimated. Majority consensus is represented onto a composite tree (Figure 4 B) that represents uncertainly related to the phylogenetic hypothesis about branching model of divergent evolution. The branches supported by less than 50% of the bootstrap interactions were collapsed and represented as unresolved. These unresolved ancestral relations occur when an internal node of a cladogram  Since the tree obtained is additive, the length of the branches implies an amount of evolutionary change. Both resolved and consensus tree were used to plot lineage to time trajectories (Nee et al. 1995), which are bivariate graphs that represent accumulation of species number against branching times. This plot depicts the pattern of diversification throughout the evolutionary time (which is defined from the distance between tree branches from root to tips). While this method function is mainly for exploratory analyses, their plots are useful to explore the relationship between diversification and extinction ratios, because concave exponential curves (linear under logarithmic transformation) is expected under constant diversification rates (Nee et al. 1995). In this case, we present the number of untransformed scale, since no significant differences were observed between the two methods and the raw frequencies are easier to interpret (see Figure 4).
However, one of the most known pitfalls of NJ is that character identity is loss due to the use of a "secondary" pairwise distance matrix. So, shared evolutionary novelties (synaphomorphies) or independent evolution (homoplasy) in a set of characters cannot be evaluated directly. On the contrary, with the Maximum Parsimony (MP) method the phylogenetic tree that minimizes the total number of character state changes is to be preferred, and perhaps for this reason is the most common method of phylogenetic reconstruction in general and the more widely applied in archaeological research (see section 1.2). Unlike distance-based methods, a result of a maximum parsimony search could be one or more most parsimonious trees with similar number of steps. Also, character mapping onto the tree allows reconstruct ancestral states, as well as estimate the degree of homoplasy in each one. For this reason, the main objective in using this method is to compare with the results obtained through previously used procedures (NJ) and to evaluate the degree of phylogenetic signal (as well as the homoplasy) contained in the shape. As we prefer a method that allows us to use metric and shape continuous characters as such, we select a method of phylogenetic reconstruction implemented for Catalano et al. (2010, also Goloboff & Catalano 2011 in the TNT program (Goloboff et al. 2008;Goloboff & Catalano 2016). Different recent research shows the good performance of parsimony reconstruction in the study of shape evolution using aligned landmark coordinates as characters (Catalano & Torres 2017;Catalano et al. 2010;Goloboff & Catalano 2011).
In this method, a set of landmarks is considered a configuration and is equivalent to one character. In this case the total set of landmark points could be considered the same configuration and used in the phylogenetic reconstruction. However, as was observed in experimental studies (Catalano et al. 2015;2017), the robustness of the results increases as the number of configurations increases. For this reason, two configurations or sets of landmarks were isolated for this analysis: one defining the blade and other for the stem (see File 6 of supplementary material). A previous work aimed at assessing the modularity of southern Patagonia projectile point designs has shown the blade and the stem function as independent modules (González-José & Charlin 2012). Even more important, the division of these two sections of the morphology allows us to evaluate the existence of a phylogenetic signal, as well as the degree of homoplasy in these configurations of shape independently. As mentioned above, we believe that this may allow to evaluate (in particular in the blade), the factor of convergence between classes, due in part to the life history of projectile points. To generate two separate modules, the landmarks number 9 and 17 were left aside, since they were in the middle of the two configurations (see Figure 3 A). In the implementation of phylogenetic method with sets of landmarks considered directly as characters, the homoplasy is calculated as the difference between the observed and the minimum possible displacements of the landmarks, and therefore the minimum possible displacement for each character (landmark configurations) is calculated (see also Klindemberg 2010).
For the search of MP trees we use 1000 independent heuristic searches of Wagner trees and successive rearrangements with tree bisection and reconnection (TBR), which try all possible re-connections between edges of a tree in order to reduce tree length. Also a maximum of 100 of the best trees was stored in each run, and suboptimal trees were discarded only when a new tree with smaller number of steps was found. In this case, reliability of the results has been evaluated by means of 1000 replications of symmetric resampling (see Goloboff et al. 2003). Then ancestral character states of blade and stem configuration was reconstructed and degree of homoplasy for each one was computed.
To polarize both NJ and MP trees and give direction to the branching pattern, the mean shape of Fishtail points from Southern Patagonia was selected as out-group.

Phylogenetic shape space
The resulting tree of NJ was used for phylogenetic morphospace reconstruction (Sidlauskas 2008). Since we start from a tree constructed from shape coordinates, its representation within morphospace (determined by the principal component resulting from the OTUs, see File 4 of supplementary material) serves to represent the evolutionary trajectories of morphological change within space in a total way. As we explain in the section 1.2, this allows us to visualize which clades were more diversified in shape and which were not. To assess the relationship between morphological diversification and spatial distance the average coordinates for each group was estimated and a distance matrix between them was generated. In parallel, the distance between branches of the tree without the out-group was calculated. Both distance matrices were correlated by means of the Mantel test. The p-value for the observed correlation was estimated by 10.000 permutations. More information about the analytic procedure can be found in the R script file 1 of supplementary material.
For OTU searching and phylogenetic analysis the package R 3.5.0 (R developed core team 2015) was used. Also, permutational MANOVA and the visualization of deformations by Thin-Plate-Spline were obtained with Past 2.14 (Hammer et al. 2001).
Supplementary material composed by an R script for main analysis steps (File 1), the aligned coordinates (File 2), mean shape for each group (OTUs, File 3), three main components of shape variation (File 4), mean spatial coordinates for each OTU (File 5) and TNT file for Maximum Parsimony searching with two shape configurations (File 6), are available online.

General shape trends
Twelve groups were recognized first with the gap statistics and by K-means clustering on aligned shape coordinates ( Figure 5). Permutational MANOVA with 10.000 bootstrap replicas and Bonferroni correction (F=280.8 p=0.0001) suggest that all shape means are statistically different at α=0.05 level. Between-pair comparisons are all significant at a corrected p-value of p=0.007.
First two principal components (PCs) showing main shape variations and the distribution of groups can be observed in Figure 5. They explain 80% of overall variation in projectile point shape. The first PC (PC1 58%) shows projectile points with shorter and wider blades and bigger stem areas in the positive scores whereas longer and narrow blades with smaller stems are located in the negative scores. The second axis (PC2 22%) shows wider blades with smaller stems in negative scores and narrow blades with longer stems in positive ones.
Some overlapping is observed in the 95% ellipses of the PC plot, but such overlap is partly related to the projection of the multidimensional group boundaries (where the search of the groups was carried out) onto a bidimensional space which we use here to represent the general variation trends.

Tree reconstruction
NJ gives one fully resolved tree (Figure 4 A) with bootstrap support values at the bottom of each node. Results suggest two big clades (if taking into account a basal split with 82% bootstrap support). Clades with elongate shapes (5 A, 1) have better resolution than clades mainly composed by wider and rounded blades (Figure 4 A, 2). These differences in support also can be seen in the majority consensus tree at the right of Figure 4 B. In the fully resolved tree we depicted the distance to the root on colour. More supported branches appear to be more distant to the root or more derived than the others, suggesting more evolutionary time. Two lineages through time plots were estimated for each tree (resolved "A" and consensus "B" tree in Figure 4). In the first case, the plot suggests high diversification at the beginning followed by a high rate averaging the mean evolutionary time. In the second case, as low support branches were collapsed (values lower of 50% bootstrap), high diversification is depicted only at the beginning of the tree with a sudden fall for the remained evolutionary time, suggesting lower diversification rates.
A single tree was also obtained by the MP method ( Figure 6). One interesting result is that this tree is very similar to that obtained by NJ (see Discussion). The only OTU in a different position is G5, which appears in an early node than in the NJ tree but in the same clade. Also group support is better than in NJ tree but in general those nodes with better support in the NJ tree also have high values in MP.
Likewise, the reconstruction of the configuration for blade and stem shows two main tendencies: one less derived closest to the out-group of robust forms, with more rounded blades with expanded stems (in particular stem neck) and another more derived towards more contracted stems and expanded blades with more lateral compressed stems, which also defines G4 that comprises barbed and expanded blade shapes, an evolutionary novelty not shared with G8 or other members in any clade (autopomorphy, see discussion).
Overall results and resampling supports suggest that blade shape carry with phylogenetic information. However the overall homoplasy for blade configuration is 3.29, or 2.42 times greater than the stem (1.36). This pattern supports the idea of more independent change in the blade configuration that could be related to different factors, as life history (see Discussion).

Morphospace occupation
Tree distribution onto shape space suggests gradual displacement of diversification from left to right (Figure 7). The two main clades occupied different areas in the shape space. Projectile points characterized by blade contraction (depicted in blue colours) and stem expansion (depicted in red colours) are related to earlier evolved shapes, while the ones with expanded blades and contracted stems are linked to more derived ones (more distant to the root) (Figure 7).
Projectile points characterized by blade contraction (i.e. shorter and wider blades depicted in blue colours) and stem expansion (i.e. bigger stems depicted in red colours) are related to earlier evolved shapes (like G1, G3, G6, G10, G11 and G12), while the ones with expanded blades (i.e. longer and narrow blades in red colours) and contracted stems (i.e. smaller ones in blue colours) are linked to more derived ones (more distant to the root, like G4, G5, G8 and G9 ) The phylogeny displacement within the morphological space does not show overlap or tangle between the branches, suggesting a gradual shift from one morphospace sector to another, where different traits of total shape evolves.
Finally, the Mantel correlation between cophenetic distance and mean spatial coordinates of each group yields not significant results (r= -0.25, p>0.05). These results contrast with previously observed patterns , where a correlation between both distances was observed, which may be due to different factors (see Discussion). (2018) Figure 4A). Branch colour shows distance to the root. Thin plate spline reconstructed shapes represent deformation between out-group and tips.

Discussion
The results obtained support previous hypothesis on the existence of a phylogenetic signal in the lithic projectile points from continental Patagonia. General pattern suggests that projectile points with smaller stems and elongated blades evolve later, which is consistent with the existing archaeological information, which relates this type of points as belonging to release systems by means of bows (Banegas et al. 2014;Bettinger & Eerkens 1999;Bird 1988;Ratto 1994;Shott 1993;. Increasing sample size and the statistical based definition of projectile point classes allowed us to model lineage diversity through time and explore its relation with the shape evolution. Indeed branching pattern of phylogeny within the shape space is consistent with this idea, and suggests that both sectors of the morphospace were occupied at different times and gradually filled, although in a relatively homogeneous way. We believe that changes in morphospace occupation are related with a functional dimension, showing the evolution of different weapon systems along the Holocene. As Lyman et al. (2008) observed, the increase of the total variation is compatible with the appearance of new technologies, where the design space (sensu Stankiewicz 2000) changes, whenever the functional requirements also changed. (2018)  On the other hand, the bootstrap consensus tree of NJ points out to the existence of homoplasy in some clades, in particular those closest to the tree root. It is feasible that this homoplasy is related to the fact that some of the OTUs belong to a morphological continuum linked to the allometric trajectories followed by the designs throughout their life history, as was suggested previously as a source for potential convergence in projectile point shapes .

Journal of Lithic Studies
It is also interesting to note that both methods result in similar topologies, although the support values for the tree obtained by MP are higher. These results are in agreement with what was observed by Catalano & Torres (2016) and we also observe in previous research (Cardillo & Charlin 2016: 265). The authors carried out comparative studies of phylogenetic reconstruction with different methods in 41 shape datasets and observed that MP and NJ generated trees with similar topologies in all cases (Catalano & Torres 2016).
This method also allowed reconstructing the expected morphologies for the different nodes of the tree. Although both configurations have relevant synapomorphies, estimation of overall homoplasy indicates that on the blade is 2.41 times higher than on the stem, a pattern also recognized by other methods . As previously mentioned, the case study of morphological change trajectories for three classes of stemmed points commonly differentiated in Late Holocene southern Patagonia (IV and V Fell-, Magallanesor Bird-types and Ona type sensu Bird 1938Bird , 1946Bird , 1988, indicated that there was a high potential for convergence throughout life history. In this way, although they differ in size and shape at the beginning of their lifespan, these point types tend to converge in shape as the rejuvenation process advances . For this reason, it is expected that some classes show different degrees of reactivation and allometric trajectories. Also is possible that these differences are due to modifications that occurred during the life history and that do not have a hereditary component while others do. This issue and the way to deal with are still in discussion from Gould (1985) "Ontogeny and Phylogeny" book to the present. However, take into account these factors (in particular heritable changes in development time or heterochrony) need further research. Independent measures of reduction could be obtained for each projectile point and then map the degree of reduction intensity detected on the resolved phylogeny. It is to be expected that if the use-life generates convergence or parallelism between classes, differences between each clade that are not present in the immediate common ancestor (and that can then be confused as an evolutionary novelty, an apomorphy) will be observed. Taking into account different configurations, however, can be used as a first step to estimate the degree of homoplasy in this character.
Considering the pattern of diversification observed through the reconstruction of the number of lineages in relation to the evolutionary time of the tree, two hypotheses can be proposed: a first scenario of initial cladogenesis followed by a relatively high rate of diversification and a later fall-off when averaging the evolutionary time (resolved tree Figure  4 A) or a rapid initial diversification followed by a sudden descent (consensus tree Figure 4  B). This second scenario regards the uncertainty in the resolution of the tree, linked to the potential convergence and homoplasy of the OTUs. Considering homoplasy or not, both patterns are consistent with what was previously observed, where the largest diversification event was located close to the root ). In relation to this possible pattern, it is interesting to note that it is probably the most commonly observed in biological phylogenies. Hughes et al. (2013) suggest that the greatest biological variety (disparity) tends to occur early in the process of diversification, possibly related to the emergence of key innovations that are then shared by all members of the descendant clade. In technological terms, it is compatible with the generation of a basic set of different shape features (like the combination of more rounded blades with expanded stems or sharp blades with contracted stems), which in turn, could play as structuring factors in the subsequent diversification process.
With respect to the relationship between space and phylogeny, the new results are not in agreement with the previously observed spatial pattern of two main clades in northern and southern Patagonia . We believe that one of the main causes is that the larger spatial coverage reduces the contrast between the extreme north and south of Patagonia, where the largest differences seem to be located, as other lines of evidence have also suggested. However, we believe that this model of isolation by distance, as it was previously defined (Cardillo & Charlin 2006), should be put into contrast again considering other sectors of space (such as Tierra del Fuego, in insular Patagonia), or using spatially based methods for the definition of morphological groups. These issues "shape" our upcoming research agenda.

Conclusions
The aim of this paper was to perform a detailed analysis of shape diversification pattern and morphospace occupation for stemmed projectile points of late Middle and Late Holocene of Argentinean Patagonia by means of statistical and cladistic methods.
Results suggest both NJ and MP methods successfully detect a phylogenetic signal in projectile point shape data. The observed pattern support previous results that projectile points with longer and narrow blades and smaller stems evolved later in Patagonia and occupy a different sector of morphospace that could be related to the emergence of different technical systems like bow and arrow. Also, basic elements of shape configurations seem to have evolved more or less quickly, what was later manifested as a relative reduction in the rate of morphological innovation of divergent lineages that occupied the shape design space in a relatively homogeneous way. However, the relationship of shape divergence with spatial distance is not more supported in actual results. Finally, results support the idea that the average blade shape could carry a clear phylogenetic signal although it is expected that it contains greater independent changes related to life history.
Future analyses need to explore in more detail the role of design, performance and life history in the evolution of different projectile point classes. This could be done including results of regression between allometric trajectories into phylogenetic models. (2018)