Including autapomorphies is important for paleontological tip-dating with clocklike data, but not with non-clock data

Tip-dating, where fossils are included as dated terminal taxa in Bayesian dating inference, is an increasingly popular method. Data for these studies often come from morphological character matrices originally developed for non-dated, and usually parsimony, analyses. In parsimony, only shared derived characters (synapomorphies) provide grouping information, so many character matrices have an ascertainment bias: they omit autapomorphies (unique derived character states), which are considered uninformative. There has been no study of the effect of this ascertainment bias in tip-dating, but autapomorphies can be informative in model-based inference. We expected that excluding autapomorphies would shorten the morphological branchlengths of terminal branches, and thus bias downwards the time branchlengths inferred in tip-dating. We tested for this effect using a matrix for Carboniferous-Permian eureptiles where all autapomorphies had been deliberately coded. Surprisingly, date estimates are virtually unchanged when autapomorphies are excluded, although we find large changes in morphological rate estimates and small effects on topological and dating confidence. We hypothesized that the puzzling lack of effect on dating was caused by the non-clock nature of the eureptile data. We confirm this explanation by simulating strict clock and non-clock datasets, showing that autapomorphy exclusion biases dating only for the clocklike case. A theoretical solution to ascertainment bias is computing the ascertainment bias correction (Mkparsinf), but we explore this correction in detail, and show that it is computationally impractical for typical datasets with many character states and taxa. Therefore we recommend that palaeontologists collect autapomorphies whenever possible when assembling character matrices.

Starting tree for Beast2 analyses. Obtaining starting trees that would match all of the "hard" constraints mandated by the uniform tip-date priors proved difficult, so initial analyses were conducted using a random starting tree, with all of the uniform date constraints converted to "soft" normal distributions. This was effected by setting the "convert_to_normal" field to "yes" for each tip-date in the BEASTmasteR Excel settings file, an option created to avoid tedious manual re-setting of the prior for each tip-date in BEAUTi. Once this analysis had burned in, the last sampled tree served as a valid starting tree under the hard constraints.
Implementation of Mkv and Mk parsinf . We implemented Mkv and Mk parsinf in BEASTmasteR with functions (num_unobservable_patterns_ParsInf, list_unobservable_patterns_ParsInf, make_table_num_unobservable_patterns) that list all unobservable site patterns under these corrections as dummy characters. The dummy characters are included in the Beast2 XML using the "ascertained," "excludefrom," and "excludeto" options of the "data" XML tag.
MCMC runs. The relatively small size of the morphological dataset and the large uncertainties in some tip-dates necessitated substantial sampling to achieve burn-in and sufficient mixing. However, the small size of the data also sped calculation, so each analysis was run for 10 9 generations, with each run taking approximately 12 hours on a 2015 Mac Pro, for the Mk and Mkv analyses. However, the Mk parsinf ascertainment bias correction requires calculating the likelihood over all unobservable site patterns. For the  dataset , this adds up to 52 additional site patterns for the 2-state characters, and 1953 additional patterns for the 3-state characters (calculated with the function num_unobservable_patterns_ParsInf; see Appendix 1). An Mk parsinf run of 10 9 generations took approximately 15.25 days on the same machine.
Trees and parameters were sampled every 500,000 generations, resulting in 2000 saved samples. Burn-in was achieved well before this, but to be conservative, the first 25% of the samples were discarded for calculating the summary tree and parameter 95% highest posterior densities (HPDs).
Simulation analysis. In order to test our explanation of why excluding autapomorphies did not influence dating inference in the empirical eureptiles dataset, we constructed two simulation that would dramatically show the difference between a clocklike and non-clocklike dataset.
Two character datasets were simulated on the phylogeny, using sim.char from the R package ape: 1. "strict clock" --simulate 1000 binary characters under an Mk model. The rate was set to be low enough (0.05) that many characters (577/1000) are invariant.
2. "non-clock" --the time branchlengths on the true tree were randomly reshuffled once. The 1000 characters were then simulated on this new tree (by chance, again 577/1000 were invariant). The tip dates are still from the true tree.
Filtering the character datasets according to the possible ascertainment biases produced: -two "alldata" datasets of 1000 characters of all types -two "variable" datasets of 423 variable characters (both had 423, by chance) -two "informative" datasets of 190 (strict clock) and 352 (no clock) characters Twelve Beast2 inference runs (BDSS, lognormal relaxed clock, gamma-distributed sitewise rate variation with 4 sites -the same settings as for the empirical eureptile analyss) were then done, using each applicable combination of Mk, Mkv, and MkParsInf on each of these datasets.
The names of the folders in Supplemental Files containing these analyses are as follows.
On the strict clock dataset: 01_Mk_on_strict_clock_alldata  02_Mk_on_strict_clock_variable  03_Mk_on_strict_clock_informative  04_Mkv_on_strict_clock_variable  05_Mkv_on_strict_clock_informative   5   06_MkParsInf_on_strict_clock_informative   On the non-clock dataset:   07_Mk_on_no_clock_alldata  08_Mk_on_no_clock_variable  09_Mk_on_no_clock_informative  10_Mkv_on_no_clock_variable  11_Mkv_on_no_clock_informative  12_MkParsInf_on_no_clock_informative Each run took 100 million generations, sampling every 50,000. The trace files were inspected for convergence in Tracer (all converged easily), MCC trees were generated with TreeAnnotator, and the results were compared with custom R scripts (Supplemental Files; directory _06_TreeSim).

Supplemental Results
Morphological matrix. The data matrix was 84.64% complete, or 82.76% complete when autapomorphies were excluded (Supplemental Table S3).
Parsimony analysis. As expected, parsimony analyses of the two datasets (Figs S1-S10) yielded an identical strict consensus tree, with bootstrap values very similar to that of   . They reported finding 4 trees with length (TL) of 252 parsimonyinformative steps, CI=0.4405, and RI=0.6493. Here, the TNT analysis yielded found 2 trees with TL=253, CI=0.395, and RI=0.619. This may be due to a slight difference in the TreeBase dataset, or calculating the statistics on the 50% majority rule consensus tree (Müller & Reisz) versus the strict consensus tree (this study), which is slightly less resolved. As expected, MP trees had identical topology between the full and no-autapomorphies datasets, but some terminal branches were longer in the full dataset (compare Fig. S1 vs. S6; the whole tree was also longer: TL=295, i.e. 42 steps longer, corresponding to adding 42 autapomorphies). The most dramatic example of this was Procolophonidae, for which the terminal branchlength changed from 9 to 22 steps (Fig. S6).
Test for correlation between time and character change. Figures S11 and S12 show the results of regression of tip age on the number of parsimony character steps above the root. Neither dataset shows impressive correlation with time, probably in part because of the uncertainty in some OTU time ranges, and because many of the sampled taxa come from about the same time period, 315-300 Ma. Nonetheless, including autapomorphies does make the relationship somewhat more positive (without autapomorphies: slope=0.08, R 2 =0.015, P=0.56; with autapomorphies: slope=0.17, R 2 =0.067, P=0.21).
If the main purpose of this study were a "best possible" dating effort, the observed lack of correlation would be a reason for concern (although, in tip-dating, a substantial part of the signal informing the dating inference is from the dates of the tips, not only the morphological clock model), along with the non-uniform sampling of OTUs in time. However, the primary goal of the present study is to compare the impact of autapomorphy-related decisions on otherwise identical datasets.
Changes in estimates of morphological relaxed clock model parameters. Although the difference in clock rates between autapomorphy-including and autapomorphy-excluding runs is sufficiently obvious that no test is needed (the average clock rate estimate is 10 times higher when autapomorphies are excluded), some readers may prefer that the tests be presented. Tests were run on post-burnin MCMC samples of the clock mean rate parameter (N=1500 postburnin samples of clock rate, which are uncorrelated in this case, due to only 1500 samples being taken from a run of 10 9 ). Testing the null hypothesis of no difference in means between the samples of clock rates yielded P<2.2e-16 for all variants, using t-tests and Wilcoxon signed rank tests, on both untransformed and log-transformed rate samples. All tests were two-sided and non-paired. The standard deviation of the relaxed clock rate across branches is also higher in the no-autapomorphies dataset (Table 1, all P<2.2e-16, same tests).
Comparison of dates and HPDs across all 10 runs. On this dataset, results using PaleobioDB dates closely paralleled those using best-practices dates (e.g., Fig. 1e), although dating uncertainty was higher for some tips and nodes. Figure S13 contains the pairwise correlation plots between the mean node ages for the shared nodes found in the summary trees of all 10 analyses (5 analyses under the best-practices tip dates, plus the same 5 analyses under the PaleobioDB tip dates). Figure S14 contains the same for the 95% HPD widths.

Supplemental Discussion
The most dramatic difference in parameter inference on the eureptile dataset was the order-ofmagnitude difference in the mean rate parameter of the relaxed morphological clock model between analyses that include, or exclude, autapomorphies. The reason for this difference seems clear. In morphological clock analyses, autapomorphic characters are interpreted as lowrate characters. In parsimony terms, each autapomorphic state only occurred once and thus comprises only one change on the tree. In probabilistic models, there is a chance for multiple changes, but the expectation is nevertheless low. Some parsimony-informative characters, on the other hand, might have the minimum number of changes, if they exhibit no homoplasy. However, in real datasets, homoplasy is common, and this will translate to higher rates in a model-based framework. The inclusion of 42 low-rate characters in the data matrix thus shifts the distribution of rates downward, as well as the overall mean.
The effect of this change in rates on inferred dates, date uncertainty, and posterior probabilities were modest in the studied dataset. However, they are predictable consequences of the models employed in tip-dating, and thus important to understand. The increase in clock rate 7 under a no-autapomorphies dataset detectably decreased bipartition posterior probabilities, and increased dating uncertainty. The increase in rates has these effects because increasing rates increases the probability of multiple convergent character changes, and thus increases the probability of alternate topologies that fit the character data almost equally well. This translates into increased uncertainty.  Figure S1. Strict consensus tree from MP analysis of the autapomorphyexcluding dataset. Branchlengths, and the numbers on each branch, represent the number of parsimony steps as estimated by TNT. Figure S2. Strict consensus tree from MP analysis of the autapomorphyexcluding dataset. The numbers on each branch represent the number of synapomorphies as estimated by TNT. Figure S3. Strict consensus tree from MP analysis of the autapomorphyexcluding dataset. The numbers on each branch represent the Bremer support (decay index) as estimated by TNT's aquickie.run script. Figure 4. Strict consensus tree from MP analysis of the autapomorphy-excluding dataset. The numbers on each branch represent the relative Bremer support as estimated by TNT's aquickie.run script. Figure S5. Strict consensus tree from MP analysis of the autapomorphyexcluding dataset. The numbers on each branch represent the bootstrap support frequency (out of 100 trees) as estimated by TNT's aquickie.run script. Figure S6. Strict consensus tree from MP analysis of the autapomorphy-including dataset. Branchlengths, and the numbers on each branch, represent the number of parsimony steps as estimated by TNT. Figure S7. Strict consensus tree from MP analysis of the autapomorphy-including dataset. The numbers on each branch represent the number of synapomorphies as estimated by TNT. Figure S8. Strict consensus tree from MP analysis of the autapomorphy-including dataset. The numbers on each branch represent the Bremer support (decay index) as estimated by TNT's aquickie.run script. Figure S9. Strict consensus tree from MP analysis of the autapomorphy-including dataset. The numbers on each branch represent the relative Bremer support as estimated by TNT's aquickie.run script. Figure S10. Strict consensus tree from MP analysis of the autapomorphyincluding dataset. The numbers on each branch represent the bootstrap support frequency (out of 100 trees) as estimated by TNT's aquickie.run script. Figure S11. Linear regression of tip height (mean of the prior distribution for that tip) against number of character steps above the root on the MP tree, for the autapomorphyexcluding dataset. Figure S12. Linear regression of tip height (mean of the prior distribution for that tip) against number of character steps above the root on the MP tree, for the autapomorphyincluding dataset.

Supplemental
Supplemental Figure S13. Pairwise correlation plots between the mean node ages for the shared nodes found in the summary trees of all 10 analyses. Figure S14. Pairwise correlation plots between the node age HPD widths (95% highest posterior densities) for the shared nodes found in the summary trees of all 10 analyses. shows the same analysis using preliminary dates obtained from PaleobioDB, done in order to compare the effects of a "naïve" dating analysis with dates taken directly from a database, versus an expert review.

Supplemental Tables
Supplemental Table S1. Date ranges for Operational Taxonomic Units (OTUs) in  (eureptiles and outgroups), following the "best practices" recommended by Parham et al. (2012). Table S2. Date ranges for Operational Taxonomic Units (OTUs) in  (eureptiles and outgroups) used in a preliminary "practice". Dates derived from FossilWorks/the Paleobiology Database (PaleobioDB). Table S3. Summary statistics on morphology data matrix from . Calculated with BEASTmasteR functions. Table S4. Comparison of summary statistics from the five Beast2 runs using tip dates derived from the Paleobiology Database. Table S5. Number of patterns that are unobservable in the Mk parsinf model, which excludes parsimony-uninformative characters (invariant and autapomorphic characters). The ascertainment-bias correction requires that the likelihood be calculated for each of these patterns, which obviously becomes problematic when there are millions of such patterns. The calculations here assume unordered characters. This table duplicates Table 2, but is included in Supplemental Material for comparison. For ordered characters, see Supplemental Table S8.

Supplemental
Supplemental Table S6. Total number of possible patterns, determined by (# of states) (# of taxa) . "Inf" is not literally infinity, it just means the result exceeds R's numerical limit.
Supplemental Table S7. Fraction of the total number of possible patterns that are unobservable under Mk parsinf ascertainment bias correction. The total number of possible patterns is determined by (# of states) (# of taxa) . For ordered characters, see Supplemental Table S10.
Supplemental Table S8. Number of patterns that are unobservable in the Mk parsinf model, assuming ordered characters. For unordered characters, see Supplemental Table S5. Table S9. Total number of possible patterns, determined by (# of states)^(# of taxa). "Inf" is not literally infinity, it just means the result exceeds R's numerical limit.

Supplemental
Supplemental Table S10. Fraction of the total number of possible patterns that are unobservable under Mk parsinf ascertainment bias correction, assuming ordered characters. For unordered characters, see Supplemental Table S7.
Supplemental Table S11. Comparison of summary statistics from the six Beast2 runs on a simulated "strict clock" dataset. Table S12. Comparison of summary statistics from the six Beast2 runs on a simulated "no clock" dataset.               The potential oldest seymouriamorph is Utgenia from the Upper Carboniferous-Lower Permian of Russia (Klembara & Ruta 2004). The age of this taxon is not well-constrained, so we conservatively give it a maximum age of Upper Pennsylvanian, whose lower boundary is dated to 307 ± 0.1 Ma in the 2012 GTS. The youngest seymouriamorph is Karpinskiosaurus (=Kotlassia) from Russia (Klembara 2011), whose stratigraphic range extends to near or the end of the Permian Period (Benton et al. 2004), which is dated to 251.9 ± 0.02 Ma by Burgess et al. (2014).

Caseidae 302 257
The recently described Eocasea is the oldest-known caseid and Virgilian/Stephanian in age (Reisz & Fröbisch 2014). The upper boundary of the Stephanian is between 302-300 Ma in the 2012 GTS. The youngest well-dated caseid is Ennatosaurus from Russia, which is Tatarian (Maddin et al. 2008). Reisz et al. (2011) Ruthenosaurus and Euromycter could be as young as early Lopingian, but new radioisotopic ages from the nearby correlative Lodève Basin make an Artinksian-Kungurian age more likely (Michel et al. 2015). The Tartarian is 269-257 Ma according to the 2012 GTS.
Millerettidae 265.5 251.9 The oldest millerettid is Broomia from the Tapinocephalus Assemblage Zone of southern Africa ; the minimum age of this biostratigraphic zone is constrained by U-Pb dates (Rubidge et al. 2013;Day et al. 2015), but its maximum age is not. Therefore, we are conservative in using the lower boundary of the Capitanian as a maximum age for this taxon; this is dated to 265.1 ± 0.4 Ma in the 2012 GTS. The youngest millerettids occur in the Dicynodon Assemblage Zone (Rubidge 2005;Cisneros et al. 2008), the upper boundary of which is thought to coincide with the Permian-Triassic boundary. This is dated to 251.9 ± 0.02 Ma (Burgess et al. 2014).
Procolophonidae 260.3 201.6 The maximum age for Procolophonidae is difficult to constrain, because fossils are either phylogenetically well-constrained but poorly dated (Pintosaurus), or poorly-constrained phylogenetically but relatively-well dated (Spondylolestes and Kinelia) (Cisneros 2008). Nonetheless, these data collectively suggest a Late Permian origin for procolophonids, which is consistent with the Late Permian appearance of their sister clade Owenettidae. Thus, we conservatively use the base of the Wuchiapingian as the maximum age of Procolophonidae, which is dated to 259.9 ± 0.4 Ma in the 2012 GTS. Procolophonids went extinct at or near the end-Triassic mass extinction, which has been dated to 201.6 Ma (Blackburn et al. 2013).

290
Romeria is from the Archer City Formation of Texas (Hook 1989;Hentz 1989), which is equivalent to the Moran through Santa Ana Branch Shale formations to the south (Hentz 1989). This series of units is upper Asselian to Sakmarian in age (Wardlaw 2005), which spans 297 Ma to 290.1 ± 0.1 Ma based on the 2012 GTS.

Protocaptorhinus 288 282
This taxon is known from the Petrolia Formation of Texas; Hook (1989) listed occurrences from the Nocona and Waggoner Ranch formations, but these have not been described and we do not consider their stratigraphic range here. The Petrolia is equivalent to the Elm Creek, Jagger Bend, and Valera Shale formations further south (Hentz 1989), which are middle Artinskian in age (Wardlaw 2005). Based on the 2012 GTS, this gives an age range of approximately 288 to 282 Ma.

295.1
All known specimens of this taxon are from a single locality (Camp Quarry) in the Arroyo del Agua Formation (Cutler Group) of New Mexico (Lucas et al. 2005a,b). Although the placement of the Carboniferous-Permian boundary in this section is controversal, all recent authors agree that the Camp Quarry is lowermost Permian in age. We conservatively consider it somewhere within the Asselian, which spans 298.9 ± 0.2 Ma to 295.5 ± 0.4 Ma in the 2012 GTS.

276
The type of C. laticeps is from the Waggoner Ranch (=Clyde) Formation of west Texas (Heaton 1979), which correlates to the Bead Mountain through Lueders formations further south (Hentz 1989). These units are dated to the upper Artinskian through lower Kungurian (Wardlaw 2005), which is between 285 and 276 Ma in the 2012 GTS. Referred material from the Wellington Formation of Oklahoma (Heaton 1979;Modesto 1998) is thought to be lower Leonardian in age (Sawin et al. 2008), which is within the same age range. Material from the Richards Spur locality that Heaton (1979) assigned to C. laticeps has subsequently been referred to C. aguti (Bolt 1980;Modesto 1998). Captorhinus_aguti

276
C. aguti is best known by abundant remains recovered from the Richards Spur locality in Oklahoma, for which a U-Pb speleothem age of 289 ± 0.68 Ma was recently published (Woodhead et al. 2010). Potential open-system behavior (e.g., lead loss) is difficult to fully ascertain with such carbonate samples, and the duration of deposition at Richards Spur is unclear, so we prefer to take a conservative approach and consider an age range for this taxon that is larger than the analytical uncertainty for this date. Most of the material reported from Texas (Fox & Bowman 1966) has either been re-assigned to other taxa (Heaton 1979;Modesto 1998) or has not been re-evaluated since. Modesto (1998) assigned specimens from the Lueders Formation in Texas to C. aguti; Wardlaw (2005) assigned this unit to the lower Kungurian, providing a minimum age for this taxon of 276 Ma using the 2012 GTS.
Labidosaurus 279.9 271.8 Known specimens of Labidosaurus are from the lowermost portion of the undivided Clear Fork Group of west Texas (Williston 1917;Sumida 1987Sumida , 1989Modesto et al. 2007). This unit is Kungurian in age (Wardlaw 2005), but lacks more precise age constraints. The Kungurian spans 279.3 ± 0.6 to 272.3 ± 0.5 Ma in the 2012 GTS.
Labidosaurikos 291 271.8 The only known specimen of Labidosaurikos is from the Hennessey Formation of Oklahoma (Dodick & Modesto 1995). Few non-vertebrate biostratigraphic age constraints are available for this formation. Based on the vertebrate assemblage, Olson (1967Olson ( , 1970 proposed that the Hennessey was correlative with the Clear Fork Group (=Vale & Choza formations) in Texas. See the age justification of Labidosaurus above for the age range of the Clear Fork Group. In contrast, the Hennessey also shares several vertebrate taxa with the Richards Spur locality in Oklahoma, which was recently dated to 289 ± 0.68 Ma (see age justification for Captorhinus aguti). Therefore, we take the conservative approach and assign an age range to Labidosaurikos that encapsulates the age range for both potential correlations.

301
All known material of Concordia is from the Hamilton Quarry in the Calhoun Shale, Kansas (Müller & Reisz 2005). This unit is assigned to the Streptognathodus virgilicus conodont zone (Ritter 1995), which is dated to between 302.25 and 301 Ma in the 2012 GTS.

293
The type species of Protorothyris, P. archeri, is from the lower Archer City (=Moran) Formation of west Texas (Price 1937;Clark & Carroll 1973). This unit correlates to the Moran and Sedwick formations further south (Hentz 1989), which are upper Asselian to lower Sakmarian in age (Wardlaw 2005). Based on the 2012 GTS, this suggests an age range of 297-293 Ma for P. archeri. The referred species P. morani is known only from the Washington Formation of the Dunkard Group in West Virginia (Clark & Carroll 1973). This unit is poorly dated; the best age constraints come from insect and vertebrate biogeography, which suggest a broadly Asselian-Sakmarian age for the formation (Schneider et al. 2013;Lucas 2013). The closest conodont age constraint is from the Ames Limestone of the underlying Conemaugh Group, which is assigned to the Idiognathus simulator zone (Heckel et al. 2011), suggesting only that the Dunkard Group is younger than lower Gzhelian (i.e., <302 Ma based on the 2012 GTS).
Paleothyris 309.5 307.5 All known specimens were found above the Lloyd Cove/Lower Bonar Coal in the Morien Group of Nova Scotia (Carroll 1969), which is dated to the Westphalian D regional substage (Zodrow & Cleal 1985;Gibling & Bird 1994). This sub-stage has an age range of 309.5 to 307.5 Ma in the 2012 GTS.

307
The only specimen of Cephalerpeton is from the Francis Creek Shale of Mazon Creek, Illinois (Carroll & Baird 1972). Based on conodont biostratigraphy, this unit correlates to the Verdigris cyclothem in the mid-continent (Heckel 2013), which is assignable to the Neognathodus roundyi conodont zone . The age of this zone is 309 to 307 Ma in the 2012 GTS.

Anthracodromeus 309 305
This taxon is from the Upper Freeport Coal of the Allegheny Group in Ohio (Carroll & Baird 1972). The presence of the conodonts Idiognathodus and Neognathodus roundyi in the underlying Washingtonville Shale (Sturgeon & Youngquist 1949;Merrill 1972 Brouffia is known from a single specimen from the Nýřany coal measures of the Czech Republic (Carroll & Baird 1972). This unit has been dated to the Westphalian D regional sub-stage based on palynomorph biostratigraphy (Kalibová-Kaiserová 1967;Kalibová 1989;Bek 1995), which is consistent with aquatic vertebrate biostratigraphic correlations (Zajíc 2000). See Paleothyris age justification above for Westphalian D numerical age.
Hylonomus 318.5 317.5 All known specimens of Hylonomus are from the Joggins Formation of Nova Scotia (Carroll 1964), which has been dated using palynomorph biostratigraphy to the Langsettian regional stage (Dolby 1991(Dolby , 2003Calder et al. 2005). Using the 2012 GTS, this dates the formation to between 318.5 and 317.5 Ma.

271.8
Thuringothyris is represented by several specimens from the Bromacker Quarry of the Tambach Formation in Germany (Boy & Martens 1991;. This site has traditionally been dated to the Wolfcampian because it shares the vertebrate taxa Seymouria, Diadectes, and Dimetrodon with Lower Permian strata in western North America (e.g., Sumida et al. 1996Sumida et al. , 1998Berman et al. 1998Berman et al. , 2001. However, it's worth noting that the only shared species between Germany and North America is Seymouria sanjuanensis, which is only known from the poorly-dated units of the upper Cutler Groupthe Organ Rock Shale in Utah and Arroyo del Agua Formation in New Mexico (e.g., Vaughn 1966;Berman et al. 1987). Furthermore, all three genera have been reported from the Clear Fork Group in Texas (e.g., Romer 1928;Olson 1958), and both Diadectes and Dimetrodon are known from the Wellington Formation of Oklahoma (e.g., Olson 1967); both of these units could be as young as Kungurian in age (see above) and thus younger than Wolfcampian. Schneider & Werneburg (2006) assigned the Tambach Formation to the lower Artinskian because of the presence of the insect Moravamylacris kukalovae, but this biostratigraphic scheme is problematic even regionally in Europe (e.g., Michel et al. 2015). Given the uncertainties in the proposed correlations, we suggest a conservative approach, namely, that the Bromacker Quarry could be Wolfcampian or Leonardian in age (i.e., Asselian to Kungurian). As such we assign a large age range of 298.9 ± 0.2 to 272.3 ± 0.5 Ma to Thuringothyris.
Petrolacosaurus 304.3 303.6 All specimens of Petrolacosaurus are from the Rock Lake Shale of the Stanton Limestone in Kansas (Peabody 1952). This unit is assigned to the Streptognathus firmus conodont zone (Ritter 1995); in the 2012 GTS this zone is dated to between 304.3 and 303.7 ± 0.1 Ma.  (Vaughn 1955;Reisz et al. 1984). These units correlate with the Coleman Junction through lower Cleark Fork Group further south (Hentz, 1989) which are assigned an upper Sakmarian through Kungurian age (Wardlaw 2005). Based on the 2012 GTS, this gives an age range of 292 to 272.3 ± 0.5 Ma.