Total DNA was extracted from specimens preserved in 70–95% ethanol using the QIAamp DNA Micro Kit (Qiagen, Germany), following the manufacturer's protocol. Six markers were sequenced for this study, including three nuclear genes (one protein-coding gene, H3, and two rRNA genes, 18S and 28S) and three mitochondrial genes (one protein-coding gene, CO1, and two rRNA genes, 12S and 16S). Each gene was amplified using the primers listed in Table 1. PCR reactions were carried out in a 30 μl final volume using the following conditions: 10x reaction buffer with MgCl₂, 3 μl; dNTP mix (6.6 mM), 3 μl; primers (10 μM), 1.5 μl; H₂O, 19.3 μl; Sigma Red Taq DNA polymerase, 0.70 μl, and DNA template, 1 μl. The cycling protocol included an initial denaturation step for four minutes at 94 °C, 30 cycles of 30 s at 94 °C, 30 s at the appropriate annealing temperature and 60 s at 72 °C. Final extension followed for 10 min at 72 °C. Purification and cycle-sequencing reactions were performed at the Genoscope (Evry, France). Both DNA strands were sequenced for all PCR products. Sequences were edited and assembled using Sequencher 4.7 (Gene Codes Corporation, Ann Arbor, MI, USA).

Two different datasets were analyzed in this study: a dataset containing only 18S sequences and a dataset including six genes. The nine outgroup taxa used for the analyses of both datasets were chosen using the two following criteria: (1) complete sequences available for all the six selected markers; and (2) no evidence of reverse strand bias in the mitochondrial sequences. They include two Crustacea, two Euchelicerata, two Hexapoda and three Myriapoda (Table 2).

The first dataset was constructed by using all nearly complete 18S sequences available for Pycnogonida in the GenBank/EMBL/DDBJ nucleotide databases. Several sequences were characterized by high levels of missing data or degenerate nucleotides, suggesting problems during PCR amplification or chromatograms of bad quality (double peaks, weak signal, etc.). The more the amount of missing data is important in a sequence, the more the quality of the rest of the sequence is doubtful and can have an impact on phylogenetic reconstruction. We used therefore two successive criteria for selecting high-quality 18S sequences: a sequence length superior to 1200 nucleotides and less than 10% of missing data and/or gaps in the final alignment. Our alignment includes 115 sequences and represents all the 10 families of pycnogonids currently recognized.

The second dataset was constructed to compare the phylogenetic signals between six markers (12S, 16S, 18S, 28S, CO1 and H3). To limit the impact of missing data, we used only taxa for which all the six genes were sequenced for more than 90% of the total length. Our final alignment includes 43 sequences, and represents eight pycnogonid families (Table 2).

Alignments were performed manually with Se-Al v2.0a11 [18]. All regions in the alignments involving ambiguity for the position of the gaps were excluded from the analyses to avoid erroneous hypotheses of primary homology. The alignments of rRNA genes were found to be problematic for the most divergent outgroup taxa. Therefore, non-ambiguous regions were first determined by aligning ingroup species, and this framework was then used to align outgroup taxa. The alignments are available upon request from the authors.

Separate and combined analyses were carried out using Bayesian inference (BI) and maximum likelihood (ML). The model of molecular evolution was selected under jModelTest 0.1.1 [19] using Akaike and Bayesian Information Criteria (AIC and BIC). The best-fit model was GTR+I+G for all markers. Bayesian posterior probabilities (PP) were calculated with MrBayes 3.1.2 [20] using four independent Markov chains run for 10,000,000 Metropolis-coupled MCMC generations, with tree sampling every 100 generations and a burn-in of 50,000 trees. The node robustness was also estimated after 1000 replicates by ML bootstrap analyses under RAxML v7.0.4 (BP) [21]. The parameters of the GTR+I+G model were estimated using different partitions under MrBayes and RAxML, corresponding to each of the markers used and each of the three positions for protein-coding genes (CO1 and H3).

The strand bias in nucleotide composition was analyzed at third codon positions of CO1 sequences using the approach detailed in Hassanin et al. [16]. As similar trends were found for two- and four-fold degenerate sites, we used a simplified method where all third positions were analyzed together. The frequencies of complementary nucleotides were compared and tested by skewness: AT skew = [A − T]/[A + T] and CG skew = [C − G]/[C + G]; AT and CG skews were statistically significant if the null hypothesis of symmetry was rejected at a confidence level of 0.05.

The analyses were also carried out for each of the 13 protein-coding genes of three complete mitochondrial genomes: two species of Pycnogonida, A. bituberculata (NC_009724; [22]) and N. gracile (NC_008572; [17]), and the horseshoe crab L. polyphemus (NC_003057; [23]).

3.1.1 Dataset construction

3.1.2 Phylogenetic results

In the Bayesian tree (Fig. 2), the class Pycnogonida constitutes a well-supported monophyletic group (PP = 1, BP = 100). Three published sequences present a discordant phylogenetic position, as they result in the apparent polyphyly of three distinct genera and families: Anoplodactylus (DQ389912; [1]), Callipallene (AY210808; [25]) and, Colossendeis (AF005440; [26]). For instance, the sequence of Colossendeis AF005440 is found within a clade including Ascorhynchidae, Callipallenidae and Nymphonidae families. Three hypotheses can be advanced to explain this result: real polyphyly, taxa misidentification or carryover contamination. All the nine other sequences of Colossendeis are robustly enclosed in the family Colossendeidae with the sequences of three other genera (Decolopoda, Hedgepia and Rhopalorhynchus). The fact that these sequences were produced by three independent teams [1,12, and this study] supports their authenticity, and strongly questions that of AF005440 sequence, which may be therefore interpreted as taxa misidentification or carryover contamination. The same argumentation can be made for the sequences of Anoplodactylus DQ389912 and Callipallene AY210808.

3.2.1 Phylogenetic information of markers

3.2.2 Separate analyses

3.2.3 Combined analyses

Taking into account sequencing errors detected in published H3 sequences, and the difficulties encountered during PCR amplification of H3 in our laboratory, we choose to present only the results of combined analyses without H3. However, Bayesian and ML analyses including H3 produced similar results, but the family Ascorhynchidae was not found monophyletic (data not shown). The Bayesian and ML analyses of the total data matrix produce similar topologies. Comparisons with the analyses of mitochondrial and nuclear datasets allow us to identify the most reliable nodes (Table 4 and Fig. 3). All three analyses (total, mtDNA, and nuDNA) strongly support the monophyly of Chelicerata and that of several other higher taxa, such as Colossendeidae, Endeidae, Pallenopsidae, Phoxichilidiidae, Pycnogonida and Pycnogonidae (PP = 1, BP = 91–100).

3.3.1 Analysis of CO1

3.3.2 Analysis of complete mitochondrial genomes

The nucleotide composition of third codon positions was analyzed for each protein-coding gene of three mitochondrial genomes (Fig. 4): Limulus polyphemus [23], which has a gene order identical to that deduced for the common ancestor of Chelicerata [16], and two pycnogonids, i.e., Achelia bituberculata [22], belonging to the family Ammotheidae sensu stricto, and Nymphon gracile [17], belonging to the family Nymphonidae.

Adequate taxonomic sampling is a fundamental prerequisite in phylogenetic reconstruction. Particularly important is the choice of outgroup taxa, because incorrectly rooted trees may result in erroneous interpretations of phylogenetic relationships. The inclusion of at least three outgroups, which have a paraphyletic position with respect to the ingroup, is therefore highly recommended, because it allows the monophyly of the ingroup to be tested, reduces the chances of inappropriate outgroup selection, and it may be used to break up the branch between ingroup and outgroup, thus reducing long-branch attraction problem [29,30]. In case of substantial acceleration of evolutionary rates in the sister taxon of the ingroup, more distantly related, but less divergent, outgroups may also be used to prevent spurious placement of the root [31].

On the basis of 18S sequences, Nakamura et al. [12] have concluded that Ascorhynchus and Eurycyde are early offshoots of sea spiders. However, they used only one outgroup for rooting their tree (Euchelicerata, represented by five species). We included two additional outgroups in our analyses: Myriapoda and Pancrustacea. Both 18S and multigene analyses show a different pattern for basal relationships, since the families Pycnogonidae and Colossendeidae appear as early offshoots, a topology which is in agreement with the results of Arango and Wheeler [1]. It is straightforward that the tree provided by Nakamura et al. [12] was incorrectly rooted, due to the long-branch separating pycnogonids from euchelicerates.

How reliable are the sequences available in nucleotide databases? A portion of nucleotide sequences is known to contain systematic sequencing errors [32], or to be generated from misidentified taxa [33] or due to carry-over DNA contamination [27]. The identification of arthropods may be especially problematic, particularly in the case of pycnogonids, which are poorly studied and for which morphological characters are not easily discerned.

Our analysis of all available 18S sequences of Pycnogonida indicate that two published sequences may be misidentified at the genus level: AF005440 (Colossendeis) and AY210808 (Callipallene). Despite the fact that these sequences were not included in previous pycnogonid phylogenies, their integration into further analyses would be at the origin of a misinterpretation of phylogenetic results. Unfortunately, the first sequence cannot be linked to a voucher specimen, which prevents further studies to determine its taxonomic status. If keeping voucher specimens has always been the ideal approach for systematics, it was not a standard practice in the beginnings of molecular systematics. From this point of view, the Barcode of Life Data system (BOLD) represents a real improvement for future phylogenetic studies. DNA barcoding is a DNA-based identification system, which is founded on the mitochondrial gene CO1. It focuses on the assembly of reference libraries of barcode sequences for voucher specimens with authoritative taxonomic identifications [34]. This strategy would usually make it possible to determine whether the organism from which the DNA was extracted had been correctly identified.

DNA contamination is also problematic for phylogenetic reconstruction. The use of chimeric sequences can have dramatic effects, as it can lead to robust but incorrect conclusions, or can lead to lack of resolution by increasing the levels of phylogenetic incongruence amongst sites [27,35]. For instance, we pointed out that one of the three fragments of the 18S sequence DQ389912 of A. lentus is a carry-over DNA contamination by Ascorhynchus. This kind of contaminations is very difficult to detect. It should be standard practice to analyze the segments separately where concatenated sequences are combined from multiple PCR products. The deposit of chromatograms in sequence databases would dissipate such cases of uncertainty.

Missing data have long been considered as a source of error in phylogenetic analyses, because the inclusion of taxa with incomplete data can lead to a lack of resolution [36,37]. Simulation and empirical studies suggest, however, that missing data are not problematic for phylogenetic inferences if the number of characters is sufficient to provide a robust signal [38,39].

Two kinds of missing data have to be distinguished: partial DNA sequences and incomplete taxa sampling for one or several markers. Incorporating low quality sequences in our datasets introduces noise which may have two different consequences: the lack of signal can lead to a lower phylogenetic resolution, and the introduction of sequencing errors can lead to erroneous phylogenetic inferences. Our strategy was to minimize missing data in our analyses in order to provide equivalent phylogenetic information among taxa, a necessary condition for comparing the signals of different markers. This approach allowed us to study the contribution of each marker to the phylogenetic results, and to identify conflicting signals between genes. By contrast, Arango and Wheeler [1] included several taxa with a large proportion of missing data. For example, Rhynchothorax australis (family Rhynchothoracidae) was only sequenced for H3 (327 nt) and a small part of the 28S gene (280 nt), which represent less than 10% of the total number of characters. Although Rhynchothorax was found associated with Ammotheidae sensu stricto (Fig. 1), Arango and Wheeler [1] did not provide any support value. Our ML analyses indicate that R. australis belongs to a large clade composed of all pycnogonid families except Pycnogonidae and Colossendeidae (BP = 76), but there is no signal for supporting a relationship with Ammotheidae sensu stricto (data not shown). Using an 18S sequence of R. mediterraneus, Nakamura et al. [12] have obtained a different position for the genus Rhynchothorax: it was found to be allied with Pycnogonidae (BP = 86), hypothesis here confirmed by our ML analyses (BP = 91, data not shown).

There are reasons to think that there may be a problem with the sequences of Rhynchothorax provided by Arango and Wheeler [1]. Firstly, only two genes on six have been sequenced, which indicates a poor DNA conservation for this specimen. Secondly, absence of contamination cannot be guaranteed because the sequences are the result of short (327 nt and 280 nt) and unique amplification products. Interestingly, a sister-group relationship between Pycnogonidae and Rhynchothoracidae is corroborated by the presence of a unique pair of genital pores in females of these families [40]. Considering correct the identification of R. mediterraneus [12], we are in favor of a sister–group relationship between Pycnogonidae and Rhynchothoracidae.

Two categories of multiple sequence alignment can be distinguished: the traditional approach, and the direct optimization approach [41]. In the traditional approach, homology is defined prior to tree reconstruction, whereas alignment and tree search are performed simultaneously in the direct optimization approach implemented in the program POY [42].

Comparisons of rRNA sequences often show important variations in the length of loop segments of secondary structure. The alignments of rRNA loops are therefore characterized by many indels and most of them have uncertain positions [43]. Indeed, all of our rRNA datasets contain high percentages of ambiguously aligned sites (63.4, 48.4, 16.4 and 51% for 12S, 16S, 18S, and 28S, respectively; Table 3). This issue is particularly problematic for the mitochondrial rRNA genes (12S and 16S), because they evolve much faster than the nuclear ones (18S and 28S). We first aligned each marker, focusing on ingroup taxa. Integrating outgroup taxa increased uncertain positions and the strict application of the primary homology criterion to all taxa would have constrained us to discard 82.9 and 66.7% of sites for 12S and 16S, respectively (data not shown). Given these difficulties, the use of mitochondrial rRNA genes for deciphering basal relationships within Pycnogonida is questionable.

Under direct optimization, ambiguous characters are not removed. According to Gatesy et al. [44], the exclusion of these ambiguously aligned sites is subjective and associated with a loss of parsimony-informative characters. As pointed by Simmons [45], positions that would be considered ambiguously aligned with the similarity criterion, and therefore excluded because of violation of primary homology, may be considered unambiguously aligned using the DO criterion.

Several studies have concluded that the dynamic approach is less accurate than the traditional approach [46,47]. Direct optimization does not offer a realistic model of indel evolution and the way in which POY handles rRNA loops remains problematic. Moreover, it has been demonstrated that the removal of ambiguous regions can lead to better trees [48]. Because we consider that character homology assessment is an essential step in phylogenetic analyses, we recommend the exclusion of ambiguous regions.

Another advantage of traditional alignment is that it can reveal sequencing errors in published data, such as a shift of the reading frame in protein coding gene sequences. This strategy has the advantage of checking each gene matrix prior to proceeding to phylogenetic inferences and allowed us to identify major sequencing errors in sequences extracted from nucleotide databases. Indeed, several H3 sequences published by Arango and Wheeler [1] are affected by a reading-frame shift. Because they are not identified when performing a direct optimization approach, such sequencing errors contradict the primary homology criterion and can lead to misleading conclusions.

A typical metazoan mitochondrial genome is a circular and double-stranded DNA molecule. Most mitochondrial genomes present a strand asymmetry: one strand is rich in Adenine and Cytosine, whereas the other is rich in Thymine and Guanine [16,49] The strand with positive AT and CG skews can be defined as the “positive strand”, whereas the other strand can be defined as the “negative strand” [5].

Several independent reversals of asymmetric mutational constraints have been demonstrated during the evolution of Metazoa. The most spectacular cases include the echinoderm Florometra, the mollusk Katharina and several arthropods, such as the cephalocarid Hutchinsoniella, the copepod Tigriopus and the hemipterans Aleurodicus and Trialeurodes [5,16]. The study of Echinodermata has shown that asymmetric mutational constraints can be reversed through two different mechanisms: (i) inversion of the control region, which results in a global reversal, and (ii) gene inversion, which results in a local reversal [16]. Interestingly, two main groups of Chelicerata are characterized by global reversals of asymmetric mutational constraints in their mitochondrial genome, including all sequenced species of scorpions and Opisthothelae spiders [5,50]. It has been therefore suggested that two independent inversions of the control region occurred during the evolution of Chelicerata: one in the ancestor of scorpions, and another one in the ancestor of opisthothele spiders [5].

In this study, the strand bias was described by skewness (AT and CG skews) at third codon positions of CO1. The results revealed a strong heterogeneity in nucleotide composition within sea spiders (Fig. 3): six unrelated species of Ammotheidae sensu stricto and Nymphonidae present a reversed nucleotide composition, which suggests that several reversals of asymmetric mutational constraints have occurred independently within Pycnogonida. We suggest that the reversals of base composition observed in the CO1 gene of these species are the consequences of independent genomic inversions, concerning either CO1 alone, a genomic fragment including CO1 and other genes, or the control region. With time, substitutions accumulated in CO1 sequences, which has led to the establishment of several independent strand bias reversals. The analyses of the two complete mitochondrial genomes of Pycnogonida available in the databases confirm our hypothesis (Fig. 4). Both genomes differ from Limulus polyphemus, a species that presents the same organization as the ancestor of Chelicerata [16]. In Nymphon gracile, the genomic fragment including the three protein-coding genes Nd2, CO1 and CO2 is inverted [17] and, logically, the strand compositional bias of this fragment is reversed with respect to that of Limulus. The case of Achelia bituberculata is more delicate to interpret because the values of strand bias are not significant, i.e. A ≈ T and C ≈ G. These results suggest that the transposition of tRNA-Q between the 12S gene and the control region [22] was a recent event accompanied by the inversion of the control region, thus inducing a global reversal of asymmetric mutational constraints.

The analysis of base composition in the CO1 gene revealed a strong heterogeneity between species of Achelia: A. assimilis exhibits significant positive AT and CG skews, whereas A. hispida has significant negative values of strand bias (Fig. 3). These results suggest multiple changes in asymmetric mutational constraints during the recent evolution of Achelia, which may have been caused by several independent genomic rearrangements. This hypothesis needs to be confirmed by sequencing the complete mitochondrial genome of A. hispida and A. assimilis.

Several previous studies have shown that reversals of strand-specific bias can be strongly misleading for phylogenetic inference, since taxa with reverse nucleotide composition tend to group together [5,16]. To prevent the misleading impact of reversals of asymmetric mutational constraints on phylogenetic reconstruction, we did not used outgroup taxa showing negative CG and AT skews, such as Onychophora, Araneae and Scorpiones. As discussed above, several changes in asymmetric mutational constraints have been demonstrated in Pycnogonida. Despite some differences between phylogenetic analyses (total, mtDNA and nuDNA), there is no important conflict in the position of species characterized by strand bias reversals. This result therefore suggests that reversal events are too recent in the evolution of pycnogonids to erase the phylogenetic signal among closely related species. Comparison of the mitochondrial genome organization of various pycnogonid genera in a phylogenetic framework will be essential for a better understanding of the molecular evolution of the group and for a detailed analysis of mitochondrial features.

In an attempt to resolve pycnogonid relationships, we sequenced all the six markers currently used in arthropod phylogeny: CO1, 12S, 16S, H3, 18S and 28S. Estimations of parsimony-informative sites for each gene highlight the poor efficiency of 18S and 28S genes to answer this question, with respectively 4.4 and 5.7% of parsimony-informative sites for the ingroup (Table 3). Other markers show higher levels of phylogenetic information (38–55%). As revealed by separate analyses, the results of H3 are divergent from those found with other markers (Table 4). The origin of these conflicts may be linked to DNA contamination and sequencing errors, because several uncommon changes have been detected in the sequences of H3, including shifts of translation reading-frame or changes in highly conserved amino acids. Such sequencing problems are expected to be particularly frequent during PCR amplification of H3 due to its small size (less than 400 nt) and the use of degenerated primers that amplify across Arthropoda. Consequently, the H3 fragment was excluded from our combined analyses and we do not recommend this marker for future molecular studies on Pycnogonida and more generally on Arthropoda.

The analyses based on mtDNA, nuDNA, and total datasets allowed us to compare the phylogenetic signals provided by two independent evolutionary units: the mitochondrial genome, which is maternally inherited, and the nuclear genome, which has a maternal and paternal inheritance. We consider that nodes supported by both mtDNA and nuDNA are the best indicators of “true relationships” among species. Using this repeatability criterion, we identified 17 reliable nodes in the tree of Pycnogonida (Fig. 3).

Previous studies have provided discordant conclusions for basal relationships within Pycnogonida (Fig. 1). On the one hand, Arango and Wheeler [1] have proposed a major dichotomy separating Austrodecidae, Colossendeidae, Pycnogonidae and Phoxichilidium from all other Pycnogonida. On the other hand, Nakamura et al. [12] have obtained a basal position for Ascorhynchus. Our analyses instead support an early divergence of the family Austrodecidae, and the existence of a large clade including all other families except Colossendeidae, Pycnogonidae and Rhynchothoracidae. Although not strongly supported, these results are probably reliable, because they are recovered by both mtDNA and nuDNA analyses (Table 4 and Appendix A). Unfortunately, the position of Rhynchothoracidae as sister-group of Pycnogonidae was only supported by 18S sequences, and that of Phoxichilidium could not be tested here due to lack of data. It is worth noting that the family Austrodecidae presents significant morphological autapomorphies: a slender and annulated proboscis; a biradiate mouth [1,51]. However, characters supporting its early divergence with respect to other pycnogonid lineages have not been demonstrated.

Nakamura et al. [12] removed Ascorhynchus and Eurycyde from the family Ammotheidae and included them into their own family, Ascorhynchidae. Surprisingly, their analyses show a paraphyletic pattern for both Ascorhynchidae and Ascorhynchus at the base of the pycnogonid tree. However, as explained above, this pattern is due to misrooting, the consequence of which is that all basal relationships within Pycnogonida were erroneously interpreted in Nakamura et al. [12]. According to Arango [52], several morphological characters support the monophyly of Ascorhynchidae: presence of terminal claws; multiple spine rows in ovigers; absence of auxiliary claws and heel spines in walking legs. Our phylogenetic inferences did not produce a strong signal for the monophyly of Ascorhynchidae. However, Ascorhynchus and Eurycyde are distantly related to Ammotheidae, suggesting that they should be excluded from this family.

The monophyly of Ammotheidae sensu stricto is retrieved in the combined analysis. Nevertheless, the matter is not completely settled, since mtDNA and nuDNA analyses instead indicate paraphyly of the group. Monophyly of Ammotheidae sensu stricto was not recovered in Arango and Wheeler [1] due to the inclusion of Rhynchothorax. In Nakamura et al. [12], Ammotheidae sensu stricto is not supported because of the sister-group relationship between Tanystylum and Anoplodactylus (Fig. 1). Ammotheidae is constituted by 395 species in 28 genera [1], so adding taxa in future analyses will be necessary to better understand relationships within this family and determine its validity. Furthermore, as suggested by Bamber [51], only the addition of distinct ammotheid genera such as Bathyzetes, Calypsopycnon or Heterofragilia can clarify the validity of both Ammotheidae and Ascorhynchidae.

Paraphyly of Callipallenidae was obtained in all separate and combined analyses, supporting the reliability of this result. Contrary to Arango and Wheeler [1], we found Nymphonidae to be monophyletic. However, the taxon involved in the polyphyly of Nymphonidae in their study, Nymphon floridanum, was not integrated in our datasets due to our sequence selection criteria. The clade uniting Callipallenidae and Nymphonidae is retrieved as monophyletic in all analyses. The long-branch leading to the clade was found in all trees (obtained from separate and combined datasets), suggesting close affinities among members of these two families and their recent diversification. This result was also obtained in Arango and Wheeler [1] and Nakamura et al. [12]. In view of this, their grouping into a single family should be considered as a better taxonomic choice.

The first occurrence of Pycnogonida in the fossil record is that of Cambropycnogon klausmuelleri, a larva from the Upper Cambrian of Sweden [53]. The origin and diversification of the extant pycnogonid groups is more enigmatic, and two periods have been proposed in the literature: Paleozoic or Mesozoic.

A Paleozoic origin of extant pycnogonids is supported by the morphological analyses of Siveter et al. [54] and Arango and Wheeler [1], who concluded that there were close relationships between extant families and Paleozoic fossils, such as Haliestes from the Silurian of England (425 Mya), and Palaeothea and Palaeopantopus from the Lower Devonian of Germany (390 Mya).

More recently, Charbonnier et al. [55] discovered three fossils from the Middle Jurassic of France (160 Mya), which have been assigned to three different extant families: Palaeopycnogonides (Ammotheidae), Colossopantopodus (Colossendeidae) and Palaeoendeis (Endeidae). They suggested that the evolution of pycnogonids proceeded in two steps, with an initial Paleozoic diversification, followed by a Mesozoic radiation that gave rise to the extant families.

Our molecular analyses tend to support the hypothesis of Charbonnier et al. [55]. Indeed, all of our phylogenetic trees show a long branch between the separation of Pycnogonida from other groups of Arthropoda and the common ancestor of extant sea spiders (Figs. 2 and 3). This result suggests that a very long period of time occurred between the origin of Pycnogonida, and the subsequent diversification that led to extant families.