Introduction
Characterizing morphologic variation within fossil groups is fundamental for determining taxonomic boundaries and interpreting evolution of extinct taxa (George Reference George1956; Gingerich Reference Gingerich1985), yet challenges remain in capturing the full range of shape variability within species (Allmon Reference Allmon2013). Paleontological species cannot rely on the biologic species concept (Mayr Reference Mayr1957, Reference Mayr1996) and thus are determined using the morphological species concept (George Reference George1956; Gingerich Reference Gingerich1985; Allmon Reference Allmon2013). Commonly, species are classified using qualitative and ultimately subjective morphologic metrics that may overemphasize the significance of shape difference or fail to capture the full range of form variation within a given species’ population (Vogt et al. Reference Vogt, Bartolomaeus and Giribet2009; Allmon Reference Allmon2013). Temporal and spatial shape variability within a fossil group further complicates characterizing shape variation (Allmon Reference Allmon2013). Finding objective methods to record the spectrum of morphologic variation within extinct taxa is particularly valuable for biostratigraphically useful groups such as conodonts, because subjective morphologic criteria commonly lead to disagreements that affect stratigraphic correlations (Girard and Renaud Reference Girard and Renaud2011).
Since their discovery 160 years ago by Pander (Reference Pander1865), conodonts have become highly important biostratigraphic markers for Paleozoic strata (Donoghue et al. Reference Donoghue, Purnell and Aldridge1998; Sweet and Donoghue Reference Sweet and Donoghue2001), but debates persist surrounding conodont taxonomy (Aldridge and Purnell Reference Aldridge and Purnell1996; Barrick et al. Reference Barrick, Lambert, Heckel and Boardman2004; Brown et al. Reference Brown, Rexroad and Zimmerman2016). Certain conodonts exhibited cosmopolitan distributions and evolving morphologies that provide an excellent record of morphologic change throughout the Paleozoic and Triassic (Sweet Reference Sweet1988; Sweet and Donoghue Reference Sweet and Donoghue2001) and an evolutionarily significant history of early chordate evolution (Purnell Reference Purnell1995; Aldridge and Purnell Reference Aldridge and Purnell1996; Donoghue et al. Reference Donoghue, Purnell and Aldridge1998; Donoghue et al. Reference Donoghue, Sansom and Downs2006). Nevertheless, disagreements regarding morphologic variation within and among species persist and fuel continual debate regarding evolutionary interpretations and biostratigraphic utility of certain conodont taxa (Barrick et al. Reference Barrick, Boardman and Heckel1996, Reference Barrick, Lambert, Heckel and Boardman2004; Rexroad et al. Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown et al. Reference Brown, Rexroad and Zimmerman2016).
To help resolve these debates, conodont workers have employed quantitative methods to characterize certain conodont morphologies. Many researchers used Fourier analysis of outlines to test species groupings (Klapper and Foster Reference Klapper and Foster1986, Reference Klapper and Foster1993; Sloan Reference Sloan2003; Girard et al. Reference Girard, Renaud and Serayet2004b, Reference Girard, Renaud and Feist2007). These studies validated the use of outline methods for characterizing conodont morphology and demonstrated that species boundaries are commonly gradational. Others demonstrated that patterns of morphologic change can be characterized quantitatively through time with outline methods (Murphy and Cebecioglu Reference Murphy and Cebecioglu1987; Renaud and Girard Reference Renaud and Girard1999; Girard et al. Reference Girard, Renaud and Korn2004a; Roopnarine et al. Reference Roopnarine, Murphy and Buening2004; Jones and Purnell Reference Jones and Purnell2007). These results provide a useful foundation for the validity of conodont morphometric work, but all studies lacked a direct evaluation of how morphometric results would affect evolutionary interpretations and biozonations. We concur that assessing the influence of morphometric results on previous biozonations is highly valuable for such a biostratigraphically useful fossil group.
Here, we build on previous morphometric conodont work by conducting the first landmark-based geometric morphometric (GM) analysis (Bookstein Reference Bookstein1991; Zelditch et al. Reference Zelditch, Swiderski, Sheets and Fink2004) of the biostratigraphically useful conodont genus Neognathodus and by evaluating how our results compare with evolutionary interpretations and previous biozonations. In addition, we conduct the first statistical assessment of Neognathodus evolutionary models in the Desmoinesian using methods developed by Hunt (Reference Hunt2006) and used by many other researchers (Novack-Gottshall Reference Novack-Gottshall2008; Jones Reference Jones2009; Piras et al. Reference Piras, Sansalone, Marcolini, Tuveri, Arca and Kotsakis2012; Van Bocxlaer and Hunt Reference Van Bocxlaer and Hunt2013). GM analysis allows us to both quantify the continuous shape variation of P1 elements of Neognathodus morphotypes without placing specimens into different groups a priori and test whether established morphotype groups are reliably distinct from one another. Hunt’s (Reference Hunt2006) approach also allows us to minimize a priori assumptions, as the method does not assume a null hypothesis for any given evolutionary mode and does not favor any one evolutionary model a priori over others. Instead, it allows us to test all three modes simultaneously and compare each with equal weight based on the numeric data (Hunt Reference Hunt2006).
Evolution of Neognathodus species through time is particularly well documented, as the species are used to correlate Desmoinesian (Middle Pennsylvanian) strata in North America, most commonly in the Illinois Basin (Fig. 1); (Merrill Reference Merrill1972, 1975a,Reference Merrillb; Grayson et al. Reference Grayson, Trice and Westergaard1985; Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2013, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009). Neognathodus species are defined as populations of morphotypes, and morphotypes are classified by the morphology of the P1 (platform) element (Merrill Reference Merrill1972, 1975a,Reference Merrillb, Reference Merrill1999). The species name corresponds to the predominant morphotype in the population (Merrill Reference Merrill1972, Reference Merrill1999). Interpretations of Neognathodus evolution and biozonation were summarized using the Neognathodus Index (NI), a weighted average of the morphotype population (i.e., the species) in a given stratigraphic unit (Brown et al. Reference Brown, Rexroad and Zimmerman2016). The NI relies on qualitative and subjective groupings of morphotypes (Brown and Rexroad Reference Brown and Rexroad2009). Furthermore, Merrill (Reference Merrill1972, 1975a,Reference Merrillb, 1999) documented a gradational transition in shape between morphotype groups. Thus, classifying specimens into separate morphotype groups may not effectively capture the range of transitional shape variation.
Figure 1 Summary of previous Neognathodus work from the Desmoinesian of the Illinois Basin. A, Extent of the Illinois Basin during the Desmoinesian (Tri-State Committee 2001). B, Chart of European and North American stages in the Pennsylvanian; adapted from Cohen et al. (Reference Cohen, Finney, Gibbard and Fan2013). C, Stratigraphic nomenclature of the Desmoinesian (Middle–Late Pennsylvanian) in Indiana. Stratigraphic units examined in this study are highlighted in gray. Stratigraphy is adapted from Shaver et al. (Reference Shaver, Ault, Burger, Carr, Droste, Eggert, Gray, Harper, Hasenmueller, Hasenmueller, Horowitz, Hutchison, Keith, Keller, Patton, Rexroad and Weir1986), Tri-State Committee (2001), and Swezey (Reference Swezey2009). Lithologic unit thicknesses are not to scale. Fm, formation; Gr., group; LS, limestone; and SH, shale. D, Histograms showing the distribution of Neognathodus morphotypes for each stratigraphic interval. Height of histogram bars are in percent, and all bars within a histogram total 100%. Morphotypes are numbered from oldest to youngest along the x-axis (1=N. bassleri, 2=N. bothrops, 3=N. medadultimus, 4=medexultimus, 5=N. roundyi, 6=N. dilatus, and 7=N. metanodosus) The study that published each distribution is shown to the right in part V. E, Published Neognathodus Index (NI) value, named biozone for each interval, and references to published works for each unit. The NI is the weighted average of each morphotype distribution.
Our study addresses three questions and their subsequent hypotheses. (1) Are the GM groups of Neognathodus congruent or incongruent with published morphotype groups of N. bassleri (Harris and Hollingsworth Reference Harris and Hollingsworth1933); N. bothrops Merrill, Reference Merrill1972; N. medadultimus Merrill, Reference Merrill1972; N. medexultimus Merrill, Reference Merrill1972; N. roundyi (Gunnell Reference Gunnell1931); N. dilatus (Stauffer and Plummer Reference Stauffer and Plummer1932); and N. metanodosus Merrill, Reference Merrill1975b? (2) How do statistical GM interpretations of morphologic change through the Desmoinesian compare with previous evolutionary interpretations? (3) Do GM results significantly affect previous biozonations of Illinois Basin units during the Desmoinesian? GM results will allow us to capture the full spectrum of shape variability within and among morphotype groups, provide a quantitative record of morphologic change through the extinction of the genus at the end of the Desmoinesian (Merrill Reference Merrill1975a; Merrill and Grayson Reference Merrill and Grayson1989; Boardman et al. Reference Boardman, Heckel, Barrick, Nestall and Peppers1990), and supply the first quantitative evaluation of how GM-based biozonation compares with previous biozonation patterns. Due to the gradual transition in shape between morphotype groups (Merrill Reference Merrill1972, 1975a,Reference Merrillb, Reference Merrill1999), we hypothesize GM groups will not be entirely congruent with previous morphotype groups. Previous evolutionary interpretations and biozones are based on the population-based NI; thus, we hypothesize that GM morphologic change through the Desmoinesian will support previous interpretations of evolution, and GM biozones will be similar in structure to published biozonations.
Certain considerations must be addressed regarding taxonomic boundaries and temporal shape variation. Our methods assume that an outline-based morphometric analysis will effectively characterize shape variability (Klapper and Foster Reference Klapper and Foster1993; Girard et al. Reference Girard, Renaud and Korn2004a, Reference Girard, Renaud and Feist2007) within and among morphotype groups and between Neognathodus populations through time. We also judge that morphologic changes of the Neognathodus P1 element are sufficiently documented through the Desmoinesian in the Illinois Basin (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009) to serve as a valid and established biostratigraphic proxy for comparing and statistically correlating with new GM patterns of Neognathodus morphologic change through time.
Quantitative characterizations of Neognathodus morphologies could offer valuable numeric evidence to document how conodont evolution was influenced by environmental change, especially during the global Desmoinesian glacial–interglacial iterations (Heckel Reference Heckel1991; Falcon-Lang et al. Reference Falcon-Lang, Heckel, DiMichele, Blake, Easterday, Eble, Elrick, Gastaldo, Greb, Martino, Nelson, Pfefferkorn, Phillips and Rosscoe2011). The Desmoinesian was dominated by an icehouse regime, as continental glaciations repeatedly occurred in southern Gondwana (Veevers and Powell Reference Veevers and Powell1987; Joachimski et al. Reference Joachimski, von Bitter and Buggisch2006). These glacial repetitions altered climate conditions and sea levels in the tropical regions including the Illinois Basin (Poulsen et al. Reference Poulsen, Pollard, Montanez and Rowley2007; Tabor and Poulsen Reference Tabor and Poulsen2008). During glacial intervals, the tropical climate was drier and more seasonal, whereas during interglacial intervals, the tropical climate was more humid to subhumid (Falcon-Lang and DiMichele Reference Falcon-Lang and DiMichele2010). An intense glacial phase and regression followed by a warming and transgression likely occurred at the Desmoinesian–Missourian boundary (Falcon-Lang et al. Reference Falcon-Lang, Heckel, DiMichele, Blake, Easterday, Eble, Elrick, Gastaldo, Greb, Martino, Nelson, Pfefferkorn, Phillips and Rosscoe2011), and the extinction of multiple conodont groups at that boundary, including Neognathodus (Merrill Reference Merrill1975a; Merrill and Grayson Reference Merrill and Grayson1989; Boardman et al. Reference Boardman, Heckel, Barrick, Nestall and Peppers1990), could have resulted from these major climatic changes (Falcon-Lang et al. Reference Falcon-Lang, Heckel, DiMichele, Blake, Easterday, Eble, Elrick, Gastaldo, Greb, Martino, Nelson, Pfefferkorn, Phillips and Rosscoe2011).
Materials and Methods
We photographed 390 specimens of Neognathodus P1 elements for GM analysis. Our project examined Neognathodus specimens published in six prior studies (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009). Specimens were spatially sampled from locations in eastern Illinois, northwestern Kentucky, and southwestern Indiana, and stratigraphically sampled from the following Desmoinesian units: Perth and Holland Limestone Members of the Staunton Formation, Mecca Quarry Shale and Velpen Limestone Members of the Linton Formation, Alum Cave and Providence Limestone Members of the Dugger Formation, and the West Franklin Limestone Member of the Shelburn Formation (Fig. 2). The Mecca Quarry Shale and Velpen Limestone Members were combined and sampled as one interval, because previous conodont sampling did not differentiate between the two members (Rexroad et al. Reference Rexroad, Wade, Merrill, Brown and Pagett2001). For all units, the depositional environment is interpreted as a nearshore, deltaically influenced marine setting (Ferm et al. Reference Ferm, Horne, Swinchatt and Whaley1971; Trask and Palmer Reference Trask and Palmer1986; Rexroad et al. Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown et al. Reference Brown, Rexroad and Zimmerman2016). Thus, any Neognathodus sample variation was likely not affected by changing paleoenvironments.
Figure 2 Map showing sampling localities and stratigraphic units of the 390 photographed specimens used in geometric morphometric (GM) analysis. Unit names are arranged in stratigraphic order, and dashed gray lines separate members of different stratigraphic groups. Detailed stratigraphic information is shown in Fig. 1. Outcrop data are sourced from Gray et al. (Reference Gray, Ault, Keller and Harper2010) for Indiana and from Noger et al. (Reference Noger, McDowell, Grabowski and Moore1988) for Kentucky. Detailed, statewide outcrop maps were not available for Illinois. Generalized paleogeographic reconstructions of deltas and rivers during the Desmoinesian are superimposed and adapted from Jacobson (Reference Jacobson2000). LS, limestone; Mbr, member; SH, shale.
All specimens are housed in the Indiana University Paleontology Collection in the Department of Earth and Atmospheric Sciences at Indiana University, Bloomington, Indiana (catalog numbers in Supplementary Table 1). We selected specimens with a completely preserved and unbroken P1 element for photography. We photographed Neognathodus P1 elements with a digital microscope at 20× magnification at 1024 × 1660-pixel resolution. Neognathodus morphotypes are identified from an oral view of the P1 element, and we photographed and used the same view for GM analysis.
The 390 specimens comprised two different subsets, because we conducted two separate analyses. One subset was used to analyze morphotype group affinities, and we sampled specimens across stratigraphic boundaries to acquire 30 specimens per morphotype group. The other specimen subset was used to analyze Neognathodus morphologic change through the Desmoinesian, and we sampled 30 specimens per stratigraphic interval. Using two different specimen subsets allowed us to form a consistent sample size for each separate analysis, thereby avoiding the influence of sample size differences on results.
■
Sample Selection
To test morphotype group affinities among N. bassleri, N. bothrops, N. medadultimus, N. medexultimus, N. roundyi, N. dilatus, and N. metanodosus we compiled data tables from previous studies (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009) and then counted the number of times each morphotype group occurred in every stratigraphic interval. For example, the N. bassleri morphotype group has 110 occurrences in the Perth Limestone Member, six occurrences in the Mecca Quarry Shale and Velpen Limestone Members, and four occurrences in the Alum Cave Limestone Member (Supplementary Table 2). We then totaled the number of occurrences for each morphotype group and used that sum to convert the raw occurrence data to percent occurrence data. For example, there were 120 N. bassleri occurrences across all units, which demonstrates that 92% of all N. bassleri specimens occur in the Perth Limestone Member (110/120=0.92), 1.5% occur in the Mecca Quarry and Velpen Limestone Members (6/120=0.015), and 1.0% occur in the Alum Cave Limestone Member (4/120=0.010). To replicate this stratigraphic morphotype distribution in our GM analysis, we selected 30 specimens from each of the seven morphotypes (total of 210 specimens), and then multiplied these percentages by 30 to determine the number of specimens to select and photograph from each stratigraphic interval (Table 1). For example, to acquire our 30-specimen sample for N. bassleri, we selected 27 N. bassleri specimens from the Perth Limestone (92% × 30 ≈ 27), two from the Mecca Quarry Shale and Velpen Limestone (1.5% × 30 ≈ 2), and one from the Alum Cave Limestone (1.0% × 30 ≈ 1). We chose a sample size of 30 for each subset, because it provides statistically useful and significant results without requiring an excessively high amount of specimen and data processing.
Table 1 Table showing the stratigraphic distribution of Neognathodus specimens selected to test morphotype group affinities. Specimens were chosen between multiple stratigraphic boundaries. For example, to select the 30-specimen sample of N. bassleri, 27 N. bassleri specimens were chosen from the Perth Limestone Member, two were chosen from the Mecca Quarry Shale and Velpen Limestone Members, and one was chosen from the Alum Cave Limestone Member. LS, limestone, Mbr, member; SH, shale. The number of specimens selected from each unit is proportional to the total number of morphotype group occurrences in that unit.
To characterize Neognathodus morphologic change through the Desmoinesian in the Illinois Basin, we selected a representative, random sample of 30 Neognathodus specimens from each stratigraphic interval (total of 180 specimens). We used published occurrence data to select randomly, with replacement, 30 Neognathodus specimens from each stratigraphic interval (Supplementary Table 2) (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009). For example, one sampling run of the Perth Limestone Member selected 15 N. bassleri and 15 N. bothrops specimens. We repeated random sampling 5000 times and then recorded the number of specimens in each morphotype group for each run. After randomly choosing 30 specimens 5000 times, we averaged the results of all runs to produce a summarized mean list of 30 specimens from each of the six stratigraphic intervals (30 × 6=180 total specimens). This final mean list provided a representative random sample of 30 Neognathodus specimens to select and photograph from each stratigraphic interval (Table 2). For example, the final 30-specimen sample for the Perth Limestone Member consisted of 10 N. bassleri specimens, 14 N. bothrops specimens, 5 N. medadultimus specimens, and 1 N. dilatus specimen. We again chose a sample size of 30 for consistency with sample sizes in prior analyses.
Table 2 Table showing the morphotype group counts of the specimens selected to analyze Neognathodus morphologic change through the Desmoinesian. Specimens were chosen from multiple morphotype groups within each stratigraphic interval. This selection of 30 specimens from each stratigraphic interval represents a random sample of the population of Neognathodus within that stratigraphic unit. For example, the 30-specimen random sample from the Perth Limestone Member consisted of 10 N. bassleri specimens, 14 N. bothrops specimens, and 1 N. dilatus specimen. LS, limestone, Mbr, member; SH, shale. The number of specimens to select from each morphotype group was determined using the average of 5000 random resampling runs.
Geometric Morphometric Data and Processing
Our project followed standard landmark-based GM techniques that use landmarks (x, y points on an image) to quantitatively characterize shape (Bookstein Reference Bookstein1991; Zelditch et al. Reference Zelditch, Swiderski, Sheets and Fink2004). The GM method requires that an equal number of landmarks be placed in the same order for each specimen’s image. To meet this requirement, we rotated and vertically reflected photographs to a consistent orientation, with the outer parapet facing upward and the posterior tip of the P1 element facing to the right (e.g., see Fig. 3A). We removed the background of each image and a digital outline of 241 landmarks was placed to characterize the two-dimensional shape of each Neognathodus P1 element (Fig. 3A). Outline points were automatically placed at equal intervals using the freeware program tpsDig (Rohlf Reference Rohlf2016a), which greatly reduced the likelihood of subjective landmark variability. We configured the coordinate data as a TPS file with tpsUtil (Rohlf Reference Rohlf2016b) and then imported the TPS landmark file into Mathematica (Wolfram Reference Wolfram2016) software using the freeware extension Geometric Morphometrics for Mathematica (Polly Reference Polly2016).
Figure 3 Diagram showing the methodology and results of morphotype group affinities. A, Photograph of a N. bassleri P1 element with a fully placed outline and annotated Neognathodus morphology. B, Thirty superimposed outlines of N. bassleri. Stars show similar locations in panels A and B. C, Scatter plot of PC 1 and PC 2 values for each specimen used to analyze morphotype group affinities. Every specimen is represented by a number corresponding to its morphotype group. The black dots correspond to the mean PC 1 and PC 2 coordinates of each morphotype group. Photographs of the P1 element for each morphotype group are also shown. PC signs are arbitrary, and the origin coordinates (0,0) represent the mean of PC 1 and PC 2, respectively. There is a gradual progression from the oldest morphotype (N. bassleri=1) on the right to the youngest morphotype (N. metanodosus=7) on the left.
Our configuration of 241 landmarks captured detailed P1 element shape variation while still using minimal computational power. The outline traced both the outer and inner parapets of the P1 element, from the tip of the platform to the anterior section where both parapets meet with the carina. We placed one additional landmark at the posterior end of the carina to characterize carinal length (Fig. 3A). We did not outline the anterior edge of the carina, because it is usually broken and the morphology is indistinct and not significant among morphotypes (Merrill Reference Merrill1972, Reference Merrill1975b). We placed only one landmark on the posterior end of the carina, because (1) adding more points did not significantly affect later statistical results and (2) using only one point provided clear shape comparison of the P1 element outlines in subsequent visual analysis. We determined that 241 landmarks captured P1 shape variation better than an outline with a lower number of landmarks, because P1 shape changes were not as well visualized with fewer landmarks and statistical results were not stable below about 200 landmarks. Increasing the number of landmarks did not statistically change results and resulted in more computational resources and time for subsequent GM calculations.
Per common GM practices (Bookstein Reference Bookstein1991; Zelditch et al. Reference Zelditch, Swiderski, Sheets and Fink2004), we superimposed the landmark coordinates from the images using Procrustes analysis and then calculated principal component (PC) scores. Procrustes analysis, described in detail by Rohlf and Slice (Reference Rohlf and Slice1990), first calculates the centroid of all landmarks, then centers that point on the origin of a Cartesian coordinate frame, scales the shape to a standard size, and finally rotates the shape to minimize distance variability (Fig. 3A). Following Procrustes superimposition, we performed principal components analysis (PCA) for each specimen outline to calculate PC scores. The GM data set is multidimensional; thus, PC scores act as convenient and consolidated shape variables to use in subsequent calculations and statistical tests. PC 1 shows the largest percent of shape variance, while each subsequent PC represents a progressively smaller amount of variance.
We created a scatter plot of PC 1 and PC 2 scores for both sets of specimens to visualize shape change in a PC morphospace. For the morphotype specimen set, we calculated the mean shape of each section in the PC scatter plot and displayed it as Procrustes coordinates in the form of thin plate splines (TPS). The TPS morphospace diagram shows average and extreme morphologies and illustrates which morphologic changes are described by the first two PC axes. The TPS calculations follow the methods of Hammer and Harper (Reference Hammer and Harper2008), and we performed the mathematics using the function ‘PrincipalComponentsOfShape’ in Geometric Morphometrics for Mathematica (Polly Reference Polly2016).
Quantitative Evaluations
To evaluate shape similarity among both morphotype groups and stratigraphic units, we performed a multivariate analysis of variance (MANOVA) hypothesis test examining mean PC differences using the function ‘ShapeMANOVA’ in Geometric Morphometrics for Mathematica (Polly Reference Polly2016), which performs a nonparametric randomization test for differences in group means. The null hypothesis of this test states there is no difference among the mean PC scores for each individual group of specimens, while the alternate hypothesis states there is a statistically significant difference among the mean PC scores for each individual group. This test is two tailed, uses all of the nonzero PC scores for each group, and compares each group to all other groups. We also performed nonparametric randomization to evaluate the significance of the difference in means among individual groups. This process randomized the groups and recalculated the difference 5000 times to assess significance. To account for the likely variability in Neognathodus P1 element shapes, we chose a significance level (alpha) of 0.10. This alpha also helps minimize type II statistical error (i.e., false-negative). Thus, a p-value greater than 0.10 would allow the null hypothesis, or no difference among mean shape, to be confidently rejected.
To evaluate how GM results might change morphotype diversity, we compared the number of GM groups with the number of morphotype groups at each sampling location. GM analysis coalesced morphotype groups, so the total number of distinct GM groups was less than the total number of morphotype groups. Therefore, we subtracted the number of GM groups from the number of morphotype groups at all previously published sampling localities of Brown et al. (Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016), Rexroad et al. (Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001), and Brown and Rexroad (Reference Brown and Rexroad2009). To summarize these results, we calculated the mean and median difference for each stratigraphic interval.
To assess morphologic changes through the Desmoinesian, we visually compared mean PC 1 values for each stratigraphic interval with published trends of the NI (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009). We calculated the mean PC 1 value for each stratigraphic interval and paired the mean PC 1 value with a stratigraphic depth to form a generalized time series showing average PC 1–based shape change through the Desmoinesian. We plotted the PC 1 time series with published NI trends of Brown et al. (Reference Brown, Rexroad and Zimmerman2016) to visualize and compare the two patterns of changing Neognathodus morphology.
To quantitatively characterize Neognathodus P1 element shape change through the Desmoinesian, we conducted a maximum-likelihood evaluation of three evolutionary models following methods developed by Hunt (Reference Hunt2006). The three models tested were an unbiased random walk (URW) with no directional component, a general random walk (GRW) that has a directional component, and stasis. This method fits the mean and variance of observed distribution of first differenced PC 1 values to their expectation under each of the three models (McKinney Reference McKinney1990). Following Hunt (Reference Hunt2006), small-sample unbiased Akaike information criterion (Akaike Reference Akaike1974) and Akaike weights (Anderson et al. Reference Anderson, Burnham and Thompson2000) were used to select the model that best fits the observed PC 1 data. In terms of the evolutionary models, the mean determines the directionality, or tendency to change, and the variance determines the volatility of changes around the mean. The procedure requires an estimated age value, trait mean, and trait sampling error (i.e., variance) for each interval. We used age estimates for each stratigraphic interval from Tri-State Committee (2001) and Swezey (Reference Swezey2009), calculated the mean and variance of PC 1 for each stratigraphic interval (i.e., no pooled variance), and performed the analysis using the function ‘ThreeModelTest’ in the Phylogenetics for Mathematica package (Polly Reference Polly2018).
To quantitatively evaluate the relationship of GM morphologic patterns with the previous methods of interpreting evolution using NI values, we conducted a multivariate least-squares regression of shape using the function ‘ShapeRegress’ in Geometric Morphometrics for Mathematica (Polly Reference Polly2016), which is a fully multivariate nonparametric randomization test for shape variables. The predictor value was the NI, and the response variable was each nonzero PC value of the mean shape for each stratigraphic interval. We determined the mean shape of all 30 specimens in each interval. All six mean outlines were Procrustes superimposed. We then calculated the PC scores of the averaged Procrustes-superimposed shape and regressed the PC scores onto the input, or NI, variable. To determine significance, we calculated the R-squared (R 2) values for the regression and the scores randomized with respect to the input variable, the NI, 10,000 times.
To assess how GM results affected previous biozonations, we determined GM-based biozones. First, we used the compiled data tables from prior studies (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009) to determine the percent abundance of every distinct GM group in each stratigraphic interval. If a GM group consisted of two previously separate morphotype groups, we combined the counts of both groups to find the total percent abundance for the GM group. We then assigned the GM biozone to the predominate GM group in each stratigraphic unit.
Results
■
Morphotype Group Affinities
Calculations of all PC scores for the 210 Neognathodus specimens selected for testing morphotype group affinities show that PC 1 accounts for 45.7% of total shape variance, PC 2 accounts for 18.9%, and subsequent PCs progressively decrease in percent variance. MANOVA tests used all PCs and compared each morphotype with all other morphotypes, but for visual clarity, only PC 1 is shown in figures.
Plotting PC 1 and PC 2 values for each of the 210 Neognathodus specimens and the PC 1 and PC 2 mean of each morphotype group illustrates there is considerable shape overlap and similarity among morphotypes (Fig. 3C). Note that the signs of the PC axes are arbitrary and the x-axis and y-axis represent the mean of PC 1 and PC 2, respectively. The most noticeable point pattern follows PC 1, as the clusters for each morphotype group gradually progress from the oldest morphotype, N. bassleri (group 1), on the far right to the youngest morphotype, N. metanodosus (group 7), on the far left. Also, the morphotype pairs N. medadultimus/N. medexultimus (groups 3 and 4), and N. dilatus/N. metanodosus (groups 6 and 7) plot extremely close together, indicating each pair has highly similar average shapes.
The TPS morphospace diagram visually shows that PC 1 seems to capture the overall width of the P1 element, while PC 2 appears to capture the degree of parapet curvature (Fig. 4). First, it is important to note that PC 1 and PC 2 are statistically uncorrelated, thus the shape change across the diagram should not be interpreted diagonally. The vertical distance between P1 element landmarks increases from left to right along PC 1 for each individual row. The left TPS diagrams in the middle and bottom rows show the narrowest P1 elements. This region of the morphospace is primarily occupied by N. dilatus and N. metanodosus, the morphotypes with the narrowest P1 elements (Fig. 3C). The upper two diagrams in the middle and right columns are primarily composed of N. roundyi specimens that exhibit a similar distinct flare of the outer parapet (Fig. 4).
Figure 4 Morphospace diagram illustrating the average shape of Neognathodus P1 elements across PC 1 and PC 2 as thin plate splines. PC 1 and PC 2 are statistically uncorrelated, and shape change across the diagram should not be interpreted diagonally. This morphospace diagram spans a wider range of shape space then observed in Fig. 3 to better visualize and identify shape variation across PC 1 and PC 2. Note P1 element width increases from left to right along each row, and curvature of the parapets decreases from top to bottom along each column.
MANOVA tests of all PC means among Neognathodus morphotype groups provide p-values (Fig. 5) that support distinct groups among certain morphotypes and indistinct groups between others. Statistically distinct groups (i.e., p-values<0.10 when compared with all other morphotype groups) include N. bassleri, N. bothrops, and N. roundyi. These statistically significant results, combined with the visually separate placement on morphospace plots, provide evidence that these groups are distinct in shape. Statistically indistinct groups (p-values>0.10) include the morphotype pairs N. medadultimus/N. medexultimus and N. dilatus/N. metanodosus. The high p-values, in addition to the extremely close placement on morphospace plots, indicate these morphotype pairs are not distinct, and there is no meaningful, quantitative shape difference captured between N. medadultimus and N. medexultimus, and between N. dilatus and N. metanodosus.
Figure 5 Diagram of PC 1 distributions for the 30-specimen sample of each analyzed Neognathodus morphotype group. Photographs of the P1 element for each morphotype are also shown. The p-values from MANOVA tests of all PC means are plotted on the x-axis. Dark gray intervals indicate morphotypes with statistically similar mean shapes (p-value>0.10). MANOVA tests used all PCs and compared each morphotype with all other morphotypes. For visual simplicity, only PC 1 distributions of large p-values (p-value>0.10) and of certain low p-values (p-value<0.10) are shown. This analysis shows the morphotype pairs of N. medadultimus and N. medexultimus and N. dilatus and N. metanodosus do not exhibit a statistically significant difference in average P1 element shape.
Coalescing these statistically similar morphotype pairs decreases morphotype diversity in many of the previously published Desmoinesian sampling localities of Brown et al. (Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016), Rexroad et al. (Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001), and Brown and Rexroad (Reference Brown and Rexroad2009) across Indiana, Kentucky, and Illinois (Table 3). The only member to show no diversity loss (mean difference of 0) is the Perth Limestone Member. This member only contained N. medadultimus and older morphotypes, which remained as distinct groups in the GM analysis. All other members show a diversity loss up to one or two morphotype groups at certain sampling localities.
Table 3 Table summarizing, for each stratigraphic unit, how geometric morphometric (GM) group diversity results compare with previous diversity results based on morphotype group counts. For each stratigraphic interval, the median and mean value for the number of morphotype groups and the number of GM groups is shown, as is the median and mean difference between the two. Localities for each stratigraphic interval are plotted in Fig. 2. The stratigraphic units least affected by diversity loss are the Holland and Perth Limestone Members.
Morphologic Change and Biozonation through the Desmoinesian
Calculations of all PC scores for the 180 Neognathodus specimens selected for characterizing morphologic change through the Desmoinesian demonstrate that PC 1 accounts for 38.0% of total shape variance, PC 2 accounts for 24.9%, and subsequent PCs progressively decrease in percent variance. MANOVA tests used all PCs and compared each morphotype with all other morphotypes, but for visual clarity, only PC 1 is shown in figures.
Plotting PC 1 and PC 2 values for each of the 180 Neognathodus specimens illustrates there is considerable shape similarity between stratigraphic interval samples (Fig. 6). The most noticeable pattern follows PC 1 as the point clusters for each stratigraphic interval quickly progress from the oldest unit, Perth Limestone Member, on the right, to youngest unit, West Franklin Limestone Member, on the left. The means for the Mecca Quarry Shale/Velpen Limestone Members and the Alum Cave and Providence Limestone Members plot relatively close together, indicating all three units have similar average shapes.
Figure 6 Scatter plot of PC 1 and PC 2 values for each specimen used to analyze Neognathodus morphologic change through time. Every specimen is represented by a letter corresponding to its stratigraphic interval. The black dots correspond to the mean PC 1 and PC 2 coordinates of the average shape of each stratigraphic interval sample. PC signs are arbitrary, and the origin coordinates (0,0) represent the mean of PC 1 and PC 2, respectively. There is a gradual progression from the oldest stratigraphic interval sample (Perth Limestone Member=A) on the right to the youngest sample (West Franklin Limestone Member=F) on the left. LS, limestone; Mbr, member; SH, shale.
MANOVA tests of all PC means between Neognathodus samples from each stratigraphic interval show p-values (Fig. 7C) that exhibit periods of both statistically distinct and similar mean shape through the Desmoinesian. There are statistically significant p-values (p-value<0.10) from the Perth Limestone Member to the Holland Limestone Member and to the Mecca Quarry Shale/Velpen Limestone Members. These statistically significant results, combined with the visually separate placement on the PC plot, provide evidence that the mean shape for each successive interval is distinct. In contrast, statistically insignificant p-values (p-value>0.10) are shown from the Mecca Quarry Shale/Velpen Limestone Members to the Alum Cave and Providence Limestone Members. The high p-values, as well as the extremely close placement on the PC plot, indicate that the mean shapes in these intervals are not distinct. Statistically significant p-values and separate PC means are again observed from the Providence Limestone Member to the West Franklin Limestone Member, indicating mean shape is distinct between the two.
Figure 7 Diagram comparing geometric morphometric (GM) and published patterns of Neognathodus morphologic change through the Desmoinesian in the Illinois Basin. A, Summary chart showing the relative stratigraphic placement of each analyzed interval. B, Diagram of PC 1 distributions for the 30-specimen sample of each analyzed stratigraphic interval through the Desmoinesian. C, The p-values from MANOVA tests of all PC means are plotted along the y-axis. Dark gray intervals indicate stratigraphic intervals with statistically similar mean shapes (p-value>0.10) and are interpreted as morphologic stability. MANOVA tests used all PCs and compared each stratigraphic interval with all intervals. For visual simplicity, only PC 1 distributions of large p-values (p-value>0.10) and of certain low p-values (p-value<0.10) are shown. D, Graph comparing the new GM interpretations to the published patterns of Neognathodus shape change, which are represented by Neognathodus Index (NI) values (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009). The dashed line labeled “GM” uses the bottom scale of “PC1,” while the solid line labeled “Literature” uses the top scale of “NI (Literature) Values.” Despite minor differences in magnitude, the overall trends (i.e., morphologic change and stability) are relatively similar. Gr., group; LS, limestone; SH, shale.
A visual comparison of PC 1 averages and NI values for each stratigraphic interval shows a first-order similarity between both time series, and although certain differences are noticeable (Fig. 7D), additional statistical tests provide a strong correlation between the two (Fig. 8). The most apparent difference occurs where the NI value shows a large change from the Perth to the Holland Limestone Members, but the PC 1 mean exhibits relatively less change. Multi- and univariate regressions show primarily strong correlation strengths between PC means and NI values for each stratigraphic interval (Fig. 8). The multivariate regression of all nonzero PC means onto the NI resulted in an R 2 value of 0.67, meaning that 67% of the total shape variance in the PC means is predictable from the correlation to the NI. In contrast, the univariate regression of PC 1 means onto the NI increased the R 2 to 0.91.
Figure 8 A univariate least-squares regression of shape for all mean PC 1 values onto all Neognathodus Index (NI) values. The NI values represent the previously published patterns of Neognathodus morphologic change through the Desmoinesian (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Brown and Weinrick1996, Reference Rexroad, Brown, Devera and Suman1998, Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009). Regressing the NI onto only PC 1 provides an R 2 of 91%, while using all of the PCs decreases the R 2 to 67%. LS, limestone; Mbr, member; SH, shale.
The maximum-likelihood tests for evolutionary modes resulted in almost equally strong support for a URW with no directional component and a GRW with a directional component (Table 4). Akaike weights indicate that the model with the highest support is URW at 0.539, but GRW also has substantial support at 0.405, while stasis has the least support at 0.055.
Table 4 Table summarizing results of maximum-likelihood evolutionary model fitting for patterns of Neognathodus morphologic change through the Desmoinesian in the Illinois Basin. Thirty PC 1 values from each of the six stratigraphic intervals were used to assess evolutionary model support for three models, an unbiased random walk (URW) with no directional component, a general random walk (GRW) that has a directional component, and a stasis model (Stasis). Terminology and methods follow Hunt (Reference Hunt2006). K is the number of parameters in each model; AICc is the small-sample, unbiased Akaike information criterion; σstep is the estimated step variance for the URW and GRW models; μstep and θ are the estimated step means for the GRW and Stasis models, respectively; ω is the estimated step variance for the Stasis model; and Akaike weights show model support. Variance was not combined across stratigraphic intervals due to significant variance differences. The σstep parameter is rounded to zero but is actually just very small (<0.00000000001). URW (0.539) and GRW (0.405) models provide the strongest support (in bold), and Stasis (0.055) supplies the least support.
Determination of GM biozones resulted in the creation of four distinct biozonations through the Desmoinesian for all examined stratigraphic intervals (Fig. 9B). The GM biozones for the stratigraphic units are assigned as follows: Perth Limestone Member lies in the N. bothrops Zone; Holland Limestone Member resides in the N. medadultimus and N. medexultimus Zone; Mecca Quarry Shale, Velpen Limestone, Alum Cave Limestone, and Providence Limestone Member lie in the N. roundyi Zone; and the West Franklin Limestone Member resides in the N. dilatus and N. metanodosus Zone.
Figure 9 Comparison of published Neognathodus biozones with new geometric morphometric (GM) results. A, Percent abundance of each Neognathodus morphotype group through the Desmoinesian and named biozone for each stratigraphic interval. B, Percent abundance of each Neognathodus GM group and the resulting biozone name for each interval. Although the structures and names of the biozones are identical, the literature and GM biozonations were determined using significantly different methodologies. The biozones for the literature data were determined and named using the Neognathodus Index. In contrast, the GM biozones were determined with statistical tests of GM-based shape and named using the percent abundances of GM groups. The GM biozones (right side) use our synonymized morphotype pairs of N. medadultimus/N. medexultimus to N. medadultimus and N. dilatus/N. metanodosus to N. dilatus. Gr., group; LS, limestone; SH, shale.
Discussion
■
Implications for Neognathodus Taxonomy
Variability and transitional forms between Neognathodus morphotypes were recognized by Merrill (Reference Merrill1972, 1975a,Reference Merrillb, 1999), and our GM analysis alleviates the problem of documenting and classifying gradual morphologic transitions between morphotype groups. GM analysis does not place specimens into separate bins, and results illustrate and record the entire gradual spectrum of shape variability among morphotypes (Fig. 4). Morphotype groups transition along PC 1, indicating the shape traits used for morphotype group classification are best captured by PC 1 (Fig. 3C). Subsequent PCs are valuable for capturing additional aspects of shape variability, and all nonzero PCs should be used. For example, PC 2 differentiates N. roundyi more effectively than PC 1, because PC 2 captures the unique curvature of the outer parapet. Using all PCs provides full quantitative characterization of these gradual morphology shifts and allows for thorough and objective statistical testing of shape similarity.
Analysis of morphotype statistical tests (Fig. 5) and the morphotype PC plot (Fig. 3C) shows that GM groups are not entirely congruent with previous morphotype groups. This result supports our initial hypothesis of some dissimilarity among GM and morphotype groups. The statistically distinct GM groups align with and support previously established classifications of N. bassleri (Harris and Hollingsworth Reference Harris and Hollingsworth1933), N. bothrops Merrill, Reference Merrill1972, and N. roundyi (Gunnell Reference Gunnell1931). In contrast, the statistically indistinct pairs of GM groups are incongruent with and do not support literature designations of N. medadultimus Merrill, Reference Merrill1972/N. medexultimus Merrill, Reference Merrill1972 and N. dilatus (Stauffer and Plummer Reference Stauffer and Plummer1932)/N. metanodosus Merrill, Reference Merrill1975b. Coalescing indistinct morphotype groups decreases the number of unique groups (i.e., diversity) across most sampled localities (Table 3). Despite the apparent loss of any diversity, GM analysis still records the full spectrum of morphologic variability for all photographed specimens in the form of objective and quantitative PC distributions.
In both indistinct morphotype pairs, we judge that GM statistical results provide sufficient support to officially synonymize each morphotype pair. The key morphologic difference between N. medadultimus and N.medexultimus is that the outer posterior parapet of N. medexultimus is relatively shorter and shows slightly more anterior fusion with the carina (Merrill Reference Merrill1972). The primary difference between N. dilatus and N. metanodosus is that the inner posterior parapet in N. metanodosus is relatively shorter and exhibits slightly more anterior fusion with the carina (Merrill Reference Merrill1975b). Both the outer and inner posterior parapet length, shape, curvature, and position relative to the carina are captured effectively with our GM outlines, and the difference between the two pairs remains statistically indistinct. Therefore, we propose N. medadultimus and N. medexultimus be synonymized as N. medadultimus, and N. dilatus and N. metanodosus be synonymized as N. dilatus (synonymy naming prioritizes the first-named morphotype within each pair). Detailed systematics for the synonymy are presented in the Systematic Paleontology section of the Supplementary Material.
Regardless of the differences between the numbers of GM and morphotype groups, the observed overlap among morphotype specimens in PC morphospace (Fig. 3C) supports the published concept that Neognathodus species consist of populations of morphotypes (Merrill Reference Merrill1972, Reference Merrill1999). GM results exhibit a gradual progression from one morphotype group to the next (Fig. 3C). Even if morphotype groups were coalesced following GM group results, the concept of Neognathodus species as populations of morphotypes (Merrill Reference Merrill1972, Reference Merrill1999) would remain valid. Overall, the quantitative characterization of form variation, along with the objective verification of Neognathodus species concepts without reliance on species-level taxonomy, further validate the use of a landmarked-based GM approach to conodont morphometric research. Our work adds additional support to the research that shows outline-based morphometric analyses effectively characterize shape variability within a given conodont genus (Klapper and Foster Reference Klapper and Foster1993; Girard et al. Reference Girard, Renaud and Korn2004a, Reference Girard, Renaud and Feist2007). Given this validation, it is vital to evaluate how GM analysis describes Neognathodus shape change through time and how GM group results affect previous interpretations of evolution and biostratigraphic correlations.
Evaluating Neognathodus Evolution
Previously, the NI was used to interpret Neognathodus evolution through the Desmoinesian in the Illinois Basin (Brown et al. Reference Brown, Rexroad and Zimmerman2016). As mentioned earlier, the NI is a weighted average of the morphotype population (i.e., the species) in a given stratigraphic unit (Brown and Rexroad Reference Brown and Rexroad2009). Using the NI, Brown et al. (Reference Brown, Rexroad and Zimmerman2016) summarized the following morphologic patterns: (1) change from the Perth to the Holland Limestone Members of the Staunton Formation and to the Mecca Quarry Shale and Velpen Limestone Members of the Linton Formation, (2) stability from the Mecca Quarry Shale and Velpen Limestone Members of the Linton Formation to the Alum Cave and Providence Limestone Members of the Dugger Formation, and (3) change from Providence Limestone Member of the Dugger Formation to the West Franklin Limestone Member of the Shelburn Formation (Fig. 7D).
Statistical tests and PC plots of Neognathodus samples from each examined stratigraphic interval support our hypothesis that GM patterns of Neognathodus P1 element shape change are generally similar to established patterns of morphologic change based on the NI. We interpret the following patterns through the Desmoinesian: (1) morphologic change occurs from the Perth to the Holland Limestone Members of the Staunton Formation and to the Mecca Quarry Shale and Velpen Limestone Members of the Linton Formation, (2) relative morphologic stability persists from the Mecca Quarry Shale and Velpen Limestone Members of the Linton Formation to the Alum Cave and Providence Limestone Members of the Dugger Formation, and (3) morphologic change again occurs from the Providence Limestone Member of the Dugger Formation to the West Franklin Limestone Member of the Shelburn Formation (Fig. 7). These morphologic patterns provide a fully quantitative characterization of Neognathodus shape change progressing into extinction of the genus at the end of the Desmoinesian (Merrill Reference Merrill1975a; Merrill and Grayson Reference Merrill and Grayson1989; Boardman et al. Reference Boardman, Heckel, Barrick, Nestall and Peppers1990), and they also align with a summary of published NI trends (Brown et al. Reference Brown, Rexroad and Zimmerman2016). Furthermore, GM patterns of P1 element shape change through the Desmoinesian retained a strong statistical correlation with the NI (Fig. 8). These morphologic patterns should be considered with additional evidence from lithologic information, previous interpretations of the depositional environment, and results of maximum-likelihood tests of evolutionary modes.
Overall, patterns of Neognathodus P1 element morphology through the Desmoinesian appear to be relatively independent of the lithology, sampling location, and depositional setting. The regional paleoenvironment remained relatively consistent through the Desmoinesian in the Illinois Basin. During this time, the basin was located along tropical latitudes and likely maintained a tropical to humid subtropical climate that was subject to wet and dry seasons (Cecil Reference Cecil1990). The lithologies for each examined member consist primarily of discontinuous limestone and interbedded shale (Shaver et al. Reference Shaver, Ault, Burger, Carr, Droste, Eggert, Gray, Harper, Hasenmueller, Hasenmueller, Horowitz, Hutchison, Keith, Keller, Patton, Rexroad and Weir1986), and the depositional setting is interpreted as a nearshore marine environment with deltaic influence (Ferm et al. Reference Ferm, Horne, Swinchatt and Whaley1971; Trask and Palmer Reference Trask and Palmer1986; Rexroad et al. Reference Rexroad, Wade, Merrill, Brown and Pagett2001; Brown and Rexroad Reference Brown and Rexroad2009; Brown et al. Reference Brown, Rexroad and Zimmerman2016). A specific example of the independence between morphologic pattern and lithologies is evident between the Perth and Holland Limestone Members. Both consist of gray, fossiliferous limestone and belong to the same formation (Shaver et al. Reference Shaver, Ault, Burger, Carr, Droste, Eggert, Gray, Harper, Hasenmueller, Hasenmueller, Horowitz, Hutchison, Keith, Keller, Patton, Rexroad and Weir1986), but morphologic change is interpreted between the two units. In contrast, morphologic stability is interpreted between the clastic-rich Mecca Quarry Shale and Velpen Limestone Members of the Linton Formation and the non-clastic Alum Cave and Providence Limestone Members of the Dugger Formation. These examples indicate that temporal patterns of Neognathodus P1 element morphology were not likely driven by localized environmental factors and do not obviously represent ecophenotypic variation. Instead, we evoke morphologic evolution as the major driver affecting P1 element morphology through the Desmoinesian.
Maximum-likelihood tests of evolutionary modes indicate the most likely models of P1 element morphologic evolution are a URW or a GRW (Table 4). Hunt (Reference Hunt2006) states that only an evolutionary mode with an Akaike weight greater than 0.90 should be treated confidently as the most likely model. As such, we recognize our results show equivocal support for both URW (0.539) and GRW (0.405) modes. By the standards set by Hunt (Reference Hunt2006), the mean for GRW is relatively small (μstep=−0.033), indicating that if GRW was the true model, then the directionality of the steps was relatively small. We also interpret the relatively low Akaike weight of the stasis model (0.055) to be significantly smaller than the other mode weights, and we support that stasis is an unlikely model for Neognathodus PC 1 evolution.
An important consideration is that restricting the analysis to only PC 1 values can potentially bias results toward the GRW model, because PC 1 shows the axis of greatest variance and likely shows a directional pattern in multivariate space (Bookstein Reference Bookstein2013). Our results show statistically similar PC 1 means between certain stratigraphic intervals, and this suggests a directional bias is not inherent with our PC 1 data. Furthermore, Jones (Reference Jones2009) performed the same model tests (Hunt Reference Hunt2006) with PC 1 data from conodont P1 elements and demonstrated GRW was not overly favored.
Overall, our work supplies the first statistical assessment of Neognathodus evolutionary models in the Desmoinesian. Previous interpretations of Neognathodus evolution through the Desmoinesian provided highly useful generalized patterns of morphologic change based on the NI (Brown et al. Reference Brown, Rexroad and Zimmerman2016), and our GM work builds on that foundation by providing statistically significant, maximum-likelihood test results that support URW and GRW modes. Hunt’s (Reference Hunt2006) methods to interpret modes of evolution with maximum likelihood are applied by many researchers (Novack-Gottshall Reference Novack-Gottshall2008; Jones Reference Jones2009; Piras et al. Reference Piras, Sansalone, Marcolini, Tuveri, Arca and Kotsakis2012; Van Bocxlaer and Hunt Reference Van Bocxlaer and Hunt2013), and our results contribute to this work of documenting evolutionary patterns of significant fossil taxa. In sum, using our GM results with Hunt’s (Reference Hunt2006) evolution tests provides a fully quantitative means of characterizing conodont morphology and evolutionary partners without a priori classification of morphotypes or assumption of evolutionary modes.
Assessing Biostratigraphic Utility
Neognathodus biozones have been used to correlate strata in the Illinois Basin throughthe Desmoinesian, and the NI was used to determine the Neognathodus biozone of each stratigraphic interval examined in this study (Brown et al. Reference Brown, Rexroad and Zimmerman2016). The morphotype name closest to the NI value represents the biozone name (1=N. bassleri, 2=N. bothrops, 3=N. medadultimus, 4=medexultimus, 5=N. roundyi, 6=N. dilatus, and 7=N. metanodosus). For example, the Perth Limestone Member has an NI of 2.0, thus the member resides in the N. bothrops biozone (Rexroad et al. Reference Rexroad, Brown, Devera and Suman1998). Biozonations for the remaining units are as follows: the Holland Limestone Member belongs to the N. medadultimus Zone (Rexroad et al. Reference Rexroad, Brown and Weinrick1996); the Mecca Quarry Shale, Velpen Limestone, Alum Cave Limestone, and Providence Limestone Members all fall in the N. roundyi Zone (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991, Reference Brown, Rexroad and Zimmerman2016; Rexroad et al. Reference Rexroad, Wade, Merrill, Brown and Pagett2001); and the West Franklin Limestone Member belongs to the N. dilatus Zone (Brown and Rexroad Reference Brown and Rexroad2009) (Fig. 1C). In sum, previous, NI-based work shows these Illinois Basin units reside in four Neognathodus biozones.
Our GM-based biozonations support the previous structure of four separate biozones through the Desmoinesian in the Illinois Basin units (Fig. 9) and our hypothesis that GM analysis would yield four similar biozonations. GM results coalesced the morphotype group pairs of N. medadultimus/N. medexultimus and N. dilatus/N. metanodosus (Fig. 5), but the structure of four distinct biozones remains intact (Fig. 9). Also, our GM-based biozonations show that stratigraphic intervals with statistically indistinct shapes lie within the same biozone. For example, N. roundyi is the published biozone for the Mecca Quarry Shale and Velpen Limestone Members (Rexroad et al. Reference Rexroad, Wade, Merrill, Brown and Pagett2001), the Alum Cave Limestone Member (Brown et al. Reference Brown, Rexroad and Zimmerman2016), and the Providence Limestone Member (Fig. 1) (Brown et al. Reference Brown, Rexroad, Eggert and Horowitz1991). GM results show all of these units are also statistically indistinct in mean shape (Fig. 7). In addition, the stratigraphic intervals with statistically distinct mean shapes each reside in a different biozones, with the Perth Limestone Member in the N. bothrops Zone (Rexroad et al. Reference Rexroad, Brown, Devera and Suman1998), the Holland Limestone Member in the N. medadultimus Zone (Rexroad et al. Reference Rexroad, Brown and Weinrick1996), and the West Franklin Limestone Member in the N. dilatus Zone (Brown and Rexroad Reference Brown and Rexroad2009).
The identical biozone structure between previous NI work and our GM analysis illustrates that Neognathodus-based biostratigraphic correlations would not change between methods. Thus, previous Neognathodus-based correlations both within the Illinois Basin and around North America would remain intact. Examples within the Illinois Basin include the correlation of the Alum Cave Limestone Member in Indiana with the St David Limestone Member of the Carbondale Formation in Illinois (Brown et al. Reference Brown, Rexroad and Zimmerman2016) and the correlation of the West Franklin Limestone Member of Indiana to the Lonsdale Limestone Member of the Shelburn Formation in northwestern Illinois (Brown and Rexroad Reference Brown and Rexroad2009). Additional examples across North America include correlating units from Ohio (Merrill Reference Merrill1972), Texas (Grayson et al. Reference Grayson, Trice and Westergaard1985), and New Mexico (Brown et al. Reference Brown, Rexroad and Zimmerman2013) to units within the Illinois Basin.
Our GM-based biozonations also provide the first direct evaluation of how morphometric results compare with previous biozonations. Prior studies have coalesced certain conodont species that were useful for biostratigraphy and zonation (Renaud and Girard Reference Renaud and Girard1999; Girard et al. Reference Girard, Renaud and Korn2004a, Reference Girard, Renaud and Feist2007; Roopnarine et al. Reference Roopnarine, Murphy and Buening2004), but these studies have not evaluated how morphometric-based biozones can be compared with previous, qualitative-based biozonations. Our GM analysis shows that objective and statistically robust morphologic results can be used to construct biozonations without relying on previous morphotype-level taxonomy. The structural similarity between previous NI-based and GM-based biozonations showcases that determining GM-based biozones is neither redundant nor unwarranted, as this comparison validates the use of landmark-based GM work for constructing viable biozonations and for subsequent stratigraphic correlations.
Conclusions
Our study provides the first morphometric work on the conodont genus Neognathodus, supplies the first statistical analysis of Neognathodus evolution in the Desmoinesian, and provides a broad-scale methodology to quantitatively test morphotype designations, interpret statistically robust morphologic change through time, and construct valid biozones. GM analysis fully differentiates N. bassleri, N. bothrops, and N. roundyi from all other Neognathodus morphotype groups and supports previous morphotype designations (Gunnell Reference Gunnell1931; Harris and Hollingsworth Reference Harris and Hollingsworth1933; Merrill Reference Merrill1972). In contrast, there is no statistically significant difference between the species pairs of N. medadultimus/N. medexultimus and between N. dilatus/N. metanodosus, suggesting these GM group pairs are incongruent with and do not support established species designations (Stauffer and Plummer Reference Stauffer and Plummer1932; Merrill Reference Merrill1972, Reference Merrill1975b). As such, we synonymize these pairs to N. medadultimus and N. dilatus, respectively. Coalescing taxonomic boundaries does not adversely affect patterns of morphologic change or biozonation structure. GM analysis of morphology patterns through time still shows periods of morphologic change and stability through the Desmoinesian that support previously published patterns (Brown and Rexroad Reference Brown and Rexroad2009; Brown et al. Reference Brown, Rexroad and Zimmerman2016). Maximum-likelihood statistical tests of evolutionary modes indicate the most probable models of P1 element morphologic evolution are a URW (i.e., a Brownian motion random walk) or a GRW (i.e., directional change). Finally, our GM-based biozonations corroborate the previous structure of four distinct biozones in the Illinois Basin units through the Desmoinesian.
Although this study location is limited to the Illinois Basin, our quantitative methodology can be broadly applied to other geologic settings and additional conodont genera. Similarly, our work contributes an example of how GM results can be used with Hunt’s (Reference Hunt2006) maximum-likelihood tests to fully quantify conodont morphology and evolutionary patterns without categorizing morphotypes or assuming a null evolutionary model a priori. Previous morphometric outline analyses conducted on other conodont genera, such as the Late Devonian genus Palmatolepis, also exhibit a spectrum of transitional forms between species (Klapper and Foster Reference Klapper and Foster1993; Girard et al. Reference Girard, Renaud and Serayet2004b, Reference Girard, Renaud and Feist2007). Although separated by more than 50 million years, both Palmatolepis and Neognathodus show similar gradational species boundaries. This indicates that it is valid and useful to pursue a standardized morphometric methodology of photography, outlining, and statistical testing to describe shape variation among many different conodont genera. In addition, the genera Idiognathodus and Streptognathodus are also used for Middle Pennsylvanian correlation, and species in both genera exhibit gradual shape variation between species morphologies (Barrick and Boardman Reference Barrick and Boardman1989; Heckel Reference Heckel1991; Lambert Reference Lambert1992; Barrick et al. Reference Barrick, Boardman and Heckel1996, Reference Barrick, Lambert, Heckel and Boardman2004). Thus, the genera are prone to the same inherently subjective taxonomic debate that affects Neognathodus zonation. Our quantitative approach to characterizing the range and significance of shape variation within each genus would form objective groups of species and may ultimately provide a robust framework with biostratigraphic utility in which to develop data-driven hypotheses for this significant chordate group.
Acknowledgments
We would like to thank L. Brown, professor emeritus in the Department of Geology at Lake Superior State University, and P. Novack-Gottshall, assistant professor in the Department of Biological Sciences at Benedictine University, for providing thorough and valuable reviews that significantly enhanced the quality of this manuscript. This research was supported in part by the Galloway-Perry-Horowitz Fellowship of Indiana University.
Supplementary Material
Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.1fm52vv
