Sauropod Trackway Morphometrics: An Exploratory Study Using Highway A16 Excavation at the Courtedoux-Tchâfouè Track Site (Late Jurassic, NW Switzerland)

Sciscio, Lara; Belvedere, Matteo; Meyer, Christian A.; Marty, Daniel

doi:10.3389/feart.2022.805442

ORIGINAL RESEARCH article

Front. Earth Sci., 17 March 2022

Sec. Paleontology

Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.805442

This article is part of the Research Topic Technological Frontiers in Dinosaur Science Mark a New Age of Opportunity for Early Career Researchers View all 9 articles

Sauropod Trackway Morphometrics: An Exploratory Study Using Highway A16 Excavation at the Courtedoux-Tchâfouè Track Site (Late Jurassic, NW Switzerland)

Lara Sciscio^1,2*^†

Matteo Belvedere³^†

Christian A. Meyer⁴^†

Daniel Marty⁵

¹JURASSICA Museum, Porrentruy, Switzerland
²Department of Geosciences, University of Fribourg, Fribourg, Switzerland
³Dipartimento di Scienze della Terra, Università di Firenze, Firenze, Italy
⁴Department of Environmental Sciences, University of Basel, Basel, Switzerland
⁵Natural History Museum Basel, Basel, Switzerland

Ichnogeneric classification of sauropod trackways is determined using qualitative and quantitative descriptions of morphological parameters. More recently, the validity of several of these parameters has been called into question (e.g., trackway gauge). This paper aims to test traditional and more novel landmark-based geometric morphometric (GM) analysis to describe sauropod tracks and trackways. The Courtedoux-Tchâfouè (TCH) tracksite, in NW Switzerland, has been used as a test site because it hosted several Late Jurassic sauropod track-bearing levels, over a short time period (∼40 ka), creating a time capsule of sauropod activity. Exploratory statistical analyses suggest that the TCH trackway dataset can be partitioned into clusters based on trackway parameters (e.g., width of pes angulation pattern/pes length) that demonstrate differences between trackways. Nevertheless, clustering reflected a moderate degree of intercluster similarity and a continuum in TCH trackway morphologies. The GM analysis of TCH pes impressions, from several trackways, indicates a similar morphotype and did not significantly differentiate pes tracks from various stratigraphic levels. Currently, the results indicate that neither linear nor landmark-based geometric morphometric methods strongly segregate tracks at TCH nor show a time-control (i.e., via stratigraphic level) over the short interval studied. The methods reinforce that sauropod pes impression shape is conservative and that combining morphological methods to be inclusive of sauropod trackway data is important. The TCH site demonstrates that variability may reflect morphological continuums and behavioural factors and does not unambiguously differentiate unique sauropod trackmakers. More exploration on the minimum number of trackway parameters as predictors of sauropod ichnotaxonomy is needed. Integrating complementary techniques will provide more insight and practical means of identification and delineation of sauropod tracks and their trackways.

1 Introduction

The Courtedoux-Tchâfouè track site (TCH) is one of several that bear sauropod tracks from the Late Jurassic (Kimmeridgian) tidal-flat deposits of the Reuchenette Formation of NW Switzerland (Ajoie district, Canton Jura). Documentation of the TCH track site, and several others, was completed by the industrious work of the Palaeontology A16 project. This team curated and preserved information from track-bearing palaeosurfaces under threat during the construction of the ‘Transjurane Highway’ (Swiss federal Highway A16; ∼2002–2011; Figure 1A). To date, the 254 sauropod trackways documented represent the largest repository of Late Jurassic sauropod trackways, globally. Furthermore, the ∼18,500 m² of documented track-bearing palaeosurfaces act as a narrow window into the tidal flats, associated with the Jura carbonate platform, that were periodically exposed and trampled during this time. A wealth of biogeographic information in addition to sources of ichnological, palaeontological, and palaeoclimatic evidence is hosted in these Highway A16 track sites (e.g., Meyer and Thüring, 2003; Jank et al., 2006a; Jank et al., 2006b; Marty et al., 2007; Marty, 2008; Marty et al., 2010; Schwarz et al., 2020).

FIGURE 1

FIGURE 1. (A) Courtedoux-Tchâfouè tracksite (TCH) tracksite (NW Switzerland; 47°24′13.7″N 7°01′21.6″E). (B) Lithostratigraphy of the Reuchenette Formation with emphasis on the biolaminate horizons bearing sauropod tracks/trackways (level numbers in boxes). The TCH track-bearing biolaminate levels of this study are shown to the right and numbered sequentially (levels 1000–1070 and indicated with green track symbols). (C) Section of trackway 1055-S4 from Level 1055.

The TCH tracksite (Figure 1A) lies ∼6 km southwest of Porrentruy in Canton Jura (NW Switzerland) and is documented in the catalogues du patrimoine paléontologique jurassien—A16 (Paratte et al., 2017). At TCH, the limestone sequence contains several discrete bio-laminated tidal horizons (level 1000; Figure 1B) bearing sedimentary structures, such as desiccation cracks, ripple marks, invertebrate traces, evidence for microbial mats, and last, but not least, tracks/trackways left by sauropod and theropod trackmakers (Figure 1C). The track-bearing bio-laminated layers (“levels” herein) cover a period of approximately 20–40 ka (Marty et al., 2007), fall within the mutabilis zone, and are numbered sequentially (levels 1000—1070; Figure 1B). Several studies have focused on the tridactyl tracks and trackways from these layers (e.g., Marty et al., 2007; Marty et al., 2017; Castanera et al., 2018), and of particular interest is TCH level 1000’s large theropod trackway, Megalosauripus transjuranicus (Razzolini et al., 2017).

Marty et al. (2007) first noted that the TCH sauropods tracks were similar in morphology and trackway pattern, independent of track size, with a pronounced heteropody and narrow gauges; with Marty et al. (2010) expanding on one well-preserved trackway from level 1055 (Trackway 1055-S4; Figure 1C). In general, TCH sauropod pes impressions are subtriangular/sub-oval in shape (longer than wide), with the greatest width recorded in the anterior third of the pes impression. Pes digital and claw impressions may or may not be present. Pedal rotation relative to the trackway midline is positive (i.e., outwardly rotated). Manual impressions have slight variations in shape between sites and trackways reported, from semi-circular, ellipsoidal to crescent (depending on the degree of overprinting by the pes) but are always wider than long. Manus placement is slightly exterior to the pes track, with manual supination angles shown as positive (outwardly rotated). The degree of manual supination is varied (<88°; Marty et al., 2012), and always recorded as greater than that displayed by the pes impression. TCH tracks were assigned to Parabrontopodus (Lockley et al., 1994; Marty et al., 2007, 2010). Parabrontopodus represents trackways with pronounced heteropody, narrow gauge, pes tracks of medium to large size (∼40–90 cm in pes length, longer than wide), rotated outwards from the trackway midline, and manus tracks, semicircular, with high outward rotation relative to pes tracks (Lockley et al., 1994; Marty et al., 2010). The type ichnogenus and ichnospecies, Parabrontopodus mcintoshi, was described from the Upper Jurassic Morrison Formation (United States) and is commonly believed to have been made by diplodocid sauropods (Lockley et al., 1994; Meyer and Thüring, 2003; Henderson, 2006; Marty et al., 2007; Marty et al., 2010). Despite the scarce Swiss sauropod body fossil record, the fragmentary remains of several non-neosauropod eusauropods (Amanzia greppini) are known from the lower Reuchenette Formation (early Kimmeridgian age; Meyer and Thüring, 2003; Jank et al., 2006a; Schwarz et al., 2020) stratigraphically older than the TCH site. The current A. greppini material suggests small, adult sauropods < ∼10 m in body length but with a ∼35% size difference between the largest and smallest individuals in the group (Schwarz et al., 2020).

Altogether, the current ichnotaxonomy of the Reuchenette Formation’s sauropod tracks places them, broadly, as Parabrontopodus isp. (e.g., TCH and Courtedoux tracksites; Marty et al., 2003; Marty et al., 2010; Marty et al., 2013), Brontopodus isp. (e.g., Chevenez–Combe Ronde, Courtedoux–Pommerat sites; Marty et al., 2003; Marty et al., 2010), or Breviparopus/Breviparopus-like (e.g., Lommiswil and surrounds; Meyer, 1990; 1993). More recently there is some debate as to whether Parabrontopodus and Breviparopus are synonymous (Marty et al., 2010; Castanera et al., 2014; Xing et al., 2017). Trackways of Breviparopus are narrow gauge, with an oval to bell-shaped pes impression, outwardly rotated, intersecting the trackway midline and the manus impressions located further from the midline than the pes (Dutuit and Ouazzou, 1980; Marty et al., 2010). The manus and pes morphology is like Parabrontopodus and shows pronounced heteropody. A Breviparopus/Parabrontopodus plexus has been suggested by Santos et al. (2009) because of the gross morphological similarities between the two trackway types. The lack of a formal erection of Breviparopus, despite detailed description by Dutuit and Ouazzou (1980), compounds the problem.

Defining sauropod ichnogenera such as Breviparopus, Parabrontopodus, or Brontopodus ispp., has relied on track shape and trackway configuration, and the latter heavily skewed towards trackway gauge (narrow, medium, wide) despite inherent issues with this metric (see Marty, 2008; Castanera et al., 2016; Meyer et al., 2018, 2021). In recent years, qualitative descriptions in ichnology have been increasingly influenced (generally) by geometric morphometrics to quantitatively assess morphology (e.g., Rodrigues and dos Santos, 2004; Lallensack et al., 2016; Abrahams et al., 2020). As such, the main aim of this study is to assess whether TCH track-bearing surfaces (“Levels 1000-1070”) all bear the same ichnogenus, and to examine if there is variation in morphology of the TCH Parabrontopodus isp. using linear-based and geometric morphometrics methods. The TCH site with its density of trackways and trackway-bearing levels, through a narrow timeframe, forms the perfect testing ground. The expected outcome establishes that any morphological variation is within the variation error of the population and that all tracks group together with low variability. Our alternate hypothesis is that the TCH track-bearing surfaces show strong variation in pes morphology and trackway morphometrics. The latter’s variations are time-transgressive across track-bearing levels and reflect different sauropod ichnogenera and sauropod faunas utilizing the area. Importantly, if morphological distinct ichnological groups are time transgressive at TCH this could reflect the likely mixing of sauropod faunas.

2 Method and Material

Standard track and trackway acronyms and measurements were used following Marty, (2008; Figure 2). Acronyms used here and in-text: PL: pes length; PW: pes width; α or av. pes. rot: pes rotation; ML: manus length; MW: manus width; β or av. manus sup.: defines manus rotation or supination; PS: pes stride length; PPA: pes pace angulation; WAP: width of the pes angulation pattern; MS: manus stride length; MPA: manus pace angulation; WAM: width of the manus angulation pattern; digW: digit width. Important to note that glenoacetabular distance (GAD), the distance between the center of the hip joint and the centre of the shoulder joint, is measured and not calculated. The ratios utilized include IPS/IMS (index of pes size to index of manus size; Marty, 2008: [PL×PW]/[ML×MW]^0.5), WAP/PL, and WAP/WAM. The WAP/PL ratio has been used as an estimate of trackway gauge, after Marty (2008): for narrow- (≤1), medium- (1–1.2), and wide- (≥1.2) gauged trackways.

FIGURE 2

FIGURE 2. (A) Trackway parameters (adapted from Marty, 2008). (B) Landmark placement on the pes. Acronyms used here and in text: PL: pes length; PW: pes width; α: pes rotation; ML: manus length; MW: manus width; β: manus rotation; PS: pes stride length; PPA: pes pace angulation; WAP: width of the pes angulation pattern; MS: manus stride length; MPA: manus pace angulation; WAM: width of the manus angulation pattern; digL: digit length.

To estimate the trackmakers’ size, the measured glenoacetabular distance (GAD) and hip height (h) are commonly used to define trackways (see discussion in Stevens et al., 2016). Sauropod hip height (h) has been estimated via three main methods utilizing pes length: h = 4×PL (Alexander, 1976), h = 5.9×PL (Thulborn, 1990), and h = 4.586×PL (González Riga, 2011). The latter is based on the articulated hind limb reconstruction of a titanosaur (González Riga, 2011). Ultimately, the use of the conventional h = 4×PL (Alexander, 1976) was adopted for all subsequent analyses, and this has been previously supported by Tschopp et al. (2015). However, González Riga (2011) equation did not lend to higher variability in the dataset analysis. Altogether, GAD and h rely on several assumptions that cannot be directly tested with the current body fossil record (see Schwarz et al., 2007 regarding sauropod limb size estimates based on preserved A. greppini humerus cartilage capsule), and while h is an independent means to predict size, it is dependent on pes length (Tschopp et al., 2015; Stevens et al., 2016).

All measurements were taken in the field and form part of the Palaentology A16 (PALA16) TCH dataset, now hosted/managed by the Jurassica Museum, with selected photographs of TCH pes used for geometric morphometric methods. Linear-based measurements (traditional morphometrics) were normalized and analyzed via correlation analysis, cluster analysis (Xu and Tian, 2015; Tibshirani et al., 2001), and principal component analysis (PCA). Trackway data was standardized and trackways with missing values were excluded for correlation and clustering analyses. Listwise deletion of incomplete trackway data (i.e., those with missing values from a given trackway) was undertaken so as not to compromise further statistical analysis (Kang, 2013), and missing data fields are a result of incomplete step cycles/missing impressions from the trackway. Estimation of the optimal number of clusters was determined using the Gap statistic (run at B = 500), Elbow, and Silhouette methods (Supplementary Figure S2; Supplementary Material). This was done to determine the most parsimonious k, beyond which increases in the value would be insignificant, although it should be stressed that the choice of k is not stable. In the resulting HK- and K-means clustering, the nstart number (the number of configurations) was set at 50 to generate 50 configurations that then provide an averaged centroid. This was conducted in R version 4.0.4 (2021-02-15) (R Core Team, 2021) with figures produced using the package ggplot2 (Wickham, 2016).

To define trackways, parameters that reference the trackway gauge, track heteropody, and size of the trackmaker have been utilized. This stems from their repeated use in literature as key parameters to study sauropod trackways and to assign specific ichnogenera (particularly Brontopodus, Breviparopus, and Parabrontopodus ispp.). These measurements also use a combination of track and trackway data, precluding the need to use “high-quality”/high-preservation grade tracks (as in Belvedere and Farlow, 2016; Gatesy and Falkingham, 2017; Marchetti et al., 2019; 2020; Falkingham and Gatesy, 2020), but necessitating a trackway. Specifically selected for cluster and PCA analysis are variables: body size (h) using Alexander (1976) equation (h = 4×PL), heteropody (IPS/IMS), manual supination, and trackway pattern as defined by the MPA, PPA, and ratios of WAP/PL (gauge; Marty, 2008), and WAP/WAM (distance of manus to trackway midline vs. that of the pes). In isolating these key variables, further study may provide more insight into their success in distinguishing ichnogenera. Other variables more specifically designating the speed of the animal such as pace, progression, and stride, and further ratios and calculations from them were not included so as not to be compounded by the relative speed of travel.

Cluster Analysis was followed by Principal Component Analysis (PCA), using the factoextra R packages, on the previously selected variables. Given several of these variables are highly correlated (see correlation matrices), the PCA is used to reduce variables and redundancy because of correlation. Altogether the biplot(s) produced here describe variances in and correlations between the variables. Distances between trackway points/ichnogenera indicate differences between trackway data. Conversely, vector length describes the degree to which that variable explains the distribution of the data, and the angle between vectors approximate their correlations. Higher cos² values (square cosine) indicate variables/trackways are well represented in the principal components. The factor map shows the cos² of each variable and illustrates the quality of representation of that variable. The scree plot (Supplementary Figure S3; Supplementary Material) illustrates the variance in the data (%; referenced as dimensions). Variables and their contributions to clustering were extracted using the get_pca_var () function in R. Data were compared to that presented in Marty et al. (2010) for Brontopodus from the Swiss sites of Chevenez—Combe Ronde and Courtedoux—Pommerat (CRO, PMM, respectively), Parabrontopodus from Courtedoux –Tchâfouè (level 1055, trackway S4; Figure 1C) and the Moroccan Iouaridène track site Deio Lav (TCH, Deio Lav, respectively), and Breviparopus from the Swiss site Courtedoux—Sur Combe Ronde and the Moroccan Iouaridène site Deio (SCR, Deio, respectively).

Currently, geometric morphometrics is an approach not typically used in ichnology but is slowly gaining popularity (Moratalla et al., 1988; Rodrigues and dos Santos, 2004). It allows the analysis of information on spatial covariation between anatomical landmarks (Rohlf and Marcus, 1993). The placement of landmarks (i.e., landmark configuration) assists in defining information on the shape of the track while excluding size. This makes an ideal means to assess tracks in an unbiased way, discounting size as a criterion and a means to understand shape and the sources of morphological variation. For this geometric morphometric analysis, 10 landmarks were digitally placed onto 28 pes tracks with higher morphological preservation (i.e., distinguishable digits, clear track margins, etc. Figure 2B) from the various levels at the TCH site (Figure 1B). Methodologies were followed according to Rodrigues and dos Santos (2004), Webster and Sheets (2010), Viscosi and Cardini (2011), and Savriama (2018). All photographs used were perpendicular to tracks.

Landmarks were placed to define the overall shape of the track (Figure 2B): Landmarks (1, 3, 5) reflect digit I, II, and III apices, respectively (excluding claw impressions); landmarks (2, 4, 6) mark the central point of the hypices between digits I^II, III, and IV, respectively; landmark (7, 9) mark the width of the track but also, respectively, the base of digit IV and I, landmark (8) is placed at the “heel” of the track, perpendicular to landmark (1), and landmark (10) marks the base of digit I, if an imaginary line projected from landmark (1) dissected the digit width-wise. Variability in the hypices of the tracks is recognized and therefore the utility of those landmarks is taken into consideration. Landmark digitization utilized the freeware TPS Series, whereby TPS files were generated from images using tpsUtil and imported into tpsDig2 (Rohlf, 2013, 2015) where all landmarks were digitized in the same order. Analysis of shape variation was then completed in the freeware program MorphoJ (Klingenberg, 2011) and landmarks were superimposed by Generalized Procrustes Analysis (Rohlf and Slice, 1990). Procrustes ANOVA, multivariate regression analysis, Principal Component Analysis (PCA), and Canonical Variate Analysis (CVA) of Procrustes shape coordinates were utilized to assess variation in shape on different levels. Independently, a second test using a curve (semilandmarks) was used, tracing the outline of the track from the hypex of digit III to the hypex of digit I (from landmarks 6—2) creating 33 points in two dimensions, and the same procedures, as outlined above, were followed. The curved semilandmarks were converted using tpsUtil 1.8 (Rohlf, 2013, 2015).

3 Results

Qualitatively, TCH footprint morphology (as generalized across stratigraphic levels) records oval pes and semi-circular manus, with the pes and manus impressions moderately sized (av. PL:39.9 cm; av. PW 34.8 cm; av. MW: 23.7 cm; av. ML: 14 cm; Table 1). Pedal rotation and manual supination are outwards relative to the trackway midline with that of the manus impression greater (av. 41.8°, max. 68°) than the pes (av. 20.1°, max. 49°; Table 1). An outlier in terms of track size is present on layer 1060 (hereafter referred to as level 1052, trackway 1052S1:UTw), with a PL and MW of 101.7 and 74.3 cm, respectively; and the smallest tracks were also recorded on the same level 1060 with average PL and PW of 27 and 18 cm, respectively. Trackway S1 from level 1052 (1052S1:UTw) is an underprint of trackway 1060S5, from the overlying level 1060, but as no trackway measurements were made of 1060S5, the former’s have been has been utilized in the current analysis. The WAP/PL for the various levels ranged from narrow-gauge (averaged trackways from levels 1003,1070), medium-gauged (levels 1000, 1014, 1015, 1020, 1025, 1030, 1055) to wide-gauged trackways (levels 1052, 1060). The speed of the averaged trackways for all levels was 2.0–3.5 km/h (walking/trotting) with level 1070 showing a speed of 10.2 km/h (running) and lacking manus impressions.

TABLE 1

TABLE 1. Tchâfouè (TCH) track site mean measurements. Measurements in cm unless specified and acronyms as per text.

3.1 Traditional Morphometrics

Linear data utilized for the traditional morphometric analysis is presented in Table 1 (see also Supplementary Table S1). From the 72 trackways available at the TCH site, across all stratigraphic levels, we have utilized the complete data from levels 1000, 1003, 1010, 1014, 1015, 1020, 1025, 1030, 1052, 1055, and 1060 (Figure 1B) for correlation, cluster, and principal component analysis. Levels 1000 and 1015 having the highest density of trackways [e.g., layers 1000 (n = 26), 1003 (n = 1), 1010 (n = 1), 1014 (n = 1), 1015 (n = 10), 1020 (n = 6), 1025 (n = 1), 1030 (n = 5), 1052 (n = 1), 1055 (n = 7), 1060 (n = 8), and 1070 (n = 2)]. Box and whisker plots (Figure 3) show the degree of variability within and between track-bearing levels as measured by the thickness of the boxes and the skewness of the median (as well as the length of whiskers and associated outliers). There were no stratigraphically well-defined trends reflected in the box and whisker plots.

FIGURE 3

FIGURE 3. (A–L) Box and whisker plots of all TCH track site data by track and trackway parameters for each stratigraphic track level (Level 1000-1070). Acronyms as per text.

3.1.1 Correlation Matrices of Track and Trackway Parameters

Correlation matrices illustrate the relationships between multiple TCH trackway variables (Figure 4; Supplementary Figure S1). Correlation analysis used the Pearson parametric correlation test and the Spearman rank correlation analysis. The main findings for TCH trackway data are summarized as follows: the PL and MW, as expected, have positive linear relationships with pes/manus pace, pes/manus stride, speed, and correspondingly larger GAD and h (Figure 4A; Supplementary Figure S1). Neither pes nor manus size shows a relationship to the degree of manus or pes rotation during locomotion.

FIGURE 4

FIGURE 4. (A,B) Correlation matrices for TCH track/trackway parameters showing pair-wise scatter plots and correlation coefficients; p-values (<0.001, <0.01, <0.05, 0.1, 1) correspond to symbols [“***”, “**”, “*”, “.”, “”]). Acronyms used as per text. Av.pes.rot: average pes rotation; Av.manus.sup.manus sup.: manus supination.

Body size, as reflected in h, has a highly statistically significant negative relationship with MPA and WAP/PL (with p values < 0.001), but a significant positive relationship with the PPA (p < 0.001) and PS (p < 0.0001). Several of these relationships are likely constrained by the calculation of h. GAD is strongly negatively correlated with WAP/PL gauge ratio (r = −0.52, p =<0.001) and positively correlated with WAM, WAP and PPA (Figure 4B).

Manus supination is weakly negatively correlated with PL (r = −0.36, p = <0.05) and MW (r = -0.27, p =<0.05; Supplementary Figure S1). There is no relationship between pes rotation and manus supination nor with either manus/pes rotation and WAP/PL (Figure 4A). Manus supination is negatively correlated with GAD (p = <0.01), h (p = 0.05), WAP (p = <0.01), PS (p = <0.01), and speed (p = 0.05). Manus stride shows a highly statistically significantly negative correlation with manus supination (p = <0.001). The average pes rotation, contrastingly, is only negatively correlated with the MPA (p=<0.01) and speed (p = 0.05), and otherwise shows no significant relationships.

MPA and MS are weakly negatively correlated (r = −0.60, p=<0.001; r = −0.36, p=<0.05, respectively; Figure 4), with no relationship between PPA and MPA (r = −0.062). Both MPA and WAM are strongly negatively correlated (r = −0.85, p = <0.001) in TCH trackways. As is expected, the WAM and MW are strongly positively correlated (r = 0.83, p = <0.001). A strong positive linear correlation of WAP with PS (r = 0.65, p = <0.001), and the width of the manus angulation pattern (WAM; r = 0.53, p = <0.001) is noted.

Speed (Alexander, 1976) is significantly positively correlated with PPA and MPA (p values < 0.05), and PS is positively correlated with WAP and WAM (p values < 0.0001). Trackway gauge ratios of WAP/PL and WAM/MW were weakly negatively correlated (r = −0.35, p =<0.05), and WAP/PL is significantly negatively correlated with both PPA and PS (p = <0.001), and significantly positively correlated with MPA (p < 0.001).

3.1.2 Cluster Analysis

Partitioning of the TCH site data via K-Means and HK-Means analysis was conducted using selected variables: manus supination, body size estimate via hip height (h; Alexander, 1976), heteropody (IPS/IMS), and trackway gauge (WAP/PL, WAP/WAM, PPA, MPA; Figure 5). Hybrid hierarchical clustering and the K-Means methods generated a clustering solution showing the relationship between levels and their trackways (Figures 5A,B). The most parsimonious number of clusters, k, via the Silhouette and Elbow Method, was 3, with cluster sizes of 12, 19, and 9 trackways, and a total variance of 48.8%. Assessment of whether the dataset contains meaningful clusters was achieved by utilizing the Hopkins statistic to determine the spatial randomness, and with an H-value of 0.61 (using get_cluster_tendency func in factoextra on TCH variables), the TCH dataset is clusterable. Substituting hip height calculated using González Riga (2011) didn’t alter the clustering or the total variance and substituting for the measured GAD did not significantly alter the clustering (k = 3; 11, 19, 10 cluster sizes) but lowered the total variance by 2% (45.5%).

FIGURE 5

FIGURE 5. (A) HK-Means cluster analysis plot of all TCH trackway data using k = 3 via the Silhouette and Elbow Method. (B) Dendrogram of HK-means. (C) HK-means cluster analysis plot of TCH data and ichnogenera from Marty et al. (2010) with k = 4.

These steps were repeated with the incorporation of data taken from other Swiss and Moroccan tracksites (Marty et al., 2010; Table 2), specifically: Breviparopus (SCR, Deio), Brontopodus (CRO, PMM), and Parabrontopodus (TCH, Deio Lav). HK-means clustering grouped this dataset into 4 clusters (14, 17, 18, and 4 trackways) with a total variance of 59.8% (Figure 5C). The HK-means clustering places Breviparopus and Parabrontopodus with the majority of the TCH dataset and distinctly segregates Brontopodus trackways and two TCH trackways (1003S1, and 1052-S1:UTw) (Figure 5C).

TABLE 2

TABLE 2. Ichnogenera for comparison with the TCH site. Data taken from Marty et al. (2010) with acronyms as per text.

3.1.3 Principal Component Analysis

PCA was employed to explore TCH trackway pattern variation, and individual trackways were plotted using the selected variables discussed (Figure 6). Figure 6A illustrates the spread of TCH trackways across PC1 and PC2, with PC1 explaining 49.4% of the total variation. PC1 was impacted, primarily, by WAP/WAM and WAP/PL (trackway gauge parameters), and less so by MPA and h (Figure 6A). WAP/WAM, MPA, and WAP/PL show negative loadings on PC1 whereas PPA and h are positively loaded. Gauge parameters WAP/PL and WAP/WAM are positively correlated, define partially overlapping and equal vector lengths, and increases in these variables are associated with increases in MPA. PC2 is defined by the manus supination angle (strong positive loading on PC2) and by heteropody (as measured by IPS/IMS). Manus supination is negatively correlated with IPS/IMS along PC2 (Figure 6A). The trackways on the same side of a given variable have a high value for this variable, e.g., trackways on level 1060 have higher average manus supination values and WAP/PL ratios than most of the level 1000 and 1020 trackways.

FIGURE 6

FIGURE 6. (A) Biplot analysis for trackway similarity. (B) TCH data overlaid with that of Marty et al. (2010). Highly correlated trackway components show vectors pointed in the same direction. The maps to the right of the biplots reference the quality of representation of each variable to the PCs (i.e., the squared cosine, cos², value). Higher cos² values (darker blue) for a variable indicates its better representation for that pC. Correlated trackway components and trackway samples are in the same quadrant.

Figure 6B shows the distribution of the TCH trackways with the additional Swiss and Moroccan ichnogenera from Marty et al. 2010; Table 2. PC1 explains 34.7% of the total variation and is impacted by both trackway gauge ratios (WAP/WAM and WAP/PL) and PPA (Figure 6A). PC2 is defined by hip height, manus supination, and MPA. Clusters 1, 2, and 4 are spread along PC1 and thus reference changes in trackway gauge (WAP/WAM, WAP/PL), PPA, and heteropody. Furthermore, trackway gauge rations are primarily driving the grouping of certain trackways, i.e., Brontopodus isp. and trackways of TCH cluster 1, that are influenced by WAP/PL. TCH cluster 1 is more strongly influenced by MPA, whereas clusters 2 and 4 (containing Breviparopus and Parabrontopodus ispp.) are characterized by heteropody, narrow-gauge (opposing WAP/PL), and PPA. Less acute manus supination values affect several trackways from cluster 4 and cluster 1 along the PC 2 (Figure 6B) and less so clusters 2 and 4; but this variable is less significant than PPA, MPA, WAP/WAM, and WAP/PL.

3.2 Geometric Morphometrics

Landmark coordinates were superimposed by generalized Procrustes analysis (Rohlf and Slice, 1990) to generate Procrustes coordinates, and deviation around consensus coordinates is shown in Figure 7.

FIGURE 7

FIGURE 7. (A) Procrustes coordinates with large blue dots (centroids), labelled 1–10, representing the consensus coordinates and the deviation around them (small dots). (B) Curve landmarks (landmarks 11–33).

3.2.1 Size and Shape Variation

In using 10 landmarks (Figure 7A), the difference in centroid size (used as a proxy of size) between TCH pes specimens were tested by Procrustes ANOVA and was found not to be significant (p = 0.46). The ANOVA results show that the individual track variation is not significantly larger than the error. In testing for allometry, the influence of size on shape, multivariate regression of the Procrustes coordinates onto log centroid size was used (permutation test with 10,000 rounds; Figure 8A). This shows an insignificant, but positive linear correlation (p = 0.22) between these parameters indicating no constraint on shape by size (no allometry). Thus, 4.74% of the shape variation can be explained by size and as such, we cannot reject the null hypothesis (i.e., that all tracks are Parabrontopodus isp.). Changes were noted with the heel region and the projections of digit I and digit III (Figure 8B).

FIGURE 8

FIGURE 8. TCH geometric morphometric plots. (A) Multivariate regression of the TCH pes shape onto size (as logarithm of centroid size) classified by Level. (B) Deformation grid reflecting shape changes (blue dots represent landmarks labelled 1–10). (C) Principal components plot of PC1 vs. PC2 of TCH pes tracks pooled by Level (stratigraphic level). (D) Shape changes are visualised by the transformation grids. (E) PC2 vs PC3 of pes tracks from TCH site pooled by Level. (F) Shape changes are visualised by the transformation grids. The black numbered circles (centroids) indicate the average landmark coordinates of the mean shape of the sample; the sticks indicate the changes in their relative positions and are emphasized by the green wireframe. Numbers reference landmarks.

3.2.2 Geometric Morphometric Principal Component Analysis

To quantify and assess shape variation in the dataset, PCA of the Procrustes superimposed coordinates from the 28 tracks used across the stratigraphic levels was implemented. The PCA produced 16 components of which the first 4 components accounted for 78% of the total variation.

Figure 8C shows the PC1 axis, accounting for 40% of the shape variation, doesn’t segregate any pes impressions by level. The PC1 axis is strongly loaded or shows the maximum variation, in variables that are related to the pes hypices and base of digit IV, i.e., specifically, the hypices between digits I^II and II^III (landmarks 2, 4), and landmark 7 which also contributes to the measurement of the width of the pes (Figures 8C,D). This demonstrates that with a minor increase in extension of digit I (directing digit I anteriorly; see also PC2, Figure 8C) there is a decrease in the extension of digit III as its orientation shifts more laterally. Conversely, the base of digit IV moves upwards and medially. PC2 (18.9% of the variance) exhibits variation between landmarks 1, 4, 10 which are related to the extent of digits I and II and the total width of digit I. The latter variables control the loading on the second axis (PC2) and reflect longer digit extensions and impressions of digit I. For Figure 8C, two specimens (S16 and S13) from level 1000 form outliers to the left- and right-hand sides of the graph along PC1. Both these specimens display expulsion rims around the digit and heel regions which may play a factor in their distribution and loading as this may emphasize (e.g., deeper digit impressions making more marked hypices, S16) or de-emphasis (e.g., infill from an expulsion rim, S13 digit I) variable features. PC3 (Figures 8E,F) explains 10.9% of the variation, and this axis is loaded by the position of all the digit tips (landmarks 1, 3, 5, and 10) and the maximum width of digit I.

3.2.3 Canonical Variate Analysis

CVA is a tool to assess distinction between groups by assessing the differences of between-group and within-group variance (Gardner et al., 2006; Klingenberg, 2011). At TCH, pes shape differences for each level were scattered on the first two canonical variate axes (CV1 and CV2) and show a single group spread across CV1. The first two canonical variates (Figure 9) accounted for 75.3% of the total variation of the TCH level pes impressions. Along the first axis (CV1 = 49.7% of among-group variance) the main cluster is subdivided along a gradient by level that would suggest slight differences in shape at distinct stratigraphic levels. Pairwise differences in mean shapes are then tested for statistical significance by using permutation tests (10,000 rounds) with Mahalanobis distances and Procrustes distances (Klingenberg, 2011). Nearly all the Mahalanobis distances obtained by pairwise comparisons of the 7 levels were not significantly different (p > 0.0001–p > 0.01), ranging from 11.9 (1060–TCH1368, p = 0.1) to 3.2 (samples from 1020 to 1000, p = 0.0003). Procrustes distances also did not show any significant differences between several levels (permutation 10,000 rounds in MorphoJ: p = > 0.17), ranging from 1.0 (1060 and 1030) to 0.17 (1020 and 1030). Shape changes along CV1 axis were found with landmarks 1, 2, 9, and 10, and along the CV2 axis with 1, 2, 5, 6, 7, 8, 9, and 10 (Figure 9). The former is associated with digit I width and the apex of digits I, II, as well as the outer width of the track (landmark 9). These appear to reflect gradational changes along-track shape space. Two samples (TCH1368, TCH4265, from levels 1060 and 1030, respectively) were outliers.

FIGURE 9

FIGURE 9. (A) Scatter plot showing the variation in shape of pes tracks of ?Parabrontopus isp. from TCH site along the first two canonical variate (CV1 and CV2) axes. There is a gradient of separation of track shape by site level across CV1. (B) Shape changes, illustrated by the deformation grids, from overall mean along CV1 and CV2 axes in positive directions. The numbered circles indicate the landmark locations in the mean shape of the sample; the sticks indicate the changes in the relative positions of the landmarks.

3.3 Curved Landmarks

In creating a curve (Figure 7B), 33 landmarks were converted for the 28 pes impressions used. Thisrepresents a low pes sample size/landmark ratio, as often happens with paleontological datasets, and as such interpretation of the results should be cautious. The difference in centroid size was tested by Procrustes ANOVA (p = 0.0012) and was found to be significant. In assessing allometry (Figures 10A,B), regression analysis onto log centroid size was used (permutation test with 10,000 rounds) and shows a statistically insignificant, but positive linear correlation (p = 0.09) between these parameters indicating no constraint on shape by size. A pooled regression shows that size explains 17.3% of the shape variance and the regression is significant (p = 0.0001).

FIGURE 10

FIGURE 10. (A,B) Multivariate regression of the TCH pes shape, via curving landmarks, onto size (as the logarithm of centroid size) classified by Level, and reflecting shape changes. (C,D) Principal component analysis plot of PC1 vs. PC2 and shape changes shown on transformation grids. (E) Canonical variate analysis scatter plot showing the variation in shape of pes tracks along the first two canonical variate (CV1 and CV2) axes.

The PCA produced 27 components of which the first 16 accounted for 97% of the dataset’s total variation, with the first 4 components accounting for 73% of the total variation. The PC1 and PC2 expressed 37.6 and 18.788% of the total variation, respectively (Figure 10C). The PC1 axis displays the variation related to the external margin of the pes and “heel” region, i.e., specifically, the landmarks 16–24, and landmark 8 (Figures 10C,D) whereas PC2 displaces variation between landmarks 3, 4, and 5 in the digit apex and hypices linked to extension and splay of digits II and III and the heel region. These PCs show shape variation related to an elongation of the pes and compression of the width of the heel region. In the CVA, grouping by level, the first two canonical variates accounted for 64.3% of the total variation (CV1 = 36.8%, CV2 = 27.5%), and tracks were not clustered into distinct groups that would suggest differences in shape at distinct levels. Nearly all the Mahalanobis distances obtained by pairwise comparisons of the 7 levels were not significantly different (permutation 10,000 rounds in MorphoJ: p > 0.0001 and p > 0.01).

4 Discussion

Sauropod tracks are found, often, forming significant trackways that can assist in characterizing the tracks and trackway patterns, and, potentially, the trackmaker’s locomotive capabilities (e.g., Henderson, 2006). As Late Jurassic sauropod pes impressions are visually similar and this is reinforced by the conservative pes skeletal record, it limits their utility in distinguishing ichnogenera. Thus, both track and trackway parameters have been and must be used in comprehensive descriptions. Quantitative trackway parameters (e.g., trackway gauge) and qualitative/preservational characteristics (e.g., presence of manus digit I claw impressions) have been typically used to describe trackways and new ichnogenera (e.g., Day et al., 2002; Santos et al., 2009; Torcida Fernández-Baldor et al., 2021). In general, much emphasis has been placed on the trackway gauge (i.e., wide- vs. narrow-gauges; Farlow et al., 1989, Farlow, 1992; Lockley et al., 1994, see discussion in Moratalla, 2009) for identifying different sauropod ichnogenera. The increasing awareness that a single individual can vary its gauge, from narrow to wide, within the same trackway (e.g., Meyer et al., 2018) illustrates this metric’s unreliability for ichnotaxonomical attribution and suggests some sauropod ichnogenera need revision (viz Xing et al., 2017). Nevertheless, in this study, we have applied repeatedly used track/trackway parameters (i.e., heteropody, trackway gauge, manus supination) to assess their utility within a larger dataset, and their suitability/distinctiveness in defining sauropod trackways. Concurrently, the pes geometric morphometric data characterized shape and size differences of these TCH sauropod tracks.

4.1 Traditional Morphometrics

Overall, the TCH box and whisker plots (Figure 3) illustrate that across stratigraphical levels 1000-1070 (through time), there is little variability in certain trackway parameters (e.g., PL, MW, speed, or hip height) with well-constrained boxes and median positions showing no/little skew. However, track level 1060 (and by extension underprint of level 1052) is an outlier with smaller and very large pes and manus impressions. Conversely, other parameters are not well constrained, and boxes for manus and pes rotation reflect variability within and between track levels. Measures of trackway gauge, such as WAP/PL, WAM/MW, MPA, together with PS, and MS show similar degrees of variability within and between track levels as evidenced by both median and box positions and thicknesses with skewed whiskers and outliers. It was interesting to note the difference displayed by measures of pes and manus pace angulation for TCH trackways (Figures 3G,H). PPA shows consistency (well-constrained boxes) and a slight decline across stratigraphic levels, as measured by decreasing medians, while the measure for the manus (MPA) shows no trends, is highly variable within/between levels, and juxtaposes the PPA. This suggests underlying factors, such as manus supination or locomotive behavior, may influence MPA more strongly than PPA. Altogether, trackway parameter variation can be a function of several overlapping factors, ranging from substrate conditions, changes in velocity, to (unmeasurable) behavioral characteristics.

The correlation matrices indicate the significance of relationships between parameters (Figure 4, Supplementary Figure S1). Body size estimates established that GAD was positively, linearly correlated with manus/pes size, PS, PPA, WAP, and WAM. This was to be expected as longer (and by inference larger) animals (i.e., longer GAD) would also, hypothetically, have incrementally larger feet/pace/stride. However, it was strongly negatively correlated with the WAP/PL-gauge ratio, suggesting trackway gauge and body size are not associated. Body size via estimated hip height (h) demonstrated a highly statistically significant negative relationship with the MPA (p < 0.001) and manus supination (p = 0.1) illustrating a likely disconnect from this estimate and the manus orientation/forelimb. Conversely, h shares a significant positive relationship with the PPA (p < 0.001) implying greater h lends to larger PPA. Gauge ratio and h show a statistically significant negative relationship (p < 0.001) but this can be influenced by how both gauge and h are built.

Generalized trackway pattern relationships show that the larger WAP corresponds to acute manus supination and negligible pes rotation (Figure 4, Supplementary Figure S1). Statistically significant, positive linear correlation of PS with WAP and WAM (p=<0.0001) shows that with increasing stride lengths, the WAP and WAM will be larger (or, potentially, display more variation). Both MPA and WAM are strongly negatively correlated (p=<0.001) in trackways. This suggests that increases in the WAM correspond with decreases in the MPA, but neither has a significant effect on manus supination (Figure 4B). Although there is not a strong relationship between WAM, MPA, and manus supination, larger WAPs are related to more acute manus supination angles (p=<0.05) and correspond to wider WAMs. Previously, MPA has been weakly negatively correlated with MW in sauropodomorph trackways (Farlow et al., 2020), and within the TCH dataset a strong negative correlation (r = −0.60, p=<0.001) between MPA and MW is noted. Less linear trackways have been related to low MPA and increasing manus size (cf. Farlow et al., 2020). The strong negative correlation between MPA and WAM suggests that with increases in WAM there are corresponding decreases in the MPA. Although there is no relationship between WAM, MPA, and manus supination, larger WPAs are related to more acute manus supination angles and correspond to wider WAMs.

Manus supination, although highly variable within and between TCH trackways, was suggested by Stevens et al. (2016), together with its morphology and track placement, to be a diagnostic trait. Previously, Lallensack et al. (2018) and Farlow et al. (2020) found that manus supination decreases with increases in manus size, but TCH manus supination is not affected by pes or manus size nor by the rotation of the pes (Figure 4, Supplementary Figure S1). Manus supination is increasingly affected by the size of the animal (as approximated by the GAD; p < 0.05), and the degree of supination decreases as GAD, WAP, and PS increase. Manus supination is only significantly affected by the MS, with increasing manus stride lengths promoting decreases in supination angle. In Lallensack et al. (2018) atypical trackway (A9; Figure 6A), longer strides and increased speed reduce manus supination angles. TCH manus supination is likely driven by a combination of locomotion speed, manus stride length, and/or behaviour given that there are not significant changes in averaged relative speeds (Figure 3; Supplementary Figure S1) for the trackways studied.

Variation in speed and corresponding manus supination variability has been previously considered a function of the active use of the forelimb for propulsion (Ishigaki and Matsumoto, 2009; Xing et al., 2015b; Lallensack et al., 2018 and refs therein). Collectively, data suggests that the manus undergoes rotation from a posterolateral to more anterior position with increasing locomotive speed and/or MS length. The variability in manus position and rotation along a trackway appears reactive to momentary and adaptive changes in speed and/or MS that are likely concurrent with shifts in substrate rheology/behavior, etc. Generally, columnar sauropod forelimbs are considered weight-bearing structures built to minimize excessive movement, and anterior manus rotation could be achieved by angled forelimb posture (at least for wide-gauge titanosaur trackmakers; Bonnan, 2003). However, Bonnan et al. (2013) indicate that thick articular cartilage in sauropods, while reducing joint stresses, may have allowed some range of forelimb movement during locomotion, and this is also discussed by Xing et al. (2015b) in relation to constraining manus placement in trackways. Importantly, in the study area, a cartilaginous articular capsule is known from A. greppini forelimb (Schwarz et al., 2007), covering up to one-third of the distal side of the humerus and as thick as ∼3–5 cm, thus reinforcing the possibility for forelimb mobility.

Sauropod locomotion based on selected Highway A16 trackways was characterized by Stevens et al. (2016) as having shorter strides relative to pes length, and although gaits are unknown, they have suggested that the sauropod manus may be ahead of the pes by two or more steps. While the current use of the TCH dataset cannot assist in the discussion around gait, the data shows that PS and PL are strongly positively correlated (p=<0.001). This suggests that increasing stride lengths are common in trackways with larger PL. While PS and speed (via Alexander, 1976) are, expectedly, significantly positively correlated, speed also shows a positive correlation with PPA and MPA, suggesting sequential increases in pace angulation and the reduced linearity of a trackway. Romano et al. (2007) have proposed greater pace angulation corresponds with smaller trackway gauges (i.e., narrow-gauge types via their gauge ratio, PTR). Trackway gauge has been linked to body size, morphology, and behavior of the trackmaker (Wilson and Carrano, 1999) with TCH trackway gauge approximated by the WAP/PL ratio. TCH data suggest a significant negative relationship between gauge (WAP/PL) and pes angulation (PPA, p=<0.001) but a significant positive relationship with MPA (p=<0.001). Furthermore, speed and PPA are significantly positively correlated (p=<0.001) implying greater speeds result in wider (less acute) PPA despite no relationship between speed and WAP/PL. No relationship between speed and gauge could be determined regardless of Romano et al. (2007) suggesting a relationship between stride, speed, and their gauge ratio (PTR). Castanera et al. (2012) demonstrated that decreasing speeds could correspond to both narrow- and wide-gauged trackways.

4.1.1 Cluster and Principal Component Analysis

Cluster analysis, a partitioning technique, grouped (clustered) TCH trackways with higher similarity to one another (i.e., high intracluster similarity) based on their trackway parameters. Similarity is measured with respect to the mean value of each cluster and recovered groupings use the non-Euclidean K-Means method, considered a robust way to discover the accurate, original groupings in a large dataset (Vavrek, 2016). However, while the data is clusterable the total variance of the clustering is moderate (48.8%) indicating that there is similarity between cluster groups. The HK-Means (Figures 5A,B) grouped TCH data into three clusters and clustering was dominated by trackway pattern specifiers (PPA, MPA, and WAP/PL). Clustering is not stratigraphically defined, i.e., clusters show mixing between different levels, but there is a degree of partitioning that accounts for some grouping based on stratigraphy. For example, levels higher in the stratigraphy are associated.

To explore clustering and trackway pattern variation further, PCA analysis (Figure 6A) was utilized. PC1 is a measure of trackway pattern and reflected trackway gauge (WAP/PL), the distance of manus vs. the pes to trackway midline (WAP/WAM), and manus pace angulation (MPA). PC2 shows the negative relationship between manus supination and heteropody. For instance, the level 1060 trackways grouping is governed by higher trackway gauge ratios (WAP/PL and WAP/WAM: wide gauge), higher manus supination angles, and are opposingly less dictated by heteropody (IPS/IMS) and PPA. Although Level 1060 records the smallest footprint sizes (PL, MW) of the TCH dataset (av. 27 cm PL), and smaller h, its heteropody ratio is the same as the average level 1000 trackways (Table 1), and thus size is not a determining factor. However, it is worth considering if ontogenetic variation is partitioning these trackways. Trackway ontogenetic variability has been previously suggested by Lockley et al. (2002) but behavioral factors may also play a role in trackway pattern (i.e., such a turning during locomotion; see Swiss examples in Meyer, 1993; Marty, 2008; Stevens et al., 2016, fig. 13.7; and globally Meyer et al., 1994; Mezga and Bajraktarević, 1999; Wright, 2005; Ishigaki and Matsumoto, 2009; Xing et al., 2015a; Goodell et al., 2021; Lockley et al., 2021). The high WAP/PL, WAP/WAM and manus supination for smaller pedal and manual sizes has been noted previously at the Swiss Chevenez-Combe Ronde (CRO) track site where behavioral changes are considered a determining factor (Marty, 2008). It is noteworthy to see a repeated portioning of manus supination, WAP/WAM, and WAP/PL for smaller pedal impressions sizes. Figure 6A shows that WAP/PL is negatively correlated with hip height, and this negative relationship is statistically significant (p = <0.001), as discussed above and in Section 3.1.1 (Figure 4).

PPA, h and trackway gauge are negatively correlated along PC1 (Figure 6A), and trackways from levels 1000, 1015, 1030, and 1055 show narrower trackway gauges and opposingly greater h than level 1060, for instance. PPA and h show similarly positively correlated vectors and vector lengths indicating concurrent increases, but this is not reflected by associated increases in the MPA. Thus, there is a disconnect between MPA and PPA, as shown in Figure 3, that can be important in defining a trackway and the role the manus plays in trackway pattern. However, high variability in manus placement/supination, as discussed, is commonly noted and this has a significant bearing on MPA and locomotion (Ishigaki and Matsumoto, 2009; Lallensack et al., 2018 and refs therein). In the PCA-biplot (Figure 6A), the manus supination vector is divergent from that of IPS/IMS (and speed indices when included in the PCA, and to a lesser extent PPA and h) implying their negative correlation. Therefore, larger supination angles are reflected by lower speeds as has been suggested herein.

The incorporation of Swiss and Moroccan trackways from Marty et al. 2010 (Table 2) provides a means to place ichnogeneric classification within the TCH clusters and PCA-biplot (Figure 6B). The addition of these 6 well-defined trackways shows a more apparent subdivision within the data. Tracks assigned to Brontopodus isp. from the Swiss sites Chevenez–Combe Ronde (CRO) and Courtedoux–Pommerat (PMM), and TCH trackway 1052S1:UTw are an outlier cluster (cluster 3) while the others (Clusters 1, 2, and 4) show variation along the PC1 axis which is controlled by trackway pattern and heteropody. Cluster 3 (Figure 6B) separates from clusters 1, 2, and 4 along PC2, and is driven by h, trackway gauge (WAP/PL, WAP/WAM), with lower values in manus supination, PPA, and heteropody. Given that cluster 3’s Brontopodus isp. is defined as a wide-gauge trackway, the placement of these specimens and the TCH trackway, reinforce the driver WAP/PL as appropriate for trackway gauge and highlight the potential role of PPA in wide-gauge trackways. Additionally, Breviparopus isp. from Courtedoux-Sur Combe Ronde (SCR; Switzerland) and Deio (Morocco) plot within the cluster group 2 with Parabrontopodus isp. within clusters 2 and 4 (Figure 6B). These are governed by PPA, IPS/IMS, manus supination, and MPA, and opposing the trackway gauge parameter in WAP/PL, with both Breviparopus specimens tracking to the outermost extent of the grouping as most affected by PPA and heteropody. The latter is a key in the definition of Breviparopus which typically shows pronounced heteropody. In Figure 6B, PC2 is governed by h as a contributing variable to the variation seen along the axis. Hip height and manus supination are negatively correlated. Interestingly, clustering is not significantly affected by h although it contributes to the variation seen along PC2 (Figure 6B).

Breviparopus and Parabrontopodus isp. in the PCA-biplot (Figure 6B), and the overlapping TCH trackways, reflect natural clustering along PC1 governed by trackway pattern and heteropody. These trackways likely reflect a morphological continuum and as Santos et al. (2009) suggested, a Breviparopus/Parabrontopodus plexus or grouping into Parabrontopodus; although more trackway data is needed. Altogether, the cluster (Figure 5) and PCA analyses (Figure 6) show more variability than the GM method alone, as is likely a function of the limitations of GM data (i.e., a limited number of morphologically “well-preserved” individual pes tracks, per level and their specifiers—landmarks subject to placement error, etc.). This likely means that there are similar trackway “shapes” occurring within the Parabrontopodus “ichnomorph” (i.e., clusters, 1, 2, and 4 shown in Figure 6B) at this level of analysis. The observed variation within the groupings is likely a result of preservation, substrate conditions, or individualistic variability that cannot be directly accounted for herein. Finally, the PCAs did not show that the trackways were clearly separated by stratigraphic level and thus no one trackway pattern or level is time-specific at TCH.

4.2 Geometric Morphometric Analyses

Describing size and shape variation in the pes of TCH was done to determine any notable morphological differences between stratigraphic levels irrespective of trackway. GM looked at a pes impression subset of TCH trackways, based on their morphological fidelity, from the various levels. As pes impressions are commonly diagnostic of ichnotaxa but can be confounded by size, preservation, etc., the GM analysis was useful in refining morphological arguments, and to assist in comparing pes impressions, at TCH, through time.

Landmark-based analysis of TCH pes impressions established that size and shape are not correlated. Although testing allometry may be made more robust by incorporating greater variation in sizes, increasing the number of samples within the dataset, or by incorporating “poorly preserved” pes impressions, and concurrently reducing the number of landmarks to include a wider range of track preservation types. Regardless, herein there are no statistically significant size-related morphological trends recorded in the tested TCH stratigraphic levels.

The CVA results, however, indicated that there are some gradational differences in shape, through time (as grouped by level), although these are not statistically significant. Thus, it is likely that the TCH represents one [pes] morphotype across all levels and through time, but which may also reflect slight (statistically insignificant) morphological changes across time. It was assumed that given the disconnect between shape and size, that CVA results would assist in defining and differentiating any grouping in shapes linked to stratigraphic levels. Rather, these results reinforce that pes impression shape is (relatively) conservative through time and complement the traditional morphometric PCA results showing gradational spread in TCH trackways (Figure 6B). This would also strengthen the need to combine pes morphometric analysis with other morphological methods and to be inclusive of trackway data—particularly relating to the manus which can be more diagnostic of specific ichnogenera. Further work will be needed to determine if landmark-based analysis will assist in refining manus morphological differences between ichnogenera. In our current study, manus morphology did not lend well to geometric morphometric work because of their generally simpler shape (i.e., fewer landmarks can be accurately placed) and their, often, poor preservation relative to the pes. Thus, our sample size of manus and pes for comparison would not have been robust. However, this work may be performed on Iberian sauropod manus tracks (Castanera et al., 2016) through time as the record and preservation quality of many of those tracks is superior to those seen at TCH.

The GM-PCA results (Figures 8C,E, 10C,E) demonstrated that the most significant variability of the TCH tracks is a function of the digit impressions with variation in the digit and hypices extension and location, laterally. The least variation is shown in the “heel” landmark between track levels. Substrate variability likely dictates the variation seen in TCH digit (and potentially ungual) impressions. It is hypothesized that sauropod pes employed plantar extension in “softer” substrates (cf. Hall et al., 2016) and would likely result in the PC1 loading shown here (Figures 8C–F), whereas the heel and width of the track remain a constant and is controlled by pes anatomical morphology (i.e., presence of a singular fatty pad). A trial of including the curved landmarks didn’t seem to improve the results—but may help with a larger dataset to explain GM more accurately. The landmark-based description of pes impression morphology is not a precise technique to separate sauropod pes ichnogenera, at least within the TCH site, and the most significant changes were noted with the heel region and the projections of digit I and digit III.

5 Conclusion

The TCH track site was useful in testing quantitative methodological approaches in understanding sauropod tracks and trackways because of the expansive dataset. Generally, quantitatively, TCH track-bearing surfaces (levels) show little significant variation in pes and trackway parameters. The TCH track site appears to reflect the same ichnogenus (Parabrontopodus isp.; Marty et al., 2007, 2010; bar those aligned with Brontopodus) with expected variation in trackway measurements likely a function of substrate, individual behavior, and relating to ichnomorphological spread within Parabrontopodus. Trackway variability can be explained via cluster and PCA analysis more concisely than via the GM of pes tracks alone. However, both methods did outline a spread in TCH trackway and track shape, illustrating a morphological continuum in Parabrontopodus across the levels studied. The incorporation of additional trackways with ichnogeneric assignments (Brontopodus, Breviparopus, and Parabrontopodus ispp.) established overlapping clusters in TCH PC space with both Breviparopus and Parabrontopodus and highlights the ongoing debate around Breviparopus–Parabrontopodus plexus. A larger dataset of known ichnogenera is needed for firmer resolution, in addition to the formal erection (or otherwise) of Breviparopus.

Interestingly, the results illustrated a strong correlation between the manus supination and the forelimb stride length; this fact highlights and reinforces the hypothesis that manus supination is related to behavioral aspects of the trackmaker (e.g., speed, gait, adapting to substrate/topography; see Lallensack et al., 2018). Considering this, manus supination should be more carefully assessed, and may be discarded, as a proxy in identification. Nonetheless, a purely quantitative approach to ichnotaxonomy should also be taken cautiously as some morphological differences can be better described qualitatively, e.g., manus morphology and features (Santos et al., 2009).

The methods used here may allow a more focused driver of qualitative descriptions and data collection of (large) bulk samples because it can define clusters of tracks with morphological similarities within the same ichnotaxon or among different ichnotaxa. Additionally, the data noted that trackway pattern may have potential in assessing ontogenetic variability, but this needs to be constrained, expanded, and described with data relating to small pedal/manual sizes. However, an explanation of this variability may be speculated as relating to the proportions of the limbs to body length in smaller/juvenile sauropods relative to larger/adult sauropods. Finally, while the trackways and tracks are relatively conservative, our results cannot rule out that several sauropod genera generated the trackways and inhabited the area throughout this time span.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

LS and MB are supported by the Swiss National Science Foundation (SNSF Grant 200021_192036).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor is currently co-organizing a Research Topic with one of the authors MB, and confirms the absence of any other collaboration.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Palaeontological team PALA16 (Paléontologie A16) for careful and detailed excavation, documentation, curation, and cataloging of the Courtedoux-Tchâfouè track site and trackway materials. JURASSICA Museum scientific team for assistance, discussion, and access to resources. Beneficial comments, revisions, and the time of the two reviewers, AK and ML, who assisted in improving this manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.805442/full#supplementary-material

References

Abrahams, M., Bordy, E. M., and Knoll, F. (2021). Hidden for One Hundred Years: A Diverse Theropod ichnoassemblage and Cross-Sectional Tracks From the Historic Early Jurassic Tsikoane ichnosite (Clarens Formation, Northern Lesotho, Southern Africa). Hist. Biol. 33 (10), 2504–2519. doi:10.1080/08912963.2020.1810681