Estimating the evolutionary rates in mosasauroids and plesiosaurs: discussion of niche occupation in Late Cretaceous seas

Bayesian inference

Both mosasauroid and plesiosaur phylogenetic datasets were analyzed using the protocol discussed in Madzia & Cau (2017), integrating the morphological data matrices with absolute ages of the least inclusive stratigraphic range including each terminal unit. The Sampled Ancestor Fossilized Birth Death Skyline Model (SAFBD) of Gavryushkina et al. (2014) and Gavryushkina et al. (2017) implemented in BEAST 2.4.4. (Drummond et al., 2012; Bouckaert et al., 2014) was used as tree model. Since the character matrices did not include autapomorphies of the sampled taxa, the Lewis’s (2001) model was conditioned to variable characters only using the implementation included in BEAST 2.4.4. Stratigraphic information for the mosasauroid and plesiosaur taxa was taken from the literature (Polcyn et al., 2014; Fischer et al., 2017; respectively), and converted to geochronological ages. Stratigraphic data and age constraints for each terminal were obtained from the Paleobiology Database (http://paleobiodb.org/), checked against the International Chronostratigraphic Chart (v2019/05; http://stratigraphy.org/), and included as uniform prior for tip-dating (Supplemental Information I).

The impact of using (or omitting) age priors incorporating stratigraphic uncertainty in tip-dating has only recently been addressed (Barido-Sottani et al., 2019; Cau, 2019). Note that in their Bayesian analysis of Mosasauroidea, Madzia & Cau (2017) used a punctiform age prior for each terminal taxon (i.e., the mean value of the shortest age range encompassing the stratigraphic uncertainty), thus they did not incorporate age uncertainty in tree reconstruction. Such a strategy may arbitrarily set the age of several taxa sharing the same stratigraphic uncertainty to an identical value, thus enforcing a strictly cladogenetic pattern for their relationships even under the SAFBD model, and biasing tree reconstruction favoring longer ghost lineages. Furthermore, punctiform tip-dating priors may lead to inflated divergence rates for taxonomic units scored from multiple non-contemporary specimens (see Cau, 2019).

The following protocol was used for both mosasauroid and plesiosaur datasets. Each BEAST analysis involved 3 replicate runs (with different random starting trees and random number seeds). Each of the 3 replicate runs involved 10 million steps with sampling every 1,000 generations, with a burnin of 4 million steps. This protocol is similar to the one followed by Simões et al. (2017) but used an additional independent run for each analysis (i.e., three instead of two) and set a more conservative burnin (40% instead of 25%). We used Tracer 1.5 (Rambaut & Drummond, 2009) to determine whether the runs reached stationary phase, and to assess convergence of the independent runs. The post-burnin parameter and tree samples were retained for the analysis and concatenated using LogCombiner in the BEAST package. Estimates (mean and 95% highest posterior density) for all numerical parameters were generated using Tracer 1.5 (Rambaut & Drummond, 2009). We used the MCCT to reconstruct the cladogenetic events (median age) and to infer the divergence rate (the amount of morphological change per branch relative to the whole topology) for both clades. Note that the absolute rate values are inversely related to the sample size (i.e., rate value in a branch is proportional to the probability of sampling each state transition of the clade history in that branch); thus, direct comparisons between the mosasauroid and plesiosaur rate values is meaningless. Given the rate distribution inferred along the MCCT in the two clades, we here define ‘high rates’ all those values equal or higher than the value at the 75 percentile in each rate distribution.

The mosasauroid matrix was acquired from Simões et al. (2017), which is a recent version of the dataset first introduced in Bell’s (1993) PhD thesis and published by Bell (1997). The recent version of Simões et al. (2017) was further modified with respect to some taxon names (as in Madzia & Cau, 2017). In sum, the dataset consists of 44 operational taxonomic units (OTUs) analyzed using 125 characters (see Supplemental Information II for BEAST file, Supplemental Information III for NEXUS file, and Supplemental Information IV for character list).

The phylogenetic assessment of Plesiosauria was performed using a modified version of the dataset first assembled by Benson & Druckenmiller (2014). We first took a recent version of that dataset, published by Madzia, Sachs & Lindgren (2019), and updated it based on personal observations and recently published literature, to include representatives of distinctive plesiosaur clades. The changes include: modifications to the scores of Brancasaurus brancai and ‘Gronausaurus wegneri’ as in Sachs, Hornung & Kear (2016); addition of Lagenanectes richterae from Sachs, Hornung & Kear (2017); addition of Nakonanectes bradti, Albertonectes vanderveldei, Aristonectes quiriquinensis, Elasmosaurus platyurus, ‘Hydralmosaurus serpentinus’, Mauisaurus haasti, ‘Libonectes’ atlasense, Terminonatator ponteixensis, Tuarangisaurus keyesi, Zarafasaura oceanis, Kawanectes lafquenianum, and Vegasaurus molyi from Serratos, Druckenmiller & Benson (2017); addition of Neusticosaurus pusillus and Nothosaurus marchicus, and modifications to the scores of Yunguisaurus liae and Pistosaurus OTUs as in Wintrich et al. (2017); addition of Acostasaurus pavachoquensis, ‘Kronosaurus’ boyacensis, and Sachicasaurus vitae from Páramo-Fonseca, Benavides-Cabra & Gutiérrez (2018), with amended scores for A. pavachoquensis and S. vitae as in Páramo-Fonseca, Benavides-Cabra & Gutiérrez (2019); modifications to the character scores of Thililua longicollis and addition of Eopolycotylus rankini, Manemergus anguirostris, Dolichorhynchops tropicensis, Georgiasaurus penzensis, Dolichorhynchops sp. (specimen ROM 29010), Dolichorhynchops herschelensis, Sulcusuchus erraini, and Mauriciosaurus fernandezi following Fischer et al. (2018), with amended scores for Trinacromerum bentonianum, Dolichorhynchops osborni, Dolichorhynchops bonneri, Mauriciosaurus fernandezi, and Polycotylus latipinnis as in Morgan & O’Keefe (2019); addition of Styxosaurus snowii from Sachs, Lindgren & Kear (2018); and modifications to the scores of Kronosaurus queenslandicus and Stenorhynchosaurus munozi following Holland (2018) and Páramo-Fonseca, Benavides-Cabra & Gutiérrez (2019).

It is essential to note that although the elasmosaurid phylogenetic relationships were a subject of several recent papers (e.g.,  Otero, 2016; Sachs, Hornung & Kear, 2016; O’Gorman et al., 2017; Sachs & Kear, 2017; Serratos, Druckenmiller & Benson, 2017; Sachs, Lindgren & Kear, 2018; O’Gorman et al., 2019), interpretations of morphologies observed in some elasmosaurid specimens differ between these studies. See, for example, conflicting scores for ‘Libonectes’ atlasense in Sachs & Kear (2017) and Serratos, Druckenmiller & Benson (2017), and for Styxosaurus snowii in Serratos, Druckenmiller & Benson (2017) and Sachs, Lindgren & Kear (2018). We have not studied these taxa in person; as such, we had to choose between scores provided in other publications. We decided to adopt those scores that derive from more recent studies in which the taxa were assessed based on direct observations. For that reason, ‘Libonectes’ atlasense is here scored as in Serratos, Druckenmiller & Benson (2017) and Styxosaurus snowii as in Sachs, Lindgren & Kear (2018). Naturally, such differences in interpretations of character states might have an impact on inferred evolutionary rates. Nevertheless, our decisions should not have any impact on the findings of the present study as Elasmosauridae is of marginal importance here.

Finally, we have also modified several character states of Anguanax zignoi based on personal observations of the type specimen (3: 0→?; 4: [12]→?; 121: 0→?; 137: 0→?; 150: 1→?; 207: [01]→?; 270: 1→0) and re-scored Megacephalosaurus eulerti for character 27 (1→0). This score has been already advocated by Madzia, Sachs & Lindgren (2019: p. 1208) but the character was erroneously scored as ‘1’ rather than ‘0’ in that study.

Additionally, we have modified several character definitions in the character list of Benson & Druckenmiller (2014):

Character 25. The character description was changed from “Maxilla, posterior extent of maxillary tooth row” to “Maxilla and dentary, posterior extent of maxillary tooth row”; after Serratos, Druckenmiller & Benson (2017).

Character 138. As noted by Madzia, Sachs & Lindgren (2019), the current state definitions for character 138 are problematic because they do not cover all plesiosaurs. In the original character list of Benson & Druckenmiller (2014), state ‘0’ was defined as codable for taxa with 12–17 maxillary teeth, state ‘1’ for taxa with 20–25 maxillary teeth, and state ‘2’ for taxa with more than 28 maxillary teeth. However, the brachauchenine pliosaurid Megacephalosaurus eulerti was shown to possess 18 teeth in the right and 19 in the left maxilla, thus falling between states ‘0’ and ‘1’. Two options were considered for M. eulerti: to score it as ‘0’, extending the state to cover taxa with 12–19 maxillary teeth, and as ‘1’, extending the state to cover taxa with 18–25 teeth in their maxillae. Madzia, Sachs & Lindgren (2019) used both these options and explored the effects of such settings. Considering that the last brachauchenines have reduced numbers of teeth in their jaws (Madzia, Sachs & Lindgren, 2019), scoring these taxa in the same way as their older relatives (that fall near the upper boundary of state ‘1’) might hinder the inference of some potential phylogenetic signal. Therefore, in this study, state ‘0’ covers taxa with 12–19 maxillary teeth, state ‘1’ covers taxa with 20–27 maxillary teeth (note that the upper boundary was extended to eliminate the gap between states ‘1’ and ‘2’), and state ‘2’ covers taxa with at least 28 maxillary teeth.

Character 139. State ‘2’ (“intermediate between states 0 and 1, with a flattened labial surface, but this surface [is] not substantially expanded anteroposteriorly [= subtrihedral]”) was added after Benson et al. (2013). Even though Benson et al. (2013) described the state ‘2’ as “intermediate”, the morphological transition from “round or sub-rounded” (‘0’) to “sub-triangular [= trihedral]” (‘1’) does not need to have the appearance of ‘2’. Later, Serratos, Druckenmiller & Benson (2017) used the dataset of Benson & Druckenmiller (2014) to infer the interrelationships of elasmosaurids and modified character 139 to include another new state (‘2’): “suboval”. However, in their data matrix, this state was scored as ‘3’. In this study, the state ‘2’ is equivalent to state ‘2’ of Benson et al. (2013), and state ‘3’ follows the new “suboval” state introduced by Serratos, Druckenmiller & Benson (2017). Nevertheless, the perception of what is “sub-rounded” (‘0’) and what “suboval” (3) may be partially dependent on subjective criteria. As such, future larger-scale phylogenetic studies of Plesiosauria should probably quantify the difference (for example, using the ‘width-to-length ratio’ [WLR] of Madzia (2016) or similarly defined parameter). Due to the lack of apparent transitional nature of particular character states, this character should stay unordered in parsimony analyses.

Character 153. State ‘3’ (“even longer, corresponding to the ‘can’ shaped morphology of Otero et al. (2016a) or ‘elongate’ morphology of O’Keefe & Hiller (2006)”) was added after Serratos, Druckenmiller & Benson (2017). We did not modify the state definition though, again, its quantification would enable to maintain an unambiguous application.

Character 248. The character description was changed from “Propodials, angle between long axes of epipodial facets in dorsal view” to “Humerus, angle between long axes of epipodial facets in dorsal view”; after Serratos, Druckenmiller & Benson (2017).

The full dataset consists of 125 OTUs analyzed using 270 characters (see Supplemental Information V for BEAST file, Supplemental information VI for NEXUS file, and Supplemental information VII for character list). This is the largest dataset used in phylogenetic assessment of Plesiosauria published to date.

Reconstruction of geographic origin

A detailed assessment of the historical biogeography of mosasauroids and plesiosaurs is beyond the scope of our study; however, for the sake of discussion related to possible competitive interactions of mosasauroids and plesiosaurs, we provide a brief insight into possible geographic origins for mosasauroid and plesiosaur clades. We took the MCCTs resulting from the Bayesian analyses of both datasets, loaded them to Mesquite 3.61 (Maddison & Maddison, 2019) and created a new character matrix with a single new character scored for each OTU—the continent of discovery. Such demarcation only very roughly corresponds to the real ancient biogeographic settings (such as Cretaceous epicontinental seas); however, for the purposes of our study such highly simplified approach seems to be sufficient as the ancestral areas in the time interval that is of special interest for us (e.g., the Western Interior Seaway [covered under ‘North America’] in the Turonian; see Figs. 1 and 2) do not seem to comprise water bodies with geographically isolated populations.

The geographic origins for particular clades were reconstructed using the ‘Trace Character History’ option and with ‘Parsimony Ancestral States’ as the ancestral state reconstruction method. The resulting ancestral areas were mapped on Figs. 1 and 2.

Rates of morphological evolution

The rate of evolutionary divergence inferred among the branches of the mosasauroid MCCT (both internodes and terminal leaves; Fig. 1) ranges between 0.0093 and 0.0800 changes per branch (median value = 0.0289, 25–75 percentile range: 0.0195–0.0367). During their evolutionary history mosasauroids experienced intervals of elevated rates of divergence. High rates are especially apparent at the origin of Russellosaurina (∼100–98 Mya; rate = 0.0588), the clade consisting of plioplatecarpines and tylosaurines (∼98–93 Mya; rate = 0.0472), the base of Tylosaurinae (∼93–87 Mya, rate = 0.0674), and near the base of Mosasaurinae (∼97–94 Mya; rate = 0.0577). An increase in rates of morphological evolution is also apparent at the base of the node comprising Mosasaurus and Plotosaurus (∼77–75 Mya; rate = 0.0588).

Only three plesiosaur lineages might have affected mosasauroid evolution or been affected by interactions with them. These include the polycotylid leptocleidians, elasmosaurids, and brachauchenine pliosaurids. However, with their elongated necks and long and pointed teeth, elasmosaurids were substantially distinct from all known members of Mosasauroidea, thus occupying dissimilar niches (e.g., Massare, 1987). As such, additional comparisons will be limited only to the brachauchenine and polycotylid plesiosaurs (see Fig. 2; due to the large size of the MCCT of Plesiosauria, the full tree is provided as the Supplemental Information VIII).

The origin of the least inclusive clade formed by Russellosaurina and Mosasaurinae (+ Halisaurinae) occurred ∼100 Mya, and the splitting of the clade uniting plioplatecarpines and tylosaurines ∼93 Mya; thus covering the last few millions of years of the brachauchenine evolution. In turn, during the splitting of the smallest clade comprising globidensins and mosasaurins (∼81 Mya) and the origin of the Mosasaurus node (∼75 Mya), pliosaurids were probably already extinct and polycotylids rare (e.g., Sato, 2005; Salgado, Parras & Gasparini, 2007; O’Gorman & Gasparini, 2013; Novas et al., 2015; Fischer et al., 2018; Morgan & O’Keefe, 2019).

The last increase in the rates of morphological evolution within Pliosauridae occurred at the base of the clade formed by Luskhan itilensis, Stenorhynchosaurus munozi, and Sachicasaurus vitae (∼142–139 Mya) and within that clade. Last brachauchenines of the early Late Cretaceous (Cenomanian-Turonian) have not experienced significant oscillations in their evolutionary rates, indicating potential evolutionary conservatism among the last pliosaurids.

In contrast, polycotylids experienced elevated rates in the ‘middle’ Cretaceous members of the clade.

(See Supplemental Information IX for resulting ‘log’ files from the analyses of both datasets).

Bayesian analysis of plesiosaur phylogenetic relationships

The phylogenetic relationships of Mosasauroidea inferred by means of Bayesian analysis of the dataset presented herein have already been thoroughly discussed in Simões et al. (2017) and Madzia & Cau (2017). Owing to the fact that the overall tree topology of the mosasauroid MCCT is the same as in Madzia & Cau (2017), with minor differences in the mosasauroid outgroup and in the ‘Prognathodon’ and ‘Mosasaurus’ grouping, detailed discussion of the results is provided only for the plesiosaur tree (see Supplemental Information VIII for the full tree).

Inference of plesiosaur interrelationships through Bayesian analysis, including estimates of rates of divergence for some clades, has previously been published as well (Cau & Fanti, 2016); however, the matrix in that study was significantly reduced to include only 39 OTUs. Pliosaurids have been represented as in Benson et al. (2013), the study that Cau & Fanti (2016) took the dataset from, with the only addition being their newly-established pliosaurid taxon Anguanax zignoi. The representatives of other major clades (rhomaleosaurids and plesiosauroids) were mostly excluded. Anyway, the original dataset has been substantially expanded and modified since 2016 (see some of the recent additions and modifications in ‘Methods’). Furthermore, there are two substantial methodological differences between the Bayesian inference of Cau & Fanti (2016) and the one performed in the present study:

1. In the present study, the age uncertainty of each taxon is incorporated in the analysis: age prior of all fossil taxa is defined as a uniform range including the absolute age limits of the shortest stratigraphic range unambiguously including any taxon. This approach differs from that followed in Cau & Fanti (2016), who used for each tip a fixed age prior defined arbitrarily by the mean value of the stratigraphic uncertainty of each taxon.

2. The tree model used here discriminates between anagenetic and cladogenetic patterns of evolution; therefore, it may test whether some of the taxonomic units that are included actually form anagenetic sequences. The analysis in Cau & Fanti (2016) was run on a previous version of BEAST which did not implement the Sampled Ancestor Fossilized Birth Death Skyline Model (Gavryushkina et al., 2014), and thus was a priori constrained to reconstruct exclusively cladogenetic frameworks.

As discussed by Cau (2019), the results of the analysis of Cau & Fanti (2016) may thus be methodologically biased in potentially inflating the extent of the inferred ghost lineages and also in pre-dating the ages of the divergence events.

The present study includes the first Bayesian analysis of the full dataset, with most representatives of all major plesiosaur clades, and including the aforementioned modifications.

The overall tree topology broadly agrees with those reconstructed through more recent parsimony analyses (see, e.g., Ketchum & Benson, 2010; Benson et al., 2013; Benson & Druckenmiller, 2014; Fischer et al., 2015; Cau & Fanti, 2016; Otero, 2016; Sachs, Hornung & Kear, 2016; Fischer et al., 2017; O’Gorman et al., 2017; Sachs, Hornung & Kear, 2017; Serratos, Druckenmiller & Benson, 2017; Fischer et al., 2018; O’Gorman, Gasparini & Spalletti, 2018; Páramo-Fonseca, Benavides-Cabra & Gutiérrez, 2018; Sachs, Lindgren & Kear, 2018; Madzia, Sachs & Lindgren, 2019; Morgan & O’Keefe, 2019). Plesiosauria (posterior probability [pp ] = 1; node origin estimated at ∼241 Mya) basally branches into Rhomaleosauridae (pp = 0.96; ∼215 Mya) and Neoplesiosauria (pp = 0.88; ∼215 Mya), consisting of Pliosauridae (pp = 1; ∼206 Mya) and Plesiosauroidea (pp = 0.89; ∼210 Mya). Within the pliosaurid branch, two nodes have been named—Thalassophonea (pp = 0.59; ∼174 Mya) and Brachaucheninae (pp = 0.98; ∼152 Mya). Interestingly, contrary to recent studies assessing the phylogenetic relationships of pliosaurid plesiosaurs by means of parsimony analyses (e.g., Fischer et al., 2015; Fischer et al., 2017; O’Gorman, Gasparini & Spalletti, 2018; Páramo-Fonseca, Benavides-Cabra & Gutiérrez, 2018; Madzia, Sachs & Lindgren, 2019), the monophyly of Pliosaurus may be considered supported (pp = 0.83). All ‘major’ plesiosauroid sub-clades appear to be very well supported by our Bayesian analyses (with possible exception to the basal branching of Elasmosauridae; see below). The clade Cryptoclidia (pp = 1; ∼180 Mya) consists of Cryptoclididae (pp = 1; ∼176 Mya) and Xenopsaria (pp = 1; ∼158 Mya), which, then, comprises the clade Leptocleidia (pp = 1; ∼145 Mya)—including Leptocleididae (pp = 0.9; ∼140 Mya) and Polycotylidae (pp = 1; ∼119 Mya)—and its closest relatives (Brancasaurus and ‘Gronausaurus’), and the clade Elasmosauridae (pp = 0.44; ∼144 Mya). However, the low pp value for the base of Elasmosauridae might be due to the problematic placement of Lagenanectes richterae which might be a non-elasmosaurid xenopsarian (D. Madzia, unpublished results). The pp value of the more ‘traditional’ elasmosaurid grouping (that is, exclusive of L. richterae) is very high (pp = 0.99).

Estimates of evolutionary rates and potential biases

Comparing mosasauroid and plesiosaur traits

Tooth crown morphologies and trophic guilds

Driven by “a lack of precision in correlating tooth type and preferred prey” in Mesozoic marine amniotes, Massare (1987) designed seven trophic guilds based on tooth crown morphologies: (I) crush, (II) crunch, (III) smash, (IV) pierce I, (V) pierce II, (VI) general, and (VII) cut (Massare, 1987: 130–131), and provided characteristics of particular crown morphologies used in determining assignment to guilds (Massare, 1987: 132, Table 3). Her classification has been further modified by Hornung & Reich (2015) who divided the ‘cut’ guild into two categories, ‘cut I’ and ‘cut II’, to distinguish between taxa with labiolingually expanded crown with cutting edges (‘cut I’) from those with strongly compressed and blade-shaped teeth (‘cut II’) (see Fig. 3). Even though the guild classification proposed by Massare (1987) has not been accepted universally (Buchy, 2010), recent quantitative studies using teeth of marine amniotes to evaluate feeding ecologies validated such system (e.g., Foffa et al., 2018).

The Late Cretaceous brachauchenines and polycotylids can be categorized relatively easily within that system. With their conical and slightly curved teeth, both these clades belong to the non-carinate/unserrated ‘general’/‘smash’/‘pierce I’ guilds, though polycotylids usually possess less robust teeth leaning further towards the ‘pierce I’ part of the system (Fig. 3). The guild assignment of the mosasauroids is problematic due to the apparent and widespread pseudoheterodonty in some clades. Specifically, the shape of the mosasauroid teeth often differs depending on the positions of the teeth in jaws (see, e.g., LeBlanc, Caldwell & Bardet, 2012; Otero et al., 2016b; Madzia, 2019). For example, the derived mosasaurine mosasaurid Mosasaurus lemonnieri, which is known from a number of well-preserved specimens with reasonably complete jaws, shows anterior-to-middle tooth crowns with subcircular/ovoid cross-sections and carinae variably developed either on a short apical segment of the tooth crowns (in IRSNB R 366, 368, 369) or along their entire apicobasal length (serrated only in the somewhat problematic specimen IRSNB R 377). Its posteriorly positioned teeth, in turn, tend to be increasingly labiolingually compressed (Madzia, 2019). Further differences are appearing through the ontogeny as larger (presumably older) individuals appear to show more robust teeth. Such morphological diversity of its tooth crowns makes the taxon occupying a wide field of the ternary graph (from ‘pierce I’/‘pierce II’ through the ‘cut I’/‘cut II’ guilds). In tylosaurines, in turn, the teeth occupy the ‘smash’ to ‘general’/‘pierce II’ field, and mosasaurins show dental morphologies indicative of the ‘crunch’/‘general’ to ‘cut I’/‘cut II’ guilds (Schulp et al., 2006; Ross, 2009; Hornung & Reich, 2015). Therefore, many mosasauroids could likely occupy the same trophic levels of generalists and represent the same trophic guilds as robust-toothed short-necked plesiosaurs.

It is also essential to note here that even though most mosasauroids possess carinate teeth (which distinguishes them from the teeth of the last brachauchenines as well as polycotylids), their distal carinae are often displaced labially (see Fig. 4), especially in the teeth from the anterior to middle section of the jaws. In such cases, the distal carinae resemble (and are often less pronounced than) the apicobasal ridges observable in plesiosaur teeth (see Fig. 4 and the distribution of apicobasal ridges in brachauchenine teeth in Madzia (2016), discussion of outer enamel structural elements in pliosaurid teeth by Zverkov et al. (2018), and assessment of ridge evolution and function in marine amniotes by McCurry et al. (2019)) and certainly does not play the role of a cutting edge (as in labiolingually strongly compressed tooth crowns).

The record of interactions between mosasauroids and plesiosaurs

It is beyond doubt that sympatric mosasauroids and plesiosaurs interacted at the individual scale. Direct evidence, however, is scarce. Everhart (2004) published on a partial plesiosaur specimen preserved as the stomach contents of an 8.8-meter-long adult of Tylosaurus proriger. It was discovered in 1918 in the lower Campanian strata of the Niobrara Formation, near Twin Butte Creek, Logan County, Kansas, and first mentioned by Sternberg (1922) but it has not been properly described until 2004. Everhart (2004) admitted a poor state of preservation of the plesiosaur remains but suggested that they most likely belong to the polycotylid Dolichorhynchops osborni. Einarsson et al. (2010), in turn, described a plesiosaur propodial of latest early Campanian age, discovered at the Åsen locality, Kristianstad Basin, Sweden. Though incomplete, the specimen was identified as an indeterminate polycotylid. It possesses a distinctive bite mark that was interpreted by Einarsson et al. (2010) to be caused by a large mosasaurine comparable to Dollosaurus. Nevertheless, none of these finds could have been unequivocally inferred as representing either predatory or scavenging behavior. Further finds, of interest with respect to early mosasauroid diversification patterns, include mosasauroid and brachauchenine specimens that overlap in time and space (see specimens discussed in Martin & Stewart, 1977; Polcyn et al., 2008; Schumacher, 2011; Kear et al., 2014; Everhart, 2016). These discoveries, however, do not show any evidence of direct interactions between the members of the two clades.

Concluding remarks

We provide the first estimates of the evolutionary rates for mosasauroid and plesiosaur clades and use the results as a baseline for discussion and comparisons of traits that might have had some impact on the shape of the niche structures in Late Cretaceous seas. Owing to the known stratigraphic distribution of the mosasauroid and plesiosaur lineages, only three plesiosaur clades might have competed with mosasauroids; the elasmosaurids, brachauchenine pliosaurids, and polycotylid leptocleidians. However, considering the overall body plans of the taxa belonging to these groups, and their tooth crown morphologies, which are key indicators of dietary niche partitioning (e.g., Massare, 1987; Schulp et al., 2013; Hornung & Reich, 2015; Foffa et al., 2018), we suggest that only brachauchenines and polycotylids might have been possible niche competitors of mosasauroids.

With respect to the possible competitive interactions between mosasauroids, brachauchenines, and polycotylids, three hypotheses were considered:

1. First mosasauroids diversified following competition with plesiosaurs.

2. At least some mosasauroids competed with contemporary plesiosaurs or seized the opportunity and occupied their niches when they were in demise or became extinct.

3. The fates of plesiosaurs and mosasauroids were independent of each other (no suggested competitive interactions between mosasauroids and plesiosaurs).

Having the results of our Bayesian analyses in mind, we have focused on several traits related to niche occupation. Specifically, we have compared the body size and swimming abilities of mosasauroids and plesiosaurs, discussed the thermoregulation and metabolic rates in these groups, and considered their reproductive strategies and early life history. Available evidence shows that the earliest mosasauroids differed from plesiosaurs in all these traits. Earliest Late Cretaceous (Cenomanian–Turonian) plesiosaurs were large generalists, excellent swimmers, giving birth to a single large progeny, and with metabolic rates higher than those of contemporary mosasauroids. In turn, the first mosasauroids were small, semi-aquatic animals that gave birth to multiple offspring. These marked differences, however, began to blur in the Turonian. Around the time the last brachauchenines died out (middle Turonian; though see Madzia (2016) for possible younger occurrences) or immediately after the demise of the clade, mosasauroids experienced high evolutionary rates; they evolved first larger-sized taxa and apparently also first good swimmers (e.g., Bell, Barnes & Polcyn, 2013; Polcyn et al., 2014). No evidence suggests substantial changes in the metabolic rates or reproductive strategies of larger-sized mosasauroids. In fact, preliminary reports seem to indicate that even larger mosasauroid taxa gave birth to multiple progenies (e.g., Bell & Sheldon, 2004). Still, we refrain from drawing any far-going conclusions based on the available data as the evolutionary meaning of the discussed differences is unknown.

Considering the record of direct mosasauroid-plesiosaur interactions, the fates of plesiosaurs and mosasauroids were probably not directly independent of each other. At least two mosasauroid groups containing large taxa, the tylosaurines and advanced mosasaurines, clearly interacted with polycotylids (Everhart, 2004; Einarsson et al., 2010). However, the assertion that the extinction of the last brachauchenines had been accelerated by the diversification of the tylosaurine russellosaurinans or that tylosaurines competed for and eventually ‘usurped’ the niche previously occupied by pliosaurid plesiosaurs would be speculative. The early radiation of polycotylids during the Early/Late Cretaceous transitional interval apparently produced a burst of disparity in the clade (Fischer et al., 2018), which could have had an impact on brachauchenines as well. Fischer et al. (2018) noted that the increased disparity was not an aftermath of the extinction of ichthyosaurs and pliosaurids, and that the vanishing of the high disparity in polycotylids during and after the Turonian is consistent with a model of ‘early experimentation/late constraint’. Considering the suspiciously thalassophonean-like body plans in the mid-Cretaceous polycotylid plesiosaurs, such as Edgarosaurus muddi, Mauriciosaurus fernandezi, and Plesiopleurodon wellesi, that co-occurred with the last pliosaurids (at a time pliosaurids show low evolutionary rates), it is possible, nevertheless, that competition with polycotylids could have contributed to the extinction of Pliosauridae.

The potential role of mosasauroids in re-shaping the early Late Cretaceous marine environments is unclear. However, we speculate that the demise of brachauchenines and decrease in both, phylogenetic diversity and disparity in polycotylids around the time mosasauroids experienced high evolutionary rates, might have resulted in that some mosasauroids, such as tylosaurines, seized the opportunity and inhabited the niche previously occupied by robust-toothed short-necked plesiosaurs. However, further clade- and trait-specific studies are necessary in order to elucidate the patterns of niche occupation in Late Cretaceous seas.