Dental casts of 160 Greek subjects (80 males, 80 females) were scanned by a structured-light scanner. The upper and lower right first molar occlusal surface 3D meshes were processed using geometric morphometric methods. A total of 265 and 274 curve and surface sliding semilandmarks were placed on the upper and lower molar surfaces, respectively. Principal component analysis and partial least square analysis were performed to assess shape parameters. Molars tended to vary between an elongated and a more square form. The first two principal components (PCs), comprising almost 1/3 of molar shape variation, were related to mesiodistal-buccolingual ratios and relative cusp position. Distal cusps displayed the greatest shape variability. Molars of males were larger than those of females (2.8 and 3.2% for upper and lower molars respectively), but no shape dimorphism was observed. Upper and lower molar sizes were significantly correlated (r(2) = 0.689). Allometry was observed for both teeth. Larger lower molars were associated with shorter cusps, expansion of the distal cusp, and constriction of the mesial cusps (predicted variance 3.25%). Upper molars displayed weaker allometry (predicted variance 1.59%). Upper and lower molar shape covariation proved significant (RV = 17.26%, P < 0.0001). The main parameter of molar covariation in partial least square axis 1, contributing to 30% of total covariation, was cusp height, in contrast to the primary variability traits exhibited by PC1 and PC2. The aim of this study was to evaluate shape variation and covariation, including allometry and sexual dimorphism, of maxillary and mandibular first permanent molar occlusal surfaces.