New insights into the reading of Paleozoic plant fossil record discontinuities

Studying the discontinuity patterns of Paleozoic vascular plants provides a global vision of these key events from the multivariate methods viewpoint. Non-metric multidimensional scaling, detrended correspondence analysis and cluster analysis have been employed together with a set of diversity and abundance measures and an evaluation of the geologic constraints from the plant fossil record data. The results reveal four clear significant discontinuities in terms of taxonomic composition and record representativeness during the early-middle Devonian, Devonian–Carboniferous, Mississippian–Pennsylvanian and early-late Permian. Due to the controversial character of the plant fossil record data and the effect of mass extinction events, the results can be explained in taxonomic turnover and ecological reorganisation terms which emphasise the crucial role of the geologic constrains in paleobiological inference.


Introduction
It is universally accepted that the Paleozoic is probably the most important geologic era in plant life history. During the Paleozoic, the first land plants appeared (Becker and Marin 2009;Steemans et al. 2009), the early evolutionary trends of the tracheophytes evolved (Niklas and Smocovitis 1983;Thomas and Spicer 1986), the mean key evolutionary innovations emerged, with the exception of the flower origin, and the most important extinction events happened which have doubtlessly channelled the plant history originating during the evolutionary process Thomas 1994, 1999;Foote 2000Foote , 2003Jablonski 2005). All these general aspects about the origin and evolution of plants during the Paleozoic have been classically studied from several methodological viewpoints of paleobiological inferences, generally a series of analyses of evolution rates and several abundance and percentage methods in the 1980s (Niklas et al. 1980Niklas 1982;Knoll 1984;Knoll and Niklas 1987).
Nevertheless in recent years, many of the previously accepted considerations on the paleobiological inferences and mass extinction events of vascular plants have led to considerable debate, especially concerning the possibility of studying plant evolution in accordance with the quality of the fossil record data (Benton et al. 2000;Holland 2000;Kidwell and Holland 2002;Crampton et al. 2003;Smith and McGowan 2005) and the consequences of the mean mass extinction events (Raup and Sepkoski 1984) on plant life (Wing 2004;McElwain and Punyasena 2007;McElwain et al. , 2009. In this way, and as a result of the record's nature, new methodological approaches are constantly being considered in inference through the fossil record in order to provide the best possible information about past life (Alroy 2000;Sumrall et al. 2001;Wilson 2003;Vermeij and Herbert 2004;Wang and Bush 2008). As the fossil record merely represents a sample of past life on this planet, paleontologists who seek answers to crucial questions in paleobiology, such as the measure of diversity and the effects of mass extinction events, have learned to deal with the imperfections of the fossil record by means of many different methodologies.
Traditionally, multivariate methods, as well as ordination and clustering analyses, have been used to infer geographical patterns and ecological vegetation gradients, and even today they are widely used (Gavilan 2005;Ezcurra et al. 2008;Collier 2009;Matthews et al. 2009;Yu and Ehrenfeld 2010), and paleobotanical sciences are not an exception. In the last few years, the paleofloristic and paleogeographical works of Cleal, based on multivariate analyses, have improved our knowledge of Paleozoic flora, particularly the Carboniferous plant fossil record (Cleal 2007(Cleal , 2008a(Cleal , 2008bHilton and Cleal 2007). Following this methodological philosophy, the aims of the present paper are to (1) produce a global assessment of Paleozoic plant life using multivariate statistics, (2) emphasise the heterogeneity of the plant fossil record over time and (3) understand the consequences of discontinuities in the history of Paleozoic plant life in terms of taxonomic composition by taking into account the geologic constraints in the plant fossil record interpretation.

Dataset and timescale
A dataset of plant fossil families was compiled from the several classic monographic studies that represent a welldocumented plant life history (Cleal 1993a(Cleal , 1993bStewart and Rothwell 1993;Taylor and Taylor 1993). This large fossil dataset has been recently contrasted to provide a global evolutionary view of all plant history. Detailed information about the several data sources, the type of fossil consulted and the methodological bias relating to paleobiological inference, e.g. sampling bias (Alroy 2003;Hunter and Donovan 2005), taphonomic and geologic constrains (Behrensmeyer and Kidwell 1985;Behrensmeyer et al. 2000;Caracuel et al. 2000), considerations related to the rock-record bias (Smith et al. 2001;Crampton et al. 2003), which are intrinsic in the paleobiological inference based on the multivariate methods presented below, appear specifically in Cascales-Miñana et al. (2010). All the morphotaxa were included in the analyses, and a similar approach to that employed recently in similar methodological works in diversity analyses of British Westphalian-aged macroflora was adopted (Cleal 2005(Cleal , 2007. The systematic compilation of taxa in the construction of the dataset was undertaken by following a three-step procedure. First, taxonomy was standardised in accordance with recently published works (Kvaček et al. 2006;Willis et al. 2007;Jud et al. 2008;Wagner and Á lvarez-Vázquez 2008;Cascales-Miñana et al. 2010) using the taxonomic system of Cleal (1993aCleal ( , 1993b in The Fossil Record 2. Second, the data were codified in a primary double-entry matrix from floristic lists by placing taxa in rows and time intervals in columns in agreement with the methodology exposed by Digby and Kempton (1987) and Kovach (1988). Recent articles have used the binary nomenclature to codify data (Hilton and Cleal 2007;Cleal 2008b;Coiffard et al. 2008). According to this methodology, each box of the matrix takes the value of one or zero when taxa are present or absent, respectively, in a given time unit. Third, nine time units of the geological timescale were used in concordance with the absolute ages of Gradstein and Ogg (2004). The first eight time units correspond to Paleozoic time intervals, whereas the last time unit is a Mesozoic time interval in order to observe the impact on the record of the events occurred during the Paleozoic -Mesozoic transition. The time units have been named according to the International Commission on Stratigraphy with the time-stratigraphical terms, and the discussion has been written accordingly (Table 1).

Processing multivariate methods
Initially, four groups of primary analytical techniques were used in this study: principal coordinates analysis, non-metric multidimensional scaling (NMDS) analysis, correspondence analysis and detrended correspondence analysis (DCA). Nevertheless, preliminary results reveal that all these ordination methods present important differences in terms of robustness and methodological bias-related problems, particularly the principal coordinates and the correspondence analyses which exaggerate the differences between groupings (Davis 1990;Hammer et al. 2001;Podani and Miklos 2002;Hammer and Harper 2006). Detailed information about the preliminary results of these methods is given in Figure 1(A) -(C) and Tables S1 -S3. Accordingly, two ordination methods were chosen mainly for the discussion of the results: NMDS and DCA. An agglomerative hierarchical cluster analysis was also performed. All these methods were performed with the PAST software package (Hammer et al. 2001).

Non-metric multidimensional scaling
In general terms, ordination methods are essentially operations in a community data matrix (or taxa by sample matrix). In this analysis, the data matrix shows taxa as columns and time units as rows. Here, the ordination and the classification methods are used to interpret patterns in plant fossil composition terms. In agreement with these observations, I opted for NMDS that can also be applied to data of a binary nature (Shepard 1962;Kruskal 1964). On the other hand, NMDS offers the advantage of not requiring data linearity. This approach does not intentionally take absolute distances into account and is more robust than the eigen-analysis methods (Kenkel and Orlóci 1986;Taguchi and Oono 2005;Buja et al. 2008).
I performed the NMDS method using two different similarity measures: Jaccard's similarity coefficient (Jaccard 1908(Jaccard , 1912 and Raup -Crick's probabilistic similarity measure (Raup and Crick 1979). Jaccard's index was chosen because it uses a range from zero to one, is independent of double absence, emphasises presence instead of absence and is not vulnerable to sample-size differences Mapples 1987, 1989; Tables S1-S4. For abbreviations, see Table 1. Magurran 1988;Mapples and Archer 1988;Hengeveld 1990). In order to contrast the ordination results obtained with Jaccard's index, I used Raup -Crick coefficients because it is well documented that the 'Monte Carlo' randomisation procedure provides more robust results (Costeur et al. 2004;Cleal 2007Cleal , 2008bMaridet et al. 2007;Costeur and Legendre 2008). The implemented algorithm was based on the approach developed by Taguchi and Oono (2005).

Detrended correspondence analysis
A second approach based on the ordination methods was performed using the DCA. DCA is a relatively simple method where output is rescaled and detrended to correct the parabolic arch (known as the 'horseshoe effect'), where extreme points are compressed together at the ends of the zone, which represents the real primary axis (Hill and Gauch 1980;Hammer and Harper 2006). On the other hand, this ordination analysis has the advantage of not making assumptions on data normality or data linearity (Jongman et al. 1995;Shi 1995). Ordination axes are extracted by an iterative method (Braak 1995). For details, see Table S4.

Cluster analysis
To elucidate the robustness of the discontinuities shown by ordination methods, complementary cluster analyses have been performed using four similarity measures. I have used the following indices: Jaccard, Simpson, Dice (also known as the Sørensen Index) and Raup -Crick (Jaccard 1912;Dice 1945;Sørensen 1948;Simpson 1949;Raup and Crick 1979). In short, in concordance with Travouillon et al. (2007), I have used three descriptive binary indices together with a probabilistic measure to avoid any possible methodological bias in the interpretation of the results. In line with Kovach (1988), clustering was performed using the unweighted pair group method with Arithmetic mean.

Abundance and diversity analyses
In order to elucidate the possible role of the geologic and biological processes in the data interpretation and to provide a synthetic vision of the several floras present in the time units analysed, the data were analysed based on a set of abundance measures. In this study, two complementary abundance analyses were attempted. First, a relative abundance measure of the representativeness of each type of structural plant regarding the total record of each group of vascular plants during the Paleozoic, which is known as abundance analysis Type I. This analysis was calculated from the expression A k ¼ (T ck /N c ) £ 100, where k is the time unit, T ck corresponds to the number of the vascular plant type c taxa presenting a record during this interval and N c is the total record of this particular group of vascular plants during the Paleozoic.
Subsequently, the richness of each vascular plant type was measured in all the taxa of each time unit. In this paper, this parameter represents abundance analysis Type II. The richness percentage (R) was calculated by the expression R k ¼ (T ck /T k ) £ 100, where T k is the total number of taxa with a record in this time unit by adding all the vascular plant types considered as a single total (rhyniophytes, lycophytes, horsetails, ferns, progymnosperms, pteridosperms and gymnosperms). Relative abundance values are expressed as percentages and calculated for each analysed time unit.
Finally, in order to show the global patterns of changes which occurred in the plant life history during the Paleozoic, a global diversity analysis has been performed by employing absolute taxonomic counts; it has been shown in a single representation that includes all the considered plant groups as a measure of paleofloristic diversity. The raw data are available in Table S5.

Sample bias evaluation
In agreement with Krug and Patzkowsky (2004), and having taken into account that there are many inherent sources of bias in measuring diversity over time, including the variations in sampling intensity and the rock-record restrictions (Miller and Foote 1996;Foote 2000;Alroy et al. 2001;Krug and Patzkowsky 2004), I have attempted to illustrate these biases by compiling the data for this study. With this purpose in mind, a second database has been compiled. The data were extracted from the Paleobiology Database (www.paleodb.org; Alroy 2000).
The compiled dataset is formed by 435 plant fossil occurrences from 148 collections at the family level between the Silurian and the Triassic, and by considering only the paleobotany research group's paleontological data. To avoid variations in sampling not only across the environmental gradients (Sepkoski and Miller 1985;Smith et al. 2001), but also among the geographic regions (Maurer 1999;Fisher et al. 2010;Ricklefs 2010), the database configuration has considered all the collections present in Western Laurasia, Eastern Laurasia, Western Godwana, Eastern Godwana and Eastern Paleotethys, together with all the available environmental zones. Detailed information about the major contributors of this dataset are given in Table S6.
From the dataset, the following parameters were calculated for each time unit: the number of terrestrial rock formations with the presence of plant fossil occurrences, the number of plant fossil occurrences and the number of plant fossil collections. Subsequently, the following relations between the raw data for each time unit were calculated: the number of occurrences per formation, the number of collections per formation and the number of occurrences per collection.

Report of the ordination and agglomerative methods
The results of the ordination methods using the NMDS based on Jaccard and Raup-Crick coefficients are shown in Figure 1(D),(E). First, both similarity measures seem to present parallel patterns. Second, although the level of resolution and the differences between the Paleozoic time units appear to be clearly more emphasised with Raup -Crick coefficients, the stress level expressed by the analysis is lower when the Jaccard measure was used than for the Raup-Crick similarity measure: 0.0888 as opposed to 0.1648, respectively. The major discontinuities observed in the plant fossil history from the end of Carboniferous to the Triassic can be interpreted similarly and independently of the similarity measure chosen, given the main differences appearing in the first plant history steps from the end of the Silurian to the Carboniferous extinction event. Irrespectively of the type of ordination method employed (Figure 1), the results presented herein clearly show four important discontinuities in the data structure. These discontinuities appear in the following transitions: Early-middle Devonian, Devonian -Carboniferous, Mississippian -Pennsylvanian and early-late Permian. Furthermore, and depending on the type of analysis, other discontinuities appear between the late Silurian and the early Devonian (Figure 1(D)) or at the end of the Carboniferous (Figure 1(E)). Finally, the major discontinuity observed appears at the end of the Paleozoic between the late Permian and the early Triassic periods (Figure 1(F)).
Following the results based on the multivariate analyses, the complementary cluster analyses ( Figure  2(A) -(D)) clearly reveal a deep discontinuity in the data structure in the Devonian -Carboniferous transition. This discontinuity is present in the four approaches performed. In the same way, two discontinuities from these methods can be observed between the early Devonian and middle Devonian. All the cluster analyses also show that the Carboniferous period presents an important disruption in the reading of its floristic taxonomic structure, which appears to be reflected as a major discontinuity in the Mississippian -Pennsylvanian transition. On the other hand, and with the exception of the Paleozoic -Mesozoic transition, the observed similarity patterns are in parallel if we take into account the mathematical particularities of the used similarity indices. This transition only shows a deep discontinuity in the results based on the Jaccard and Dice coefficients (Figure 2(A),(B)). In contrast, if the similarity pattern is evaluated using Simpson's index or the similarity measure of Raup -Crick, this discontinuity appears between the early and late Permian (Figure 2(C),(D)). These results are in agreement with the interpretation based on the ordination methods.

Paleofloristic diversity patterns
As regards the abundance analyses, first I note how not only the richness of each vascular plant type, but also the diversity reflected for the plant record, clearly present heterogeneous patterns (Figure 3 (Table S5). Moreover, from the data results, a trend can be observed in the richness values of progymnosperms and gymnosperms. In the former, richness decreases progressively from the end of the Devonian, when it reaches its maximum value, to the end of the Paleozoic. In the latter, gymnosperms tend to increase from the Mississippian to the Late Permian, and a third of the flora is dominated by this group. The results indicate that the well-documented extinction events in the literature did not affect the different vascular plant groups equally. This fact is reflected in the relative contribution of each group to the total flora. Thus, the early-middle Devonian transition is characterised by an important reduction in the richness of rhyniophytes and lycophytes. Nevertheless, the Carboniferous -Permian transition only lowers by 48.7% for lycophytes, whereas the richness of other groups, such as horsetails and gymnosperms, increases by 29.9 and 75.8%, respectively. On the other hand, the relative abundance values show that each group of plants reaches its maximum record depending on the type of vascular plant being considered (Table S5). The maximum record for rhyniophytes occurs in the early Devonian with 100% representativeness in relation to the total record of this group. Lycophytes, ferns and progymnosperms mainly appear in the Pennsylvanian, whereas the maximum recorded values for horsetails and gymnosperms take place during the Permian. Finally, pteridosperms are mainly represented in the Mississippian.
The diversity analyses based on the raw taxonomic counts are summarised in Figure 4. In general terms, the diversity curves show global paleofloristic patterns where it can be observed how the primitive structural plants of vascular plants underwent a setback during the transition between the early and middle Devonian, whereas the lycophytes diversified considerably. This group's maximum diversity value is registered in the Pennsylvanian. On the other hand, the horsetails reveal an important increase towards the end of the Devonian and their diversity line remains constant until the end of the Paleozoic when this group increases considerably. Ferns present two clear peaks of diversity, one in the middle Devonian and the other at the end of the Carboniferous. In the Permian, ferns underwent important regression. Regarding progymnosperms, a strictly Paleozoic plant group, its maximum diversity value is noted in the Pennsylvanian. Progymnosperms were extinct by the end of the Permian. Alternatively, seed ferns present a notable, constant increase between the middle Devonian and the Mississippian-Pennsylvanian transition where they register a maximum value to then remain constant until the Paleozoic-Mesozoic transition when their diversity line dramatically drops. Finally, gymnosperms present their first major radiation in the Carboniferous, with a diversity peak during the transition between the early and late Permian. This group decreases notably in terms of the number of families at the end of the Paleozoic.

Geologic factor estimation
Despite attempts to obtain best sampling intensities, there is considerable variability in the number of accounts per time unit (Figure 5(A)-(F)). In order to elucidate the geologic effects and the sample bias about the plant fossil data, the analyses done support the idea that the oldest time units have important limitations to be able to evaluate the diversity present in these time intervals. These geologic constraints relate with the number of terrestrial rock formations, which markedly decrease from the late Devonian towards the end of the Silurian (Figure 5(A)). However, the results do not show a direct relation between the number of rock formations and the number of fossil plant occurrences. Therefore, the results show that the Mississippian and Pennsylvanian present the highest number of occurrences ( Figure 5(B)), although these time units also offer the lowest number of formations. In contrast, and in relation to older intervals, the number of registered collections of plants reveals a pattern that runs in parallel to the pattern shown by the number of formations ( Figure 5(C)).
In terms of the richness of rock formations however, the Pennsylvanian registers the highest rates of occurrences and collections per formation ( Figure 5(D),(E)). Nevertheless, these collections do not present high-quality levels in terms of the number of plant occurrences. In this sense, the best collections are those presented in the early Devonian and early Triassic (Figure 5(F)). These results show an important geologic factor and a highly heterogeneous sampling intensity over time with an unequalled sample collection opportunity.  Table S5. For abbreviations, see Table 1.

Quality of the plant fossil record data
The analysis of the paleobiological patterns in the fossil record based on multivariate methods urgently requires the prior consideration of several aspects which directly influence this type of inference. In this sense, Darwin's works question whether the fossil record provides a clear account of evolutionary history (Darwin 1859), which also occurred more recently because of its incompleteness (Benton and Storrs 1994;Benton and Hitchin 1996). In the last 25 years, however, paleontological knowledge based on several methodological approaches has clearly improved. In paleobiology, a series of methodological problems relating to the nature of paleontological data, all of which are inherent to these methods, have led to considerable controversy as to the quality of inference through records in recent times. Consequently, these methodological problems are inevitably present in this study and, therefore, I will now comment on them briefly.
Near the end of the twentieth century, paleontologists used taxonomic counts, a pool with a number of fossil taxa (in the present study, families) from a given time unit, as a quantitative measure of diversity across geologic times (Sepkoski 1997). Unfortunately, many biases can distort the measure of biological diversity of any time unit based on this classic paleontological approach (Raup 1976a(Raup , 1976bSepkoski 1997). Some studies in the paleontological literature have used taxa compilations within a higher taxon to reflect evolutionary trends . Other studies have performed several statistical techniques that account for various biases over wide timescales (Miller and Foote 1996;Peters and Foote 2002;Foote 2003), which has sparked an important debate in recent years.
In general terms, the three main concerns relating to paleobiological inference are taphonomic bias, taxonomic problems and geologic constraints. In this sense, I wish to remark that the taphonomic and preservation bias is probably the key factor of the incompleteness of the record (Galtier 1986;Solow and Smith 1997;). The perceived imperfections of the fossil record (Darwin 1859) are true for unmineralised tissues and organisms, despite the fact that our knowledge of such organisms continues to improve (Briggs and Gall 1990). Furthermore, taphonomic biases can greatly affect the distribution of fossils and could potentially alter our perceptions of generic richness (Behrensmeyer and Kidwell 1985;Behrensmeyer et al. 2000). Several recent studies have evaluated the diversity of the Phanerozoic fossil record and have concluded that this record is strongly biased by differential preservation towards the present time (Foote 1997(Foote , 2000Alroy et al. 2001). Obviously, this fact directly relates with the probability of the sampling bias. In agreement, the dynamics of  Table S5. For abbreviations, see Table 1. taphonomic processes helps our understanding of both the evolutionary interpretation and the paleobiological inference through the fossil record data (Gomez et al. 2001(Gomez et al. , 2002. Furthermore, plant species have peculiarities that make them special groups as regards their evolution (De Renzi 1995), and they develop characteristics that mark the different forms of preservation in the fossil record (Schopf 1975), thus influencing interpretation possibilities.
On the other hand, taxonomy is a complex and dynamic area undergoing continuous reviews, and this situation plays a crucial role in the interpretation of paleontological data from a paleobiological perspective if we take into account that the fragmentary nature of plant fossil material, the variety of different preservation states and the lack of information on breeding limits make the accurate identification of paleobiological taxa difficult. For this reason, and because of the controversial nature of the taxonomic data, I have used the family level in the approach attempted herein; it not only agrees with the idea that the use of higher taxonomic levels avoids possible biases arising when interpreting organic evolution with the fossil record in use, but also provides valuable information on features of evolutionary innovation and turnover patterns (Benton 1997;Benton et al. 2000;Lane and Benton 2003). Furthermore, the taxonomic problem in paleobiological inference increases given the fact that most fossil collections are poorly preserved, fragmented or prove especially difficult for correct identifications due to taxonomic constraints (Cleal and Thomas 2010), and they are not considered. This situation clearly leads to an impaired understanding of the fossil record, which could explain the low number of plant collections described in the oldest time units analysed ( Figure 5(C)).
Finally, the third main group of paleobiological bias relates to geologic constraints and the rock-record. In general, plant fossils in ancient rocks are more likely to have eroded, melted, crushed, not collected or been misunderstood than younger fossils. In addition, ancient rocks clearly preserve less information on average than more recent rocks, and they also show results in relation to the number of collections ( Figure 5(C)).
Consequently, several approaches have documented evidence of geologic constraints, which no doubt directly affect our perception of plant life history based on the plant fossil record (Peters andFoote 2001, 2002;Smith et al. 2001;Smith and Peterson 2002;Peters 2006). In the present paper, the results presented in Figures 3 and 5 also provide two excellent pieces of evidence. The results summarised in Figure 3 clearly show that the representativeness of the lycophytes record has a positive trend between the late Silurian and the Pennsylvanian (Figure  3(B)). On the other hand, the analysis of the number of terrestrial rock formations provides direct evidence of a reduction in available rock on outcrops over time, especially from the late Devonian to the end of the Silurian (Figure 5(A)). These results suggest that there are important geologic factors that bias the sample in the oldest time unit in sample collection opportunity terms, as the number of plant fossil collections presenting the same pattern for these intervals corroborates ( Figure 5(C)). In contrast, the Permian time units present the best accounts of terrestrial rock formations; nevertheless, the number of plant fossil occurrences drops dramatically towards the Paleozoic-Mesozoic transition (Figure 5(B)). Furthermore, the abundance analyses reveal that ferns and progymnosperms also experiment this drop in their representativeness values towards this boundary ( Figure  3(D),(E)). This reduction in seed ferns and gymnosperms may also be observed in the Triassic (Figure 3(F),(G)). These results do not show such an important geological bias as the previous case, although it is in accordance with the biological context at the end of the Paleozoic when Permian extinction took place. To summarise this point, and in accordance with Montero and Wagner (2008), geologic constraints show that the plant fossil record does not depend exclusively on evolutionary processes.
On the other hand and in relation to geologic time, Foote noted that the varying duration of different time units can distort diversity counts significantly, particularly if some time units vary greatly in length Foote 2000). This problem can be minimised by using not only the finest timescale possible, but also those that are evenly divided. Ideally, the length of each time interval should be less than the duration of the taxa being counted, and the taxa should occur in multiple time intervals. Unfortunately, many plant fossil taxa are based on single occurrences and most fossil families are confined to a single time unit (the so-called singleton taxa). So in my opinion, the use of boundary crossers would be inappropriate, particularly for performing this approach to show the main floristic times in a global analysis.
The re-emerging debate provoked by these methodological problems was partially resolved by the observation that different datasets gave similar patterns of rising diversity over time (Benton et al. 2000). Benton (1995) and Benton et al. (2000) emphasised that the geological factor did not fatally damage the signal to be read in the fossil record. For this reason, and to minimise the taxonomic, sampling and rock-record biases (Raymond and Metz 1995;Foote 1997;Benton et al. 2000;Lane and Benton 2003;Wang and Bush 2008), several types of fossils were considered for each taxon in the sample (Benton 1993;Taylor et al. 2009). Furthermore, and in accordance with the methodology presented by Benton et al. (2000), I have used The Fossil Record 2 as a major source of stratigraphic data for plants' dates of origin to present a consistent and contrasted source of information in the paleobiological approach used herein. The family taxa used in the present study avoid possible bias arising from the interpretation of organic evolution using the fossil record and provide valuable information on feature of evolutionary innovation and turnover patterns (Cascales-Miñana et al. 2010). Finally, by taking into account that plant fossils are unique and direct evidence of past plant life, together with all the biases exposed herein, although no definite conclusions on this issue can be drawn from analysing only one geologic era, the results suggest that, in accordance with Krug and Patzkowsky (2004), plant discontinuities studies should attempt to remove the distorting effects of incompleteness and variation in sampling intensity if macroevolutionary and ecological processes are to be fully understood.

Discontinuities of the Paleozoic plant fossil record
The cause of extinction events and the nature of biological selectivity are central questions in paleobiology. It is not possible to deduce the cause of the several extinction events of plant life from the ordination and clustering methods. Nevertheless, their consequences can be observed in taxonomic composition terms. The negative contributions of extinctions can be seen as the elimination of many taxa. Therefore, this consequence must be reflected and evaluated from the fossil record, and must be based on multivariate analyses in terms of the discontinuities and distance between time units, as I comment hereafter.
In the last few years, some dissenting voices were raised against the concept of the mass extinction of the marine invertebrate record being applied to the plant record studies. In accordance with this idea, this challenge of the orthodoxy relating to the extinction events proclaimed several extinction episodes of plant life from the vegetation dynamics viewpoint (McElwain and Punyasena 2007;McElwain et al. , 2009 in mass ecological reorganisation terms (DiMichele et al. 2001;Cascales-Miñana et al. 2010). On the other hand, previous works have also corroborated the incredible resilience across the mass extinction events of plant families based on a complete set of paleobiological inference methods (Cascales-Miñana and Ruiz-Sánchez 2007, 2008a, 2008b, 2009) and, accordingly, the discussion of the results presented herein can also be interpreted in these terms. Subsequently, I will comment on the important discontinuities presented in the 'Results' section, as summarised in Figures 1 and 2. For instance, in the early-middle Devonian transition, the decline of the early evolutionary trends of vascular plants characterised by terminal sporangia and a dichotomous branching occurred (Cleal and Thomas 1994;Roth-Nebelsick et al. 2000), whereas the sporangia of others, such as the zosterophylls and trimerophytes, aggregated in terminal spikes (Edwards 1969;Banks and Niklas 1989) and a trifurcating branching pattern (Niklas and Banks 1985;Rothwell 1999), respectively. This could be related to the appearance of progymnosperms in middle Devonian times. Progymnosperms appear to have had a remarkable effect on plant evolution because they developed large, secondary-thickened, gymnospermlike axes and many of them achieved the arborescent habit (Fairon-Demaret and Leponce 2001;Scheckler and Galtier 2003). Note that all the great key evolutionary innovations of the Paleozoic had appeared and that plant life had explored new developmental possibilities by this point of the geologic temporal line, which probably accentuated the recession of these primitive structural plants and, consequently, the disappearance of early vascular plant taxa took place (Gensel and Andrews 1984;Doyle 1998;Gensel and Berry 2001), as the diversity curves in Figure 4 show. This discontinuity in extinction terms is reflected in the ordination methods. This recession of the primitive structural plants also clearly appears in the results of the abundance analysis where we can perceive how the percentage contribution of the rhyniophytes took place together with a significant increase in the representativeness of progymnosperms and club mosses ( Figure  3(A),(B),(E)). This result would be in agreement with the taxonomic turnover process on a macroevolutionary scale (Cascales-Miñana et al. 2010). Nevertheless, the sampling bias process cannot be ruled out because the relative abundance of the early vascular plants decreased synchronously with a considerable drop in its representativeness in these groups. This is likely because these early plants present a fossil record that relates to some latitudes (Montero and Wagner 2008), and this fact, together with the ancient rock effects (Benton et al. 2000;Tarver et al. 2007), could decisively slant the sample, as observed in the negative trend of the number of terrestrial rock formations towards the oldest time units that represent a clear reduction in available rock on outcrops across this transition ( Figure 5(A)). Accordingly, the opportunity to increase the number of sample collections lessens.
During the Devonian-Carboniferous transition, fossil plant groups, including horsetails, early seed plants and true ferns, are all acknowledged in the fossil record. Yet during the Carboniferous period, we can also see the great diversity of exclusive groups of vascular plants, such as seed ferns or pteridosperms. It is at this particular moment of geologic time that an important extinction event occurred (Godderis and Joachimski 2004). In this sense, the Late Devonian glaciation is associated with increased extinction on the Devonian -Carboniferous boundary (Raymond and Metz 2004). Nevertheless, this moment is characterised by the first vascular seed plants branching out (Figure 4), and this new structural plant seems to be a key factor to explain the discontinuity of the data structure at this point. Note also that the number of outcrop areas decreases throughout this transition ( Figure 5(A)). Nevertheless, the number of plant collections per formation increases from the middle Devonian to the Pennsylvanian ( Figure 5(E)).
Initially, Mississippian plants were predominantly a continuation of the last Devonian groups. However, they were partly displaced by arborescent lycophytes, pteridosperms and ferns, all of which became important elements in Carboniferous flora ( Figure 4); their representativeness increased as did their relative diversity in the Carboniferous time units (Figure 3(B),(D),(F)), and only a slight decrease in abundance was seen for seed ferns, probably due to a sampling bias because this fact appears together with a slight reduction in their representativeness in the Pennsylvanian (Figure 3(F)). Contrasting with this positive dynamics noted among the other vascular plant groups analysed, the relative diversity record for horsetails in the Mississippian and Pennsylvanian decreased from the late Devonian to the end of the Carboniferous period (Figure 4(C)), whereas its representativeness remained constant during this time interval. During the Carboniferous, the ordination methods presented reveal a heterogeneous taxonomic composition. This fact could relate to the appearance of new taxa (the singleton taxa) with records in a single time interval that, presumably, occupied the niches left during the last extinction event at the end of the Devonian (Cascales-Miñana et al. 2010). Furthermore, a well-documented paleoclimatic change occurred during the Carboniferous (Phillips and Peppers 1984;Phillips et al. 1985;Kosanke and Cecil 1996). So we may state that club mosses underwent a major extinction event after the drying of habitats (Taylor et al. 2009). Note also that ferns became a main group at the end of Carboniferous and an important recession is observed in its representativeness record on the Carboniferous -Permian boundary (Figure 3(D)), as the strong negative slope of its diversity curve also corroborates ( Figure 4). All these considerations, together with the first steps in the change of the dominance of vegetation types towards the gymnosperm-dominated plant ecosystems, as demonstrated by the diversity and abundance analyses (Figures 3(G) and 4), and together with true fern forms, could explain this heterogeneity nature of the record. To go about this, it is necessary to take into account that this geologic time is probably the most well-documented time of Paleozoic plant fossil history, as the analyses on the sample bias over time by geologic factors show, and where we can clearly observe the high number of occurrences and collections between the Pennsylvanian and the early Permian ( Figure  5(B),(C)).
A similarity context occurs in the Permian period in biological terms. This period marks the culmination of the first great development in land plants. Seed plants supplanted club mosses, horsetails and tree-ferns in a macroevolutionary mass turnover which is clearly reflected in the plant record (McElwain and Punyasena 2007;McElwain et al. 2009). These turnovers are based on changes in vegetation which during the Paleozoic -Mesozoic transition were mainly dominated by lycophytes, sphenophytes, filicophytes and pteridosperms, and became dominated by other vegetations such as cycads, ginkgos and conifers (Knoll 1984;Knoll and Niklas 1987;Jablonski and Sepkoski 1996;Willis and McElwain 2002). Since the Permian presents a good number of outcrop areas ( Figure 5(A)), and as the Paleozoic-Mesozoic transition collections offer relatively good quality ( Figure 5(F)), all these processes could be interpreted as a result of the environmental changes relating to the Permian extinction.
The last transitions are presented from different viewpoints given an important methodological component as their detection depends on the algorithm used. This is precisely the case of the discontinuity between the Silurian and the Devonian. This fact can be easily explained from a biological viewpoint because, from this geologic time, the great radiation of vascular plants at the early Devonian took place, and these new plant forms removed the early forms of vascular plants. This observation could relate to the high origination levels observed during the early Devonian (Cascales-Miñana and Ruiz-Sánchez 2008a; Cascales-Miñana et al. 2010), and was probably due to this fact. Nevertheless, from a geologic and methodological viewpoint, these records also present the most fragmentary record included in oldest rocks which directly affect the sample bias; therefore, this transition appears with different magnitudes depending on the data processing. In contrast, the end of the Carboniferous is a fine documented geologic time, as previously noted. In this last case, where the geologic constraints seem to be minimised, the structure of these data, more probably, shows the consequences of the extinction process that was characterised, according to the literature, by two humid periods of unequal magnitude in relation to the water regime (Phillips and Peppers 1984;Cecil et al. 1985;Phillips et al. 1985;Kosanke and Cecil 1996); the gradual drying of habitats was one of the main causes of the turnover of the structural plants at those times along with the disappearance of the giant forms of lycophytes.

Concluding remarks
This study uses a set of multivariate analyses to show the paleofloristic composition of the Paleozoic fossil record from a multi-method viewpoint, especially from ordination methods, complementary agglomerative methods, abundance analyses and an evaluation of the sample bias effect. This study provides new evidence that ordination and agglomerative methods are powerful tools for such work. However, it also emphasises the importance of using more than one type of multivariate analysis and comparing the results with other paleobiological inference methods, such as abundance analyses, measures or percentages of representativeness, together with an estimation of the geologic component in the inference from the fossil record, to obtain a better interpretation of the plant fossil record data. Nevertheless, given the controversial nature of the fossil record data, which present an inevitable set of biases, this approach will have to be evaluated in the future with new paleobotanical data and, therefore, with new methodological analyses.
Paleozoic flora is characterised by the greatest variety of anatomic forms of vascular plant life history. This study corroborates that, this geologic time interval provides most valuable information to understand the trends and patterns in the evolution of vascular plants (Cascales-Miñana et al. 2010). During the Paleozoic, three key important events in plant history took place: the first radiation of vascular plants, the appearance of the great evolutionary innovations and extinctions that have affected vascular plants the most, these being the Carboniferous and Permian extinctions. All of them are characterised by significant differences in the paleofloristic diversity composition in terms of relative abundance, which show deep structural reorganisations of the ecosystems.
This study shows that the changing composition of the Paleozoic flora is reflected in NMDS and DCA plots, and that the complementary abundance analysis clearly shows the relative importance of different plant groups over time.
The several ordination analyses presented herein reveal significant discontinuities in terms of taxonomic composition and record representativeness. Although they can be explained in taxonomic extinction terms from a classic paleontological viewpoint, the data structure can also be explained in taxonomic turnover and ecological reorganisation terms in accordance with the new recent interpretations of this topic (McElwain et al. 2009;Cascales-Miñana et al. 2010) by taking into account the crucial role of the geologic constraints.