Detecting morphological gaps in tooth outlines of a Pachyrukhinae (Hegetotheriidae, Notoungulata) lineage: systematic and palaeobiogeographical significance of the records from Northwestern Argentina

gaps in tooth outlines of a (Hegetotheriidae, Notoungulata) lineage: ABSTRACT Pachyrukhinae (Hegetotheriidae, Notoungulata) is a highly frequent clade in the Late Miocene-Pliocene outcrops of southern South America. In Argentina, two genera have been recognized for this The euhypsodont dentition of these Neogene forms creates significant difficulties when cheek are described for systematic purposes. Tremacyllus has been scarcely studied in comparison with Paedotherium , and taxonomic analyses have interpreted diagnostic features as intraspecific variations and proposed the monospecific status of the genus. Given the discussion regarding the validity of Tremacyllus species and the fact that dental elements are the most abundant remains in the fossil record, we employed a quantitative framework provided by geometric morphometrics and multivariate statistics to discriminating intra- from interspecific variability by tooth outline. We analyzed a large sample of 82 specimens and two hypotheses were tested: 1) there are morphological gaps within the analyzed sample; and 2) morphology follows a pattern of geographical variation within the sample, suitable for recognition of species. We found that morphological variability is organized into two clusters. Morphological gaps are associated with geographical patterns in the P4 and upper premolars datasets. Based on the classification of the type specimens and supported cluster structure, we recognize Tremacyllus incipiens Rovereto, 1914 as a valid taxon, endemic from western outcrops of Northwestern Argentina. Segregation between northern and southern morphologies agrees with two different palaeo-phytogeographic provinces. This approach proved to be very effective to address intra- and interspecific variation and contribute to the knowledge of available techniques to assess morphological variation. entre les morphologies du nord et du sud correspond à deux provinces paléo-phytogéographiques différentes. Cette approche s’est avérée très efficace pour aborder la variation intra- et interspécifique et contribuer à la connaissance des techniques disponibles pour évaluer la variation morphologique.


INTRODUCTION
Pachyrukhinae is a clade represented by small herbivorous mammals; they were ecomorphologically similar to leporids and rodents, with postcranial features frequently associated with fossorial habits, probably to build burrows (Reguero et al. 2010;Elissamburu et al. 2011;Sostillo et al. 2018). Along with Hegetotheriinae, both integrate the family Hegetotheriidae Ameghino, 1894 (Typotheria, Notoungulata) and represent one of the most diverse and widely distributed groups among the South American native ungulates (Cerdeño & Bond 1998;Reguero et al. 2010;Seoane et al. 2017).
So far, the occurrences of the genus Tremacyllus include outcrops from Late Miocene to Late Pliocene of Cata marca, Mendoza, La Pampa, San Luis, Córdoba, La Rioja and Buenos Aires provinces, Argentina (Cerdeño & Bond 1998;Prado et al. 1998;Tauber 2005;Tomassini 2012;Bonini 2014;Garrido et al. 2014;Tauber et al. 2014;Armella et al. 2015Armella et al. , 2016Vera & Ercoli 2018;Sostillo et al. 2018) (Fig. 1). This taxon is mainly characterized by antorbital maxillary and postorbital processes less developed than in Paedotherium; presence of a reduced sagittal crest, unlike that of Paedotherium and more similar to that of Prosotherium Ameghino, 1897 and Hegetotheriinae; mastoid bullae enormously inflated, more than in Paedotherium, projecting well beyond occipital condyles; presence of post-incisive depressions that extend beyond the level of P2 (reaching the anterior margin of P3); adult dental formula with 1/2, 0/0, 3/3, 3/3; upper and lower premolars not molariform and lower premolars notably more imbricated (less evident in upper premolars); M3 shorter or equal to M2, P2 proportionally shorter than in Paedotherium; upper premolars elliptic in cross-section; shorter mandibular symphysis; and smaller in size than Paedotherium (Cerdeño & Bond 1998;Ercoli et al. 2017;Sostillo et al. 2018). Whereas these features make Tremacyllus a very distinctive, easily distinguishable taxon, with regard to Paedotherium, the recognition of taxonomically meaningful entities below the genus level is a daunting task (Armella et al. 2015;Vera & Ercoli 2018;Sostillo et al. 2018;Armella 2019). Indeed, this issue is more evident in the case of fragmentary remains, where the simplified tooth design increases the difficulty of recognizing ranges of variability among the supposed taxa.
Regarding the taxonomic background, Tremacyllus was coined by Ameghino (1891) based on Pachyrukhos impres sus, a synonym of Tremacyllus impressus (Ameghino, 1888), from Early Pliocene outcrops of Monte Hermoso (Tomassini et al. 2013), on the Atlantic Coast of Buenos Aires Province (Fig. 1B). It and, Tremacyllus diminutus (Ameghino, 1891), constitute the first mentions of the genus. Tremacyllus intermedius Rovereto, 1914 was also described from Monte Hermoso outcrops. Ameghino (1908) assigned remains to Tremacyllus chapalmalensis Ameghino, 1908 andTremacyllus novus Ameghino, 1908 from Chapadmalal, on the Atlantic Coast of Buenos Aires Province, increasing the number of taxa in these areas. It is important to note that most of these taxa were based mainly on fragmentary remains.
The taxonomic list is also long in other regions of Argentina. Tremacyllus subdiminutus Rovereto, 1914 was created based on a fragmentary palate from Huayquerías de San Carlos, Men-doza Province (Fig. 1B). In Northwestern Argentina (NWA), Tremacyllus incipiens Rovereto, 1914 andTremacyllus latifrons Rovereto, 1914 were established from "Estratos Araucanos" cropping out in the Santa María valley, Catamarca Province (Fig. 1C). Later, Riggs & Patterson (1939) recognized two morphotype sizes among specimens from NWA, and assigned the larger to T. latifrons (synonym of T. incipiens sensu Riggs & Patterson, 1939) and the smaller to T. diminutus. Rovereto (1914) and later Cerdeño & Bond (1998) recognized many similarities between Tremacyllus impressus and T. incipiens but also subtle differences in their craniodental morphology: T. incipiens is larger than T. impressus, shows less imbricated premolars, lacks the anterior projection on p2 and presents a more excavated palate and post-incisive depressions that extend further posteriorly. They performed a broad review of known Tremacyllus species and only considered two taxa to be valid: T. impressus (synonym of T. intermedius, T. novus, T. chapalmalensis, T. diminutus)  The most recent treatment of Tremacyllus is the study of Sostillo et al. (2018), in which a large sample of these taxa from Cerro Azul Formation (Late Miocene, La Pampa Province; Fig. 1B) was analyzed . These authors interpreted the variability observed in the diagnostic characters of T. impressus and T. incipiens as intraspecific variation and proposed their synonymy. Based on these results, T. impressus is currently the only species of the genus.
The fact that Tremacyllus remains were found in several Late Miocene-Pliocene outcrops along the central and northwestern regions of Argentina suggests an interesting pattern from a palaeobiogeographic point of view. In this sense, Zetti (1972; see also Seoane et al. 2017) proposed that there is a differential distribution pattern between Paedotherium and Tremacyl lus, with the first one being more abundant in the Pampean Region of Argentina, and Tremacyllus more frequent than Paedotherium in NWA. Cerdeño & Bond (1998) suggested an interesting geographic distinction at the species level, where T. impressus and T. incipiens seem to be predominant in the Pampean Region and NWA (Santa María valley), respectively.
In the last decade, several palaeontological expeditions to NWA outcrops have yielded a significant number of new fossil remains, preliminarily identified as Tremacyllus incipiens and Tremacyllus sp. Most of these specimens were recorded from outcrops of the Andalhuala Formation (Late Miocene-Early Pliocene, Catamarca and Tucumán provinces; Bonini 2014; Armella et al. 2015Armella et al. , 2016, the India Muerta Formation (Late Miocene, Tucumán Province; Armella 2019; Armella et al. 2019) and the Las Cañas Formation (Early Pliocene, Santiago del Estero Province; Armella 2019; Armella et al. 2019Armella et al. , 2020 (Fig. 1C). Armella M. A. et al. Given the availability of a large number of Tremacyllus specimens to explore statistical techniques and the fact that dental elements are the most abundant remains in the fossil record, we employed a quantitative approach focused on discriminating intra-from interspecific variability using tooth morphology. In this regard, the combination of geometric morphometrics based on anatomical contour and multivariate statistics constitutes a powerful tool that has been successfully used for the study of morphological variation (e.g. Cassini 2013;Ercoli et al. 2017;Scarano & Vera 2018;Vera & Ercoli 2018). On the other hand, given that natural populations exhibit some degree of variability in their attributes, the sole recognition of morphological differences seems to be ambiguous to delimit species, at least as far as Tremacyllus is concerned. We base our approache on the idea that it is necessary to identify morphological gaps or discontinuities separating groups associated with distributional patterns to delimit species (Darwin 1859;Simpson 1937;Mallet 1998;Zapata & Jiménez 2012;Milla Carmona et al. 2016). Under this approach, the validity of species as taxonomic entities lies in the concept of geographical variation of the morphology (i.e., the organization of morphological variation following a geographical pattern [Wilson & Brown 1953;Drooger 1954;Thorpe 1987;Wiens & Servedio 2000;Mallet 2007]).
In this contribution, we explored the clustering structure and taxonomic significance of morphological tooth variation of Tremacyllus within the quantitative framework provided by geometric morphometrics and multivariate statistics. In particular, we included new specimens collected in the Late Miocene-Pliocene outcrops from NWA. Two hypotheses were tested: 1) that there exist morphological gaps within the analyzed sample; and 2) that morphology follows a pattern of geographical variation among them. After thorough morphometric and statistical analysis, we present a new systematic proposal for the genus Tremacyllus and palaeobiogeographic affinities between the Miocene-Pliocene areas from Argentina. Moreover, these analyses contribute to the knowledge of available techniques to assess morphological variation.

Data source
The material for this study is based on 44 specimens of Tremacyllus from NWA housed in the Museo Arqueológico Provincial "Condor Huasi" Sección Paleontología, Belén, Catamarca, Argentina (MCH-P); Colección Paleontología Vertebrados Lillo, San Miguel de Tucumán, Tucumán, Argentina (PVL); Museo Municipal "Rincón de Atacama", Termas de Río Hondo, Santiago del Estero, Argentina (MPAT); and Field Museum of Natural History, Chicago (FMNH). All of them were carefully selected from an original pool of 51 specimens taking into account that upper and lower teeth series (premolar and/or molar) needed to be preserved and fully erupted (adult stages) for the proposed analyses. In addition, to study as many samples as possible, we include 38 published specimens (e.g. Ercoli et al. 2017;Vera & Ercoli 2018;Sostillo et al. 2018) from other regions such as La Pampa, Mendoza, and Buenos Aires provinces that are housed in several institutions (see below). A complete list of fossils, collection numbers and their geographical provenance is provided in the Appendix 1.
The tooth morphology of all the studied specimens was digitized in 2D images using photographs of ventral views of skulls and dorsal views of mandibles oriented along the alveolar plane (Fig. 2). To standardize the photographs, each specimen was positioned in the same orientation, and a scale was included to account for tooth size. The mid line of the enamel perimeter of each tooth in occlusal view (Fig. 2) was used to capture its outline (cementum pattern was used when necessary and possible, i.e., when the tooth margin was broken). Adopting this criterion allowed us to increase the sample size within a standardized and comparable framework. We analyze individually the outline of each tooth (upper and lower dentition) and, following a modular design (i.e., block or structural package), two different traits pertaining to the upper and lower tooth series were quantified: premolar and molar module morphologies, each one sampled using a different set of contours (Fig. 2 L a s C u e v a s R a n g e s H u a l f í n R a n g e s A c o n q u ij a R a n g e s anatomical contour analysis In order to explore shape variation within the sample, we use elliptic Fourier analysis (EFA) (Kuhl & Giardina 1982), a procedure that allows the quantification of geometry of closed contours. By this, a given number of sine and cosine functions (harmonics) are employed to reproduce a set of target contours. In this sense, the degree of fit depends upon the number of harmonics used, more complex contours requiring a higher number of harmonics for a more accurate description (Milla Carmona et al. 2016). Each harmonic has four coefficients, which are used as shape descriptors of the contours (Appendices 2-4). For the selection of the number of harmonics needed for proper description, the most geometrically complex contours of each set were reconstructed using different numbers of harmonics. The lowest number of harmonics that deviates on average less than 0.01% of the centroid size from the best possible shape (calculated using the highest possible number of harmonics) was used to describe the shapes of each set of contours (Appendix 2). To remove the rotation, translation, and scaling differences from the studied objects (normalization) and only analyze the shape variation among them (MacLeod 2009), we performed a generalized Procrustes analysis (GPA ;Gower 1975;Rohlf & Slice 1990;Frieß & Baylac 2003). This procedure was done prior to EFA and after placing reference points on the contours (Appendix 2). A complete list of Fourier coefficient harmonics and collection numbers is provided in Appendices 3 and 4.
A principal component analysis (PCA) was made using the generated variance-covariance matrix of coefficients of all harmonics (Appendices 3; 4). Resulting principal components represent different aspects of the geometric variation present in the set of studied contours (MacLeod 2012). Only principal components representing at least 1% of the original variation were retained. Scores of specimens on principal components were used as the input variables for all subsequent statistical analyses (Fig. 2).
Both GPA, EFA, and PCA were carried out in the R environment (R Core Team 2017, v. 3.4.3) using the package Momocs v. 1.2.9 (Bonhomme et al. 2014).

statistical analysis
Morphological groups within each data set were explored using the gap statistic method (Tibshirani et al. 2001;Charrad et al. 2014) (Fig. 2). This is a standard technique for estimating the optimal number of clusters (i.e., groups) in a data set by comparing the within-cluster dispersion to the expectation under a null reference distribution with no obvious clustering structure. To this, Euclidean distances were selected as a dissimilarity measure. Three different clustering algorithms were used to determine the optimal number of groups: gap statistic non-hierarchical clustering (K-means) and Ward and centroid linkage algorithms (hierarchical clustering).
Classification of specimens in the obtained number of clusters was determined using partitioning around medoids (PAM), a nonhierarchical clustering method that is insensitive to both noise and outliers (Singh & Chauhan 2011). Moreover, the medoids are robust representations of the cluster centers, which is particularly important in the common context that many elements do not belong well to any cluster ( Van der Laan et al. 2003). This classification was used as the input categories for linear discriminant analysis (LDA) with leave-one-out cross-validation to generate the canonical axes, identify the shape variables involved in-group distinction and produce a robust reclassification of the specimens (prior probabilities were set to be equal) (Fig. 2).
The taxonomical significance of the obtained classification was tested following approaches that use gaps in morphological variation as a criterion for species delimitation (Zapata & Jiménez 2012). Under this approach, it is assumed that phenotypic variation was the result of the effect of several genes and that there was random mating among those from a geographical region, thus making intraspecific morphological variation reasonably described by a normal distribution (Sokal & Rohlf 1995;Templeton 2006;Zapata & Jiménez 2012). Hence, evidence for morphological discontinuity between hypothesized species arises in assessing the degree of overlapping (i.e., the magnitude of the morphological gap) of the probability distributions of identified groups (Wiens & Servedio 2000). This is achieved estimating tolerance regions for the distribution of each hypothesized species at an a priori non-zero frequency cut-off (i.e., certain percentage of the hypothesized species is encompassed within these regions; see Wiens & Servedio 2000) and assessing their overlap.
Univariate normal distribution of each resulting cluster was tested with a Shapiro-Wilks normality test (α = 0.05) (Fig. 2). When the p-value is less than the alpha level and the hypothesis that the sample came from a normally distributed population (null hypothesis) is rejected, we stopped the analyses because the necessary condition to use gaps in morphological as a criterion for species delimitation was not met. The overlap was estimated computing one-tailed normal tolerance region (Young 2010(Young , 2014 for the 90% probability of each group distribution of along the canonical axes resulting from LDA, with 95% confidence (γ = 0.05). The confidence level is the frequency with which estimated tolerance regions encompass a proportion ≥95% of the sampled population. An a priori frequency cut-off of β = 0.1 was selected (i.e., the tails of the distributions of hypothesized species may only overlap in an area comprising the 10% or less of their populations; otherwise, they are considered to be overlapping; see Wiens & Servedio 2000). If the overlap of tolerance regions exceeds the frequency cut-off, we stopped the analyses because the condition to support the hypothesis of a species boundary was not met.
All the analyses were performed in the R environment (R Core Team 2017, v. 3.6.3) using the packages cluster (Maechler et al. 2014 -Scheme of the methodological procedure applied in this study. The enamel perimeter of each tooth in occlusal view was digitized using EFA, analyzing individually the outline of each tooth and premolar and molar module morphologies (upper and lower dentition) using (EFA). Scores of specimens on principal components of all harmonics resulting of EFA were used as the input variables for all subsequent statistical analyses. Morphological groups within each data set were explored using the gap statistic method. Classifications determined using PAM were used as the input categories for LDA. Univariate normal distribution of each resulting cluster was tested with a Shapiro-Wilks normality test and the morphological discontinuity among clusters was achieved estimating tolerance regions between distributions. An LDA was performed in order to produce a classification of the samples based on the occurrence locality and to generate a morphospace in which those categories were maximally separated. Then, we use the previously obtained clusters and geographical patterns to perform a visual inspection by plotting the specimens in a morphospace. For more details, see Materials and methods section. Scale bars: 5 mm.

Upper dentition
Lower dentition

GeoGraphical variation analyses
Following the hypothesis of geographical variation (i.e., morphological variability is expected to be organized following a geographical pattern), the correspondence between morphological and geographical aspects was assessed (Fig. 2). To this, we generated categories corresponding to the occurrence localities of each analyzed specimen, which were grouped based on their geographic proximity: NWA, La Pampa, Mendoza, and Atlantic Coast (detailed in Figure 1B and Table 1). Although the assessment of chronostratigraphic aspects is beyond the scope of this study, we also include references about this matter (Table 1), potentially relevant regarding the temporal distribution of Tremacyllus.
An LDA was performed in order to produce a classification of the samples based on the occurrence locality and to generate a morphospace in which those categories were maximally separated. As input data for each data set, we used a specific number of variables (PCs) limited by the number of specimens of the smallest geographical category. Then, we plotted the cluster membership and geographical provenance on the mentioned morphospace to explore the relationship between these variables (Fig. 2).
To visualize the geographical distribution of the morphological pattern (i.e., clustering structure), we used the spatial datahandling capabilities of the Geographical Information System (GIS) (by QGIS 3.10 R Core Team 2017). This analysis was performed using the Kernel density estimation (KDE, Baxter et al. 1997) -based heat-map, which is one of the commonly used tools for spatial point pattern visualization and analysis (Yuan et al. 2019). Generally, a heat-map uses a warm and cold colored continuous surface to reflect the density distribution trend of point data without overlapping (Wilkinson & Friendly 2009). But, in this case, we use the heat-map to represent the geographic location of the studied specimens in each data set (i.e., tooth or module contour) by weighing the morphological variability explained for the first component of PCA.
In the Results section, we characterize each morphospace (upper and lower dentition outlines), describing the main shape changes explained for each axis (after PCA performed on resulting EFA). For practical reason, the cumulative variation in the first two components, which represent most of the total variation, are described in detail. Also, we visually explore clustering patterns in the distribution of the specimens. Then, we report the clustering structure and morphological gaps of each data set. Finally, we analyze those clusters statistically supported within a geographical context.

RESULTS
anatomical contour analysis of the upper Dentition M1 Seven harmonics were selected for the description of the tooth outlines of 36 specimens. Eight principal components (i.e.,  shape variables) resulting from EFA were retained after principal component analysis for the set of contours describing the M1 outline (96.8% of the total variation; Appendix 5). The first principal component (PC1) accounts for 53.5% of the total variation and describes differences in tooth overall shape and ectoloph outline, with almost elliptical tooth and convex ectoloph in positive extreme and rectangular-shaped tooth and relatively undulated ectoloph in negative extreme of PC1 (Fig. 3A). The second component (PC2) accounts for 18.4% of the total variation and describes changes in the labiolingual width along the mesiodistal axis of the tooth. High values of PC2 correspond to morphologies with a wider mesial face than the distal one, while low values of PC2 correspond to shapes with mesial and distal faces similar in width (Fig. 3A). There is no visual clustering pattern in the distribution of the specimens in the morphospace of the first two principal components (Fig. 3A).

M2
Eight harmonics were selected for the description of the tooth outlines of 32 s pecimens. Eight principal components were retained for the set of contours describing the M2 outline (96.7% of the total variation [Appendix 5]). The PC1 (56.2%) describes changes in the labiolingual width of the tooth and parastyle orientation, with wide tooth and parastyle labially oriented in positive extreme and narrow tooth and parastyle mesially oriented in negative extreme of PC1 (Fig. 3B). The PC2 (13.8%) mainly describes differences in ectoloph, lingual face, and parastyle. Higher values of PC2 correspond to shapes with relatively undulated ectoloph, rounded distal face, and well-marked parastyle pointing labially, while lower values of PC2 correspond to morphologies with convex ectoloph, straight distal face, and parastyle relatively defined (Fig. 3B). The morphospace of the first two principal components shows some outliers with extreme forms; however, most of the specimens present similar shapes (Fig. 3B).

M3
Eight harmonics were selected for the description of the tooth outlines of 30 specimens. Eight principal components were retained for the set of contours describing the M3 outline (96.3% of the total variation [Appendix 5]). The PC1 accounts for 47.7% of the total variation. It describes the relative change from a labiolingually wide and straight mesial face (parallel oriented to the minor axis of the tooth) and wide and rounded distal face in positive extreme to a labiolingually narrow and obliquely oriented (distolingual-mesiolabially) mesial face and small rounded distal face in negative extreme (Fig. 3C). The PC2 accounts for 21.5% of the total variation and mainly describes the shape change of the distal face, where it is well-   rounded reaching almost a lobe shape, pointing to the distal side in positive extreme, while it is small and pointing to the labial side in negative extreme of PC2 (Fig. 3C). There is no visual clustering pattern in the distribution of the specimens in the plot of the first two principal components (Fig. 3C).

Upper molars module (UMM)
Twelve harmonics were selected for the description of the tooth outlines of 30 specimens. Nine principal components were retained for the set of contours describing the UMM outline (95.7% of the total variation [Appendix 5]). The PC1 (53.4%) describes changes in the shape of the M1, degree of contact between M1 and M2, and M3 configuration with respect to the M2. High values of PC1 correspond to morphologies with M1 that show a straight lingual face and mesiolabially more projected parastyle than that of M2, short M1-2 contact and restricted to the labial region and M3 more lingually projected than the level of M2 and ectoloph of M2-3 reaching a similar level. Low values of PC1 correspond to morphologies with M1-2 parastyle similar in shape with rounded lingual face, wide M1-2 contact that occupies the middle section of the mesial face of M2 and parastyle of M3 exceeding the level of the ectoloph of M2 (Fig. 3D). The PC2 (13.5%) mainly describes changes in the size ratio of the M1-2, ectoloph morphology of M1, relative position of M2 with respect to the adjacent teeth and size ratio of the M2-3. High values of PC2 correspond to morphologies with M1 almost similar or slightly larger than M2, almost straight ectoloph of M1, M2 located more labially than M1 and M3, and M3 smaller than M2. Low values of PC2 describe morphologies with M1 smaller than M2, convex ectoloph of M1, M2 located at the same level of M1 and M3, and M3 similar in size to M2 (Fig. 3D). Although the morphospace of the first two principal components shows a few specimens with extreme forms, a clear visual clustering pattern in the distribution of the specimens is absent (Fig. 3D).

P2
Eight harmonics were selected for the description of the tooth outlines of 28 specimens. Seven principal components were retained for the set of contours (96.7% of the total variation [Appendix 5]). The PC1 (67.1%) describes differences in overall tooth shape, with elliptical P2 morphologies in positive extreme and triangular-shaped ones in negative extreme of PC1 (Fig. 4A). The PC2 (10.6%) mainly describes changes in ectoloph and mesio-lingual face. High values of PC2 correspond to shapes with markedly convex ectoloph and concave mesio-lingual face, while low values of PC2 correspond to morphologies that present an ectoloph with a small convex section and straight to convex mesio-lingual face (Fig. 4A).
In the morphospace, although most of the specimens are located in positive values of the PC1, there is no clear visual clustering pattern (Fig. 4A).

P3
Eight harmonics were selected for the description of the tooth outlines of 35 specimens. Eight principal components were retained for the set of contours (97.0% of the total variation [Appendix 5]). The PC1 (56.2%) describes changes in the labiolingual width of the tooth and lingual face extension, with narrow P3 and straight and long lingual face in positive extreme and markedly wide P3 and convex and short lingual face in negative extreme of PC1 (Fig. 4B). The PC2 (13.6%) mainly describes differences in the parastyle morphology. High values of PC2 correspond to morphologies with well-marked and acute parastyle, mesially pointing, while low values of PC2 correspond to shapes with less-developed and wide parastyle (Fig. 4B). Visually, the distribution of the specimens in the morphospace does not suggest a clear clustering pattern (Fig. 4B).

P4
Eight harmonics were selected for the description of the tooth outlines of 37 specimens. Seven principal components were retained for the set of contours (95.9% of the total variation [Appendix 5]). The PC1 (63.2%) describes changes in the parastyle and ectoloph morphology, development and orientation of the distal face and arrangement of the mesial and lingual faces. High values of PC1 correspond to shapes with wide parastyle, small and distal concave section in the ectoloph, straight distal face orientated perpendicular to the mesiodistal axis of the tooth, and lingual and mesial face slightly forming an obtuse angle. Low values of PC1 correspond to morphologies with acute parastyle pointing mesially, convex ectoloph, table 2. -Clusters and morphological gaps of the upper dentition. Principal components scores (resulting from EFA and PCA) retained for statistical procedures, clustering algorithms used to determine the optimal number of groups under gap statistic methods (K-means, Ward, and centroid indices), Shapiro-Wilk test values for each cluster (α = 0.05), and tolerance regions (γ = 0.05; β = 0.1). Note that the overlap is calculated using the upper limit of Cluster A and the lower limit of Cluster B on their probability distributions along the canonical axis resulting from the LDA. Groups without overlap are highlighted in bold font, providing statistical evidence for the existence of a morphological gap (Fig. 5). straight distal face obliquely oriented to the mesiodistal axis of the tooth and lingual and mesial face markedly forming an obtuse angle (Fig. 4C). The PC2 (12.7%) mainly describes differences in the ectoloph morphology, labiolingual width of the distal face, and arrangement of the mesial and lingual faces. High values of PC2 correspond to morphologies with markedly convex ectoloph, wide distal face, and lingual and mesial face markedly forming an obtuse angle; while lower values of PC2 correspond to shapes with concave ectoloph, narrow distal face and lingual and mesial face slightly forming an obtuse angle (Fig. 4C). There is no clear visual clustering pattern in the distribution of the specimens in the morphospace of the first two principal components (Fig. 4C).

Upper premolars module (UPM)
Twelve harmonics were selected for the description of the tooth outlines of 27 specimens. Eight principal components were retained for the set of contours describing the UPM outline (94.8% of the total variation [Appendix 5]). The PC1 (69.7%) describes changes in the degree of imbrication (overlapping) among P2-4 and the orientation of P2. High values of PC1 correspond to morphologies with a low degree of imbrication among the teeth and elliptical P2 with its major axis almost mesiodistally oriented. Low PC1 values describe morphologies with teeth highly imbricated, elliptical-shaped P2 with its major axis obliquely oriented (distolingual-mesiolabially), and almost elliptical P4 (Fig. 4D). The PC2 (9.3%) mainly describes variations in the size ratio among each tooth. Despite that, the size increase from P2 to P4, high PC2 values correspond to morphologies with P2 markedly smaller than P3 and P4, while low PC2 values describe a gradually increasing in tooth size, from P2 to P4 (Fig. 4D). The morphospace of the first two principal components shows a clustering pattern in the distribution of the specimens, where two groups are identified by the PC1 values (Fig. 4D).

clusters anD morpholoGical Gaps of the upper Dentition
The three clustering algorithms tested found that the best number of clusters in the data set of each contour analysis is two ( Table 2). The normality test performed on the distribution of specimens resulting from LDA indicates that each cluster corresponds to a normally distributed population (Table 2); except for M3 (null hypothesis rejected [ Table 2]). Tolerance regions calculated for each recognized morphological group did not overlap only for P4 and UPM contour analyses, providing statistical evidence for the existence of a morphological gap in the clustering pattern within each data set (Table 2; see Discussion for a morphological description of these groups). From the LDA, probability distributions of each cluster along the canonical axis are essentially unimodal and Gaussian (Fig. 5); however, an incipient bimodal curve is observed in one cluster of the UPM data set (Fig. 5B). On the other hand, tolerance regions overlap for the morphological groups of the other analyzed data sets (M1, M2, P2, P3 and UMM; Table 2), failing to support the hypothesis of a species boundary according to the a priori frequency cutoff.
analysis of the lower Dentition m1 Seven harmonics were selected for the description of the tooth outlines of 42 specimens. Nine principal components (i.e., shape variables) resulting from EFA were retained after principal component analysis for the set of contours describing the m1 outline (96.4% of the total variation [Appendix 5]). The PC1 (49.6%) describes variations in the development of the labial sulcus and trigonid outline, morphologies with a less-marked labial sulcus and rounded trigonid, markedly smaller than the talonid, are located in the positive extreme, while morphologies with well-marked labial sulcus and triangular-shaped trigonid, smaller than the talonid, are in negative extreme of PC1 (Fig. 6A). The PC2 (18.3%) describes changes in the labiolingual width of the tooth. High values of PC2 correspond to morphologies with a wide m1, while low values of PC2 correspond to narrow ones (Fig. 6A). There is not a visual clustering pattern in the distribution of the specimens in the morphospace (Fig. 6A).

m2
Eight harmonics were selected for the description of the tooth outlines of 35 specimens. Nine principal components were retained for the set of contours describing the m2 out- line (95.7% of the total variation [Appendix 5]). The PC1 (47.1%) describes differences in the development of the labial sulcus, trigonid outline, and relative size difference between the trigonid and talonid. High values of PC1 correspond to morphologies with well-developed labial sulcus and triangular trigonid similar in size to the talonid (Fig. 6B). Low values of PC1 correspond to morphologies with less-developed labial sulcus and rounded trigonid that is smaller than the talonid (Fig. 6B). The PC2 (17.7%) mainly describes differences in the lingual face of m2 and labiolingual width of the tooth, morphologies with almost straight lingual face and narrow tooth are observed in the positive extreme, while shapes with convex lingual face and wide tooth are in the negative extreme of PC2 (Fig. 6B). Visually, the distribution of the specimens in the morphospace does not suggest a clear clustering pattern (Fig. 6B).

m3
Nine harmonics were selected for a description of the tooth outlines of 28 specimens. Ten principal components were retained for the set of contours describing the m3 outline (96.7% of the total variation [Appendix 5]). The PC1 (49.8%) describes variations in the development of two labial sulci, the shape of the trigonid, and the development of the "third lobe" (common to all pachyrukhines). High values of PC1 correspond to morphologies with well-developed labial sulci, triangular trigonid and well-delimited "third lobe" (Fig. 6C).
Low values of PC1 correspond to morphologies with less developed labial sulci, rounded trigonid and less delimited of the "third lobe" (Fig. 6C). The PC2 (16.0%) describes changes in the shape of the trigonid, and development of the "third lobe", with shapes with triangular trigonid and less development of the "third lobe" in positive extreme, and shapes with rounded trigonid and well-differentiated "third lobe" in negative extreme (Fig. 6C). There is a partial clustering pattern identified by the first principal components values (Fig. 6C).

Lower molars module (LMM)
Thirteen harmonics were selected for the description of the tooth outlines of 26 specimens. Ten principal components were retained for the set of contours describing the LMM outline (94.6% of the total variation [Appendix 5]). The PC1 (52.4%) mainly describes changes in the imbrication (overlapping) among m1-3 (Fig. 6D). High values of PC1 correspond to morphologies with high imbrication, while low values correspond to shapes with low imbrication among the teeth (Fig. 6D). The PC2 (10.9%) describes relative changes in the development of the labial sulcus of each tooth, where morphologies with a well-marked labial sulcus are in positive extreme, while less marked labial sulcus ones are located in negative extreme (Fig. 6D). Visually, the distribution of the specimens in the morphospace does not suggest a clear clustering pattern (Fig. 6D).  Eight harmonics were selected for the description of the tooth outlines of 34 specimens. Nine principal components were retained for the set of contours describing the p2 outline (96.0% of the total variation [Appendix 5]). The PC1 (32.5%) describes differences in the morphology of the labial face and labiolingual width of the tooth, shapes with an undulate labial face with two shallow concavities ("trilobed") and narrow p2 are located in the positive extreme, while morphologies with a marked sulcus in the labial face and wide p2 (molariform morphologies) are in negative extreme of PC1 (Fig. 7A). The PC2 (26.5%) also describes changes in the morphology of the labial face and labiolingual width of the tooth. High values of PC2 correspond to morphologies with a shallow labial sulcus and wide p2, while low values of PC2 correspond to shapes with a posterior well-marked labial sulcus and an anterior shallow concavity in the labial face and narrow p2 (Fig. 7A).
There is no visual clustering pattern in the distribution of the specimens in the morphospace of the first two principal components (Fig. 7A).

p3
Seven harmonics were selected for the description of the tooth outlines of 46 specimens. Ten principal components were retained for the set of contours describing the p3 outline (97.0% of the total variation [Appendix 5]). The PC1 (39.6%) mainly describes variations in the development of the labial sulcus and trigonid outline (Fig. 7B). High values of PC1 correspond to morphologies with less-marked labial sulcus and triangular trigonid with rounded edges, while low values correspond to shapes with well-marked labial sulcus and rounded trigonid (Fig. 7B). The PC2 (19.1%) describes changes in the relative size of the talonid, where morphologies with labiolingually wider talonid than the trigonid are in positive extreme, while shapes with talonid similar in width to trigonid are located in negative extreme (Fig. 7B). Visually, the distribution of the specimens in the morphospace does not suggest a clear clustering pattern (Fig. 7B).

p4
Six harmonics were selected for the description of the tooth outline of 47 samples. Nine principal components were retained for the set of contours describing the p4 outline (97.1% of the total variation [Appendix 5]). The PC1 (42.6%) describes differences in the paraconid morphology and labiolingual width of the tooth, shapes with a rounded paraconid pointing mesially and narrow p4 are located in the positive extreme, while morphologies with acute paraconid pointing mesiolabially are in negative extreme of PC1 (Fig. 7C). The PC2 (19.7%) describes variations in the paraconid morphology and the relative size of the talonid and trigonid. High values of PC2 correspond to morphologies with an acute paraco-  nid pointing mesially and larger trigonid than the talonid, while low values of PC2 correspond to shapes with a rounded paraconid pointing mesially and smaller trigonid than the talonid (Fig. 7C). There is no visual clustering pattern in the distribution of the specimens in the morphospace of the first two principal components (Fig. 7C).

Lower premolars module (LPM)
Fourteen harmonics were selected for the description of the tooth outline of 32 specimens. Ten principal components were retained for the set of contours describing the LPM outline (93.4% of the total variation [Appendix 5]). The PC1 (56.7%) mainly describes changes in the degree of imbrication among p2-4 and the relative size of the talonid of p4, where high imbrication and larger talonid than the trigonid are observed in the morphologies located in positive extreme, while low imbrication and smaller talonid than the trigonid are observed in morphologies located in negative extreme (Fig. 7D). The PC2 (9.9%) describes changes in the relative size of the p2 and p3, and the talonid outline of p4. High values of PC2 correspond to morphologies with mesiodistally long p2, shorter p3 than the p2, and rounded talonid of p4, while low values of PC2 correspond to shapes with mesiodistally short p2, longer p3 than the p2, and triangular talonid of p4 (Fig. 7D). Visually, the distribution of the specimens in the morphospace does not suggest a clear clustering pattern (Fig. 7D).

clusters anD morpholoGical Gaps of the lower Dentition
The three clustering algorithms tested found that the best number of clusters in the data set of each contour analysis is two (Table 3). The normality test performed on the distribution of specimens resulting from LDA indicates that each cluster corresponds to a normally distributed population (Table 3). Tolerance regions calculated for each recognized morphological group did not overlap only for m3 and p3 contour analyses, providing statistical evidence for the existence of a morphological gap in the clustering pattern within each data set (Table 3). Probability distributions of each cluster along the canonical axis (resulting from the LDA) are essentially unimodal and Gaussian (Fig. 8). On the other hand, tolerance regions are overlapped for the morphological groups of the other analyzed data sets (m1, m2, p2, p4, LMM and LPM; Table 3), failing to support the hypothesis of a morphological gap according to the a priori frequency cutoff.

morpholoGical Gaps across the GeoGraphical variation
A scatterplot of the first two discriminant functions of the P4 data set (Fig. 9A) shows that the specimens from the NWA are relatively well separated by the first canonical axis (96.8% of explained variation), in positive extreme. Specimens from La Pampa, Mendoza, and Atlantic Coast a overlap in the negative extreme. The second canonical axis (2.6% of explained variation) does not achieve a good separation of specimen provenance. On the other hand, most of the specimens from NWA were previously sorted within Cluster A (morphological group), while specimens from the other regions were classified within Cluster B (Fig. 9A; Appendix 6). However, there are a few mismatches; four specimens from NWA (FMNH-P 14374, FMNH-P 14498, FMNH-P 14538, and FMNH-P 15254; see Discussion section) were included within the Cluster B and two specimens, from La Pampa (GHUNLPam 21813) and Mendoza (IANIGLA-PV 517) were included within Cluster A ( Fig. 9A; Appendix 6).
Regarding the UPM data set, the scatterplot of the first two discriminant functions allows separating the NWA specimens from the other ones ( Fig. 9B; Appendix 6). The first canonical axis (77.9% of explained variation) distinguishes specimens from La Pampa, Mendoza, and Atlantic Coast, which overlap in positive extreme, from those found in NWA, in negative extreme. In turn, the second canonical axis (22.1% of explained variation) separates specimens from La Pampa in positive values from those from Mendoza and Atlantic Coast, which overlap in the negative extreme (Fig. 9B). Specimens from NWA do not achieve a good separation by this axis (Fig. 9B). The morphological clustering pattern shows that those specimens from NWA were classified within Cluster A, except for one specimen (MPAT 727, Cluster B), while the others were included within Cluster B, except for one sample from Mendoza (MACN-Pv 8495, table 3. -Clusters and morphological gaps of the lower dentition. Principal components scores (resulting from EFA and PCA) retained for statistical procedures, clustering algorithms used to determine the optimal number of groups under gap statistic methods (K-means, Ward, and centroid indices), Shapiro-Wilk test values for each cluster (α = 0.05), and tolerance regions (γ = 0.05; β = 0.1). Note that the overlap is calculated using the upper limit of Cluster B and the lower limit of Cluster B on their probability distributions along the canonical axis resulting from the LDA. Groups without overlap are highlighted in bold font, providing statistical evidence for the existence of a morphological gap (Fig. 8).
Morphological gaps data m1 m2 m3 LMM p2 p3 p4 LPM PC scores PC1-9 PC1-9 PC1-10 PC1-10 PC1-9 PC1-10 PC1-9 PC1-10 K-means index (N clusters) Regarding the lower dentition, both m3 and p3 data set scatterplots do not to achieve an unambiguous separation of specimens based on their provenance (Fig. 9C, D; Appendix 6). Although the specimens from La Pampa of the m3 data set appear well distinguished in the positive extreme of the first canonical axis (71.3% of explained variation), the other ones overlap in both the first and second canonical axes (Fig. 9C). In turn, no correspondence between geographical and morphological distributions was observed, as none of the specimen provenance categories were representative of a single morphological cluster or vice versa (Fig. 9C, D).
Considering that P4 and UPM data sets segregate specimens in two main groups that fit well with both geographical provenance and morphological clusters, we perform an additional analysis combining both sets of contours (harmonics describing shape) and selecting the same specimen sample (n = 27). Seventeen principal components were retained that explain 96.1% of the total variation (Appendix 5). The three clustering algorithms tested found that the best number of clusters is two (Table 4). The normality test performed on the distribution of specimens resulting from LDA indicates that each cluster corresponds to a normally distributed population, unimodal and Gaussian ( Fig. 10A; Table 4). Tolerance regions calculated for each recognized morphological group did not overlap ( Fig. 10A; Table 4). From the scatterplot of the first two discriminant functions of P4 + UPM, we can see two main groups markedly separated based on different geographical categories (Fig. 10B). The first canonical axis (99.1% of explained variation) separates NWA specimens (except MPAT 727, see Discussion section) from those of La Pampa-Mendoza-Atlantic Coast (overlapip). The second canonical axis explains 0.9% of the total variation. Regarding the morphological clustering pattern, specimens from NWA were sorted within Cluster A, except one sample (MPAT 727, Cluster B) (Fig. 10B). Specimens from La Pampa, Mendoza, and Atlantic Coast were classified within Cluster B (Fig. 10B). In sum, these results indicate a strong correspondence between geographical and morphological distributions, with the specimens from NWA in Cluster A (except MPAT 727) and those from La Pampa-Mendoza-Atlantic Coast in Cluster B (Appendix 6). table 4. -Clusters and morphological gaps of the combining P4 and UPM coefficient of harmonics datasets. Principal component scores (resulting from EFA and PCA) retained for statistical procedures, clustering algorithms used to determine the optimal number of groups under gap statistic methods (Kmeans, Ward, and centroid indices), Shapiro-Wilk test values for each cluster (α = 0.05), and tolerance regions (γ = 0.05; β = 0.1). Note that the overlapping is calculated using the upper limit of Cluster A and the lower limit of Cluster B on their probability distributions along the canonical axis resulting from the LDA. Note that there are not overlapping between groups, providing statistical evidence for the existence of a morphological gap (Fig. 10A).

DISCUSSION
The morphological variability displayed by the dentition of Tremacyllus specimens is organized into two main groups within each data set. Nevertheless, these groups are separated by a statistically significant morphological discontinuity only for P4, UPM, m3, and p3 analyses. The fact that the other data sets (M1-3, P2-3, UMM, m1-2, LMM, p2, p4, and LPM) are also organized into two groups could be related to an artifact of the statistical algorithm used for determining the relevant number of clusters. Given that this procedure recognizes a minimum of two clusters (Charrad et al. 2014), if the structure of the data set does not present any clustering pattern, it would not be observable as a first output result. In this sense, the evidence for morphological discontinuity was best achieved by estimating tolerance regions and their overlap for the distribution of each hypothesized cluster (Wiens & Servedio 2000;Milla Carmona et al. 2016). Among statistically supported clusters, only the resulting groups within P4 and UPM data sets are closely related to geographical categories when specimen provenance was used as a reference for grouping (Figs 9; 10). Probably, m3 and p3 supported clusters reflect different ontogenetic stages or sexual dimorphism (see Armella 2022). Nevertheless, degree of tooth wear and stage of preservation could be involved in this observed clustering pattern (see below). Regarding Cluster A within P4 and UPM data sets, it is distinctly composed of specimens belonging to NWA geographical category, while Cluster B (P4 and UPM) is composed of specimens from La Pampa, Mendoza, and Atlantic Coast geographical categories (mismatches detailed below). In this sense, the correspondence between the morphological and geographical evidence (see Material and methods section; Wiens & Servedio 2000; Zapata & Jiménez 2012) supports the idea that there are at least two species present in the studied sample, distinguished by the P4 and UPM morphology and their geographical distribution.
Cluster A is characterized by a P4 with almost rectangular distolingual face with rounded edges, ectoloph slightly concave in the distal sector and convex near the parastyle region, and mesial face obliquely oriented (linguodistalmesiolabially), well-differentiated from the lingual one fig. 11. -Mean shapes of the two identified morphological clusters, resulting from EFA, for the P4 and UPM data sets. Note that Cluster A is strongly related to specimens from NWA, herein referred as Tremacyllus incipiens Rovereto, 1914 ( Fig. 11). Regarding the UPM of Cluster A, the premolars are less elongated in their distolingual-mesiolabial axes, clearly less imbricated than Cluster B, and the P2 is markedly smaller than P3 and P4 and the P2 of Cluster B (Fig. 11). Cluster B is characterized by a P4 with less rectangular distolingual face in comparison with the P4 of Cluster A, convex ectoloph, and mesial face more obliquely oriented (linguodistal-mesiolabially), slightly differentiated from the lingual one (Fig. 11). The overall comparison between clusters reveals a more elliptical contour of P4 in Cluster B than A. Regarding the UPM of Cluster B, the premolars are characterized by elongated main axes, stronger imbrication pattern, size decreasing gradually from P4 to P2, and P2 proportionally larger than the P2 of Cluster A (Fig. 11).
It is important to note that, although Cluster A of P4 and UPM analyses and the NWA geographical category show a strong correlation, there are few specimens that did not follow this pattern. Regarding the P4 data set, four specimens (FMNH-P 14374, FMNH-P 14498, FMNH-P 15254, and FMNH-P 14538) of 19 from NWA were classified within Cluster B (Fig. 9A). The specimens FMNH-P 14374 and FMNH-P 14498 (NWA) are badly preserved, showing broken areas mainly in the cementum of the teeth; however, these specimens have been included in the analyses since the enamel outline could be recognized (Fig. 12A-E). Probably, this preservation affects the morphology of the P4, resulting in the classification within Cluster B. Unfortunately, the fact that these specimens lack P2-3 (FMNH-P 14374) prevent their inclusion in the UPM analysis and more detailed interpretations. A particular case is observed in the specimen FMNH-P 15254, which has a P3 with two longitudinal sulci in the mesial face and a particularly short P2 (mesiodistally) with a marked concavity in the posterior region of the ectoloph (Fig. 12F). These features are not observed in other specimens and the fact that FMNH-P 15254 is only a left maxillary fragment avoids more detailed inferences. Besides, although FMNH-P 15254 dentition seems to be permanent, ontogenetic factors should be not discarded given that there are cases of Tremacyllus specimens with the absence of deciduous teeth in a clearly young individual (Cerdeño et al. 2017). For instance, GHUNLPam 21813 (La Pampa) is a juvenile, but shows a permanent dentition (P2-M1;Cerdeño et al. 2017: fig. 6). In fact, in our analysis, GHUNLPam 21813 was classified within Cluster A based on P4 morphology, as well as the specimen IANIGLA-PV 517 (Mendoza). These mismatches are here interpreted as variation related to certain preservation bias, different ontogenetic stages (e.g. age of the individual), or even stage of wear, which could explain the shape changes of the teeth. On the other hand, the specimen FMNH-P 14538 (NWA) shows central enamel bands in P3-4 and a small isolated central fossette in P2 (Fig. 12G). The presence of these traits has been interpreted as representing an early ontogenetic stage (see Cerdeño et al. 2017) and, probably, results in the classification of FMNH-P 14538 within Cluster B. When this sample was tested based on UPM morphology, it was sorted within Cluster A. This fact could be related to the fact that UPM analysis recovers variation in imbrication and size among teeth pieces, while not P4 morphology. Combining both data sets (i.e., P4 + UPM), FMNH-P 14538 was sorted within Cluster B ( Fig. 10B; Appendix 6). With respect to UPM mismatches, one specimen (MPAT 727; Fig. 12H) of 14 from NWA was sorted within the Cluster B, and two (MACN-Pv 8495 and MACN-Pv 2434) of 13 samples from southern localities were included within Cluster A (Fig. 9B). In these cases, particular preservation, ontogenetic, age or stage of wear factors are not observed. When we perform an analysis combining both P4 and UPM data sets, these mismatches disappear and only the specimen MPAT 727 maintains its status as a member of Cluster B (Fig. 10). The latter result suggests two main points that should be taken into account. Firstly, although P4 and UPM morphologies seem to be individually very informative to identify Cluster A and B, the combination of both traits is even more useful (Fig. 13). Secondly, the classification of the specimen MPAT 727 within Cluster B raises some interesting palaeobiogeographical aspects (see below).
In a taxonomic context, when examining the classification of type specimens of Tremacyllus, T. impressus (syntype MACN-A 1377; Vera & Ercoli 2018), and T. subdiminutus (holotype MACN-Pv 8494;Rovereto 1914) were sorted within the Cluster B (P4, PM, and P4 + UPM; Figs 9A, B; 10; 11). The fact that these taxa were classified in the same morphological cluster is consistent with the proposal of Vera & Ercoli (2018). In that study, the authors performed a morphological comparison of the dentition of a large sample of Tremacyllus from Mendoza with betterknown species of the genus (taking into account type material and referred specimens). Since their analyses did not support a distinction between these species, they proposed the name T. subdiminutus as a junior synonym of T. impressus. Furthermore, specimens from La Pampa analyzed here were also classified within Cluster B (P4, PM, and P4 + UPM; Figs 9A, B; 10; 11). In this sense, Sostillo et al. (2018) referred a large sample of specimens of La Pampa to Tremacyllus impressus based on the comparison with type material and referred specimens, including some samples of Mendoza. Our results agree with these previous studies highlighting that specimens from both Mendoza and La Pampa, along with specimens from the Atlantic Coast, show similar morphologies, which were statistically supported (i.e., minimal internal variation), particularly regarding P4 and UPM outlines (see above). Based on this, we consider that specimens classified within Cluster B should be assigned to Tremacyllus impressus.
Regarding other type specimens of Tremacyllus, T. incipi ens (holotype MACN-Pv 8163) and T. latifrons (holotype MACN-Pv 8157) were described by Rovereto (1914) based on specimens from Catamarca Province (Santa María valley; Fig. 1). Rovereto (1914) characterized T. incipiens by an excavated palate, elliptical P2, which is the smallest tooth, while the other pieces increase in size gradually to the first molar, M3 with well-developed cusps, a wide diastema, long lower teeth with marked labial sulci, and p3 with two cusps (i.e., roughly bilobed). It differs from T. impressus by a more excavated palate, short post incisive depressions (posteriorly to the incisive foramen), and P2 variable in shape (Rovereto 1914 Argentina (e.g. Entre Ríos Province), in which recent palaeoenvironmental studies proposed closed, wet, and forested areas developed under tropical-subtropical conditions during this period (Schmidt et al. 2020). In the context of this wide geographical distribution of Tremacyllus, we found support for a morphological segregation between southern representatives of the genus (recognized as T. impressus; from La Pampa, Mendoza, and Atlantic Coast, 38 specimens) and the northern ones (herein recognized as T. incipiens; from Catamarca and Tucumán, 42 specimens). The clustering structure shown by LDA and geographically visualization by the KDE heat-map (Fig. 14) represent the morphological variability distribution of the studied sample. These clusters have two well-differentiated morphologies in which, roughly, the specimens from NWA group in one extreme and those from the Mendoza, La Pampa, and Atlantic Coast in the opposite side ( Fig. 14A-C). This geographical pattern agrees with palaeoenvironments (Latorre et al. 1997;Barreda et al. 2007;Hynek et al. 2012;Georgieff et al. 2017;Domingo et al. 2020) and the two palaeo-phytogeographic provinces proposed by the Late Miocene-Pliocene (Barreda et al. 2007), in which Tremacyllus species fall: Neotropical province for T. incipiens (Cluster A), and Proto-Spinal Steppe province for T. impressus (Cluster B) (Fig. 14). In this sense, although the distribution of Tremacyllus specimens would be related, broadly speaking, to the presence of open and relatively dry areas, the regional palaeoenvironmental conditions were different, which probably triggered different dietary adaptations which can explain part of the tooth morphology variation here observed (see Ercoli et al. 2021). While this is of relevance for the geographical variation hypothesis and the taxonomic concerns herein addressed, it is an interesting aspect, since the latitudinal segregation could reflect, among other aspects, particular environmental conditions for each region, which probably played an important role in the establishing of these populations. Furthermore, the occurrence of more and complete specimens in middle latitudes, between both southern and northern ones (such On the other hand, a particular case is observed among Tremacyllus records from NWA. The specimen MPAT 727 shows features in the P4 and UPM that result in its inclusion within Cluster B (herein Tremacyllus impressus). This result is suggestive; however, since this is a single specimen, it is difficult to perform an integrative assessment. Moreover, a mandible from the India Muerta Formation, PVL 7559 ( Fig. 12I; see Alonso 2018) also shows an unusual position within the analysis (Fig. 9C, D) but does not seem to be very meaningful at morphological level (as it only preserves lower dentition). Considering this scenario, the scarce record of Tremacyllus in both eastern outcrops (India Muerta and Las Cañas formations [ Fig. 1C]) hampers more encompassing studies; thereby, we prefer to maintain an open taxonomic approach for MPAT 727 and PVL 7559 and avoid referring these records to a particular species of the genus.
Even though, the fact that the specimen MPAT 727 presents a different tooth morphology regarding that of western Tremacyllus records stand as an interesting point to be considered. This specimen was recorded in the Las Cañas Formation (Santiago del Estero Province [ Fig. 1C]) from beds dated as 3.73 ± 0.07 Ma (Early Pliocene; Armella et al. 2020), representing the youngest record of this genus in the region (and the easternmost along with PVL 7559). Considering that almost all of the specimens analyzed in our study are recorded in Late Miocene-Pliocene western outcrops of NWA (e.g. Santa María valley [ Fig. 1C]), explanations for this phenomenon could involve temporal or local palaeoenvironmental variations.
The main mountain systems of NWA (Cumbres Calchaquíes and Aconquija Ranges [ Fig. 1C]) began to uplift at c. 6 Ma, driving different depositional conditions in the eastern and western slopes (Sobel & Strecker 2003). At c. 3 Ma, these ranges represent a well-established orographic barrier, responsible of the rain shadow effect (Kleinert & Strecker 2001;Carrapa et al. 2012). During this lapse, the east-west precipitation gradient was gradually changing each local condition: arid at the western side and humid at the eastern side ( Armella et al. 2020). In this context, a differential distributional pattern and particular taxonomic compositions among western and eastern outcrops of NWA were proposed, at least regarding notoungulate faunas (Armella 2019;Armella et al. 2019). The differential morphology of MPAT 727 would reflect this community difference, an idea to be tested with further studies conducted on eastern outcrops. reviseD DiaGnosis. -Excavated palate; P2 elliptical with its major axis almost mesiodistally oriented, markedly proportionally smaller than P3 and P4, and the P2 of Tremacyllus impressus; P3 elliptical with its major axis obliquely oriented (distolingual-mesiolabially); P4 with almost rectangular distolingual face with rounded edges, ectoloph slightly concave in the distal section and convex near the parastyle region, straight distal face oriented perpendicular to the mesiodistal axis of the tooth and mesial face obliquely oriented (distolingual-mesiolabially) well-differentiated from lingual one, both forming an obtuse angle; upper premolars clearly less imbricated than in

CONCLUSIONS
The morphological variability displayed by the dentition of Tremacyllus specimens is organized into two statistically wellsupported clusters. Morphological gaps are evident based on P4, upper premolars (P2-4), m3, and p3 morphology. Particularly, the features of the P4 and upper premolars are markedly related to a geographical pattern. Following the hypothesis of the geographical variation, correspondence between the morphological and geographical evidence supports the idea that there are, at least, two species present in the studied sample distinguished by the P4 and upper premolars morphology. Based on the classification of type specimens of Tremacyllus and samples attributed to T. incipiens and T. latifrons within statistically supported clusters, we consider T. incipiens as a valid and endemic taxon from the Late Miocene-Pliocene western outcrops of NWA (Chiquimil and Andalhuala formations). Moreover, we formally propose the synonymy of T. incipiens and T. latifrons. The strong geographical pattern allows us even to assign herein analyzed mandibular remains as T. incipiens, while we have not found features useful to diagnostic inferences in the lower dentition. Two isolated specimens recorded in eastern outcrops of NWA (India Muerta and Las Cañas formations) show particular morphologies and geographic contexts, hampering a specific taxonomic assignment.
In this scenario, the segregation between southern representatives of the genus (recognized as T. impressus; from La Pampa, Mendoza, and Atlantic Coast) and the northern ones (herein recognized as T. incipiens; from Catamarca and Tucumán) agrees with the already proposed particular environmental conditions for each region, which probably played an important role in the establishment and morphological differentiation between these populations.
The integration of morphological criteria, biological concepts, and quantitative framework adopted in this study proved to be very effective at solving the distinction between inter-and intraspecific variation and recognizing taxa in Tremacyllus. Quantification of the morphological variation using geometric morphometrics, exploration of the grouping structure present in the data set, and quantitative assessment of the taxonomic level of variation in accordance with the chosen criteria seem to be suited to deal with the study case of these taxa. This approach is easily applicable to the occurrence of this genus in other regions, being the main restrictions the quality of preservation and sample size. In this sense, the fact that dental remains are the most abundant elements in the Tremacyllus record makes that the resulted morphological gaps could be considered, per se, strong traits for diagnostic purposes. Nevertheless, additional cranial and postcranial traits and more and betterpreserved specimens could be important for reinforcing our hypothesis. On the other hand, this general methodology can be extended to the study of other taxa, in order to give taxonomic significance to the morphological variability using a rigorous quantitative protocol.