Ethical statement. No live animals were handled in this study.
The sample examined consists of 112 gorilla specimens, including 10 Gorilla beringei beringei, 25 Gorilla beringei graueri, and 77 Gorilla gorilla gorilla. We selected only fully erupted M2s (either left or right, but no antimeres) because they provide a good general overview of the development of masticatory function in primates43. We included only molars with a slight to moderate degree of wear because in heavily worn teeth occlusal facets tend to coalesce often preventing a clear identification34. We qualitatively evaluated the level of wear based on the amount of dentine exposure and cusp removal (wear stages 1-4)44. When available we included from museum records45 information about age (subadult and adult), sex and locality (electronic supplementary material, Table S7).
Occlusal Fingerprint Analysis (OFA)
Three-dimensional (3D) digital models of teeth were post-processed using Polyworks® V12 (InnovMetric Software), a 3D metrology software. OFA consists of four consecutive steps: 1) model orientation, 2) facet identification, 3) facet area, and 4) facet inclination29. The polygonal models are oriented using a reference plane that is created along the cervical line of the tooth through the least square best-fit method. Successively, the reference plane is rotated to the xy plane obtained from the original coordinate system.
The identification of wear facets follows the numbering system originally created by Maier and Schneck46 and later modified by Kullmer and colleagues29. Facets 1, 1.1, 2, 2.1, 3 and 4 develop along the buccal slopes of protoconid, hypoconid and hypoconulid, while facets 5, 6, 7 and 8 form along the buccal slopes of the metaconid and entoconid (Fig. 3). These facets are created during the initial phase of the rhythmic chewing cycle (phase I), and they can be further divided into buccal (facets 1-4) and lingual (facets 5-8) facets37. Facets 9, 10, 11, 12 and 13 are formed during the second phase of the rhythmic chewing cycle (phase II), when the mandibular molars move out of occlusion. These facets develop along the lingual slopes of protoconid, hypoconid and hypoconulid.
During the rhythmic chewing cycle, the food bolus is initially processed by shearing, which is generated by forces parallel to the contact plane (phase I), followed by crushing between basins and cusps of molars, where the occlusal force is perpendicular to the contact plane37. After the molars reach centric occlusion (phase II), the food bolus is processed by grinding, which is the resulting action of the combination of perpendicular and parallel forces to the contact plane.
We also identify tip crush areas, which are generally formed during puncture-crushing, when food is initially pulped by a series of masticatory cycle where tooth-to-tooth contacts do not occur27. Finally, we identify three additional facets in gorilla molars: facet 5.1, facet 8.1 and facet 10.1. These facets are created by the presence of additional dental traits, which are less common in Homo. More specifically, facets 5.1 forms in mandibular molars around the area of cusp 7 (C7), and it occludes with maxillary molars in proximity of the lingual cingulum47. Facet 8.1 forms around cusp 6 (C6) and it occludes with the lingual cingulum of the maxillary molars in correspondence of the mesiolingual aspect of the protocone48. Facet 10.1 has been previously described in Neanderthal molars, and it is created by the contact between the distolingual slope of the metaconule of the maxillary molar with the mesiolingual slope of the protoconid49. Facets 5.1 and 8.1 are grouped with lingual phase I facets, while facet 10.1 is considered a phase II facet.
Once the facets have been identified, we automatically calculate the surface area using the area function available in Polyworks® V12 (InnovMetric Software). To facilitate the analysis of the occlusal contact areas we grouped phase II facets together with tip crush areas, and divided the phase I facets into buccal and lingual facets34. However, tip crush areas and phase II facets have been kept separated for the analysis of facet inclination, considering they are characterised by significant different angles40. Because larger teeth generally develop larger facets than smaller teeth, we eliminate the size factor by using relative areas only. Relative areas are calculated by dividing the absolute area of a facet with the total wear area (TWA), which is the sum of absolute area of all facets34. The facet inclination is calculated by measuring the angle between the reference, or cervical plane, with the facet plane, which is created by selecting all triangles within its perimeter and by applying the best-fir plane function of Polyworks® V12 (InnovMetric Software) (electronic supplementary material, figure S1).
Occlusal relief index (OR)
The occlusal relief index (OR) was calculated by dividing the 3D area by the two-dimensional (2D) area of the occlusal surface50. The occlusal plane, parallel to the cervical plane, was translated along the y axis until it reached the deepest point of the occlusal surface, known as the central fossa. Next, the digital model was sliced, with respect to the occlusal plane, and the 2D area was calculated. The 3D area was calculated by selecting all triangles of the polygonal models above the occlusal plane (electronic supplementary material, figure S2).
Percentage of dentine exposure (PDE) and enamel wear (PEW)
To calculate the percent of dentine exposure (PDE) is obtained by dividing the sum of total dentinal areas for each tooth by the 3D occlusal area and then multiplied by 10026. The percentage of enamel wear (PEW) is obtained by dividing the total wear area (TWA) with the 3D occlusal area and then multiplied by 100.
Statistical analysis
We employed summary statistical analyses (median and standard deviation, SD) for each variable considered in this study. Because we further divided our sample according to wear stages 1 to 444 the comparative groups became relatively small. As the presence of a small sample size prevents the assumption of a normal distribution, we used the nonparametric Mann-Whitney pairwise test, to test whether two univariate samples are taken from populations with equal medians51.
For the analysis of facet areas, we grouped together molars characterised by wear stage 2 and wear stage 3 in order to maximize the sample size. We then tested if we could combine these two groups by comparing molars wear stage 2 with molars wear stage 3 of western lowland gorillas without finding any statistically significant result (Table S8). However, for the inclination we kept the molars separated by their wear stages, because the level of wear strongly influences the values of these variables40. OR, PDE and PWE were not separated into different wear stage groups. All statistical analyses have been performed using the software PAST v4.04 (Paleontological Statistics)52.