Moving towards a better understanding of iterative evolution: an example from the late Silurian Monograptidae (Graptolithina) of the Baltic Basin

Iterative evolution has proved a difficult evolutionary phenomenon to study and interpret. Inferences of causality vary from study to study and quantitatively based phylogenetic reconstruction has never been attempted. In an effort to better understand iterative evolution we employed stratocladistics, gap analysis and disparity analysis to study the case of the Monograptidae in the aftermath of the late Silurian Cyrtograptus lundgreni extinction event. Our combination of gap analytical and stratocladistic techniques allowed us to elucidate the evolutionary relationships between the studied taxa. Based on our stratocladistic results we recommend the generic reassignment of five monograptid taxa. The stratocladistic results, in conjunction with morphological disparity analysis suggest the presence of a persistent developmental potential for the emergence of iteratively evolving characters. This persistent potential appears to be limited by extrinsic ecological constraints, which would have relaxed in the aftermath of the C. lundgreni extinction event. Our findings indicate that iterative evolution in the late Silurian Monograptidae is a product of the interaction of both intrinsic and extrinsic constraints on the acquisition of the iteratively evolving character, with the exact causality being dependent on the particular character.

morphological disparity analysis suggest the presence of a persistent developmental potential for the emergence of iteratively evolving characters. This persistent potential appears to be limited by extrinsic ecological constraints, which would have relaxed in the aftermath of the C. lundgreni extinction event. Our findings indicate that iterative evolution in the late Silurian Monograptidae is a product of the interaction of both intrinsic and extrinsic constraints on the acquisition of the iteratively evolving character, with the exact causality being dependent on the particular character.
I T E R A T I V E evolution is a uniquely informative case of evolutionary convergence. Unlike traditional convergence, iterative homoplasies occur not between separate, distinct clades, but from within the same clade. Much as in evolutionary parallelism, iteratively evolving morphologies generally share the same or very similar developmental paths to one another (Urbanek 1998). The distinction defining iteratively convergent taxa is their placement in time. Iterative evolution specifically refers to the repeated, independent development of homoplasies within the same clade at different points in time. In the most extreme cases, iterative evolution can occur with convergent taxa arising from the same long-lived parent species. This 'Groundhog Day'-like phenomenon of repeated births and deaths of the same morphological traits has been recorded in a number of different taxonomic groups including foraminifera (Coxall et al. 2007), graptolites (Urbanek et al. 2012), conodonts (Dzik 2005), ammonoids (Yacobucci 2015), turtles (Parham & Pyenson 2010) and canids (Van Valkenburgh 1991). Iterative evolution is a particularly tricky phenomenon to represent phylogenetically, as temporally convergent morphologies are difficult to reconcile using traditional morphological phylogenetics. As a result, the focus of most studies concerning iterative evolution has been limited to causality, it being attributed on a case-by-case basis both to intrinsic (behavioural, developmental or genetic) constraints (Van Valkenburgh 1991; Urbanek 1998) and to extrinsic (ecological and environmental) constraints (Jacobs et al. 1994;Parham & Pyenson 2010;Schmidt et al. 2013) as well as to the intrinsic facilitation of iterative morphologies (Xie et al. 2019).
As described by Urbanek (1998), iterative evolution is likely to be closely associated with oligophyly (the severe restriction and minimization of phylogenetic lineages and their subsequent evolutionary radiation) and, by extension, faunal turnover. Faunal turnover events are likely to incur major morphological shifts, facilitated by changes in developmental timing (Vrba 2005). If those turnovers are accompanied either by significant enough intrinsic or extrinsic constraints, they may result in the iterative evolution of specific morphologies. Disparity analyses have often been used to aid in the interpretation of these turnover recoveries, including in the case of early Silurian graptolites (Bapst et al. 2012). Particular attention is paid to the combined patterns of species diversification and morphological disparity and the implications they have on ecospace occupation (Foote 1993). In the case of heavily constrained morphology, as would be expected with iterative evolution (Rickards & Wright 2003;Parham & Pyenson 2010), species diversity should decouple from morphological disparity and increase at a much higher rate (Erwin 2007). Heavily constrained taxa should also exhibit more extreme clustering (Foote 1990). By examining the changes in disparity patterns, we may glean a greater understanding of how and why iterative evolution occurs in a given clade. However, we are still left with the question of how best to reconstruct the evolution of iteratively evolving clades.
The reconstruction of phylogenies for iteratively evolving taxa is a uniquely challenging task. There is ample opportunity for error given that many iteratively convergent taxa are best distinguished temporally, particularly those which arise from the same central lineage or species. This is the case even with convergent evolution at molecular level, which has been seen to induce false apparent phylogenetic relation (Paterson et al. 2010). This has meant that phylogenetic reconstructions of iteratively evolving taxa have had to either be based on qualitative expert opinion (Urbanek et al. 2012) or corroborated with genetic data from living clade members (Parham & Pyenson 2010). One way of resolving this issue is to employ the use of stratocladistics. Proposed by Fisher (1992), stratocladistics adds a stratigraphic character to maximum parsimony calculations, and has been shown to be more effective at producing accurate topologies than traditional maximum parsimony (Fox et al. 1999). Here stratigraphic distance between related taxa is interpreted as stratigraphic parsimony debt, and is included in the overall calculation of parsimony debt. In other words, the stratigraphic character is effectively coded as an irreversible ordered character, thus preventing older taxa from being interpreted as descendants of younger taxa. It is important to establish a framework for defining the assignment of the stratigraphic character both in reference to the stratigraphic record and with respect to the character state assignment of the involved taxa. Wagner (1995) established a system of determining first and last occurrences for fossil taxa which is used as the basis for stratigraphic ranges in stratocladistics. The definition of the stratigraphic character is harder to pin down, as changes in the number (and therefore resolution) of states should alter the effect of the character on the analysis. While the importance of the value of the stratigraphic character in stratocladistics has been explored from a theoretical standpoint (Clyde & Fisher 1997;Smith 2000;Fisher et al. 2002), the precedent set by stratocladistic analyses in practice (e.g. Leighton & Maples 2002;Pardo et al. 2008;Missiaen et al. 2013) has been not to acknowledge the assignment of the resolution of the stratigraphic character. This implies that the resolution of the stratigraphic character can be assigned arbitrarily, regardless of the quality or completeness of the fossil record. For most situations (especially at regional or global scales) the smallest logical unit should be a stratigraphic interval (e.g. period, stage, epoch, biozone, etc.) giving some basis for character state assignment (Alroy 2002). Further work has been done to better define stratigraphic resolution and character state assignment in stratocladistics by likelihood methods, but so far these have required more data than is available for most studies, particularly those on larger scales (Huelsenbeck & Rannala 2000). To improve the robustness of stratocladistic analyses, we used stratigraphic gap analysis (Paul 1982;Foote et al. 2007) for the definition of stratigraphic character resolution. The addition of this technique provides a quantitative basis for the resolution of the stratigraphic character by evaluating probabilities of given gap sizes between ancestors and descendants.
The case of Silurian monograptids provides a good test for the ability of stratocladistics and disparity analysis to represent iterative evolution. The Homerian and Gorstian ages of the Silurian are punctuated by the Mulde bioevent, a major ecological disturbance associated with significant glaciation (Jepsson & Calner 2002) which had a catastrophic impact on the diversity of shelfal fauna (Koren' & Urbanek 1994;Steeman et al. 2016;Venckut_ e-Aleksien_ e et al. 2016;Spiridonov et al. 2017a). Chief among the affected taxa were the monograptid graptolites, which drop to a single diminutive taxon, Pristiograptus dubius parvus Ulst, in the late Homerian immediately following the C. lundgreni graptolite biozone. This is generally referred to as the lundgreni extinction event (Koren' & Urbanek 1994). The monograptids which diversify in the aftermath of the extinction bear some of the same morphologies as those present before the lundgreni event (Urbanek 1998), despite not sharing those same morphologies with the proposed Pristiograptus dubius (Suess) central stem lineage (Rickards & Wright 2003), which links the two intervals together. Most obvious of these homoplasies are the incised apertural lips and (in some cases) enlarged bilateral processes of the pre-extinction Pristiograptus sardous (Gortani) and P. lodenicensis P ribyl, and of the post-extinction P. labiatus Urbanek, as well as the genera Colonograptus, Lobograptus and Neocolonograptus. This repetition of morphological traits in species originating from the same stem lineage is a clear case of iterative evolution, and has been recognized as such (Urbanek et al. 2012).
The purpose of this case study is to quantitively test and expand upon the findings of Urbanek et al. (2012) with repeatable methodology and a larger suite of taxa, integrating information on the completeness of the fossil record of given group of graptolites, and explore the underlying processes which give rise to iterative evolution. In this way we intend to produce a theoretical and methodological framework through which more potential cases of iterative evolution, in both deep time and the recent past, can be examined.

MATERIAL
The Monograptidae are found in particularly high abundance and diversity in material from the Baltic Basin of the Eastern European Platform (Fig. 1). The boreholes in Poland (Urbanek 1966;Teller 1969;Porȩbska et al. 2004) and Lithuania (Pa skevi cius 1979; Radzevi cius & Pa skevi cius 2000; Radzevi cius 2013), among others, exhibit not only a high diversity of species (~60 currently recognized) (Fig. 2), but also extremely high stratigraphic resolution (Spiridonov 2017;Spiridonov et al. 2017a). The high resolution of the stratigraphy in conjunction with the high degree of species representation makes the late Silurian Monograptidae of this region a prime candidate for stratocladistic analysis and reduces the likelihood of highly asymmetrical rates of morphological evolution confounding the results of the phylogenetic analysis (Norell & Novacek 1997).
Gap probability was estimated on the graptolite occurrences from the mid-Wenlock to lower Ludlow, more specifically between the upper part of the lundgreni and the end of the nilssoni biozones of the Vidukl_ e-61 ( Fig. 1) core section (Martma et al. 2005;Radzevi cius et al. 2014) in the depth interval between 1319.9 and 1244.7 m, a table of productive samples is available in Whittingham et al. (2020, supp. 2).

Gap analysis
In the interest of improving the robustness of our analysis, and of the stratocladistic practice as a whole, we herein justify our stratigraphic resolution on the basis of gap analysis.
Stratigraphic gap analysis (Paul 1982;Foote et al. 2007), is a method which estimates probability of occurrence of separate taxa in a temporal succession based on frequency of occurrence of productive samples between the first and the last occurrence events (thus omitting both). The average probability of occurrence of a species in a given clade, and at the given site could be estimated in following way (Foote et al. 2007): Here, Np i is the number of productive samples with a given species in a section), Nt i is the total number of samples inside the visible range of a species in a section, n is the number of species with the Np > 2 in a given section.
The estimated probability was used in the estimation of distribution of probable temporal gap sizes between ancestors and descendants ('ghost ranges ';Benton et al. 1999;Wills 1999)  uniformity in probability of occurrence of species. This is not an unfair assumption, since the Vidukl_ e-61 section represents deep water offshore environments which should experience much less perturbation in sedimentation rates than their nearshore counterparts.
To reveal the distribution of gap sizes between ancestors and descendants we used the most conservative stratigraphic scenario, in which ancestors are replaced phyletically by descendants; this could happen either gradually or by punctuated change (Fig. 3). In this case, there is no temporal overlap which will inevitably increase probability of co-occurrence of the two taxa, thus decreasing the stratigraphic debt on phylogeny. If the average probability of occurrence is approximately constant, then the probability of gaps of increasing size decreases exponentially (Paul 1982;Foote et al. 2007). In the case of ancestor-descendant transition in the illustrated scenario, two separate gaps are evaluated: the gap between the true last appearance datum (LAD) of the ancestor and its apparent LAD, and the gap between the true first appearance datum (FAD) and apparent FAD of the descendant. This will lead to the convolution of two separate exponential functions. First, we can estimate the decay constant of the exponential distribution from the average probability of occurrence per sample p, and from average sampling interval k: This constant could then be used for calculating distribution of gaps between ancestors and descendants. There are simple analytical solutions for cases where ancestor and descendant decay constants are equal to each other k ¼ k 1 ¼ k 2 . Simplifying eq. 3.1 of Akkouchi (2008), for the case of the sum of two exponentially distributed random variables with identical decay constants the probability density function is: In the case where ancestor and descendant decay constants differ k 1 6 ¼ k 2 , another function is derived. Simplifying eq. 2.3 of Akkouchi (2008), for the case of the sum of two exponentially distributed random variables with differing decay constants, the probability density function is: In both cases, due to the effects of the central limit theorem, the probability density function shows a modal (hump backed) shape. The described functions (Eqs 3, 4) could be used for determining the probable gap sizes, as well as the probability of a given gap between ancestor and descendant.
In addition to the described analytical approaches, a more flexible direct stochastic simulation approach was used in this study. The numerical model treats lineages as being composed of a series of samples (n 1 , n 2 ) which are occupied according to their measured uniform probabilities (p 1 , p 2 ). The simulation then treats distances between samples as random, truncated at '0', Gaussian variables of the average sampling interval k, and the standard deviation SD, both measured in kyr. In this way, numerical simulation accounts for convolution of the ancestor gap function, descendant gap function, and sampling distance variability function. It could account for equal as well as unequal probabilities of the occurrence of ancestors and descendants. After the simulated data are obtained, the histogram of gap durations could be evaluated in statistical terms. The simulation and analytical equation codes were written and executed in the R computational environment (R Core Team 2015) and are available in Whittingham et al. (2020, supp. 1).

Stratocladistics
In accordance with Alroy (2002), the highest available degree of biozonal resolution was used as the basis of the stratigraphic character, which is highly supported in the gap analysis (see Results). In the case of the Wenlock and Ludlow of the Baltic Basin, the highest-resolution biozonation is that defined in Urbanek et al. taxa not available from Lithuania (particularly P. sardous), no further increase of stratigraphic resolution was deemed appropriate or necessary beyond the biozone level. StrataPhy results were interpreted as anagenic in those cases where taxa fulfill Grantham's (2004) criteria for anagenetic interpretation. In addition, we interpret only those taxa which could be justifiably placed within the same species as anagenetic, thus allowing any anagenic lineage to be congruent with the lineage species concept (de Queiroz 1998(de Queiroz , 2007 Character data was acquired primarily from physical specimens in the collections of the Vilnius University Department of Geology and Mineralogy, with the remainder of the data being acquired from systematic measurements and images from the relevant systematic literature (Urbanek 1963(Urbanek , 1966Koren' 1992;Maletz 1999;Urbanek et al. 2012;Storch et al., 2014Storch et al., , 2016. As many morphological characters as possible were included, with particular care taken to include all traits relevant to taxonomic diagnosis. Care was taken to ensure that each character was unambiguously defined to avoid subjectivity in interpretation.

Bilateral process extent was defined in accordance with
Urbanek et al. (2012), with the addition of spinosity as the most extreme character state. For stratocladistic analysis the majority of morphological characters were unweighted and all were unordered. Many dimensional characters were represented by the same measurement at the three separate points on the rhabdosome, these were weighted at one-third value to reflect the interdependence of the measurements used. The character table and weight set were run through the stratocladistics software StrataPhy (Marcot & Fox 2008). This program adds stratigraphic data to traditional maximum parsimony analysis by coding stratigraphic gaps between taxa as additional parsimony debt, thereby reducing the likelihood of data being confounded by temporal convergence. Maximum parsimony analysis was conducted via heuristic search with TBR branch swapping, a maximum tree number of 1000, and 10 search repetitions.

Principal components analysis and morphological analysis
Disparity analysis was based on a binary character matrix derived from that used in the stratocladistic analysis, with the stratigraphic character removed (Whittingham et al. 2020, supp. 3). Morphological phylogenetic datasets also highly compatible with disparity analyses (Deline et al. 2018), requiring minimal adjustment beyond the simplifying of multi-state characters. Multi-state characters were split into binaries coded either as successive binaries in the case of ordered characters (where 'simpler' character states are implied to be occupied by taxa which occupy more 'complex' states), or, in the case of unordered characters, simple presence/absence binaries where each character state was treated as its own unique character. Care was taken to ensure the adequate representation of homoplasic and autapomorphic characters as outlined by Gerber (2018). The matrix was plugged into a principal component analysis (PCA) in the program PAST v3.09 (Hammer et al. 2001). The two principal components which surpassed the broken stick model of significance were then plotted against one another. Nine sub-samples were constructed to determine the relationships between the species on either side of the lundgreni extinction, and to determine the contribution of the Colonograptus and Lobograptus groups to the patterns of disparity. Those sub-samples were selected for: (1) only species which predate the lundgreni extinction; (2) only species which predate the lundgreni extinction with the exclusion of the Monograptus outgroup; (3) all species excluding members of the genera Colonograptus and Saetograptus (hereinafter regarded as the Colonograptus group); (4) all species excluding members of the genera Lobograptus, Cucullograptus, Crinitograptus, Bohemograptus and Neodiversograptus, as well as P. idoneus (hereinafter regarded as the Lobograptus group); (5) all species excluding the members of the Colonograptus and Lobograptus groups; (6) all species which appear after the lundgreni extinction; (7) all species which occur after the lundgreni extinction with the exclusion of the Lobograptus group; (8) all species which appear after the lundgreni extinction with the exclusion of members of the Colonograptus group; and (9) all species which appear after the lundgreni extinction excluding the members of both of the Colonograptus and Lobograptus groups.
Disparity analyses are a common mode of analysis of morphological diversity. Conventional ordination methods used for this purpose are principal coordinates analy- However, it is impossible to determine the contributions of individual characters in PCoA and NMDS to the ordinations produced. As these loadings are crucial for comparing the contributions of different characters to differences in species diversity in separate time intervals, it was necessary to use PCA for the purposes of this study. For the purposes of determining the validity of the use of PCA for the purposes of ordination, verifying the validity of our restructuring of multistate characters, and observing ordination with maximized explanation of variance, NMDS was conducted on a generalized Euclidean distance matrix based on the same primary dataset to compare results.
NMDS scores were additionally used to calculate average disparity through time. Even though the distributional gradients in the NMDS analysis do not preserve their original configuration in the original state space, there is no ad hoc reason not to compare their final configuration in the NMDS space using Euclidean-like measures (i.e. sums of variance). The scatter of points in the NMDS solution maximizes the optimal gradient of differences between points. Therefore, the distances between observations, even though non-metric, are meaningful, and monotonic increase in separation implies monotonically greater difference between two compared points in a new gradient.
Following the methodology of Marx & Fordyce (2015), average disparity was estimated as the sum of variances per biozone. Variance measures were chosen over rangebased measures for this estimate, as the former are more independent of sample size biases (Wills et al. 1994) and mitigation of sample size biases via rarefaction as described by Foote (1992) was considered inappropriate given the small number of species available in any given time bin.
Principal component scores from the full matrix were then tested for fit to a logistic regression with non-zero slope to determine the significance of change in morphospace occupation between the taxa before and after the lundgreni extinction. Additionally, for the purpose of illumination of the evolutionary mechanisms forcing morphological disparity, the contributions of each character were evaluated for each sub-sample. This was accomplished by producing a secondary matrix of characters which passed an arbitrary threshold loading weight of 0.2 for each significant component of each sub-sample of the dataset, as well as the full dataset. Those sub-samples with multiple significant components had their character matrices transformed such that all characters with a loading weight which passed the 0.2 threshold in any significant principal component of each dataset was given a value of 1. A Jaccard's similarity matrix was then produced in PAST v3.09 (Hammer et al. 2001) to determine the frequency of shared most valued characters (MVCs). While this similarity index is inadvisable for the direct study of discrete morphological data as it only considers shared positive values (Hopkins & St John 2018), in our case the only characters in the matrix of significance are given positive values, and inclusion of 0-valued characters would increase the influence of characters already considered to be uninformative.

Gap analysis
Based on the analysis of 15 graptolite species/subspecies which had more than two productive samples in the Vidukl_ e-61 core section it was determined that the mean probability of occurrence equals to p = 0.15. The samples in the studied interval were placed on average every 20 cm AE 10 cm. Estimated total duration of the interval studied for the frequency of gaps spans c. 3 million years (Melchin et al. 2012). Therefore, mean sampling intervals (k) translates into c. 8000 AE 4000 years. Mean duration of a graptolite lineage in the studied stratigraphic window is 58 samples or 464 kyr.
Using these parameters and applying them to Eq. 3, it was determined that the modal gap in the phyletic succession of ancestor and descendant should be~49 kyr, which is much smaller than the mean duration of a graptolite zone, and the probability of a gap that equals or exceeds 500 kyr is p = 0.00043 (Fig. 4A). The stochastic simulation showed very similar results. Mean gap size 99 kyr, median = 83 kyr, with the modal interval between 50 and 100 kyr, and a probability of having gap equal to or larger than 500 kyr of p = 0.00057 (Fig. 4B). The probable gap sizes are larger in the simulation than they appear in the analysis of analytical formula because of the variability of sampling intervals, which are explicitly incorporated in the stochastic model.

Stratocladistic phylogenies
The stratocladistic analysis produced a single island of 240 unique trees, of which 17 exhibited unique topologies. The remaining trees were distinguished from those unique topologies only in the assignment of different taxa as anagenic internal nodes. The variable assignment of internal nodes appears to be a product of taxa not being designated as ancestors, as is mentioned by Marcot & Fox (2008).
The phylogram produced by stratocladistics (Fig. 5) shows the retention of an anagenic stem lineage spanning across the Mulde event from P. dubius to P. postfrequens Urbanek et al., with no other lineages spanning the extinction. This confirms, with different methodology and an expanded taxonomic dataset, the findings of Rickards & Wright (2003) and Urbanek et al. (2012). This phylogram also shows multiple iterative originations of the bilateral process character (clades originating from P. praelodenicensis Urbanek et al., C. praedeubeli (Jaeger), L. sherrardae (Sherwin)), and of the connected thecal lip incision character (clades originating from P. sardous, P. praelodenicensis, C. praedeubeli, L. sherrardae, P. tumescens (Wood)).
The stratocladistic results also indicate a number of changes to the phyletic relationships of the studied genera. P. idoneus Koren' is shown to be polyphyletic to the remainder of Pristiograptus, showing greater affinity for the Lobograptus group. Branching off of the Pristiograptus stem lineage as early as the parvus biozone, geniculate taxa of the Ludlow (U. uncinatus (Tullberg), Heisograptus micropoma (Jaekel), Ps. dalejensis (Bou cek), Ps. antiqua Storch et al.) form a monophyletic clade, this clade is the only one to suggest a preservation gap of any size greater than a single biozone, with no known members throughout the majority of the Homerian. The genus Colonograptus is shown to be paraphyletic, giving rise to the genus Saetograptus, while the P. dubius stem lineage of Rickards & Wright (2003) is retained in its entirety. The genus Lobograptus is shown to be polyphyletic, with L. ciriffer Urbanek outlying in a clade with Cucullograptus, and with L. claudiae Koren' and L. sherrardae being the stem of the Lobograptus group as a whole. Neodiversograptus is also shown to be polyphyletic, with N. nilssoni (Barrande) and N. beklemishevi Urbanek occupying separate clades.

Principal component analysis and disparity results
PCA of the full dataset showed two principal components of significance under a broken stick model, accounting for 37.5% of the variance in the character data. The results of the full PCA indicate that the post-extinction species occupy morphospace which includes that of all but one of the pre-extinction species (the exception being Monograptus antennularius), as well as new sets of morphologies displayed by the Colonograptus and Lobograptus groups (Fig. 6A). The new morphospace is not evenly occupied, instead displaying a moderate degree of clustering distinguishing three separate loci. The members of the Colonograptus group are distinguished in terms of morphospace occupation by the most influential principal component, which explains 26.5% of the total variance, while the members of the Lobograptus group are distinguished in terms of the second most influential principal component, which explains 11% of the total variance. When tested for fit to a logistic regression with a nonzero slope between the species groups from two intervals, the difference in occupation of the morphospace on the first principal component was significant (p = 0.016), while the difference of occupation of taxa between the two intervals on the second principal component was also revealed to be significant (p < 0.001). The same   morphospatial shift via the Colonograptus and Lobograptus groups in the post-extinction interval was also seen in the NMDS results (Fig. 6B), with the slope between the two intervals also shown to be non-zero (p = 0.006 for Component 1 and p < 0.001 for Component 2). Those two NMDS components were highly explanatory, producing a Shepard plot with a stress of only 0.09. The NMDS and PCA results from the full dataset showed a high degree of correlation (PC1 vs NMDS1: r = 0.997, permutation p = 0.0001; PC2 vs NMDS2: r = 0.915, permutation p = 0.0001). Subsequent PCA of the sub-sampled data which excluded the Colonograptus and Lobograptus groups showed no change in morphospace occupation in either of the two most influential principal components (p = 0.57 for PC1 and p = 0.42 for PC2) between the pre-and post-lundgreni species (Fig. 7A), indicating that the bulk of the morphological change observed between two intervals can be explained by those two genera. This is also backed up by the NMDS results (Fig. 7B) P. d frequens .

Post-extinction species of Pristiograptus
Colonograptus group

Pre-extinction species of outgroup Monograptus
Post-extinction -group U. uncinatus H. micropoma F I G 6 . A, PCA plot for the full set of studied taxa. B, NMDS plot for the full set of studied taxa. between the intervals (p = 0.20 for Component 1 and p = 0.15 for Component 2). Those two NMDS components were moderately explanatory, producing a Shepard plot with a stress of 0.17. The NMDS and PCA results from the pruned dataset showed a high degree of correlation (PC1 vs NMDS1: r = 0.938, permutation p = 0.0001; PC2 vs NMDS2: r = 0.936, permutation p = 0.0001).
When analysed for average disparity (sums of variances) per biozone, the disparity (sums of variances) and total diversity curves are coupled (Fig. 8). Total diversity and disparity continue to increase through the Homerian and Gorstian, seeing a dramatic rise in the aftermath of the lundgreni extinction event. This is displayed by the sums of variances rising from a mean of 0.003 in the pre-extinction interval to a mean of 0.017 in the postextinction interval.
When analysed individually, the inclusive datasets of each interval showed distinctly different patterns, the prelundgreni species occupy a similar distribution to that which they exhibit in the sub-sample which spans both intervals and excludes the Colonograptus group, with distinction in morphospace shown primarily in the outgroup taxa Monograptus riccartonensis, M. antennularius and Testograptus testis (Fig. 9A). In contrast, the post-lundgreni species showed a clear separation in morphology between the members of the Colonograptus and Lobograptus groups and the remainder of the species in the interval (Fig. 9B). When those groups are removed from the post-lundgreni dataset, the distribution of species within the occupied morphospace follows a distinctly different pattern which bears closer resemblance to the distribution they exhibit in the sub-sample which spans both intervals and excludes the Colonograptus and Lobograptus groups (Fig. 9C), as is also seen in the pre-lundgreni species of Pristiograptus. Jaccard similarity comparison of the most valued characters between all datasets showed the closest similarity between the full dataset and the inclusive post-lundgreni dataset. The subsamples which excluded both groups are both most distinct from the full dataset, and also highly distinct between themselves (Table 1). Comparison between the most valued characters of the pre-and post-lundgreni sub-samples which excluded the Colonograptus group, Lobograptus group, and Monograptus outgroup show that the two groups share no MVCs, and only share Jaccard similarities of 0.26 (for the post-lundgreni sub-sample) and 0.25 (for the pre-lundgreni sub-sample) with the subsample which excludes both the Colonograptus and Lobograptus groups but encompasses both intervals. By comparison, the two samples which include the Colonograptus and Lobograptus groups share a Jaccard similarity index value of 0.67. is very unlikely. It should be kept in mind that the palaeontological examination of each new section works as an independent trial of the hypotheses on species distributions. Therefore, each new search diminishes probability of finding an artefactual gap between ancestor and descendant. For example, the probability of finding gap of a given size (e.g. 0.5 myr) or larger, assuming that we have five geological sections with fossil record quality and the rate and variability of sampling as in Vidukl_ e-61 core, will be p = 6 9 10 À17 . In the general case, the probability of observing gap of duration T ≥ t could be estimated using following expression:

DISCUSSION
Here, j is the number of examined sections, and the terms in brackets are the key parameters of the simulation for each stratigraphic section i. It could be seen that even having low sampling probabilities (p 1 , p 2 ) and large sampling interval (k), the overall probability will be quickly diminished by increasing the number of examined sections (j). This is a probable reason why graptolites, which are known from many sections from each time interval, are so robust in their consistency of temporal distributions. Therefore, the temporal succession of graptolite species which is observed in the Wenlock to Ludlow time interval should be very informative for the reconstruction of phylogenies, since gaps (or 'ghost ranges') that are much larger than a single standard graptolite zone are highly unlikely given the known properties of their fossil record.
The gap analytical results give us a very high degree of confidence in the stratocladistic results, as only one stratigraphic gap in the ingroup clade exceeds a length of one biozone. That long gap being that which spans the majority of the Homerian between P. d. parvus and the U. uncinatus -H. micropoma clade. The likelihood of this gap being a product of sampling bias is negligible based on the results of our gap analysis, indicating that there

Legend:
F I G . 9 . Principal component plots of the two most explanatory components A, full set of pre-extinction taxa. B, full set of postextinction taxa. C, set of post-extinction taxa which excludes the Colonograptus and Lobograptus groups.
should be other factors influencing this anomalous result. The most immediately obvious candidate for the cause of this gap is geography. As our analysis focuses solely on taxa which occur in the Baltic Basin, there is a distinct possibility that the 'ghost' ancestor of the U. uncinatus -H. micropoma clade did not disperse into the Baltic Basin. This is a much more likely explanation than Urbanek's (1998) invocation of the 'Lazarus Effect' and disagrees with Lenz's (1994) proposed alignment of Ps. dalejensis with Colonograptus. Our results follow more closely with Rickards & Wright's (2002) interpretation of the evolution of U. uncinatus, but with the species originating from the P. dubius stem lineage instead of M. flemingii. The other most likely explanation would be extreme asymmetry in the rate of evolution of the hooded thecae. Evidence from bryozoans, considered to be similar to graptolites both in development and ecology (Urbanek 2003), suggests that morphological change may occur on timescales shorter than are available even in an exceptionally high-resolution fossil record (Cheetham 1986), though the trend in graptolites appears to be more towards a gradual progression of morphologies (Urbanek 1966;Maletz et al. 2019). Both of these explanations highlight the pitfalls associated with the use of stratocladistics. Evolutionary refugia (sensu Rickards & Wright 2002) and asymmetrical rates of evolution can bias stratocladistic analyses and result in false topologies or falsely assigned ancestor-descendant relationships, as outlined by Benton et al. (1999). However, these are factors that also bias both maximum parsimony (Maddison 2006) and molecular phylogenetics (Goldberg et al. 2011). As such, we consider a stratocladistic method to be more appropriate for reconstructing clades experiencing iterative evolution than the use of traditional maximum parsimony, as stratocladistics successfully differentiates between taxa which are temporally disparate but morphologically convergent.
The reconstruction of the tree with temporally convergent taxa kept as separate clades shows the utility of stratocladistics as a means of reconstructing iterative evolution. The stratocladistic results show exact support for the phylogram of Urbanek et al. (2012), with iteratively evolving bilaterally lobate taxa stemming from the central P. dubius stem lineage. The support of an anagenic stem P. dubius lineage also means that the iterative evolution described by Rickards & Wright (2003) within that lineage is intraspecific, with all members falling under the species P. dubius under a lineage species model (sensu de Queiroz 1998de Queiroz , 2007. This follows even the most stringent and restrictive interpretations of anagenesis in morphospecies (Strotz & Allen 2013). And as virtually all Baltic Basin species of Pristiograptus are included in this study, it is unlikely that a lack of completeness would have biased this result. The use of a lineage-based species definition also means that all previous subspecies of P. dubius not included in the reconstructed stem lineage (e.g. P. labiatus, P. magnus) should be considered to be their own species. The persistence of long-ranging lineages also necessitates the interpretation of direct ancestor-descendent relationships and, by extension, the use of T A B L E 1 . Jaccard similarity matrix of most valued characters (MVCs) for each subsample of the studied monograptid taxa. Subsamples are numbered as follows: 1, only species which predate the lundgreni extinction; 2, only species which predate the lundgreni extinction with the exclusion of the Monograptus outgroup; 3, all species excluding members of the genera Colonograptus and Saetograptus (i.e. the Colonograptus group); 4, all species excluding members of the genera Lobograptus, Cucullograptus, Crinitograptus, Bohemograptus and Neodiversograptus, as well as P. idoneus (i.e. the Lobograptus group); 5, all species excluding the members of the Colonograptus and Lobograptus groups; 6, all species which appear after the lundgreni extinction; 7, all species which occur after the lundgreni extinction with the exclusion of the Lobograptus group; 8, all species which appear after the lundgreni extinction with the exclusion of members of the Colonograptus group; 9, all species which appear after the lundgreni extinction excluding the members of both of the Colonograptus and Lobograptus groups. The absence of common MVCs of the pre-and post-extinction groups which exclude the Colonograptus group, Lobograptus group, and Monograptus outgroup (subsamples 2 and 9) indicates a major shift in the direction of evolution between the two intervals.
a model of budding cladogenesis. It should not be particularly surprising to find these ancestor-descendent pairs in the Baltic Basin graptolite data, as the high temporal resolution and preservation rate greatly increase the likelihood of finding direct ancestors (Foote 1996). Bifurcating cladogenesis is inferred only for cases where sister taxa originate simultaneously, or where multiple branching events occur in short succession, as ancestry cannot be reliably determined (Wills 1999). Our phylogram also suggests that a number of Wenlock and Ludlow monograptid genera (Neodiversograptus, Monograptus, Lobograptus) are polyphyletic. These polyphylies may be remedied by the reassignment of the aberrant taxa 'Monograptus' insperatus, Lobograptus? sherrardae, Lobograptus? claudiae, Lobograptus progenitor and Lobograptus cirrifer, to new genera. There is support for generic reassignment in the morphology of these taxa. 'Monograptus' insperatus fits in with the description of Colonograptus (P ribyl, 1943) in having large bilateral processes and a conical sicula. Lobograptus? claudiae and Lobograptus? sherrardae are both remarkably robust compared to Lobograptus spp., the broader thecae and shorter thecal interspaces align them with Bohemograptus per Urbanek's (1970) revision of the genus. Lobograptus cirrifer aligns with Cucullograptus based on the overlapping asymmetry of its bilateral processes. This is a characteristic not shared with any Lobograptus spp. and is present in every Cucullograptus spp. Lobograptus progenitor fits the diagnosis of Neodiversograptus (Urbanek, 1963) in both dimensions and thecal apertural morphology. The only major discrepancy is in the lack of a dorsal apertural spine on the sicula. Our recommendation for the reassignment of N. progenitor would, therefore, make it an aberrant member of the genus, but would solve the polyphyly in Neodiversograptus and would align the species with more morphologically alike taxa. The reassignments of L.? sherrardae and L.? claudiae would extend the range of Bohemograptus back into the praedeubeli biozone, placing Bohemograptus and P. idoneus at the base of the Lobograptus clade. This confirms the speculations of Koren ' & Urbanek (1994), though not all of their phylogenetic interpretations are supported. Both their placement of C. deubeli as a sister to the remainder of the Colonograptus clade and their placement of 'M'. insperatus outside of the clade are refuted by the stratocladistic results. In service of addressing both the issues of polyphyly and taxonomic consistency, we recommend the following generic reassignments: Colonograptus insperatus (formerly 'Monograptus'), Bohemograptus claudiae (formerly Lobograptus?), Bohemograptus sherrardae (formerly Lobograptus?), Neodiversograptus progenitor (formerly Lobograptus), and Cucullograptus cirrifer (formerly Lobograptus). As the purpose of this paper is not to revise the taxonomy, we still refer to the above taxa by their original generic assignments, but recommend reassignment pending more thorough morphological assessment.
The reconstructed evolutionary succession of the Colonograptus and Lobograptus clades shows a distinct transition between astogenetic modes. The growth of these species is adherent to Urbanek & Uchmanski's (1990) interpretation of a morphogen gradient for the determination of changes in character expression with astogeny. For example, in each of C. praedeubeli, C. deubeli, C. ludensis and C. roemeri, thecal morphology is clearly split between astogenetic 'zones of changes' and 'zones of repetition' (Urbanek 2003), the former being decreasingly lobate and the latter expressing no thecal ornamentation. This is exaggerated in C. colonus and S. chimaera, which show ontogenetic hypermorphy in the form of bilateral spines in the earliest thecae, but still reduce to a consistent morphology unadorned by ornamentation later in astogeny. In both Colonograptus and Saetograptus there are also cases of fixation of the proximal morphology (C. gerhardi and C. insperatus in the former and S. leintwardinensis and S. fritschi in the latter), resulting in monomorphic thecae throughout astogeny. The same fixation occurs in the evolution of the Lobograptus clade, with Bohemograptus and Crinitograptus maintaining bimorphic thecae, and Lobograptus experiencing fixation of the proximal morphology. This fixation represents a change in developmental timing, and a rare case of astogenetic heterochrony in the vein of that seen in bryozoans (Anstey 1987;Pachut et al. 1991).
The Colonograptus and Lobograptus groups also represent the only significant morphological expansion between the two study intervals. The results of the disparity analysis of the full dataset clearly show the retention and expansion of morphospace from the pre-extinction interval to the post-extinction interval via the Colonograptus and Lobograptus groups. This fits the classical model of ecological and morphological expansion by surviving clades in the wake of an extinction event, wherein both disparity and diversity are seen to increase (Foote 1993). Expansion of diversity and disparity can be facilitated by the relaxation of competitive selection (Schluter 2000;Buress et al. 2018), as has been observed in semionotid fishes, another group which exhibits repeated evolution of specific morphologies (McCune 1990). Given the paucity of graptoloid taxa with similar morphologies in the early post-extinction interval, and the ecological significance of those morphologies (Rigby & Rickards 1989;Underwood 1993), it is logical to assume a relaxation of competition in this environment, allowing for the expansion of diversity of and disparity. This increase in morphological disparity is limited to the Colonograptus and Lobograptus groups, as illustrated by the morphospatial stasis observed both when those two groups are removed from analysis (Fig. 6). These results fit Koren ' & Urbanek's (1994) notion of morphological conservatism in their Pristiograptus group, with limited expansion via the Colonograptus and Lobograptus groups. This morphological conservatism in the remainder of the group gives another piece of evidence for the presence of a persistent stem lineage accompanied by iterative development of temporally convergent characters as described by Urbanek et al. (2012), as the pre-and post-extinction samples which excluded the Colonograptus and Lobograptus groups occupy a nearidentical morphospace. The differentiated clustering observed in the total and post-extinction species subsets supports the action of constraint on the Monograptidae as a whole to a limited set of general morphological layouts (Foote 1990). However, the coupling of average disparity and total diversity in the aftermath of the lundgreni event ( Fig. 9) indicates that the constraint was so not strong as to prevent morphological expansion (Erwin 2007), but rather funnel it in three distinct directions.
These general patterns of morphospace conservation and expansion are shared between both PCA and NMDS ordination methods, indicating comparable validity of the two methods in this particular case. The convergence of the results between PCA and NMDS from the same base data adds to the findings of Romano et al. (2017) that discrete and continuous character data from the same taxa will also converge to the same results. The high similarity of the results of both analyses verify the restructuring of ordered and unordered multistate characters for the PCA, as well as the use of NMDS results for total disparity through time.
Though morphology is generally conserved in the non-Colonograptus and Lobograptus group member species, the internal structure of that morphospace (i.e. what specifically controls the distinctions between taxa) does shift in the wake of the lundgreni event. This is shown in that a majority of characters which best describe the most significant ordination are not shared between each of the two intervals. Additionally, only a small number of the most valued characters in each sub-sample (2 pre-extinction and 0 post-extinction) are characters which are considered to be directly affected by changes in development (thecal shape characters sensu Urbanek 2003), with emphasis instead being on size characters. Frequent fluctuations in the direction and strength of morphological evolution are understood to be commonplace in modern (Siepielski et al. 2009) and fossil (Nehm 2005) biota, and variability is especially high for size characters in evolutionary static clades (Hunt 2007). As such it is unlikely that these characters would contribute significantly to the morphological expansion observed in the wake of the lundgreni extinction, and this is borne out in the morphological conservatism observed in the non-Colonograptus and Lobograptus group taxa, evidenced by the small number of novel morphological traits (eg. dorsally oriented siculae, contorted thecal lobes, spinose thecae) as compared to the diverse range of apomorphies seen in Colonograptus and Lobograptus. As even these heavily constrained taxa do not display complete stasis, our results concur with the view of intrinsic constraints as a non-factor in the maintenance of morphological stasis (Eldredge et al. 2005;Nehm 2005). By contrast, a majority of characters distinguishing the Colonograptus and Lobograptus groups from the remainder of the species are likely to be directly affected by changes in developmental timing, these primarily being the presence, extremity, and fixation of bilateral process and spines. These are the same characters associated with iterative evolution in the Monograptidae. As very little is understood about graptolite ecology outside of the functional morphology their general colonial structure (Rigby & Rickards 1989) and inferences made from their closest living relatives, the Pterobranchia (Maletz & Bates 2017), none of the studied characters could be inferred as being directly related to changes in nichespace occupation. The limited understanding of graptolite functional morphology potentially biases the interpretation of the causality of the iterative evolution in monograptids away from ecological forcing mechanisms. However, these findings suggest some change in the directionality of evolution, and by extension ecology, after the lundgreni event in conjunction with the increase in total morphospace occupation of the group after the lundgreni extinction event.
The rapid increase in diversity and disparity in the praedeubeli Biozone, and their subsequent expansion could be explained by extrinsic factors like increased possibilities of allopatric speciation events due to cyclic changes in the global sea level. Cyclic perturbations in the sea level gives rise to episodes of vicariance (when sea level is low) and dispersal (when it is high) speciation events in the oceans ( Moreover, the conservatism in developmentally indicative characters in the non-Colonograptus and Lobograptus groups falls in line with Gould & Lewontin's (1979) model of developmental constraint. The directional evolution of developmental characters in the Colonograptus and Lobograptus groups could provide a good case for developmental drive (Arthur 2001), which operates under the same conceptual framework. This, in conjunction with the patterns shown by the average disparity curve, suggests that the shift in morphospace in the aftermath of the lundgreni event is a singular breakthrough of developmental drive during C. praedeubeli Biozone.
The rapid increase in abundance and geographic expansion of species in the late Homerian probably enhanced the production and fixation of new morphological forms due to spatial sorting and drift in the absence of strong selection (sensu Spiridonov et al. 2015;Lehman & Miikkulainen 2015). While no new examples of the iterative bilobate morphology are recorded until P. labiatus, and later Neocolonograptus, appear in the upper Ludfordian after another major monograptid extinction event at the end of the Neocucullograptus kozlowskii graptolite biozone, the lip incisions associated with bilobate thecae are observed in P. tumescens, indicating the persistent intrinsic potential of the iterative morphology. It should be noted that it is possible that the development of the iteratively evolving bilobate thecae could also be extrinsically constrained by competitive ecological factors, given the high diversity of bilobate monograptids co-occurring with P. tumescens and members of the P. dubius stem lineage. The action of competitive constraint as a limiting factor would also help to explain why the expansion of morphology past the apparent barrier of the development of bilobate thecae occurs only in aftermath of a major extinction, where competitive restrictions should be more relaxed (Vermeij 1987;Dietl et al. 2004). The post-extinction relaxation in ecological constraint makes a strong candidate for the allowance of an expansion of morphology (Oyston et al. 2015), in our case a breakthrough in the development of bilobate thecae. It is then likely that causation of iterative evolution is rooted in the interaction of both intrinsic developmental drivers and extrinsic competitive constraints, at least in the case of the late Silurian Monograptidae.
It is important to note that without detailed knowledge of monograptid functional morphology it is hard to determine the exact effects of ecological constraint on morphological expansion. It is our view that iterative evolution is a phenomenon with causality highly reliant on the specific characters which define it. This is supported by the causes of iterative evolution events in ammonites (Jacobs et al. 1994) and turtles (Parham & Pyenson 2010) having been linked to extrinsic ecological factors, as the iteratively evolving characters in both cases are closely linked to the functional morphology and ecology of the studied taxa. In the same vein, the case of Xie et al.'s (2019) repeated losses of pelvic hindfins in stickleback shows intrinsic genetic facilitation of morphological change, in much the same way as we observe in late Silurian monograptids, but only in conjunction with the relaxation of extrinsic ecological constraints, indicating that an ecological factor (in our case suspected to be competitive interaction) may also be at play in the Silurian Monograptidae.

CONCLUSIONS
The combination of gap analysis-informed stratocladistics and disparity analysis indicates that developmental driving of morphological expansion and the action of ecological constraint were the key factors in the iterative evolution observed in the late Silurian Monograptidae. Developmental driving is supported by the paucity of developmentally relevant characters distinguishing species in the stem lineage producing the iterative characters, and the developmental relevance of those iterative characters. Ecological constraint, and the relaxation thereof, are indicated by the synchronous rapid increase of diversity and disparity in the praedeubeli biozone. These two factors represent the action of both intrinsic and extrinsic constraints affecting iterative evolution.
The completeness estimates show, under realistic conditions, that large gaps (larger than a standard graptolite zone) are highly unlikely for an analysed clade. Therefore, temporal distribution of graptolite taxa should be accounted for in searching the optimal topology of an evolutionary tree. The produced phylogeny and observed conservation of morphospace supports reconstruction of the group with the iterative evolution of taxa bearing bilobate thecae from the central P. dubius stem lineage. Based on that phylogeny we recommend the reassignment of five monograptid taxa to new genera (Co. insperatus, B. claudiae, B. sherrardae, N. progenitor, Cu. cirrifer) pending more detailed morphological analysis.