The Gobi Desert hosts an extensive record of Cretaceous terrestrial ecosystems with some of the richest dinosaur sites in the world. However, many sites lack temporal information such as pollen and zircons, while limited exposures limit lithostratigraphic work. Biostratigraphy offers a potential solution. Oftentimes the most abundant diagnosable fossils from key Upper Cretaceous localities are the protoceratopsids. However, intra- and inter-specific variation and interrelationships within this clade from multiple localities have not been studied in detail, currently limiting their biostratigraphic use.

Here, we examined protoceratopsids from eight localities producing multiple specimens from four sedimentary basins. A character matrix was constructed to capture the full range of potential variation that may be taxonomically significant in their distribution, and included ontogenetic and individual variations. Character states and polarity were determined based on outgroup conditions and ontogenetic variations. For each taxon from each locality, the character states were plotted against the respective specimen size to effectively allow comparison among specimens of equivalent ontogenetic stages, given similar adult sizes. We recovered a total of four protoceratopsid species from the Gobi localities.

Taxonomic revision and the inferred character evolution merit a reassessment of the temporal relationships of the localities hosting the protoceratopsids. Specimens from the classic Djadokhta Formation localities Bayan Zag and Tögrögiin Shiree cannot be reliably distinguished and belong to a single species, Protoceratops andrewsi, as previously noted using geometric morphometric analyses and ontogenetically conserved traits. This shared presence of a single species may support the coeval, rather than successive, deposition of the two localities. When the evolutionary trend is coupled with co-occurrence data, the inferred biostratigraphic correlation disagrees with previous assessment of temporal ordering of localities, while cladogenesis and multiple dispersal events are inferred.

To test whether the sequence of protoceratopsid character evolution match the overall faunal biostratigraphic signal, quantitative biochronological analyses were conducted using the Unitary Association Method (UAM). Occurrence data were updated with the taxonomic revision of protoceratopsids from various localities. The UAM only recovered the overall signal of aeolian redbeds with small animals (Djadokhta and Baruungoyot formations) versus fluvial deposits with large animals (Nemegt Formation). This is consistent when only the protoceratopsid-bearing localities are accounted for and when aquatic taxa are omitted, thereby detecting no temporal resolution in the typical Djadokhta and Baruungoyot formation localities.