Palaeontology meets metacommunity ecology: The Maastricthian dinosaur fossil record of North America as a case study

Documenting the patterns and potential associated processes of ancient biotas has always been a central challenge in palaeontology. Over the last decades, intense debate has focused on the organisation of dinosaur–dominated communities, yet no general consensus has been reached on how these communities were organised in a spatial context and if primarily affected by abiotic or biotic agents. Here, we used analytical routines typically applied in metacommunity ecology to provide novel insights into dinosaurian distributions across the latest Cretaceous of North America. To do this, we combined fossil occurrences with functional, phylogenetic and palaeoenvironmental modelling, and adopted the perspective that more reasonable conclusions on palaeoecological reconstructions can be gained from studies that consider the organisation of biotas along ecological gradients at multiple spatial scales. Our results showed that dinosaurs were restricted in range to different parts of the Hell Creek Formation, prompting the recognition of discrete and compartmentalised faunal areas during the Maastrichtian at fine-grained scales, whereas taxa ranges formed quasi–nested groups when combining data from various geological formations across the Western Interior of North America. Although groups of dinosaurs had coincident range boundaries, their communities responded to multiple ecologically–important gradients when compensating for differences in sampling effort. Metacommunity structures of both ornithischians and theropods were correlated with climatic barriers and potential trophic relationships between herbivores and carnivores, thereby suggesting that dinosaurian faunas were shaped by physiological constraints and a combination of bottom-up and top-down forces across multiple spatial grains and extents.

Methods
Dinosaur occurrences for the Maastrichtian of North America were retrieved from the Palaeobiology Database <PaleoDB.org> on May 2020, using the taxon name ‘Dinosauria’ and a time span of 72.1 – 66.0 Ma. Critically, although studies on modern community associations are limited to relatively brief periods of sampling time, fossil assemblages are windows on the faunas of ancient worlds occurring within particular chronostratigraphic units (Benson et al. 2018). Although this coarse temporal resolution will undoubtedly confound the data (which is addressed in detail below), it would be problematic to subdivide the time bins further, not least because only a handful of fossil assemblages are sufficiently informative to provide confident community-level estimates so far (Vavrek & Larsson 2010). Additionally, due to an insufficient amount of comparative data within high–resolution time bins (Dean et al. 2020) and the inherent errors in radiometric dating (Gates et al. 2010), the creation of a more tightly constrained correlative window is presently impractical. Here, we only retained occurrences belonging to Ornithischia and Theropoda since these two clades were the most diverse and abundant non–avian dinosaur groups in the latest Cretaceous of North America (Brusatte et al. 2015). Generic–level identifications were used in our study, and all avian taxa were excluded when delineating community types to keep our data more comparable to previous works (e.g. Vavrek & Larsson 2010; Dean et al. 2020). While birds are phylogenetically part of the dinosaurian clade, the different habits and habitats of latest Cretaceous Avialae (either diving or volant taxa) separates these faunas enough from ground-dwelling dinosaurs to justify their functional distinction in the context of the communities modelled here (see Heino et al. 2015b for an example on present-day biotas). Although the value of generic taxonomic ranks in community analyses has been debated, palaeontologists have used generic–level clades to investigate distributional patterns and variation in community composition of fossil taxa (e.g. Vavrek & Larsson 2010; Chiarenza et al. 2019; Dean et al. 2020). Indeed, generic–level identifications are preferred over species taxonomic ranks in dinosaur palaeobiology studies as most dinosaur genera (c. 87%) are easily diagnosed and monospecific (Weishampel et al. 2004; Mannion et al. 2012). Moreover, genus-level and species–level diversity patterns generally appear to track each other for Mesozoic tetrapods (Barrett et al. 2009), and genera are more taxonomically stable than species for many groups (Robeck et al. 2000). Here, however, taxa with unclear genus identification were discarded (i.e. we did not incorporate ‘cryptic’ diversity represented by taxonomically undiagnostic fossil remains that potentially represent distinct taxa, nor we did infer ghost lineages based on phylogenetic diversity estimates; Barrett et al. 2009; Mannion et al. 2011), and so were collections lacking formational assignment. If questionable ages appeared (e.g. ages notably deviating from ages of other collections from the same formation), they were either revised or excluded. These data are an up–to–date record of North American dinosaur faunas and therefore incorporate new Late Cretaceous fossils discovered over the past few years. Overall, our pruned dataset comprised 43 dinosaur genera, and consisted of 11 formations across the WIB of North America and 17 well–sampled locations across the Hell Creek landscape.

Palaeoclimatic general circulation model. In this study, we used palaeoclimatic model outputs (here, near-surface [1.5 m] mean annual temperature (TempMean), near-surface [1.5 m] annual temperature standard deviation (TempSDann), annual average precipitation (PrecMean) and annual precipitation standard deviation (PrecSDann)) from the fully coupled atmosphere-ocean GCM HadCM3L v. 4.5 Atmospheric–Ocean General Circulation Model (Valdes et al. 2017). More specifically, we followed the nomenclature of Valdes et al. (2017) and applied the HadCM3BL–M2.1aE version of the model. The conditions of the model simulations for the Maastrichtian consist of an atmospheric CO2 concentration of 1120 ppmv, which is within the range of uncertainty provided by the recent proxy pCO2 reconstructions of Foster et al. (2017). The model simulations were run for a total of 1422 years, and the climate variables used in our analyses were an annual average of the last 30 years of these simulations. HadCM3L has contributed to the Coupled Mode Intercomparison Project experiments demonstrating skill when it comes to reproducing present-day climates (Collins et al. 2001; Valdes et al. 2017) and has also been used for an array of different palaeoclimate evaluations during the Eocene (Lunt et al. 2012), the Oligocene (Li et al. 2018) and the Miocene (Bradshaw et al. 2012). Detailed information on this palaeoclimatic model, including large–scale circulation (and associated energy and momentum fluxes) and temporal fluctuations, as well as the impacts of fine-scale orographic features on climate signals, are available elsewhere (e.g. Lunt et al. 2016; Chiarenza et al. 2019).

Palaeogeographical digital elevation models (DEMs). The Maastrichtian palaeogeography for this study is that of Scotese & Wright (2018), which has been compiled as a palaeo-digital elevation model to facilitate grid-based analyses. In brief, these maps were created from publicly available stratigraphic literature, supplemented by fieldwork, including lithology, palaeoenvironmental information and broad-scale facies identification. For large–scale analyses, these palaeogeographies were upscaled to the palaeoclimatic model resolution (3.75° x 2.5°). This means that topographic and bathymetric information was broadly conserved, as it was resolved at a lower resolution (see Chiarenza et al. 2019 for a similar approach).

Functional features. Each dinosaur taxon was classified into several functional guilds based on body mass (very small, small, medium, large and very large), locomotor mode (bipeds, facultative bipeds –capable of both quadrupedal and bipedal motion– and quadrupeds) and trophic habits (carnivores, omnivores and herbivores, and for the latter, low and high browsers).

Body mass is perhaps the single most important and meaningful functional trait for animals, as it ultimately affects many aspects of their biology including metabolic rates, mechanical constraints, ecological performance and lifestyle strategies related to feeding, locomotion and reproduction (Loeuille & Loreau 2006; Iossa et al. 2008). Here, we used body mass estimates (very small ≤ 10 kg; 10 kg < small ≤ 100 kg; 100 kg < medium ≤ 1000 kg; 1000 kg < large ≤ 10000 kg; very large > 10000 kg; Noto & Grossman 2010) based on adult representatives from the comprehensive dataset of Benson et al. (2014), which provides a wide list of dinosaur taxa using the scaling relationship of limb bone robustness (stylopodial circumference; Campione & Evans 2012). To obtain a more comprehensive understanding of body mass distributions in our data, we further applied an inflection point criterion based on the Barry & Hartigan (1993) product partition model with Markov chain Monte Carlo (MCMC). More specifically, this algorithm used the posterior probability of changes over 10000 MCMC iterations, excluding the first 1000 as burn in, to distinguish among different body mass categories in the latest Cretaceous dinosaurs of North America. Interestingly, this Bayesian analysis roughly identified most of the original body mass categories used in our study, with each category broadly representing an order of magnitude (García–Girón et al. 2020b, appendix S1, fig. S1).

Trophic habits refer to the food processing strategies and diet of an animal, and it generally includes three primary categories, i.e. carnivores, herbivores and omnivores. Further subdivisions depend on the biological knowledge of the morphology (e.g. teeth morphology and skull) and behaviour of the study organismal group. Here, we assigned herbivores to categories of browse height rather than plant type due to the virtually unknown nature of plant preferences in dinosaurs. More specifically, we roughly assigned a simple maximum browsing limit (low ≤ 2 m; high > 2 m) based on characters such as limb length and neck posture using Noto & Grossman (2010) and Mallon et al. (2013).

We further divided locomotor mode into two major categories: quadrupeds and bipeds. For those taxa with intermediate axial and limb morphologies in proportions between those of bipeds and obligate quadrupeds (e.g. Hadrosauridae), we included an additional locomotor division, i.e. facultative bipeds (see Noto & Grossman, 2010 for a similar approach). For the following analyses, we applied the mixed–variables coefficient of distance (i.e. a generalisation of Gower’s distance; Pavoine et al. 2009) to extract a functional distance matrix, which described the functional differences between all taxon pairs based on body mass, trophic habits and locomotor mode (e.g. Heino & Tolonen 2017). Thereafter, the pairwise output values for the functional distance matrix were synthesised into separate axes using principal coordinate analysis (PCO) and following Duarte et al. (2012).

See the main text for References.

揭示古生物群落的分布模式及其潜在关联过程，始终是古生物学领域的核心挑战之一。近数十年来，学界围绕恐龙主导的古群落结构展开了激烈讨论，但关于这些群落在空间维度上的组织方式，以及其主要受非生物还是生物因子调控，至今尚未达成广泛共识。本研究采用集合群落生态学（metacommunity ecology）中常用的分析范式，对北美晚白垩世最晚期的恐龙类群分布模式展开全新解析。为此，我们将化石产出记录与功能性状、系统发育及古环境建模相结合，并秉持如下研究视角：只有同时考虑生物群落在多空间尺度下沿生态梯度的组织模式，才能为古生态重建提供更为合理的结论。

研究结果显示，恐龙类群的分布范围局限于地狱溪组（Hell Creek Formation）的不同区域，这表明在马斯特里赫特阶（Maastrichtian）的精细空间尺度下，存在多个独立且分隔的动物群分区；而当整合北美西部内陆区不同地质组的数据时，类群的分布范围则形成了准嵌套群落结构。尽管部分恐龙类群的分布边界存在重合，但在校正采样强度差异后，其所在群落仍对多个具有生态重要性的梯度响应显著。鸟臀目（Ornithischia）与兽脚亚目（Theropoda）的集合群落结构均与气候屏障及植食性与肉食性类群间的潜在营养关系显著相关，这表明恐龙动物群的形成受到生理限制，以及多空间粒度与空间范围下上行调控与下行调控共同作用的塑造。

## 方法
本研究于2020年5月从古生物学数据库（Palaeobiology Database，<PaleoDB.org>）获取北美马斯特里赫特阶的恐龙类群产出记录，检索分类单元为‘恐龙总目（Dinosauria）’，时间跨度为72.1~66.0 Ma。值得注意的是，尽管现代群落关联研究的采样时间跨度相对有限，但化石组合是窥探特定年代地层单元内古生物群面貌的窗口（Benson等，2018）。尽管这种较粗的时间分辨率无疑会对数据造成干扰（下文将详细讨论该问题），但进一步细分时间间隔仍存在诸多问题，其中最关键的一点是，截至目前仅有少量化石组合拥有足够信息，能够可靠地开展群落水平的估算（Vavrek & Larsson，2010）。此外，由于高分辨率时间间隔内的比较数据不足（Dean等，2020），且放射性定年本身存在固有误差（Gates等，2010），目前难以构建约束更为严格的关联时间窗口。

本研究仅保留隶属于鸟臀目（Ornithischia）与兽脚亚目（Theropoda）的产出记录，因为这两个演化支是北美晚白垩世最晚期最为多样且丰富的非鸟类恐龙类群（Brusatte等，2015）。本研究采用属级分类鉴定，在划定群落类型时剔除所有鸟类类群，以确保研究数据与既往研究（如Vavrek & Larsson，2010；Dean等，2020）具有更好的可比性。尽管从系统发育角度而言，鸟类属于恐龙总目的演化支，但晚白垩世最晚期的今鸟亚纲（Avialae）类群（包括潜水或飞行类群）的习性与栖息环境均与陆生恐龙存在显著差异，因此在本研究构建的群落模型中，有理由将其与陆生恐龙在功能上区分开来（可参考Heino等，2015b对现代生物群的相关研究案例）。

尽管学界对属级分类单元在群落分析中的应用价值仍存在争议，但古生物学家已广泛采用属级演化支来研究化石类群的分布模式与群落组成变异（如Vavrek & Larsson，2010；Chiarenza等，2019；Dean等，2020）。事实上，在恐龙古生物学研究中，属级鉴定相较于物种级鉴定更受青睐，因为约87%的恐龙属均可通过形态特征明确界定，且多数为单型属（Weishampel等，2004；Mannion等，2012）。此外，中生代四足动物的属级与物种级多样性模式通常具有一致性（Barrett等，2009），且多数类群的属级分类单元比物种级更为稳定（Robeck等，2000）。但本研究剔除了属级鉴定不明确的类群（即未纳入由无法通过分类学界定的化石遗骸所代表的‘隐蔽’多样性——这类遗骸可能代表独立类群，也未基于系统发育多样性估算推断幽灵支系；Barrett等，2009；Mannion等，2011），同时剔除了缺乏组级归属的化石产出记录。若出现存疑的年代数据（如与同一地层组内其他产出记录的年代偏差显著），则要么对其进行修正，要么直接剔除。本数据集为北美恐龙动物群的最新记录，纳入了过去数年新发现的晚白垩世化石。最终，经过筛选的数据集共包含43个恐龙属，涵盖北美西部内陆盆地（Western Interior Basin, WIB）的11个地层组，以及地狱溪地貌区的17个采样充分的点位。

### 古气候通用环流模型
本研究采用耦合大气-海洋通用环流模型（General Circulation Model, GCM）HadCM3L v4.5的古气候模拟输出结果，具体变量包括：近地表（1.5m处）年平均气温（TempMean）、近地表（1.5m处）年气温标准差（TempSDann）、年平均降水量（PrecMean）以及年降水量标准差（PrecSDann）（Valdes等，2017）。具体而言，本研究遵循Valdes等（2017）的命名规范，采用该模型的HadCM3BL–M2.1aE版本。马斯特里赫特阶的模型模拟设置为大气CO₂浓度1120 ppmv，该数值处于Foster等（2017）基于代用指标重建的pCO₂不确定性区间范围内。模型模拟总时长为1422年，本研究分析所用的气候变量为模拟最后30年的年平均值。HadCM3L曾参与耦合模式比较计划（Coupled Model Intercomparison Project）试验，其在重现现代气候方面表现优异（Collins等，2001；Valdes等，2017），同时也被用于始新世（Lunt等，2012）、渐新世（Li等，2018）以及中新世（Bradshaw等，2012）的多类古气候评估研究。该古气候模型的详细信息，包括大尺度环流（及相关能量与动量通量）、时间波动，以及精细尺度地形特征对气候信号的影响，均可参考其他文献（如Lunt等，2016；Chiarenza等，2019）。

### 古地理数字高程模型（DEMs）
本研究采用的马斯特里赫特阶古地理数据来自Scotese & Wright（2018），该数据已被整理为古数字高程模型（Palaeo-Digital Elevation Models, DEMs），以支持基于格网的分析。简言之，该古地理图基于公开出版的地层学文献编制，并结合了野外调查数据，包括岩性、古环境信息以及大尺度相带识别结果。针对大尺度分析，我们将古地理数据重采样至古气候模型的分辨率（3.75°×2.5°）。这一操作可在较低分辨率下最大限度保留地形与水深信息（类似方法可参考Chiarenza等，2019）。

### 功能性状
本研究基于体重、运动方式与食性，将每个恐龙类群划分为多个功能群：体重等级分为极小型、小型、中型、大型与极大型；运动方式分为两足行走、兼性两足行走（即可同时采用四足与两足运动方式）与四足行走；食性分为肉食性、杂食性与植食性，其中植食性类群进一步划分为低食性与高食性。

体重或许是动物最为关键且具有生物学意义的功能性状，因为它最终会影响动物的诸多生命特征，包括代谢速率、机械约束、生态表现以及与摄食、运动和繁殖相关的生活史策略（Loeuille & Loreau，2006；Iossa等，2008）。本研究基于Benson等（2014）的综合数据集，采用肢体骨骼粗壮度（肢骨周径；Campione & Evans，2012）的缩放关系估算成年个体体重，并按照如下等级划分：极小型≤10kg；10kg<小型≤100kg；100kg<中型≤1000kg；1000kg<大型≤10000kg；极大型>10000kg（Noto & Grossman，2010）。为更全面地解析研究数据中的体重分布模式，本研究进一步采用基于Barry & Hartigan（1993）乘积分割模型的拐点判定准则，并结合马尔可夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）方法。具体而言，该算法通过10000次MCMC迭代得到的后验概率变化（舍弃前1000次迭代作为预烧期），对北美晚白垩世最晚期恐龙的不同体重等级进行区分。有趣的是，该贝叶斯分析大致识别出了本研究中使用的多数体重等级，每个等级大致对应一个数量级（García-Girón等，2020b，附录S1，图S1）。

食性指动物的食物处理策略与饮食类型，通常分为肉食性、植食性与杂食性三大类。进一步的细分则依赖于研究类群的形态（如牙齿形态与头骨结构）与行为学相关生物学知识。由于恐龙的植物偏好几乎无从得知，本研究将植食性类群按取食高度而非植物类型进行分类。具体而言，我们基于肢体长度与颈部姿态等特征，参考Noto & Grossman（2010）与Mallon等（2013）的方法，粗略设定了两个取食高度等级：低食性≤2m，高食性>2m。

我们进一步将运动方式划分为两大类别：四足行走与两足行走。对于那些兼具两足与专性四足动物轴向及肢体形态比例的类群（如鸭嘴龙科（Hadrosauridae）），我们增设了兼性两足行走这一运动类型（类似方法可参考Noto & Grossman，2010）。在后续分析中，我们采用混合变量距离系数（即Gower距离的推广形式；Pavoine等，2009）构建功能距离矩阵，该矩阵基于体重、食性与运动方式描述所有类群对之间的功能差异（如Heino & Tolonen，2017）。随后，我们按照Duarte等（2012）的方法，采用主坐标分析（Principal Coordinate Analysis, PCO）将功能距离矩阵的两两配对值合成为独立轴。

参考文献详见正文。