We define fossil occurrence as the presence of a particular taxon at a particular location and time based on its geologic strata, building on the definition by Butler and Barrett (2008). Our global occurrence database for Diplodocoidea is based on data from online repositories (namely, the Paleobiology Database, or PBDB; Uhen et al. 2023) and the literature. We started by downloading a compilation of all recorded Diplodocoidea body fossil occurrences from the PBDB (Suppl. material 1: table S2). Then, we individually verified and cross-referenced occurrences with the latest literature, with particular attention to taxonomy and litho- and chronostratigraphic context, e.g., maximum and minimum ages, early and late intervals, lithological unit, geological formation. We also conducted a comprehensive literature review, allowing us to add new and previously missing occurrences to our dataset. Finally, we applied the following occurrence inclusion criteria based on De Celis et al. (2020): (1) valid genera and/or species, as defined by the latest published literature; (2) taxonomically justified sensuDe Celis et al. (2020), i.e., taxa that are assigned to a species based on diagnostic characters rather than overall morphological similarity and/or similar temporal/geographical distribution; (3) unique genera per stratigraphic level and collection site; (4) assignment to genus/clade level; (5) described/reported in published literature. Ootaxa and ichnotaxa occurrences were excluded. Taxonomically verifiable occurrences with missing data for key variables were kept but not used in analyses, unless they could be estimated (e.g., approximating GPS coordinates based on nearest reported city/locality).

Our dataset includes occurrences published until June 2025 and totals 495 entries (466 of which were used in our analyses) corresponding to 43 genera (‘diplodocoidea_fullinfo’ in “Database S1.xls” Suppl. material 1). Among the 466 occurrences considered herein (Suppl. material 1: table S3; also ‘diplodocoidea_dataset’ in “Database S1.xls” Suppl. material 1), 68.24% derive from localities constrained to a single chronostratigraphic stage, whereas 31.76% originate from geological units with poorer temporal resolution, spanning more than one chronostratigraphic stage. Detailed information on the temporal ranges of the genera and clades considered in this study is provided in Suppl. material 1: tables S11, S12.

The Morrison Formation data was compiled starting with occurrences listed in Foster (2020) and Tschopp et al. (2022), and then complemented with PBDB occurrences that were verifiable but not included in these publications. The oldest rebbachisaurid occurrence considered here is Xenoposeidon proneneukos Taylor & Naish, 2007 from the Ashdown Formation (Late Berriasian–Early Valanginian) of England (Taylor and Naish 2007; Taylor 2018). We decided not to include the occurrence corresponding to Maraapunisaurus fragillimus (Carpenter 2018) in our datasets. Despite the fact that it could be the oldest record of a Rebbachisauridae, this taxon is based on a middle-posterior dorsal neural arch initially attributed to Amphicoelias fragillimus (Cope 1878) that has been lost, and its reinterpretation by Carpenter (2018) as a new rebbachisaurid genus Maraapunisaurus is just supported on the drawing of Cope (1878). Although Mannion et al. (2021) tentatively agree with Carpenter (2018) that Maraapunisaurus fragillimus differs from Amphicoelias altus and might represent a rebbachisaurid, the validity of the taxa itself is debatable (Lerzo et al. 2024, 2025; Van Der Linden et al. 2025).

The detailed reasoning for the inclusion or exclusion of each verified occurrence in our analyses, as well as notes pertaining to updates and/or changes to the data (e.g., taxonomic updates), can be found in the Suppl. material 1: “Database S1.xls”, tab ‘diplodocoidea_fullinfo’, under the “occurrence comments” column.

The resolution of time bins can vary depending on several factors, such as the length of time under consideration, the type of data being analyzed, and the need to generate sample sizes large enough to yield statistically significant results. Despite being commonly used, considering stage-length time bins at both global and regional scale may obscure events that are registered within a single stage and can introduce bias, as longer bins tend to show higher diversities (Rasmussen et al. 2019). This can be addressed by using time bins of equal duration, but this approach also has drawbacks of its own, e.g., reducing the number of data points available for statistical analysis, grouping occurrences of genera that never coexisted, and still being limited by the resolution at which occurrences are dated (Gibert and Escarguel 2017; Tennant et al. 2018; Dean et al. 2020). In summary, while stage-length bins allow for finer-scale diversity analysis, equal-length bins force sample data to be evenly distributed across time intervals (Tennant et al. 2016a, 2016b). Therefore, we followed the double-pronged approach of De Celis et al. (2020) and Tennant et al. (2016a), using both stage-length bins and equal-length bins for our analyses, and thus allowing for comparisons with other studies and exploring the impact of different spatiotemporal binning strategies on diversity curves.

Approximately equal-length time bins were generated using the ‘time_bins’ function in the palaeoverse R package (Jones et al. 2024), which combines adjacent geological intervals to achieve a target bin duration. Target durations between 6 and 12 Myr were evaluated, and 9 Myr was selected because larger bins would merge the Tithonian and Berriasian stages, potentially obscuring diversity changes across the Jurassic/Cretaceous transition, whereas smaller bins were effectively equivalent to stage-length bins, except for the merger of the Bajocian and Bathonian. The resulting bins have a mean duration of 9.17 Myr, ranging from 6.62 to 11.72 Myr. We also considered two distinct datasets when calculating diversity estimates: temporally constrained (i.e., including only occurrences where the temporal range fits into a single time bin) and temporally unconstrained (i.e., including occurrences that spanned more than one adjacent time bin). This methodology allows researchers to assess the impact of including temporally uncertain occurrences on biodiversity estimates (Tennant et al. 2016a; De Celis et al. 2020).

Suppl. material 1: tables S4, S5 provide the parameters used to construct the stage-length and equal-length time bins, respectively. The four occurrence datasets generated under the temporally constrained and temporally unconstrained approaches, considering both time-binning schemes, are presented in Suppl. material 1: tables S6–S9.

We conducted our diversity estimates at genus level as many dinosaur genera are monospecific and identification of specimens to species level is challenging (Robeck et al. 2000; Smith 2001; Cleary et al. 2018, 2020). Additionally, diversity estimates at genus level closely mirror those at a species level when the genera are mainly monospecific (Benton 2008; Barrett et al. 2009; Alroy 2010a, 2010b, 2010c; Butler et al. 2011; Mannion et al. 2015). However, results obtained at different taxonomic levels may lead to different interpretations of external factors influencing biodiversity, so findings must be interpreted with caution (Wagner et al. 2007; Hendricks et al. 2014; Wiese et al. 2016; Tennant et al. 2018). In the present case, this issue may be further complicated by the fact that some flagellicaudatan genera are not monospecific, including Diplodocus, Apatosaurus, Brontosaurus, Galeamopus Tschopp, Mateus & Benson, 2015, and Dicraeosaurus.

Raw taxonomic diversity is the direct count of different taxa present at a specific taxonomic level, which we calculated at the genus level per clade and subclade for each time bin. While commonly used to predict paleodiversity (Bakker 1978; Hunt et al. 1994; Sereno 1999; Weishampel and Jianu 2000), it is considered a biased and less reliable estimator particularly due to the heterogeneous sampling of fossil occurrence data (Alroy 2010a, 2010b, 2010c; Mannion et al. 2013; Starrfelt and Liow 2016; Benson et al. 2021).

While raw taxonomic diversity provides a simple count of fossil occurrences, sample standardization allows for the estimation of fossil diversity patterns based on a statistical approach (Alroy et al. 2008; Close et al. 2018; Tennant et al. 2018). Shareholder quorum subsampling (SQS) is a sample standardization method for estimating diversity based on a statistical approach that accounts for unequal sampling across bins (Alroy 2010a, 2010b, 2014; Chao and Jost 2012). It works by standardizing samples from different time bins to equal levels of completeness, referred to as quorum or coverage, in the relative frequency distributions of taxa per time bin (Close et al. 2018). The quorum is an objective measure of diversity-sample completeness, which can be estimated using several metrics, with Good’s u being the simplest and most widely used (Good 1953; Close et al. 2018; Benson et al. 2021). In practice, the SQS algorithm works by randomly drawing fossil occurrences from a given time bin and tracking the value of Good’s u as occurrences are drawn. Once the established quorum threshold is reached, taxonomic richness is recorded. This process is repeated over multiple iterations, with the final diversity value being the median of the values from all subsampling runs. Coverage-based methods like SQS have been widely applied to estimate patterns of vertebrate diversity in the past (Cleary et al. 2018, 2020; De Celis et al. 2020; Dunne et al. 2021, 2023).

To avoid issues with small data samples, we excluded time bins with fewer than five entries, following Cleary et al. (2018) and De Celis et al. (2020). We did not calculate SQS estimates at subclade level, as essential conditions for applying this methodology were not met. Although a quorum of 0.4 has been often considered sufficient for evaluating diversity patterns in fossil occurrence data (Alroy 2010a, 2010b; Close et al. 2018), we were able to conduct all our estimates at a quorum of 0.7. More details on this methodology are available in the Suppl. material 1: section S1.2.

Sedimentary bodies that yield a high amount of information via quantity and/or quality of preservation, also known as Lagerstätten, can distort paleodiversity metrics and our interpretations through the so-called Lagerstätten-effect (Seilacher 1970; Raup 1972). The Morrison Formation (Late Jurassic, USA) exhibits a remarkable high sauropod diversity and abundance (86.63% of global Jurassic Diplodocoidea occurrences) that can act as a Lagerstätten in this context, and this diversity has been frequently interpreted as reflecting a true biological signal (e.g., Tschopp et al. 2015, 2020, 2022; Whitlock et al. 2018; Foster et al. 2019; Maidment and Muxworthy 2019; Woodruff 2019; Foster 2020; Mannion et al. 2021; Boisvert et al. 2025; Van Der Linden et al. 2025).

In order to quantitatively analyze if the obtained paleodiversity estimates are biased due to a Lagerstätten effect exerted by the Morrison Formation, a GLS (Generalized Least Squares) modeling approach was employed following previous work (e.g., Benson and Butler 2011; De Celis et al. 2021). In all models, log-transformed paleodiversity (raw generic diversity and SQS estimates) per stage was treated as the dependent variable, whereas the variables ‘Morrison’ and ‘Duration’ were included as independent variables. Stage-length data was selected for these analyses because these provided the most complete and continuous paleodiversity curves, except for the SQS constrained approach diversity estimate, which had several gaps and could not be analysed with this GLS approach. The stages duration (in millions of years) was incorporated into the GLS models to account for uneven temporal spans among stages. The ‘Morrison’ variable is a proxy to the presence or absence of the Morrison Formation diplodocoid occurrence data in each stage, and is weighted to account for the uneven distribution of these occurrences across stages relative to the total diplodocoid occurrences recorded within the formation. This weighting provides a more nuanced proxy of the potential Lagerstätten effect through time in comparison to classical presence/absence approaches. When required, temporal autocorrelation was accounted for by incorporating adequate corARMA correlation structures within the GLS framework.

Suppl. material 1: table S10 comprises the data set used for our GLS models, including all variables considered herein and their corresponding sources.

Diplodocoidea occurrences are latitudinally constrained between 60°N and 60°S (Fig. 2). Their latitudinal distribution was more restricted during the Middle Jurassic (Northern Hemisphere only; Fig. 2A) and Late Jurassic (20° to 40° latitude in both hemispheres but absent from the 20°N to 30°S belt; Fig. 2B).

The only exception to this pattern is a putative Bathonian–Callovian flagellicaudatan occurrence from the equatorial region (Fig. 2A), based on fragmentary pedal elements from Mexico’s Otlaltepec Formation (Rivera-Sylva and Espinosa-Arrubarrena 2020). However, this taxonomic assignment is currently considered uncertain (Van Der Linden et al. 2025). From the Valanginian onwards, they occupied a broader latitudinal range that included equatorial regions (20°N to 20°S; Fig. 2C, D).

Most dicraeosaurid fossil occurrences are from Gondwanan landmasses (20°S to 50°S), except for those from the Upper Jurassic Morrison Formation (34°N to 48°N; Foster 2020; Tschopp et al. 2022), and for Lingwulong shenqi Xu, Upchurch, Mannion, Barrett, Regalado-Fernandez, Mo, Ma & Liu, 2018 (Zhiluo Formation, Bathonian–Callovian, China; Xu et al. 2018; Ren et al. 2023). From the Tithonian onwards, Dicraeosauridae is constrained to the Southern Hemisphere, specifically the South American domain. Although present in both Laurasian and Gondwanan landmasses, Diplodocidae occurrences are mostly limited to the 45°N to 20°N and 35°S to 50°S latitudinal ranges, respectively, with the exception of two indeterminate diplodocids from the Salvador (Berriasian–Barremian) and Marfim (Valanginian–Hauterivian) formations of Bahia, Brazil, recently discovered in the collections of the Natural History Museum (NHM; London, United Kingdom; Bandeira et al. 2025). The bulk of the diplodocid fossil record spans the Oxfordian–Tithonian, except for a few Berriasian–Barremian Southern Hemisphere occurrences (40°S to 50°S). In contrast, Early Cretaceous rebbachisaurids (Fig. 2C) occur both within the 20°N to 30°S latitudinal belt and throughout the broader 50°N to 50°S range previously occupied by Jurassic flagellicaudatans. This pattern persists in the Middle Cretaceous in the Southern Hemisphere, whereas Northern Hemisphere occurrences appear restricted to latitudes below 25°N (Fig. 2D).

Diplodocoidea occurrences show broadly similar patterns in both the constrained and unconstrained datasets, regardless of whether stage-length or 9 Myr bins are used. In all cases, occurrence counts remain low and relatively stable from the Bathonian to the Turonian, with a pronounced Late Jurassic peak (Fig. 3).

Differences between the datasets are primarily limited to the total number of occurrences, the magnitude of the Late Jurassic peak, and the timing of the first and last diplodocoid occurrences (Fig. 3). Because the Diplodocoidea curve combines occurrences of Flagellicaudata, Rebbachisauridae, and indeterminate diplodocoids, while the Flagellicaudata curve includes Dicraeosauridae, Diplodocidae, and indeterminate flagellicaudatans, the overall diplodocoid pattern largely reflects trends observed in its constituent groups across all binning methods and datasets (Fig. 3).

All occurrence curves recover a Bathonian age for the earliest global diplodocoid occurrences (Fig. 3B–D), except for the constrained dataset using stage-length bins, which places the first diplodocoid (and oldest flagellicaudatan) occurrences in the Oxfordian (Fig. 3A). Diplodocoid occurrences peak during the Late Jurassic, spanning the Kimmeridgian–Tithonian in the constrained dataset (Fig. 3A, C) and the Oxfordian–Tithonian in the unconstrained dataset (Fig. 3B, D). This peak is driven largely by abundant diplodocoid material from the Morrison Formation (see: Foster 2020; Tschopp et al. 2022), supplemented by a smaller number of occurrences from the Iberian Peninsula (e.g., Bonaparte and Mateus 1999; Mannion et al. 2012a; Mocho et al. 2017) and from the Tendaguru Formation (e.g., Janensch 1914, 1925, 1929a, 1929b, 1935, 1961; Remes 2006, 2007, 2009). Peak magnitude varies substantially between datasets, reaching approximately 170 occurrences in the constrained dataset and over 400 in the unconstrained dataset (Suppl. material 1: tables S12, S13 and Suppl. material 1: tables S14, S15, respectively). A sharp decline occurs across the Jurassic/Cretaceous transition (Tithonian–Berriasian or Berriasian–Valanginian, depending on the binning scheme), particularly among diplodocids. Thereafter, diplodocoid occurrences remain relatively stable, while rebbachisaurids first appear in the fossil record. Consequently, all three major diplodocoid subclades co-occur until the Hauterivian in all analyses, except in the constrained dataset with stage-length bins, where rebbachisaurids first appear in the Barremian. From the Barremian onward (or Hauterivian–Barremian bin), only rebbachisaurids are recorded, with occurrences increasing gradually until the Cenomanian. In the constrained dataset, this marks the final occurrence of Diplodocoidea, implying extinction during the Cenomanian, whereas the unconstrained dataset records a subsequent decline from the Cenomanian to the Turonian, extending the diplodocoid fossil record into the Turonian (Fig. 3A, Cversus Fig. 3B, D).

Results obtained for diversity estimates are mostly consistent across time-binning procedures, but not across datasets. Here we focus on results from the unconstrained dataset using stage-length time bins, as the constrained dataset includes substantially fewer occurrences and tends to obscure broader diversity patterns that are more clearly expressed in the unconstrained data. Detailed results for raw genus-level taxonomic diversity (Suppl. material 1: tables S17–S20) and SQS diversity estimates (Suppl. material 1: tables S21–S24) of Diplodocoidea and its subclades are provided in the Suppl. material 1.

Raw genus-level diversity patterns of Diplodocoidea vary among subclades, datasets, and binning strategies, reflecting both sampling biases and the sensitivity of richness estimates to temporal aggregation (Fig. 4). Despite these differences, all approaches consistently recover two major macroevolutionary features: a pronounced Late Jurassic diversity peak and a second peak in the Cenomanian, separated by a diversity trough across the Jurassic/Cretaceous transition.

The Late Jurassic diversity peak occurs in the Kimmeridgian across all analyses, followed by a Tithonian decline (Fig. 4). It is primarily driven by flagellicaudatan diplodocoids (Diplodocidae and Dicraeosauridae), which dominate the pre-Barremian diversity in the unconstrained dataset (Fig. 4B, D), where Diplodocoidea closely track flagellicaudatan diversity until the Hauterivian–Barremian interval. The peak magnitude varies among datasets, with unconstrained analyses recovering higher diversity levels (Fig. 4B, D), but the overall pattern is seemingly consistent (Fig. 4). The 9 Myr binning approach merges earlier Jurassic stages and therefore smooths the increase from the Bathonian–Oxfordian interval, while reducing the apparent sharpness of the Kimmeridgian peak too (Fig. 4C, D). There is a marked decline in diversity following the Jurassic/Cretaceous boundary across all approaches. Stage-length binning recovers a more seemingly complex Early Cretaceous pattern, with alternating decreases and increases in diversity from the Aptian through the Barremian (Fig. 4A, B). Equal-length binning masks this fluctuation, by merging the Hauterivian and Barremian as a single interval, whilst suppressing also the evidence of low Hauterivian diversity (Fig. 4C, D). These patterns are underpinned by subclade dynamics.

Diplodocidae diversity peaks in the Kimmeridgian and declines sharply across the Jurassic/Cretaceous transition (Fig. 4). It then persists at low diversity into the earliest Early Cretaceous, disappearing by the Barremian (Fig. 4A, B) or Hauterivian–Barremian (Fig. 4C, D), depending on the binning methodology used. Dicraeosauridae also peaked in diversity during the Cretaceous but at a smaller scale (Fig. 4), and maintained a low, stable diversity into the Early Cretaceous, before disappearing from the fossil record in the Barremian (Fig. 4). Dicraeosaurids have a diversity increase in the Valanginian that may contribute to Early Cretaceous overall Diplodocoidea diversity, but it does not alter the overall declining trend of post-Tithonian Flagellicaudata diversity (Fig. 4A, B). Rebbachisauridae first appears in the Early Cretaceous (either in the Berriasian or in the Barremian; Fig. 4B–D and Fig. 4A, respectively), dominating Diplodocoidea diversity from thereafter. In the constrained dataset, they remain low in diversity until they reach a peak in the Cenomanian, after which they disappear (Fig. 4A, C). Contrastingly, the unconstrained dataset (Fig. 4B, D) shows a more gradual increase from the Hauterivian or Berriasian–Valanginian onwards, which stabilizes in the Aptian–Albian, and reaches a Cenomanian maximum before declining in the Turonian or Turonian–Coniacian. Despite the temporal difference, both datasets yield comparable peak diversity levels for Late Jurassic Diplodocidae and Cenomanian Rebbachisauridae.

The diversity patterns recovered for the Early Cretaceous seem to be the most affected by methodological differences. The equal-length time-binning approach obscures short-term dynamics by merging stages, which is particularly notable in the Hauterivian–Barremian, and can exacerbate or smooth diversity levels in combined intervals. Similarly, the use of the constrained dataset yields tendentially lower and more fragmented diversity reconstructions (namely for the Early Cretaceous), including a more marked post-Jurassic diversity drop, whereas the use of the unconstrained dataset results in smoother diversity transitions and a seemingly more sustained Early Cretaceous diversity oscillation.

In general, all analyses converge in the existence of a Late Jurassic diplodocoid diversity peak dominated by flagellicaudatans, along with a Cenomanian diversity peak driven by rebbachisaurids, indicating a shift in clade dominance through time.

For the unconstrained dataset, using stage-length time-bins, we estimated Shareholder Quorum Subsampling (SQS) values for the interval spanning from the Oxfordian to the Valanginian, and again from the Aptian to the Turonian (Fig. 5B).

The highest Jurassic SQS value was recorded in the Kimmeridgian (3.27 at a quorum level of 0.7), followed closely by the Tithonian (3.06) and the Oxfordian (3.02) at the same quorum level (Fig. 5B). SQS estimates for the Berriasian and Valanginian are similar but lower than those for the Late Jurassic, with values of 1.21 and 1.27 respectively, at a quorum of 0.7 (Fig. 5B). Due to insufficient data, no estimates could be calculated for the Hauterivian, resulting in a gap in our understanding of diversity trends between the Valanginian and the Barremian (Fig. 5B). In the Barremian, we recover an SQS value of 1.90 for Diplodocoidea diversity (Fig. 5B). SQS values peak in the Aptian, where we observe the highest value in our dataset (5.73 at quorum 0.7), followed by a decline in the Albian (Fig. 5B). Diversity estimates for the Albian, Cenomanian, and Turonian at the same quorum level are 3.37, 3.42, and 3.37, respectively (Fig. 5B). After the Turonian, no occurrences of Diplodocoidea are recorded in the fossil record.

SQS estimates for Diplodocoidea differ between datasets, primarily reflecting the limited temporal coverage achievable under the constrained approach. In the constrained dataset, SQS values could be calculated for only three intervals, the Kimmeridgian, Tithonian, and Cenomanian, primarily due to extensive bin-sharing of occurrences and a high proportion of singleton taxa (Fig. 5A, C). This pattern produces a pronounced diversity peak in the Kimmeridgian, followed by a sharp decline in the Tithonian and an isolated Cenomanian estimate. However, because adjacent intervals lack estimates, it is unclear whether the Cenomanian value reflects a genuine diversity peak or an artefact of incomplete sampling (Fig. 5A, C). In contrast, the unconstrained dataset provides more continuous and reliable estimates, with the exception of the Hauterivian and Hauterivian–Barremian intervals (Fig. 5B, D). These results indicate relatively stable diversity throughout the Late Jurassic and a major turnover across the Jurassic/Cretaceous boundary, rather than a pronounced decline within the Late Jurassic itself. Diversity reaches its highest levels in the Aptian, declines gradually through the Albian, and remains relatively stable during the Cenomanian and Turonian, after which rebbachisaurids (and therefore diplodocoids) disappear from the fossil record (Fig. 5B, D).

Comparison between approaches highlights two main uncertainties. First, Late Jurassic dynamics differ: the constrained dataset (using stage-length time-bins) suggests a sharp Kimmeridgian–Tithonian drop, whereas unconstrained estimates support sustained high diversity across the Late Jurassic (Fig. 5A and Fig. 5B, D, respectively). Second, the timing of peak diversity also changes: constrained results are influenced by isolated estimates (Fig. 5A, C), while unconstrained analyses consistently identify an Early Cretaceous (Aptian) maximum followed by decline (Fig. 5B, D). Overall, 9 Myr bin results broadly agree with these trends but smooth finer-scale variation (Fig. 5C, D and Fig. 5A, B, respectively).

Residual diagnostics confirmed that model assumptions were adequately met for the RGDc and RGDu (raw generic diversity, constrained dataset, and raw generic diversity, unconstrained dataset, respectively) models, with no significant violations of normality or homoscedasticity. Temporal autocorrelation was detected in the RGDc model and was accounted for by incorporating a corARMA(1,1) structure in the GLS analysis. The SQSu model showed a significant Shapiro-Wilk test result (p-value < 0.05), but Q-Q plot inspection suggested this was driven by two extreme tail observations rather than a systematic departure from normality, and no heteroscedasticity or autocorrelation were detected, so GLS assumptions were considered reasonably satisfied.

Both the RGDc and RGDu models (Table 1) identified the Morrison Formation as a significant positive predictor of paleodiversity (p-value < 0.001 and p-value = 0.011, respectively), while stage duration was only significant in the RGDc model (p-value < 0.001). Neither predictors, Morrison Formation nor duration, were significant in the SQSu model (p-value > 0.05 in both cases). Cox & Snell pseudo R² values (Table 1) indicate the improvement in model fit relative to the corresponding null model in each case, with higher values reflecting a greater contribution of the predictors to explaining variation in paleodiversity. The constrained raw diversity model showed the highest improvement over its null counterpart (pseudo R² = 0.848), followed by the unconstrained model (0.574), whereas the sampling-standardized model showed a substantially lower improvement (0.062). This decrease from raw approaches (RGDc and RGDu) to the sampling-standardized approach (SQSu) indicates that the explanatory contribution of the predictors is reduced once these biases are accounted for via sampling standardization.

Recent literature suggests that the geographic distribution of sauropodomorphs might be related to their physiology since they are more sensitive to temperature than other dinosaurs, which shaped their biodiversity patterns more strongly (Chiarenza et al. 2022; Dunne et al. 2023). We find that occurrences of Diplodocoidea are latitudinally constrained between 60°N and 60°S, corroborating the sauropod restrictions found by Mannion et al. (2012b), Poropat et al. (2016), Averianov and Lopatin (2019), Chiarenza et al. (2022), and Dunne et al. (2023). Namely, the marked restriction of this clade to lower latitudinal levels when compared to other non-avian dinosaur clades.

Dicraeosaurids and diplodocids show similar latitudinal distribution patterns and synchronized occurrences. In some cases, both subclades can be found in the same lithostratigraphic level, or even in the same site, e.g., Bajada Colorada Formation main site (Gallina et al. 2014, 2019) and Morrison Formation Marsh-Felch Quarry (Whitlock and Mantilla 2020), supporting the fact that they coexisted in the same terrestrial ecosystems. The mechanisms underlying this high sauropod diversity in the same areas remain poorly understood. However, previous research suggests the existence of a niche-partitioning system, showing distinct ecomorphospace distributions for different sauropod clades (Button et al. 2014; Poropat et al. 2021; Wyenberg-Henzler 2022; Filippi et al. 2024).

The currently known rebbachisaurid fossil record exhibits a lower degree of morphological diversity than that documented for flagellicaudatans (Sereno et al. 2007; Wilson and Allain 2015). Rebbachisauridae seemingly achieved a broader, less-restricted latitudinal range and therefore a wider geographical distribution than Dicraeosauridae and Diplodocidae. Somphospondyli, a clade of titanosauriform sauropods that lived from the Late Jurassic until the end of the Late Cretaceous and which display a type of pneumaticity like that of Rebbachisauridae (Wedel 2003), also appear to have reached a broad latitudinal distribution, possibly during the early Early Cretaceous (D’Emic 2012).

The discoveries of Maraapunisaurus (Carpenter 2018) and Xenoposeidon (Taylor and Naish 2007) complicate the idea of a Gondwanan origin for Rebbachisauridae, supported by Canudo et al. (2018) and Whitlock and Mantilla (2020). If their taxonomic attribution is supported, a North American origin for the clade can be suggested, with possible dispersal into Europe and later into Gondwanan landmasses during the Late Jurassic to Early Cretaceous, which aligns with the wide latitudinal range we report here for the clade (Fig. 2). This hypothesis is plausible considering: (1) a Late Jurassic faunal exchange between North America and Europe despite the narrow epicontinental sea, possibly via the “Viking Corridor” (Escaso et al. 2007; Bardet et al. 2014); and (2) the establishment of the land corridor known as the Apulian Route (Gheerbrant and Rage 2006), connecting southern Europe to northern Africa that, albeit not continuous, could have enabled the dispersion of somphospondylans and rebbachisaurids between Laurasia and Gondwana during the Early Cretaceous (Fernández-Baldor et al. 2011; Lindoso et al. 2019; Mocho et al. 2019). In contrast, Fanti et al. (2013, 2015) and Lerzo et al. (2025) suggest that South America might have been the origin of the Rebbachisauridae radiation based on their time-calibrated phylogeny, which is currently considered the most likely scenario. Nevertheless, Rebbachisauridae is recognized as a ghost lineage dating back to the Middle Jurassic (Xu et al. 2018; Averianov and Zverkov 2020) until its first unambiguous appearance in the fossil record in the Early Cretaceous, a conclusion supported by phylogenetic analyses in other studies (Whitlock 2011a; Gallina et al. 2014; Wilson and Allain 2015; Carpenter 2018; Bajpai et al. 2023; Mannion and Moore 2025; Van Der Linden et al. 2025).

Flagellicaudatans appeared while Gondwana and Laurasia were still connected (Scotese 2016), enabling them to achieve a broad geographic distribution (Upchurch et al. 2004; Mannion and Upchurch 2010; Xu et al. 2018; Dai et al. 2022; Ren et al. 2023) which we recover here for the Late Jurassic and particularly for Diplodocidae (Fig. 2). While some authors proposed that the initial radiation of Neosauropoda into diplodocoid and macronarian lineages occurred during the Middle to Late Jurassic (Wilson and Sereno 1998; Upchurch et al. 2004; Upchurch and Barrett 2005; Whitlock 2011a), others propose that these lineages may have already split by the early Middle Jurassic (Xu et al. 2018; Whitlock and Mantilla 2020; Dai et al. 2022; Salgado et al. 2022; Ren et al. 2023; Mannion and Moore 2025).

Lingwulong shenqi (Xu et al. 2018) is the only non-Morrison Formation dicraeosaurid that has been found in Laurasian landmasses. Xu et al. (2018) assigned the taxon to the Yan’an Formation of Lingwu, China, which they regard as late Toarcian–Bajocian in age; this chronology was later refined to Aalenian–Bajocian by Zhang et al. (2020). Under this interpretation, Lingwulong represents both the earliest known neosauropod and the oldest known diplodocoid. However, You et al. (2019) argued that the Lingwu dinosaur fauna instead originated from the Bathonian–Callovian Zhiluo Formation, a reassignment subsequently adopted by Bajpai et al. (2023), Ren et al. (2023), and us. Consequently, the oldest known diplodocoid occurrences are instead either the putative Flagellicaudata indet. from the Bathonian–Callovian Otlaltepec Formation (Rivera-Sylva and Espinosa-Arrubarrena 2020) or the indeterminate diplodocoid from the Callovian–Oxfordian Podosinki Formation, which may also be referrable to Flagellicaudata (Averianov and Zverkov 2020; Van Der Linden et al. 2025).

The few occurrences in our initial dataset dated from the late Middle Jurassic can be taxonomically assigned to distinct subclades of Diplodocoidea: Flagellicaudata, Diplodocidae, and Dicraeosauridae (Fig. 2). The initial radiation of diplodocoid dicraeosaurids in the Middle Jurassic of Western Asia/Europe and their subsequent dispersal to Africa and the Americas proposed by Averianov and Zverkov (2020) corroborates the biogeographic scenario we recovered, and challenges the geographic origin of the diplodocoid clade in North America, Europe, or Africa (e.g., Mannion and Upchurch 2011; Whitlock 2011a; Gallina et al. 2014). However, the phylogenetic position of Lingwulong shenqi has also recently been questioned by Mannion and Moore (2025), whose analyses recovered the taxon outside Dicraeosauridae, potentially as a stem flagellicaudatan. The authors argued that this revised placement is more compatible and parsimonious with its proposed Middle Jurassic age (Mannion and Moore 2025). If this interpretation is correct, confirmed dicraeosaurid occurrences would be restricted to the Late Jurassic onward. The broad geographic distribution of possible Middle Jurassic flagellicaudatan occurrences recovered herein may therefore indicate a Eurasian distribution of the clade during the Bathonian–Callovian (Mannion and Moore 2025). In general, and broadly, our results support: (1) a Middle Jurassic, or possibly earlier, origin of Diplodocoidea (Xu et al. 2018; Whitlock and Mantilla 2020; Gallina et al. 2022; Van Der Linden et al. 2025); and (2) an origin/initial diversification of Neosauropoda as early as the late Early Jurassic (Dai et al. 2022; Bajpai et al. 2023; Mannion and Moore 2025). Such a scenario would also be consistent with dispersal during intervals of relatively low sea level, particularly during the Bathonian or earlier (Kocsis and Scotese 2021; Dai et al. 2022; Van Der Meer et al. 2022). Notably, our paleocoordinate reconstructions for the two Middle Jurassic Dicraeosauridae occurrences in our dataset place them in Western Asia and in present-day India (Fig. 2A).

The diplodocoid fossil record reveals significant differences in the distribution of its subclades along the Jurassic/Cretaceous transition, with both Dicraeosauridae and Diplodocidae essentially disappearing from the Northern Hemisphere from the Berriasian onwards, and their occurrences becoming concentrated in areas at notably low latitudes, particularly in South America (Fig. 2). However, the South American flagellicaudatan lineages do not share a close phylogenetic relationship with the Late Jurassic faunas of North America and Iberia (Xu et al. 2018; Gallina et al. 2019; Whitlock and Mantilla 2020). According to Tschopp et al. (2015) and Xu et al. (2018), North American flagellicaudatan lineages appear to have failed to disperse and likely went extinct by the latest Late Jurassic (Fig. 2B, C). The fossil record indicates that South American flagellicaudatans persisted until the Barremian (Ibiricu et al. 2013; Windholz et al. 2023, 2024), but unlike Rebbachisauridae, these South American lineages do not appear to have dispersed northward following the Jurassic/Cretaceous transition (Fig. 2C, D).

Using the Paleobiology Database data visualization tool (https://paleobiodb.org/navigator/, PBDB; Uhen et al. 2023), we compared the global paleogeographic distributions of Dinosauria and Sauropoda during the Middle Jurassic, Late Jurassic, and Early Cretaceous. The distributions broadly overlap across all three intervals, although sauropod occurrences are generally restricted to lower paleolatitudes relative to Dinosauria as a whole, consistent with the findings of Chiarenza et al. (2022) and Dunne et al. (2023). Accordingly, the following discussion focuses on comparisons between the paleogeographic distributions of Sauropoda and Diplodocoidea. Note that the general scarcity of diplodocoids in eastern Asia, apart from Lingwulong, and its paleobiogeographic implications for the radiation of Neosauropoda have already been discussed by Xu et al. (2018), Dai et al. (2022), Ren et al. (2023), and Mannion and Moore (2025), and are therefore not considered further here.

The absence of Diplodocoidea from many Middle and Upper Jurassic formations does not necessarily reflect a true absence of the clade. In many cases, the fossil record consists only of fragmentary sauropod remains, footprints, or ootaxa that lack the diagnostic features needed to exclude diplodocoid affinities. However, several formations preserve sufficiently diagnostic sauropod taxa that appear genuinely distinct from Diplodocoidea, including the Great Oolite Group and Oxford Clay Formation (United Kingdom), Irhazer II Formation (Niger), Cañadón Asfalto Formation (Argentina), Tiourarén Formation (Niger), Walloon Coal Measures Formation (Australia), and Arcabuco Formation (Colombia). The phylogenetic positions of some key taxa from these units remain debated. Cetiosaurus Owen, 1841 is generally regarded as a non-neosauropod eusauropod (Wilson 2002; Ksepka and Norell 2010; Holwerda et al. 2021; Gomez et al. 2024; Mannion and Moore 2025), whereas Cetiosauriscus Heune, 1928 has recently been recovered within Flagellicaudata by some analyses (Mannion et al. 2019a, 2021; Mannion and Moore 2025), despite previous work placing it outside Diplodocoidea (Heathcote and Upchurch 2003; Rauhut et al. 2005; Tschopp et al. 2015). Likewise, Jobaria Sereno, Beck, Dutheil, Larsson, Lyon, Moussa, Sadleir, Sidor, Varricchio, Wilson & Wilson, 1999 and Spinophorosaurus Remes, Ortega, Fierro, Joger, Kosma, Ferrer, Idé & Maga, 2009 occupy uncertain positions near the base of Eusauropoda or Neosauropoda (Wilson 2002; Remes 2009; Carballido et al. 2011; D’Emic 2012; Nair and Salisbury 2012; Mocho et al. 2013; Xing et al. 2015; Mannion et al. 2019a; Ren et al. 2023), leaving open the possibility that some African Jurassic taxa may be more closely related to Diplodocoidea than previously recognized. Nevertheless, evidence from Argentina and central Africa suggests that the scarcity of diplodocoids in parts of Gondwana during the Middle Jurassic may represent a genuine paleobiogeographic signal (Gallina et al. 2014; Xu et al. 2018; Whitlock and Mantilla 2020; Dai et al. 2022; Gomez et al. 2024), although this interpretation has been questioned by Mannion et al. (2019a). The prevalence of non-diplodocoid eusauropods in this region, together with the presence of the Early Jurassic Central Gondwanan Desert, supports the hypothesis that geographic barriers influenced early sauropod dispersal and diversification (e.g., Bonaparte 1996; Remes 2009; Gallina et al. 2014; Xing et al. 2015; Xu et al. 2018; Mannion et al. 2019a, 2019b; Whitlock and Mantilla 2020; Poropat et al. 2021; Dai et al. 2022; Ren et al. 2023; Gomez et al. 2024; Van Der Linden et al. 2025). In contrast, the apparent absence of diplodocoids in Ethiopia, Yemen, and Colombia is more difficult to interpret because the fossil record is sparse and taxonomically poorly resolved.

By the Early Cretaceous, sauropods had achieved a nearly global distribution, including low-latitude regions that yielded few records during the Middle and Late Jurassic. Rebbachisaurid diplodocoids also became widespread but remained absent or extremely rare in three major regions: high-latitude Australia, Central and North America, and eastern Asia. In Australia, the available evidence consists largely of trackways, some of which could belong to diplodocoids, although diagnostic titanosauriform tracks are also present (Salisbury et al. 2016). In Central America and northern South America, the sparse and fragmentary fossil record prevents firm conclusions regarding diplodocoid presence (Jímenez-Moreno et al. 2012; Carballido et al. 2015; Mannion et al. 2017). In contrast, the Early Cretaceous sauropod faunas of both eastern and western North America are dominated by macronarians, particularly titanosauriforms, with additional turiasaur occurrences in the west (Tidwell et al. 1999, 2001; Wedel et al. 2000; Taylor et al. 2011; Woodruff 2012; D’Emic 2013; Britt et al. 2017; Royo-Torres et al. 2017; Yoshida et al. 2026). Although many remains are assigned only to Sauropoda indet., all diagnostically identifiable taxa belong to lineages outside Diplodocoidea, suggesting that the clade may have been genuinely absent from North America during this interval. A similar pattern is observed in eastern Asia. Despite a perhaps less sampled fossil record than North America, the region preserves diverse and diagnostic sauropod assemblages dominated by euhelopodids, titanosaurs, and other titanosauriforms (Martin et al. 1994; Allain et al. 2012; D’Emic and Foreman 2012; Averianov and Skutschas 2017; Averianov et al. 2018; Sethapanichsakul et al. 2026). The only reported diplodocoid occurrence is an undescribed and highly fragmentary specimen from Thailand whose identification remains uncertain (Shimizu et al. 2016). Consequently, there is currently no convincing evidence for a significant Early Cretaceous diplodocoid presence in eastern Asia, suggesting that the clade either failed to disperse into the region or remained exceptionally rare there.

Overall, diversity curves derived from the constrained dataset tend to obscure broader patterns that are more clearly recovered by the unconstrained dataset, particularly when stage-length bins are used. This likely reflects the fact that many occurrences, especially from the Late Jurassic, cannot be confidently assigned to a single stratigraphic stage. Restricting analyses to narrowly constrained records therefore removes a substantial amount of information and increases sensitivity to temporal binning. Consequently, we consider diversity estimates based on the unconstrained dataset and stage-length bins to provide the most reliable representation of diplodocoid diversity dynamics.

The strong correspondence between diversity peaks in the raw taxonomic metrics (Fig. 4) and the concentration of occurrences from the Morrison Formation (86.63% of all Jurassic records and 71.24% of the entire dataset; Fig. 3) suggests that a Lagerstätten effect has substantially influenced inferred raw diversity patterns. This interpretation is supported by the GLS analyses (Table 1), which identify a significant effect of the Morrison Formation on raw genus-level diversity estimates, particularly under the constrained approach. In contrast, this effect becomes weaker under the unconstrained treatment and disappears in the sampling-standardized SQSu model. This absence of a Morrison effect in SQSu is especially important, as it indicates that once sampling is standardized, the apparent Morrison-driven diversity peak is no longer recovered and, hence, SQSu results provide a more robust basis for reconstructing paleodiversity patterns in this group. Although stage duration values are identical across the three analysed approaches (RGDc, RGDu and SQSu), stage duration was significant only in the constrained raw diversity model (RGDc). This suggests that the constrained treatment of occurrences makes raw diversity estimates more sensitive to differences in stage length. This is consistent with the broader warning that geological time-bin duration can bias apparent diversity estimates (Raup 1972), although in this case this effect appears to be mediated by the exclusion of occurrences spanning several stages rather than by stage duration alone. The absence of these effects in the SQSu analyses suggests that raw metrics primarily reflect sampling artifacts rather than genuine biological signals. Together with the marked reduction in Cox and Snell pseudo-R² values from raw to subsampled approaches (Table 1), reflecting lower improvement over their corresponding null models, this supports the interpretation that much of the explanatory signal recovered in raw diversity metrics is driven by sampling structure. These results highlight the importance of accounting for temporal binning, data-treatment strategy, and uneven sampling intensity when reconstructing diversity in deep-time.

Our findings therefore provide only partial support for a possible Late Jurassic diversity maximum in Diplodocoidea (e.g., Wilson 2002; Mannion et al. 2011; Whitlock 2011a; Xu et al. 2018). Raw diversity estimates recover a pronounced Kimmeridgian diversity peak followed by a decline in the Tithonian, a pattern also observed in Diplodocidae and Dicraeosauridae (Fig. 4). However, sampling-standardized analyses indicate a different trajectory. Although the constrained SQS results may somewhat reproduce the raw diversity pattern (Fig. 5A, C), the unconstrained SQS estimates consistently recover relatively stable diversity throughout the Late Jurassic (Oxfordian–Tithonian), with the principal decline occurring at the Jurassic/Cretaceous boundary rather than within the Late Jurassic itself (Fig. 5B, D). These patterns are broadly compatible with previous hypotheses concerning the radiation of Flagellicaudata and subsequent changes in their distribution and diversity (e.g., Remes 2006; Mannion et al. 2011; Whitlock 2011a; Gallina et al. 2014; Rauhut et al. 2015; Xu et al. 2018; Whitlock and Mantilla 2020; Dai et al. 2022). Nevertheless, the sparse fossil record and poor temporal resolution between the Berriasian and Hauterivian (Fig. 3) limit our ability to determine the exact timing and magnitude of diversity change across the Jurassic/Cretaceous transition. Although the mechanisms driving these patterns remain uncertain, the unconstrained SQS results (Fig. 5B, D) are consistent with a substantial faunal turnover during this interval, as proposed by Tennant et al. (2016b, 2017) and Benton (2023a, 2023b). Similar conclusions have been drawn from broader sauropod diversity analyses, which also recover a stronger Jurassic/Cretaceous decline than previously recognized (Cashmore et al. 2020). At the same time, the contrasting outcomes produced by constrained analyses highlight the continuing uncertainty surrounding the severity of this turnover (Énay 2020; Allain et al. 2022).

Our results also support the hypothesis of a second Diplodocoidea diversity maximum during the Early Cretaceous, associated with the radiation of Rebbachisauridae (Whitlock 2011a; Haluza et al. 2012; Ibiricu et al. 2013; Lerzo et al. 2024, 2025; Bellardini et al. 2025). However, whereas raw taxonomic diversity estimates place this peak in the Cenomanian (Fig. 4B, D), the sampling-standardized analyses indicate that it occurred earlier, during the Aptian (Fig. 5B, D). Notably, Aptian SQS values exceed those recovered for the Kimmeridgian, raising the possibility that Diplodocoidea reached their greatest diversity during the Early Cretaceous rather than the Late Jurassic. This interpretation should be treated cautiously because Aptian SQS estimates could only be calculated using the unconstrained dataset, potentially inflating the magnitude of the peak. Nevertheless, the available evidence suggests that Rebbachisauridae may have attained higher diversity levels than Flagellicaudata at their respective maxima. If correct, this pattern would contrast with broader sauropod diversity trends, which indicate that sauropod diversity did not fully recover following the Jurassic/Cretaceous transition (Cashmore et al. 2020). The causes of this apparent rebbachisaurid radiation remain uncertain. As dominant megaherbivores, sauropods likely exerted substantial ecological influence on terrestrial ecosystems, shaping vegetation structure and broader community dynamics (Gordon and Prins 2019; Pringle et al. 2023). Although direct comparisons with modern herbivores are limited by the unique biology of sauropods (Bates et al. 2016), environmental changes affecting plant communities may have played an important role in their diversification. In particular, while the radiation of angiosperms does not appear to have increased overall sauropod diversity (Butler et al. 2009), it may have contributed to the ecological success of Rebbachisauridae, a possibility that remains largely unexplored. Consequently, the Early Cretaceous diversity peak appears to represent a genuine biological signal, but its underlying drivers remain unresolved.

Rebbachisauridae outlasted both Dicraeosauridae and Diplodocidae, persisting until shortly after the Turonian (Fig. 5B). In contrast to flagellicaudatans, rebbachisaurid sauropods appear to have remained abundant and diverse until their extinction (Fig. 5), suggesting that their disappearance was not preceded by a prolonged decline. This pattern has previously been noted by Coria and Salgado (2005), Bellardini et al. (2025), and Simón and Salgado (2025). Their eventual extinction may have been linked to environmental perturbations associated with the Cenomanian/Turonian boundary event, including global warming, widespread oceanic anoxia, and sea-level rise (Tennant et al. 2016b; Keller et al. 2021). These changes likely altered key terrestrial habitats and may have reduced the geographic and ecological range available to rebbachisaurids. At the same time, the diversification of Lithostrotia during the late Early Cretaceous (Upchurch et al. 2004; Zaher et al. 2011), together with the expansion of other large herbivorous dinosaur groups such as hadrosauriforms and ornithopods (Butler et al. 2009; Kobayashi et al. 2019), may have intensified ecological competition. Evidence for convergent cranial and locomotor adaptations in Rebbachisauridae and Lithostrotia suggests that these clades occupied similar ecological niches and may have exploited comparable food resources (Filippi et al. 2024). Such overlap could have promoted competitive replacement, contributing to the decline of rebbachisaurids. The subsequent diversification of several titanosaurian lineages following the reduction of rebbachisaurid diversity is consistent with this interpretation and has been proposed as evidence of ecological niche replacement (Filippi et al. 2024).

To ensure a transparent discussion of our results, it is important to highlight the major challenges encountered during our analysis and discuss possible caveats in our interpretations. One key issue is the temporal resolution of Diplodocoidea occurrences. Most occurrences from Early Cretaceous formations are poorly constrained in time. These occurrences are often assigned to multiple time bins, which leads to their exclusion from the constrained dataset. Specifically, occurrences from the Berriasian to Aptian stages are mainly found in Gondwanan landmasses, with many originating from geological units lacking precise chronostratigraphic dating. In contrast, many Late Jurassic occurrences from well-dated and time-constrained locations are included in the constrained dataset. As a result, the constrained dataset tends to have time bins with no data, particularly those from the Early Cretaceous, and is biased toward Late Jurassic occurrences, resulting in an underestimation of paleodiversity. Using an unconstrained dataset, however, may lead to an overestimation of true diversity values per time bin.

Many fossil taxa, and especially those within Rebbachisauridae, are rare and fragmentary, which introduces biases related to taxonomic resolution. Furthermore, sauropod phylogeny and taxonomy within crown diplodocoid subclades have undergone significant revisions in recent years (e.g., Harris 2006; Whitlock 2011a; Carballido et al. 2012; Mannion et al. 2012a, 2019b, 2019a; Tschopp and Mateus 2013; Mannion and Barrett 2013; Gallina et al. 2014, 2019; Fanti et al. 2015; Tschopp et al. 2015; Tschopp and Mateus 2017; Canudo et al. 2018; Xu et al. 2018; Coria et al. 2019; Whitlock and Mantilla 2020; Mannion and Moore 2025). Current literature generally accepts the following points: (1) the early evolutionary history of the group remains uncertain; (2) Rebbachisauridae is believed to have originated in South America (Whitlock 2011a; Fanti et al. 2015); (3) a ghost lineage of Rebbachisauridae is thought to span from the Middle Jurassic to the late Berriasian (Carballido et al. 2012; Fanti et al. 2015; Carpenter 2018); and (4) the “mid”-Cretaceous rebbachisaurid fauna of Africa and South America consists of closely related forms, which likely share a more recent common history before the continents were separated during the Early Cretaceous, leading to vicariant processes (e.g., Ali and Krause 2011; Whitlock 2011a).

Our global geographic distribution patterns are also affected by bias. The apparent absence of Diplodocoidea during the Middle and Late Jurassic in equatorial regions likely reflects a lack of available continental outcrops and under sampling in areas like present-day Brazil and North Africa, where the fossil record is poorly known, rather than a true absence of these taxa. Importantly, our GLS models already incorporate sampling and geographic proxies, including tetrapod-bearing collections, tetrapod-bearing formations, and available landmass area. In addition, preservational bias in the fossil record generally favors large-bodied taxa such as sauropods, which have a higher likelihood of fossilization and subsequent recovery due to their body size. Furthermore, lithogeographic proxies such as rock outcrop area and volumetric estimates of fossiliferous rock availability have been shown to correlate poorly or inconsistently with the sampling biases they are intended to represent (e.g., Dunhill 2012; Dunhill et al. 2012). At present, we consider that the most reliable way to evaluate paleogeographic sampling bias for a given clade is through increased sampling effort in undersampled regions and systematic reexamination of historical material using modern methodologies and technologies.

Finally, it is crucial to mention the underreported impact socioeconomic biases and the legacy of scientific colonialism have on the fossil record (Cisneros et al. 2022; Raja et al. 2022; Macário et al. 2025). As Raja et al. (2022) state, the heterogeneous spatial sampling of fossil material across the globe is a testament to a global power imbalance that still persists in paleontology, despite pro- inclusion, diversity, and reparation efforts conducted in the last few years. Considering all these factors, it is crucial to acknowledge the inherent biases present in the fossil record when interpreting our results.