Identifying biological affinities of Holocene northern Iberian populations through the inner structures of the upper first molars

Neolithisation was a relatively fast process that affected both the interior and coastal zones of the Iberian Peninsula, but it was also a heterogeneous process that had diverse impacts on genomic and cultural diversity. In the Late Neolithic–Chalcolithic, a change in funerary practices, cultural material and trade networks occurred, and genomic heterogeneity decreased, suggesting human mobility and genetic admixture between different Iberian populations. Dental morphology has emerged as an effective tool for understanding genomic variability and biological affinities among ancient human populations. But, surprisingly, less attention has been paid to the morphological traits of inner dental tissues in Holocene European populations and their utility for the study of population dynamics. We applied 3D geometric morphometric methods on the enamel-dentine junction (EDJ) of the first upper molars to explore the biological affinities of north-eastern Iberian Peninsula populations from the Late Neolithic–Chalcolithic to the Bronze Age. Our results show that the EDJ morphologies of the northern Iberian Peninsula populations were generally homogeneous, indicative of genetic admixture as a result of human mobility and exchange networks. However, differences in the EDJ traits in remains from the Can Sadurní site are indicative of distant biological affinities with nearby populations. Additionally, the hypocone associated dentine area and the position of the trigon dentine horns relative to each other on the occlusal surface best describe the variability found among the samples studied. This study highlights the utility of EDJ morphology as a genetic proxy in Holocene population dynamic studies when paleogenomic studies are absent.


Introduction
Agriculture reached the Iberian Peninsula (IP) around 5,600-5,500 cal BC, during what is known as the Early Neolithic (ca 5,600-4,800 cal BC) Martins et al. 2015;Oms et al. 2016;Peña-Chocarro et al. 2018;Edo et al. in press). The neolithisation process in the IP was a relatively fast event that reached both the interior and the coastal zones (Aubán et al. 2018;Oms et al. 2018;Rojo-Guerra et al. 2018). However, it was a heterogeneous process with incoming populations exerting different impacts on the various areas of the peninsula (Hervella et al. 2012;García-Martínez de Lagrán et al. 2018). This heterogeneity is reflected in cultural and funerary practices throughout the Neolithic period, and particularly in the north-eastern (NE) IP. The Middle Neolithic (or Middle Neolithic I; ca 4,100-3,500 cal BC; Oms et al. 2016) in the NE-IP was characterised by notable regionalisation, evident in the variety of funerary practices (e.g. Castany 2008; Gibaja et al. 2010;Roig et al. 2010;Subirà et al. 2016), as well as the presence of extensive exchange networks established between groups in the western Mediterranean (Blasco 1996;Villalba et al. 2011;Borrell and Bosch 2012;Bosch et al. 2012;Weller and Fíguls 2012;Edo et al. 2012;Gibaja et al. 2013;Terradas et al. 2014;Gibaja et al. 2018), which reflects the economic strength of this period (Cebrià et al. 2013). However, by the Late Neolithic-Chalcolithic period (ca 3,500-2,300 cal BC; Oms et al. 2016), this system had collapsed, leading to both cultural and socioeconomic consequences (Martin 2003;Gibaja et al. 2010;Cebrià et al. 2013;Subirà et al. 2016). Changes in funerary practices occurred in the NE-IP, ranging from inhumations of one or two individuals-either in pit burials (also known as sepulcres de fossa) or cists (stone boxes)-during the Middle Neolithic to collective burials, mostly in caves (Cebrià et al. 2013;Martin 2003;Daura et al. 2015;Gómez et al. 2008). Settlement and subsistence patterns and cultural material also underwent changes during this period (Martin 2003), suggesting the disappearance of some of the exchange networking that took place in the previous period (Oms et al. 2016). Ancient DNA (aDNA) analyses of Iberian specimens have shown signs of incoming populations in the peninsula throughout the entire Neolithic period (Hervella et al. 2012;Olalde et al. 2015;Szecsenyi-Nagy et al. 2017). These studies also evidence a decrease in genomic diversity beginning in the Late Neolithic-Chalcolithic and continuing into the Early Bronze Age (2,200-1,500 BC), suggesting human mobility and genetic admixture between different Iberian populations (Szecsenyi-Nagy et al. 2017;Olalde et al. 2018).
Although recent aDNA methods have successfully overcome the challenges DNA preservation in warm climates such as that of the Mediterranean (Olalde et al. 2019;Fernandes et al. 2020), preservation remains an issue in larger sample sizes from ancient periods. Moreover, human remains found in particular funerary contexts are often fragmented, burned and commingled, which leads to further difficulties in their genetic analysis. In such fragmentary contexts, teeth take on increase significance, as their strong genetic component (Scott and Turner 2008) makes them an effective genetic proxy, their traits remain unchanged over time, and they are abundant in the archaeological record (Hillson 2005).
Dental morphology analyses have been employed in multiple studies to interpret the phylogenetic relationships between different species in the human lineage and the variability of living human populations (e.g. Benazzi et al. 2011Benazzi et al. , 2012Delgado et al. 2019;Gómez-Robles et al. 2012;Martin et al. 2017;Martinón-Torres et al. 2006;Rathmann and Reyes-Centeno 2020;Skinner et al. 2008aSkinner et al. , 2016. They have contributed to our understanding of affinities among ancient human populations by establishing historical and biological relationships between different human groups (Cucina et al. 1999;Irish and Joel 2005;Coppa et al. 2007). Dental morphology studies of prehistoric human populations in the IP are scarce and tend to focus on non-metric dental traits (Ruiz et al. 2012;Scott et al. 2013;Ceperuelo et al. 2014;Horwath et al. 2014;Subirà et al. 2016;. Some of these studies have shown biological homogeneity in some Neolithic and Late Neolithic-Chalcolithic populations in NE-IP (Subirà et al. 2016), while others not . However, the biological affinities in the NE-IP, specifically in Late Neolithic-Chalcolithic transition populations, remain unclear and require further analysis to fully comprehend them.

Geometric morphometrics on the enamel-dentine junction as a dental morphology proxy
In recent years, high resolution imaging techniques have enabled researchers to noninvasively explore detailed tooth morphology and to analyse the inner structure of both fossil and living human dentitions (Skinner et al. 2010;Prado-Simón et al. 2012;Smith and Olejniczak 2012;Fornai et al. 2016;Hublin et al. 2017). Major aspects of crown morphology, such as the number and size of cusps, are influenced by genetic, epigenetic and environmental factors (see Townsed et al. (2012) for an overview). Although there is strong evidence of the genetic component during odontogenesis (e.g. Scott et al. 2018), environmental and developmental effects may not be underestimated to the contribution to the final phenotype (Riga et al. 2014). Dental morphogenesis is explained by the morphodynamic model (Jernvall and Jung 2000;Salazar-Ciudad and Jernvall 2002). Small changes in gene expression and signalling protein activity affecting the timing and position of the specialised inhibitors-activator signalling centres (enamel knots) during tooth morphogenesis have been shown to affect final crown morphology Jernvall 2002, 2010), and particularly with the last-forming dental traits (Riga et al. 2014;Paul et al. 2017). The distribution of these enamel knots in the inner enamel epithelium during early tooth crown formation is what ultimately shapes the folded configuration of the basement membrane (Jernvall and Jung 2000). Many features of the tooth crown originate in this membrane, which is preserved in mature teeth as the enamel-dentine junction (EDJ) prior to enamel deposition and is primarily responsible for external crown morphology (Butler 1956;Skinner et al. 2008b;Morita et al. 2014). Interestingly, the EDJ has proved a useful proxy for worn teeth in the study of taxonomic differences among some primate and Homo taxa (Skinner et al. 2008a(Skinner et al. , 2010Ortiz et al. 2017;Hanegraef et al. 2018) as well as in living human populations. This has made it possible to increase the sample size in dental morphology studies (Smith et al. 2006;Monson et al. 2020). Additionally, the EDJ has been found to vary according neutral evolution resulting from genetic drift and migration patterns, thus reflecting past human history and population dynamics (Monson et al. 2020). However, environmental and developmental effects are also likely expected to affect EDJ, as mentioned above. Three-dimensional (3D) geometric morphometric (GM) approaches have captured shape differences in morphological traits of the EDJ through landmark/semilandmark data in living humans and our closest relatives (e.g. Skinner et al. 2010;Polychronis and Christou 2013;Al-Shahrani et al. 2014;Morita et al. 2014;Fornai et al. 2015Fornai et al. , 2016Martin et al. 2017;Zanolli et al. 2018;Krenn et al. 2019). GM methods constitute an especially useful technique for describing spatial aspects of morphological variation, as well as for directly visualising shape differences (Zelditch et al. 2012;Adams et al. 2013). Surprisingly, phenotypic variability in the EDJ has scarcely been examined in any prehistoric European population after the adoption of agropastoralism (Gamarra et al. 2020;Le Luyer et al. 2016;Le Luyer and Bayle 2017).
This study makes unprecedented use of the EDJ morphology of the first upper molars to analyse biological affinities between prehistoric human populations of the NE-IP associated with socioeconomic transitions and genetic admixture. Specifically, and assuming a strong genetic component for EDJ (Monson et al. 2020), we aim to confirm the decrease in genetic variability detected during the Late Neolithic-Chalcolithic period and demonstrate the use of the EDJ as a genetic proxy in population dynamic studies of prehistoric populations.

Materials
We selected 23 upper first molars (M 1 ) from NE and central IP sites dated from the Late Neolithic-Chalcolithic to the Bronze Age, although we were only able to analyse 16 samples ( Fig. 1 and Online Resource 1). The upper first molar was selected because it is the most stable tooth in the permanent molar row (Butler 1939;Dahlberg 1945;Townsend et al. 2009;Morita et al. 2014;Paul et al. 2017). Since crown formation begins in utero (Kraus and Jordan 1965;Hillson 1996), the upper first molar is less susceptible to environmental influences than later-forming teeth in the molar row (Paul et al. 2017;Monson et al. 2020), and therefore more directly reflects the individual's underlying genotype (Paul et al. 2017). The selected teeth presented well-preserved, complete crowns, with well identified cemento-enamel junctions (CEJ). The molars showed a degree of wear equal to or less than 3 on Molnar's scale (Molnar 1971). Samples with cracks and damage on the crown dentine (see below for micro-CT acquisition) were excluded from morphological analyses. Since morphological asymmetry tends to be minor among dental antimeres Turner 1997, Gómez-Robles et al. 2007), it is a common practice in studies of inner dental crown tissues to select both right and left M 1 (e.g. Le Luyer 2016; Martin et al. 2017;Zanolli et al 2018). The teeth selected were either on the maxilla in situ or isolated. When isolated, and in order to maximise the sample size, no discrimination was made between the right and left molars. However, only the side with the highest number of samples was chosen for each site in order to ensure a single tooth per individual. When both molars were present (in the maxilla), the better-preserved tooth was selected. Details on sample composition and wear degree are given in Table 1.

Micro-CT image acquisition and image processing
The samples were scanned using two scanners housed at Centro Nacional de Investigación sobre la Evolución Humana (CENIEH) facilities, in Burgos (Spain). Most of the samples were scanned using a V|Tome|X s 240 µ-CT from GE Sensing & Inspections Technologies. Scans were performed with a 0.2-mm copper filter and a filament voltage of 150 kV, a current of 150 µA and an acquisition time of 500 ms. Output images were then processed with a median and a ring artifact filter to improve image contrast and reduce noise. The resultant output images had a voxel size ranging between 20 and 25 µm for isolated teeth and teeth included in maxillary fragments, and between 60 and 91 µm when teeth were isolated from cranium scans. Some of the samples from El Mirador cave were scanned using a Scanco Medical AG Micro-Computed Tomography 80. Scans were taken using two 0.1 mm copper filters and a filament voltage of 70 kV, a current of 114 µA and an average acquisition time of 300 ms. The resultant slice thickness was 18 µm for isolated teeth and 74 µm for teeth included in maxillary fragments.
The microCT data for each sample was imported into Avizo software (version 2020.2 trial license), where dentine Archaeological sites: 1: Cova de les Agulles, 2: Cova de Can Sadurní, 3: Cova de la Guineu, 4: Cova dels Galls Carboners, 5: Cueva de El Mirador. Generic Mapping Tools 4.5.13 (Wessel and Smith 1998) and the topographic ETOPO data set (Amante and Eakins 2009) was used to create this map Table 1 List of the samples employed in this analysis. Wear stage according to Molnar (1971) Site ID Skeletal remains M1 Wear stage Horn tips reconstructed tissue was semi-automatically segmented using a 3D voxel value histogram and grayscale values (Martin et al. 2017) with manual corrections (Olejniczak et al. 2008;Le Luyer et al. 2016;Le Luyer and Bayle 2017;Zanolli et al. 2018). Due to sample variable image resolution, and in order to avoid significant effects on differential resolution on the segmentation step, we equalled image resolution among all samples by resampling the images stack of those samples presenting the highest resolution range (20 to 25 µm) to lower resolution range (60 and 91 µm). Three-dimensional (3D) surface models of the EDJ of each tooth were generated using a constrained smoothing algorithm (Le Luyer and Bayle 2017; Zanolli et al. 2018) and imported in stl format. Lastly, the EDJ 3D surface models were edited and mirrored to left molars, when the left molar was missing or unusable. For samples presenting some degree of occlusal wear (Table 1), dentine horn apices were manually reconstructed following the interpolation approach described in Zanolli et al. (2018). Digital surfaces of the EDJ 3D surface models are available on MorphoSource (see Online Resource 1 for digital object identifiers).

Enamel-dentine junction shape analyses
Viewbox 4 dHAL software, version 4.1.0.8, was used to create a template of eight landmarks and 59 curve semilandmarks to describe the EDJ surface (Table 2 and Fig. 2), following the workflow described in Bastir et al. (2019). The occlusal surface was defined by seven anatomical landmarks corresponding to the four tips of the dentine horns and three points identified on the occlusal surface (Fornai et al. 2015). Curve semilandmarks were also used to describe the marginal ridges that connect the dentine horns and other identified landmarks on the occlusal surface (Fornai et al. 2015).
The cervical outline was also characterised along the CEJ with 24 curve semilandmarks. The initial landmark of the cervical outline was placed on the middle part of the buccal face of the crown (between the paracone and metacone) and continued mesially (Martin et al. 2017). Enough landmarks and semilandmarks were employed to ensure that the main morphological aspects of the EDJ were captured. The template was subsequently applied to the EDJ surface of each sample by means of a thin-plate spline (TPS) interpolation   Table 2 for detailed descriptions of the anatomical landmarks function, and projected onto the target. The semilandmarks were allowed to slide along the curve semilandmarks to minimise the TPS bending energy between each EDJ specimen and the template (Gunz and Mitteroecker 2013). After TPS interpolation, semilandmarks were treated as geometrically homologous points between individuals (Gunz et al. 2005), being able to be analysed in conjunction with landmark points (Le Luyer et al. 2016;Bastir et al. 2019;Sorrentino et al. 2020). Different landmark and semilandmark sets were compared in order to analyse which landmark and semilandmark configurations best discriminate the sample distribution in the morphometric space. The landmark and semilandmark sets corresponded to (1) all landmark and semilandmark configurations (EDJ+CER), (2) cervical landmarks and semilandmarks (CER) and (3) occlusal traits of the EDJ landmark and semilandmark configurations (EDJ). Shape coordinates were then imported and a Generalized Procrustes Analysis (GPA) was performed using "Morpho" package in R v. 3.5.2 (R Core Team 2020), where scale, position and orientation were standardised to minimise the overall sum of the squared distances between corresponding landmarks and semilandmarks. The size measure was represented by the natural logarithm of centroid size (lnCS).

Statistical analyses
A principal component analysis (PCA) was performed on the Procrustes coordinates for all landmark and semilandmark configurations (EDJ+CER, CER, EDJ) to explore morphological variation across the sample (see Online Resource 1 for sample identifier equivalence). A column represented by the natural logarithm of centroid size (hereafter called lnCS) was added to shape coordinates in order to perform a PCA in the form space (size + shape) for all configurations (Klingenberg 2016). We applied the broken stick model (Frontier 1976, adapted from MacArthur 1957 to estimate how many principal components (PC) to retain (Online Resource 1). According to this method, we decided to use the two first PCs for graphic representation in every case using ggplot2 R package. We calculated the shape variations associated with the extremes of the PC scores to visualise the morphological variability of the sample analysed (Morpho, Rgl and Arothron R packages; Schlager 2017; Profico et al. 2018). Since normality and homogeneity of variance were not possible to test for each site (some of the sites were represented by one individual), non-parametric statistical tests were performed. Shape variation related to size (i.e. allometry) was explored by means of Spearman correlations (rho) of the PC scores against lnCS. Site differences observed in the first two PCs were statistically tested using Kruskal-Wallis tests for all landmark and semilandmark sets in both shape and form spaces . To assess the possible error derived from digitalising the corresponding points, a repeatability test for landmark and semilandmark placement was performed. We applied the intraclass correlation coefficient (Fisher 1958) to quantify measurement error in repeated measures, which consisted of performing one-way ANOVA on repeated measurements and employing the individual as the categorical variable (Fruciano 2016). A subset sample of five individuals was randomly chosen and landmarks and semilandmarks were digitalised three times with 1 week between iterations. Our repeatability test resulted in a 94%, indicating an excellent reliability (Portney and Watkins 2009). Data were statistically analysed with R software.

Results
The shape and form-space PCA plots are represented in Figs. 3 and 4. The first two PCs in shape-space accounted for 43.4% of the total variance when the EDJ+CER configuration was used. Patterns of EDJ+CER surface morphology of individuals belonging to different sites were significantly different across the first PC (24.2%) ( Table 3). More specifically, all of the individuals sampled from the Can Sadurní site were distinct from the other Late Neolithic-Chalcolithic and Bronze Age sites, which appear to overlap in the shape-space (Fig. 3a). The corresponding shape change mainly involves the expansion of the distolingual dentine horn associated with the hypocone at the cervical outline and the relative position of the dentine horns associated with the trigon cusps on the occlusal surface (Online Resource 2). The distribution of the samples across the PCA shapespace shows that individuals from the Can Sadurní site tend to be associated with lower (negative) PC1 values, with a buccolingually elongated trigon area on the occlusal surface as a result of more lingually oriented horns associated with the protocone and relatively higher trigon horn tips, together with a narrower talonid area at the cervical outline and higher cervical height (distance from the ridges connecting horn tips to the cervical line; Monson et al. 2020) than the EDJ morphologies of individuals from other sites. Samples with higher (positive) PC1 values also tend to be associated with a mesio-distal elongation of the occlusal area derived from the distolingual extension of the cervical area associated with the hypocone, as well as inner located dentine horns associated with the trigon cusps. Although not been statistically significant (Table 3), samples associated with higher (positive) PC2 (19.2%) values also tend to be associated with broader mesio-distal occlusal surfaces derived from a broader talonid occlusal area (Online Resource 3). No significant (P>0.05) Spearman correlations were found for either PC's with lnCS. The inclusion of size information in the form-space PCA explains higher total variance for the first two PCs (62.1%), and this variable was retained in PC1 (rho= −0.99, P<0.001). However, PC1 (49.9%) did not Fig. 3 Shape space PCA plots for a) EDJ+CER, b) CER and c) EDJ configuration (description in Material and Methods section). The deformed warped surfaces (occlusal and mesial view) of enamel-dentine junction and cervical and occlusal wireframes of first upper molar are drawn at the end of the axes to illustrate the extreme morphological trends.

Fig. 4 Form space PCA plots for a) EDJ+CER, b) CER and c) EDJ configuration (description in Material and
Methods section). The deformed warped surfaces (occlusal and mesial view) of enamel-dentine junction and cervical and occlusal wireframes of first upper molar are drawn at the end of the axes to illustrate the extreme morphological trends significantly distinguish between EDJ+CER surface morphologies of individuals belonging to different sites, though PC2 was close to significance (12.2%) (Fig. 4a, Table 3). PC2 was not statistically correlated with lnCS (rho= −0.09, P=0.737). This means that size, represented here by the lnCS, did not affect shape variation between individuals from different sites. Samples from different sites distributed along the PC2 tend to be associated with the same occlusal traits as the ones described for the PC1 shape-space.
When only the cervical outline was used to detect differences among the samples (CER), the two first PCs in the shape-space accounted for a similar total variance (44.6%; Fig. 3b) as for when all landmarks (EDJ-CER) were used. Although neither PC1 (29.3%) nor PC2 (15.3%) accounted for statistical differences between sites (Table 3), the Can Sadurní individuals tend to have smaller (negative) PC1 values. The cervical morphology distributed in smaller PC1 values is associated with mesio-distally narrower cervical outlines derived from the reduction of the lingual dentine horns in the cervical area. Additionally, most of the Can Sadurní individuals presented significantly smaller cervical outlines associated with lower PC1 values (rho=0.70, P=0.003) (Fig. 5a). No significant size correlations were found for PC2 values (rho= −0.10, P=0.704). Although the samples slightly overlap in the form-space PCA (Fig. 4b) and the two first PCs accounted for a higher total variance (79.3%), none of them significantly differentiated the samples from different sites (Table 3). Individuals' distributions along the form-space PC1 were associated with the same cervical outline traits as in shape-space PC1, but with inverted negative and positive values. As for the EDJ+CER analysis, size information was retained in the form-space PC1 (72.9%) (rho= −0.99, P<0.001), but not in PC2 (6.4%) (rho= −0.023, P=0.934).
The distribution of individuals belonging to different sites overlaps in both shape and form-space PCAs when only occlusal traits (EDJ) were characterised (Figs. 3c and 4c), although in the case of the space-space PCA, this overlap was caused by two individuals from Can Sadurní. The total variance is explained by the fact that the two first PCs decreased in both PCA spaces compared to previous landmark set analyses (shape-space: 38.9%; form-space: 59.7%), and no significant differences were found in any of the PCs (Table 3). Despite the lack of significance, two Can Sadurní individuals exhibited higher (positive) PC1 values, which are associated with a reduced talon area and a mesio-distally expanded trigon occlusal area. As in the CER landmark and semilandmark configurations, PC1 was significantly correlated with lnCS in both space analyses (shape-space: rho= −0.57, P=0.020; form-space: rho= −0.96, P<0.001), while PC2 was not (shape-space: rho=0.17, P=0.512; form-space: rho= −0.18, P=0.490) (Fig. 5b).

Fig. 5
Shape variation (PC scores) related with size (lnCS) in a) CER and b) EDJ landmark and semilandmark configuration in shape-space analysis. Only significant Spearman correlations were represented

EDJ variability in the north-eastern and central Iberian Peninsula
The results obtained here suggest that EDJ morphologies of individuals from Late Neolithic-Chalcolithic sites in the NE and central IP are generally homogeneous. The EDJ morphology of most of the sites overlap in the morphospace, and although it could be affected by a bias due to a small sample, this presumably indicates phenotypic similarity and, hence, biological proximity among samples. However, individuals from the Can Sadurní site present significantly different EDJ morphologies, despite the small sample size, suggesting distant biological affinities with individuals from the same NE-IP area. Specifically, M 1 of the Can Sadurní individuals analysed here presented buccolingually larger occlusal trigon area, with horns associated with the protocone and hypocone more lingually oriented and mesio-distally narrower at the cervical area-particularly at the hypoconetogether with relative higher cervical height. This indicates that the Late Neolithic-Chalcolithic Can Sadurní population may be a biologically distinct population. This is in accordance with non-metric dental studies that show that the Late Neolithic-Chalcolithic Can Sadurní population differs from other IP and Mediterranean populations with similar chronologies (Pascual et al. 2018). Previous research has suggested that the people buried at Can Sadurní during the Late Neolithic-Chalcolithic period were part of a closed population based on the high frequency of shovel-shape incisors and enamel hypoplasia documented in them (Carrasco et al. 1989;Blasco 1993;Roca 2012). In this regard, environmental factors, such as malnutrition or systemic disease (which cause most of the enamel hypoplasia (Goodman and Armelagos 1985;Hillson and Antoine 2011), are not totally discarded for understanding Can Sadurní significant EDJ differences. Interestingly, some authors have found a link between the presence of enamel hypoplasia and a greater morphological variability in the later-forming dental traits developed during tooth morphogenesis (e.g. hypocone, extra cusps) in upper molars (Riga et al. 2014). The authors do not suggest a direct link between both traits, but state that both are the result of the same disruptive factor that influence tooth developmental process, and hence the result of more dental morphology diversity. With all the above, dental morphology analyses suggest that people continuously buried throughout the Late Neolithic-Chalcolithic in Can Sadurní might be part of an endogamic group with limited biological relationships with other human groups. But it is also possible that the Can Sadurní individuals may have been part of a descent group that occupied a specific territory, as cultural elements (Vasos de Carena Alta) in common with other pre-coastal Catalan sites have been recovered (Cova Bonica and Fou de Muntaner (Baldellou 1974); Abrics de Sota Penya I, II and IV (Marcet 1981;Petit Mendizábal 1989;Villalba et al. 1989)). Additional analyses including other archaeological sites from the surrounding area (Cova de Can Figueres (Safont and Subirà 1996); Cova de Cassimanya (Edo 1997); Cova Bonica (Baldellou 1974); Marge del Moro (Edo and Artasona 1989)), together with dental remains from the same site attributed to other periods, should shed further light on the hypothesis proposed. The homogeneity of M 1 morphology at other NE and central IP sites agrees with the decrease in genetic variability among the IP population by the Late Neolithic-Chalcolithic and Early Bronze Age, as evidenced by aDNA studies (Szecsenyi-Nagy et al. 2017;Olalde et al. 2018). The population growth during Late Neolithic-Chalcolithic in some areas of the IP, which is instantiated by the increase in the funerary archaeological record, has been suggested as the cause of mobility among communities (Villalba-Mouco et al. 2020). Although the complexity of Late Neolithic-Chalcolithic societies of the northern parts of IP is little understood, transitions to narrower exchange networks, particularly from the Ebro basin (Montes and Alday 2012), together with similar funerary practices (collective burials) and grave-goods analyses, have also been proposed for the cause of the phenotypic similarity among Late Neolithic-Chalcolithic populations from the northern IP (López-Onaindia and Subirà 2017). The exchange of goods and supplies was most probably accompanied by population movements and the admixture between them (Soriano 2017). Taken together, these events may have given rise to the decrease in genetic variability and, consequently, the homogenisation of phenotypic traits such as dental morphology. Although the results obtained here are indicative of the hypothesis posed above, such interpretations must be taken with caution, as we performed an exploratory analysis of the suitability of the method with a small sample size. The inclusion of greater number of Late Neolithic-Chalcolithic samples, especially those found along the Ebro basin and NE-IP (including the Pyrenees area), together with a multi-perspective analysis from other archaeological disciplines, will provide insight into the social complexity of the northern IP during this period.

EDJ as a genetic proxy in Holocene phenotypic affinity studies
This study also confirms the utility of 3D GM analysis of the EDJ of the first upper molars as an alternative approach in the study of biological and genetic affinities in past human populations, as shown by previous studies (e.g. Monson et al. 2020). The combination of EDJ and cervical outline traits is what significantly differentiated the samples studied here, with apparently not significant allometry component-size effect influencing shape-among the individuals analysed. The greatest variability among the samples was found in the shape of the dentine horn associated with the hypocone at the cervical outline, followed by the position of the trigon dentine horns relative to each other on the occlusal surface, both of which allowed us to differentiate one population from another. Although the cervical outline has been successfully used to distinguish between different hominin species (Benazzi et al. 2012;Bailey et al. 2014Bailey et al. , 2020 and past human populations (Gamarra et al. 2020), no differences were found between the populations studied here based on cervical traits alone. Nevertheless, the cervical outline, together with occlusal traits analysed separately, exhibited a significant allometric component, which might influence the distribution of the individuals in the morphometric space. Additionally, the high variability in the hypocone shape is also consistent with the developmental tooth formation cascade model (Jernvall and Jung 2000;Salazar-Ciudad and Jernvall 2002) and in accordance with previous studies on human populations (Polychronis and Christou 2013;Morita et al. 2014;Riga et al. 2014;Monson et al. 2020). Since the hypocone is one of the last cusps to form (Hillson 1996), and according with this model, its expression is conditioned by the available time and space left by previous reactiondiffusion signals during cusp formation and their accumulative effects.
Moreover, changes affecting the proximity of these signalling centres are reflected in the relative positions of the cusps and the dentine horn on the occlusal surface. According to this model, the lower mean distances between the main cusps on the developing crown relative to the available space and time allows for the activation of extra enamel knots, and therefore the proliferation of accessory traits on the dental crown (Jernvall and Jung 2000). Several studies have found that specific crown traits, such as relative intercuspal distances and the size of the hypocone, are correlated with Carabelli's trait expression (Romero et al. 2018;Hunter et al. 2010;Guatelli-Steinberg et al. 2013;Moormann et al. 2013;Paul et al. 2017;Scott et al. 2018). These correlations are probably the response to the same underlying morphogenetic mechanisms explained by the morphodynamic model (Moormann et al. 2013). Therefore, the expression of these dental features may covary and probably coevolved, which must be considered in phylogenetic and biological distance approaches that often treat dental characters as independent (Hunter et al. 2010;Salazar-Ciudad and Jernvall 2010;Moormann et al. 2013), such as some non-metric dental studies using discrete traits.
The morphometric analysis of the EDJ has proved an excellent method for the study of biological affinities between NE and central IP populations that is not affected by the interdependence of some discrete dental traits. This approach identifies dental variability among individuals that could not be detected by means of discrete traits on their outer surfaces, which contributes significantly to the field of dental morphology studies. Moreover, we have shown that the morphometric analysis of one dental piece, here the upper first molar, but potentially applicable to the other molar pieces, also provides insights into the relationships between these past human populations. Thus, the study of EDJ shape eliminates the need for several discrete traits from full dentition required in non-metric dental traits analyses for the study of biological affinities (e.g. Coppa et al. 2007;Rathmann and Reyes-Centeno 2020). This becomes crucial in fragmentary contexts and collective burials in which the success of the use of several non-metric dental traits depends on their presence in the sample. The analysis of EDJ has also proved useful despite moderate enamel wear, which, by contrast, may mask the expression of discrete dental traits, and thus prevent proper identification. Additionally, despite the small geographical range that is covered by the sample analysed here, the EDJ morphometric analysis have demonstrated to be sensitive enough to detect morphological differences among individuals belonging to populations that were geographically close (sites from NE area). These promising findings open the possibility to further explore populations' biological affinities, and hence population dynamics in a more regional scale.

Conclusions
The EDJ morphology of the first upper molar has revealed phenotypic similarities among Late Neolithic-Chalcolithic populations from the NE and central Iberian Peninsula, in keeping with the decrease in genetic variability documented in genetic and other dental morphology studies. The EDJ traits of the population buried at the Can Sadurní site differ significantly from those of the nearby populations compared here and the El Mirador cave site (central IP), indicating that they might be a closed group or part of descent group that occupied a specific territory with low levels of admixture. The traits that exhibit the greatest variability are those related to the hypocone area and the distribution of the trigon dentine horns on the occlusal surface.
To sum up, the morphometric analysis of the EDJ has further contributed to our understanding of past human biological affinities and has proved an effective tool for investigating population history in archaeological remains. More importantly, it provides an alternative analysis to traditional genetic studies by means of novel information of dental morphology patterns from human prehistoric populations. The EDJ therefore makes it possible to study human genetic variability while overcoming the drawbacks of aDNA analyses and the enamel dental wear that archaeological remains usually exhibit, as well as the interdependence of some dental discrete traits. Moreover, it represents a biological approach that attempts to understand the complexity of past human networks, which together with future studies in other fields will contribute to our understanding of the changes in social complexity that occurred during this time period.