Phylogeny and biogeography of the scaleless scale worm Pisione (Sigalionidae, Annelida)

Abstract Pisione is a scaleless group of small scale worms inhabiting sandy bottoms in shallow marine waters. This group was once considered rare, but now 45 described species can be characterized, among others, by their paired, segmental copulatory organs (one to multiple external pairs), which display a complexity of various accessory structures. The evolutionary significance of these unique organs was suggested in the late 1960s, but has been heavily debated since the late 1990s and remains controversial. In the present paper, we study the internal relationships within Pisione, employing combined phylogenetic analyses of both molecular and morphological data from 16 terminals of Pisione, as well as two terminals of Pisionidens, and eight additional scale worms as outgroups. Our taxon sampling covers all geographical areas where the genus has been reported, as well as most of their morphological and copulatory variability, including representatives of the “africana,” “remota,” “crassa,” and “papuensis” groups, established previously by Yamanishi. We hereby provide a first insight into the relationships of the genus, testing previously proposed hypotheses on the evolutionary significance of male copulatory structures within Pisione, while attempting to understand patterns of distribution. The phylogenetic analyses using maximum likelihood and Bayesian methods consistently recovered two large clades spanning the East Atlantic (including the Mediterranean) and the Indo‐Pacific–West Atlantic, respectively. Character optimization on our trees revealed a high degree of homoplasy in both non‐reproductive and sexual characters of Pisione, with buccal acicula found to be the sole apomorphy among the morphological features assessed herein, with none defining the biogeographical subclades within. Overall, our comparative analyses highlight the high degree of morphological variation in this widely distributed genus, rejecting previous assertions of an increasing number and complexity of copulatory structures across the genus.


| INTRODUCTION
Since the mid-nineteenth century, the placement of the small and aberrant annelid genus Pisione Grube, 1857 has been one of the trials and tribulations. Species of Pisione are unpigmented annelids, only a few millimeters in length, and with well over 50 segments. They are commonly found in sandy bottoms of shallow marine waters (Rouse & Pleijel, 2001), with one exception in freshwater (San Martín, López, & Camacho, 1998). However, their general morphology resembles various annelid groups, which accounts for the numerous suggested systematic affinities. Pisione was until recently one of four genera placed within Pisionidae nomen suppressum. A close association of this group to Aphroditiformia had long been proposed (i.e., Åkesson, 1961;Pleijel & Dahlgren, 1998). It was not, however, until 2005 that molecular and combined molecular and morphological analyses concluded that they are highly derived sigalionids (Struck, Purschke, & Halanych, 2005;Wiklund, Nygren, Pleijel, & Sundberg, 2005). A recent systematic analysis finally synonymized "Pisionidae" with Sigalionidae (Norlinder, Nygren, Wiklund, & Pleijel, 2012), a family within Aphroditiformia that includes scale-bearing annelids with compound chaetae.
Members of Pisione were once considered rare (Hartman, 1959), yet the sheer number of recent descriptions indicates that the number of species will likely continue to increase (Aguado & San Martín, 2004;Martínez et al., 2008;Rouse & Pleijel, 2001).
Based on morphological comparisons, Yamanishi (1998) suggested that Pisione evolved from a Pholoe-like ancestor. Pisione (Fig. 1) share with most sigalionids the presence of compound neurochaetae and a slender and elongated body, but lack scales (=elytra). The loss of scales is hypothesized as one of many adaptations to an interstitial lifestyle (Struck et al., 2005) and is a trait convergently shared by other interstitial scale worms including Metaxypsamma Wolf, 1986; as well as the macrofaunal Palmyra Savigny, 1818 (Watson Russell, 1989;Wiklund et al., 2005;Wolf, 1986). The interstitial lifestyle of Pisione also seems correlated with other changes, including copulation and internal fertilization. This reproductive strategy is commonplace to interstitial taxa (Giere, 2009), but unlike the external fertilization normally found in scale-bearing macrofaunal annelids (Rouse & Pleijel, 2001). Reproductive adaptations are essential for interstitial annelids, especially with limited availability of reproductive products, body size, and space limitation within their environment (Jörger, Heß, Neusser, & Schrödl, 2009;Westheide, 1984;Yamanishi, 1998). Within Pisione, the males display elaborate paired copulatory organs, which may be present in a single segment, or up to 15 or more depending on the species. Yamanishi (1998) found these male copulatory structures to be essential in the classification of the group, while emphasizing that they could be informative for understanding the evolution and even biogeography of Pisione. However, due to immaturity or seasonality, penises may be lacking from the examined collected material, and Salcedo et al. (2015) have suggested that other non-reproductive morphological characters may be systematically informative. While these structures would include characters like neurochaetae and buccal and neuroacicula, to date, no detailed study across taxa has compared the significance of these characters, nor of the copulatory organs. Yamanishi (1998) identified five groups based on a proposed evolutionary trend in male copulatory organization (Fig. 2). The simplest construct of copulatory structures was his "africana" group. According to Yamanishi, the "africana," "remota," and "crassa" groups evolved from an evolutionary progression in which accessory structures were added progressively to the copulatory organ, which consists of a thick and tapering copulatory structure adjacent to the ventral cirrus. Yamanishi (1998) further introduced two additional groups that did not fit into this proposed progression series. Still, his "papuensis" group does exhibit copulatory characters that resemble an intermediate between the "africana" and "remota" groups, whereas his "gopalai" group is characterized by the fusion of the copulatory organ stem and the parapodia, forming a bulging structure surrounded by a hood with spinous papillae. The ventral cirrus is greatly reduced in the "gopalai" group, but well developed in the "papuensis." Regardless, these evolutionary hypotheses were strictly based on observations, and thus far, never investigated by phylogenetic methods.
We here investigate the evolutionary history of the genus Pisione, from both a phylogenetic and a biogeographical perspective. We implemented combined phylogenetic analyses to investigate the character evolution within the genus, while tracing the morphological character evolution on our tree topology. Two long-standing questions were further addressed in our analytical comparisons. First, we compare the degree of homoplasy in both non-reproductive and sexual characters in order to evaluate their diagnostic value for species identification.
Second, we used comparative methods to investigate the detailed evolution of copulatory organs, testing the hypothesis of progressively increasing complexity in copulatory structures as proposed by Yamanishi (1998). And finally, we investigated the optimal distribution range for Pisione, testing for the presence of a latitudinal diversity gradient (LDG; Jablonski, Roy, & Valentine, 2006) and the preference of biogeographical hotspots (Bowen et al., 2013) based on all the records for the genus.   Table 1).

| Specimen collection
All collections were carried out between the intertidal/swash zone and 30 m depth by snorkeling or diving. Specimens were extracted from fine sand to coral/volcanic rubble and gravel sediments. Animals were anesthetized in isotonic MgCl 2 with seawater and extracted using the decantation method with a 63μm mesh (Pfannkuche & Thiel, 1988). Targeted specimens were sorted and identified to genus using a field microscope. Animals used for molecular analyses were preserved in 99% ethanol (EtOH) and stored at −20°C. Vouchers and specimens used for morphological character coding were fixed either in 3% glutaraldehyde or trialdehyde (in 0.1 mol/L cacodylate buffer with 5% sucrose), or in 2%-4% paraformaldehyde [PFA in PBS buffer; following protocols listed in (Kerbl, Bekkouche, Sterrer, & Worsaae, 2015)]. Original descriptions were used for taxonomic identification.
Full species names and detailed collected localities are given in Table 1.

| Morphological examinations
Morphological characters of all newly acquired material were examined using whole mounts prepared with glycerol on an Olympus IX70 inverted microscope mounted with an Olympus DP73 digital camera at the Marine Biological Section, University of Copenhagen, Denmark.
Most morphological character coding was possible using light microscopy (LM) techniques (Table 2).
Specimens requiring detailed examination of copulatory segments, chaetal characters, and other gross anatomy were prepared for scanning electron microscopy (SEM). Specimens prepared for SEM were first transferred to cacodylate buffer, post-fixed in 1% osmium tetroxide (in 0.1 mol/L cacodylate solution) for 1 h, and rinsed in distilled water. Specimens were dehydrated using a graded ethanol series (20%-100%) and transferred over three graded steps to 100% acetone for critical-point-drying. Critical-point-dried specimens were mounted on aluminum stubs and sputter-coated with platinum/palladium using a high-resolution fine coater (JFC-2300HR) and examined using a JEOL JSM-6335F field emission scanning electron microscope at the National History Museum of Denmark, University of Copenhagen.

| Morphological data matrix
Forty-four morphological characters ( Fig. 1; Table 2) were used to construct a morphological data matrix of 16 ingroup taxa, including their sister taxa Pisionidens Aiyar & Alikunhi, 1943, as well as several additional scale worms based on both direct observations and a review of the literature (Salcedo et al., 2015;Yamanishi, 1998). Detailed character descriptions and their states are listed in Appendix 1.
The final matrix of scored characters for all taxa included in this study is listed in Table 3.

| Molecular techniques
Newly generated sequences were deposited in GenBank with the accession numbers KY657611-KY657671 (Table 1).
Individual genes and morphological matrices were concatenated using Sequence Matrix (Vaidya, Lohman, & Meier, 2011). A total of two nested datasets were compiled based on the available information gathered during the analysis, a molecular dataset (MDS), which included 26 taxa exclusively represented by molecular data, and a combined dataset (CDS), which included the same 26 taxa from the MDS combined with 44 morphological characters for each taxon.

| Phylogenetic analyses
Individual gene datasets, 18S rDNA, 28S rDNA, 16S rDNA, and COI, as well as all concatenated gene datasets were analyzed using maximum likelihood (ML) and Bayesian methods (BA).
Maximum likelihood (ML) analyses were computed using RaxML version 7.2.8 (Stamatakis, 2006) as implemented on the CIPRES Science Gateway (Miller, Pfeiffer, & Schwartz, 2010). General time reversible (GTR) model of sequence evolution with corrections for discrete gamma distribution (GTR + Γ) was specified for each partitioned dataset, as this is the only model for molecular evolution available in RaxML. Morphological partitions were analyzed using a Markov model (Lewis, 2001). Non-parametric bootstrapping with 1,000 replicates was used to generate nodal support estimations (Felsenstein, 1985).
Bayesian analyses (BA) were performed using MrBayes version 3.2.5 (Ronquist & Huelsenbeck, 2003) as implemented on the CIPRES Gateway (Miller et al., 2010). Prior to analyses, jModelTest (Posada, 2008) was used for all multiple sequence alignments (MSA) of individual genes to infer their optimal evolutionary model estimated by the corrected Akaike information criterion (AICc). The models selected for each gene included a GTR model with gamma distribution and a proportion of invariable sites (GTR + I + Γ) for 18S rDNA and 28S rDNA.  partition was analyzed using a Mk1 model (Lewis, 2001). All individual gene datasets and concatenated MDS and CDS were run with two independent analyses using four chains (three heated and one cold).
Number of generations was set to 30 million, sampling every 1,000 generations. Burn-in was set to 10 million generations. TRACER version 1.6.0 (Rambaut & Drummond, 2007) was used to verify convergence of all the MCMC runs. Majority-rule consensus trees (50%), posterior probabilities, and branch lengths were constructed with the remaining trees after burn-in.

| Species delineation
As mentioned above in Taxon selection, species identification of terminals from the same region with similar (or lacking) morphological traits were investigated using species delineation tests. Delineations were made on those terminals of putative similar identification using ultrametric trees obtained in BEAST (see below). Only those terminals corresponding to Pisione were utilized in the two most commonly used methods of species delineation (Fontaneto, Flot, & Tang, 2015): the generalized mixed Yule coalescent model (GMYC) (Fujisawa & Barraclough, 2013), and the Poisson tree process (PTP) (Zhang, Kapli, Pavlidis, & Stamatakis, 2013). For all methods, outgroups were excluded from the analyses.

| Character evolution tracing
Character evolution was traced using parsimony methods computed with Mesquite version 3.02 (Maddison & Maddison, 2007) and MacClade version 4.0 (Maddison & Maddison, 2000). Most parsimonious reconstructions (MPR) were computed on the tree recovered from the BA of the CDS. MPR, when required, were compared using both the accelerated transformation parsimony model (ACCTRAN) and delayed transformation parsimony model (DELTRAN) available in MacClade.

| Testing relative consistency of non-reproductive versus sexual characters
Most studies on Pisione have suggested that the number and structure of male copulatory organs relative to other non-reproductive character traits show more systematic importance, providing the best basis for morphological species identification. However, a recent review by Salcedo et al. (2015) indicated that certain chaetal and parapodial features should be equally considered. All morphological characters showed a high degree of homoplasy in our character tracing (see results), supporting our decision to test whether, as previously suggested, more emphasis should be placed on copulatory structures than other character traits when performing systematic studies of Pisione. We compared the overall consistency index of all the sexual and non-reproductive morphological characters included in our analyses that exhibit variations within Pisione. All sexual and nonreproductive characters scored by Yamanishi (1998)  Characters that may be associated with groups designated by Yamanishi (1998), and scored for species included in our dataset, are listed in superscript and include: AG, "africana" group; RG, "remota" group; CG, "crassa" group; PG, "papuensis" group; GG, "gopalai" group.
T A B L E 2 (Continued) non-reproductive as explanatory characters by the number of steps.
The number of steps is considered as discrete counting of the data, for which a Poisson distribution was used as the candidate to fit the model. Boxplots were plotted, and the z-scores generated by the GLM were used to help with the interpretation.

| Evolution of copulatory organs
Increasing complexity of copulatory organs within the evolution of Pisione has been previously proposed (see Salcedo et al., 2015;Yamanishi, 1998) based on the presence of different morphotypes in different geographical areas. While our dataset does not include all known species of Pisione, it does offer an adequate framework to assess the evolution of complexity of copulatory organs. Hence, we tested whether the changes in complexity and number of copulatory organs in Pisione were better explained by body size (as in other groups); better explained based on geographical region as proposed by Yamanishi (1998); or, in opposition, phylogenetically constrained.
The number of copulatory organs per species was coded as the maximum number of copulatory segments for each species, while copulatory organ complexity was characterized by the number of accessory structures present for each of the species (characters 37-44; Table 3).
We characterized body size as the maximum body length and the Analyses were run with independent MCMC chains, which were set for 100 million generations. Sampling was set every 10,000 generations. Convergence of the reads was confirmed using Tracer version 1.5 (Rambaut & Drummond, 2007). The consensus tree based on maximum clade credibility (MCC) was obtained using TreeAnnotator version 1.6.1 with a burn-in of 20% (Drummond et al., 2012).
The hypothesis that number and complexity of copulatory organs are phylogenetically constrained was evaluated using Pagel's lambda (λ) and Blomberg's K indices with the function "phylosig" as implemented in the R package phylotools (Revell, 2012). Values close to 1 provide a high indication that the characters are phylogenetically constrained, whereas values near zero indicate that the character is highly variable between closely related species.
Alternatively, we investigated whether there was evidence of coevolution of both body size and the number and complexity of copulatory organs, while taking into account the effect of geographical areas. This was evaluated using the phylogenetic general least square methods (PGLS). The logarithm of the maximum number of copulatory structures and copulatory complexity were the response variables, treated as continuous, with a Gaussian distribution. The logarithm of maximum body length, the logarithm of the maximum number of segments, and the geographical area were selected as explanatory variables. All models were investigated using Monte Carlo Markov chains within the package MCMCglmm (Hadfield, 2010), which was run for 5.1 million generations. The first 100,000 generations were discarded as burn-in, and the thinning parameter was set to 500.
The evolution of maximum number of copulatory organs, copulatory organ complexity, maximum body length, and the maximum number of segments was estimated on the obtained ultrametric tree to visually illustrate our results. Each continuous character was traced using the function contMap (Revell, 2013) implemented in the R package phylotools (Revell, 2012).

| Analyses of geographic distribution patterns
Geographical distribution patterns were analyzed from geographic collection data for all 55 described and unidentified species of Pisione.
Latitude preferences for Pisione were inferred for tropical, subtropical, and temperate zones using the correlation between the numbers of species recorded at each 20 degrees of latitude in order to predict if there is a latitudinal gradient (LDG), where, for example, the biological diversity is increasing from the poles to the tropics (Jablonski et al., 2006). We also predicted that Pisione diversity is correlated to the Caribbean and Indo-Pacific Oceans, usually interpreted as biogeographical hotspots for several metazoan groups (Bowen et al., 2013).
We counted the number of species for latitude in steps of circa 10 degrees. Our first set of ranges started at 0 degrees to 9.9, from 10.0 to 19.9, from 20.0 to 29.9, etc., in the same manner, and following the same approach for the remaining latitudes. This was repeated and also used for the negative latitudes. Pisione preferences for the Indian Ocean, Pacific Ocean, and the Atlantic Ocean were tested using the correlation between the numbers of species recorded in each longitudinal subgroup every 40° of longitude following similar approaches to those for the latitude. The maximum diversity of Pisione was used to infer latitudinal and longitudinal preferences.
Geographical ranges with maximum diversity were estimated using spline models (Anderson, 2008;Di Domenico, Martínez, Lana, & Worsaae, 2014;Koenker, 2005;Koenker, Ng, & Portnoy, 1994) built for the 95th percentile. Models were fitted using the functions rq() and bs() ["quantreg" package for R; (Koenker, 2007;R Core Team, 2014)]. The degree of the polynomial was chosen by the AICc corrected for small samples (Burnham & Anderson, 2003;Hurvich & Tsai, 1989). The model with the smallest AICc value from a set of models with a degree of polynomial was selected, and optimal values were interpreted. Ninety-five percent bootstrap confidence intervals (Manly, 2006) were obtained for each of the 10k bootstrapped sample pairs using the polynomial degrees that were chosen for the original data.

| Phylogenetic analyses
The Bayesian and maximum likelihood analyses recovered similar topologies between the MDS and CDS analyses and are presented in Fig. 3 with both Bayesian posterior probabilities (BPP) and maximum likelihood bootstrapping (MLB) values. However, between the methods, slight differences were recovered and are described below.
The terminals representing the genera Pisione and Pisionidens formed reciprocally monophyletic clades, both with high support (Fig. 3).

| Species delineation
In the absence of diagnostic morphological penile traits or in cases of damaged specimens, species delineation tests made it possible to distinguish the least taxonomic units. While pairwise genetic distances may resolve such issues, we preferred to implement more accurate methods given the high interspecies variability within Pisione (Fontaneto et al., 2015). Species delineation techniques employed recovered a single species representing the Mediterranean-East Atlantic region, P. puzae, with representatives from both Italy and Galicia. In the West Atlantic region, P. wolfi was found throughout the Cuban coastline, and P. hartmannschroederae was recovered with large distributions throughout the West Atlantic.

| Character evolution tracing and analyses
The relationship of the members within the clade Pisionidens -Pisione was supported by five synapomorphies, traced without homoplasy (see     Overall, the level of homoplasy between sexual and nonreproductive characters was the same, reflected by comparable distribution of the number of steps across both groups of characters ( Fig. 4), of which z-scores calculate by the GLM did not yield any significant differences (z = 1.09; p = .28). Although no differences were observed between sexual and non-reproductive characters, the intercept was different from zero steps (z = 2.19; p = .028), indicating some level of homoplasy equally predicted by both characters. Tree topology based on the Bayesian analysis (BA) of the combined dataset. Only nodal support above BPP = .5 and MLB = 50 is displayed. Nodes not recovered or with low support are illustrated with a dash (-). Triangles with asterisks inside represent maximum support in all analyses, while a single asterisk (*) denotes maximum support in a specific analysis (BPP = 1.0 or MLB = 100). Boxes on branches identify apomorphies, including character number and states in parentheses. Full character coding for all terminals can be found in Table 2. (a) Schematic representation of Pisionidens, including prostomial appendages and first few segments (redrawn from Rouse & Pleijel, 2001). (b) Diagrammatic representation of Pisione, including prostomial appendages and first few segments (redrawn from Rouse & Pleijel, 2001)

| Evolution of copulatory organs
The maximum number of paired copulatory structures and their complexity, estimated as the number of accessory copulatory structures, showed low phylogenetic signal (number of copulatory organs λ = 0.00, K = 0.241; complexity of copulatory organs λ = 0.00; K = 0.252). PGLS was unable to show any relationship between maximum number of copulatory organs, complexity of copulatory organs, or any of the explanatory variables. The number of paired copulatory organs is reduced independently in two clades, including P. puzae clade and Pisione papuensis -Pisione hartmannschroederae, while increasing in Pisione sp.
A. The complexity of the copulatory structures is reduced once in P. bulbifera, but increases independently toward P. guanche and P. vestigialis.

| Geographic distribution patterns
Geographic analyses yielded a well-supported diversity gradient of Pisione, with a maximal diversity estimated at 20° latitude, with a range between −20° and 30° latitude within a 95% confidence interval. A steep decrease in diversity was present in latitudes >−20° and >30° toward both poles (Fig. 5).
The comparison of the longitudinal distribution patterns showed no significant optimal values and is not illustrated.

| Phylogenetic relationships of Pisione
Pisione was well-supported, recovered sister to Pisionidens, and with buccal aciculae as a morphological synapomorphy. We recovered two clades within Pisione supported by comparatively high nodal support values, but lacked identifiable synapomorphies. Comparably, both sexual and non-reproductive morphological characters were traced with high levels of homoplasy. Similar patterns of character evolution with numerous homoplasious characters have also been found in other lineages of interstitial annelids, including Saccocirridae (Di Domenico et al., 2014), Protodrilidae (Martínez, Di Domenico, Rouse, & Worsaae, 2015), and Nerillidae (Worsaae, 2005(Worsaae, , 1883. However, generally in these families, the major internal clades could more easily be diagnosed by unique combinations of morphological characters (Di Domenico et al., 2014;Martínez et al., 2015). an East Atlantic and Western Atlantic-Indo-Pacific clades . The similar distribution patterns between Pisione and Megadrilus might be related to the presence of a common dispersal strategy, as they both are comparatively large species with pelagic larvae. However, these similar distributions may just be due to common vicariant processes or sampling biases that have artificially produced similar distribution patterns. Yamanishi (1998) argued that reproductive characters, specifically the male copulatory organs and accessory structures, are the most important morphological diagnostic characters for describing and identifying Pisione species. However, given that male sexual maturity is seasonal within Pisione, mature males are often lacking from collections. This scenario has complicated descriptions, and those descriptions lacking males usually rely on size comparisons or on extremely divergent characters. Salcedo et al. (2015) argue that additional characters are needed in conjunction with male copulatory structures, especially those that are expressed regardless of sexually maturity or season. These characters should include type and number of neurochaetae, shape of the dorsal cirri on segment three, neuropodial lobes, and shape, size, and ornamentation of the neuro-and buccal acicula. While these character transformations rarely constitute apomorphies, they may be of systematic importance when combined, as several species of Pisione either lack aciculae (neuro-or buccal-) or may exhibit intraspecific variation within neurochaetal numbers and patterns.

| Character evolution tracing and analyses
Building atop of these previous ideas, when we compared sexual characters to that of non-reproductive characters, we find that there are no comparable differences in the degree of homoplasy between them. Both sexual and non-reproductive characters lack both clade and regional specificity. While the results do not refute Yamanishi (1998)

| Evolution of copulatory organs
While we have shown that sexual characters are equally important as non-reproductive, this is the first comparative investigation into the organization and complexity of copulatory structures in Pisione based on comparative phylogenetic methods. While attempting to address both species level identification and evolutionary and biogeographical patterns, Yamanishi (1998) explained that it was possible to group species based on their degree of complexity in male copulatory structures. While morphological phylogenetic analyses were not performed by Yamanishi (1998), five groups based on copulatory complexity were illustrated (Yamanishi, 1998; fig. 22), providing a framework for our current comparative investigations. Three of the illustrated groups were proposed to represent an evolutionary succession of increasing complexity ("africana" → "remota" → "crassa"); however, there were two groups that were characterized separately: a "papuensis" group that shared mixed similarities of the "africana" and "remota" groups and may potentially represent a hybrid or intermediate group, and a "gopalai" group that is highly modified and unlike the other groups, with fusion of several characters and lacking any protruding structures. While biogeographical patterns were not emphasized by Yamanishi (1998), he willfully excluded several of the known species at that time from his construction of his copulatory complexity progression scheme, including several members from within the Western Atlantic (including the Caribbean).
Present study included representatives of the proposed successive "africana" (P. guanche), "remota" (P. remota, P. puzae), and "crassa" (P. cf. vestigialis) groups, all having an elongated copulatory organ stem, as well as the "papuensis" (P. papuensis) group, with stem of copulatory organ broad and partially fused to the parapodia. Our comparative analyses reject the evolutionary trend of both an increase in number of copulatory organs and an increase in complexity of copulatory structures across represented species of Pisione (Fig. 3). Within the phylogeny, there were multiple instances of both increasing and reduction of complexity and number of copulatory structures. More significantly, our phylogenetic character tracing revealed a high degree of variability and homoplasy in the characters associated with copulatory and accessory structures, reflected by the paraphyly of, at least, two of Yamanishi's copulatory groups (i.e., "crassa" and "remota" groups).
Taking into account the limitation of our dataset, our results suggest that it is highly unlikely that even the addition of more taxa will change the general findings to support Yamanishi's theory of increasing complexity. Furthermore, when we consider our ancestral character reconstructions, the Pisione ancestor was estimated to bear an intermediate number of copulatory organs with a medium complexity of accessory structures. This is in direct opposition to Yamanishi's (1998) proposal that the ancestral state would represent a large bodied Pisione with simple male copulatory structure resembling more closely that of the "africana" group.
The evolution of copulatory organs within Pisione appears highly convoluted. Within our East Atlantic clade-1, P. guanche exhibits copulatory structures similar in form to that of the "africana" group, having an elongated copulatory organ proper, lacking whorls, with elongated ventral cirrus with unmodified neurochaetae (Martín et al., 1999;Yamanishi, 1998). In our analyses, this so-called primitive group (according to Yamanishi, 1998) was recovered as the sister taxon to Pisione remota -Pisione puzae clade. This clade displays copulatory organs related to the "remota" group, which bears a bidigitate process on the main stem of the copulatory organ proper (Yamanishi, 1998). When we investigated the morphology with both SEM and that of descriptions by Yamanishi (1998), a slight evolutionary progression similar to what was originally proposed by Yamanishi can be seen, especially when comparing P. guanche to Pisione remota -Pisione puzae. However, these three species are not closely related to P. vestigialis or P. wolfi, the two species included in our analyses belonging to the "crassa" group, which display highly ornamented copulatory structures. Furthermore, copulatory structures associated with the "remota" group (i.e., P. remota and P. puzae) are known to display additional accessory structures not present in the "africana" group.
However, our character reconstructions show that the homoplasious characters 42-44 are present in our Pisione representatives from both the "africana" and "remota" groups. Additional shared characters (45-46) present in the Pisione remota -Pisione puzae clade ("remota" group) are also found in terminals from clade-2 that further invalidate Yamanishi's assumptions of the "remota" group being of single origin from the "africana" group. In the Indo-Pacific-West Atlantic (clade-2), species of Pisione exhibited copulatory structures indicative of the "papuensis" group, the "crassa" group, and that of the "remota" group. Pisione cf. vestigalis and P. wolfi both share copulatory structures resembling that of the "crassa" group; however, these species are each nested within separate subclades and geographical regions, respectively, hereby also negating the monophyly of Yamanishi's "crassa" group. Pisione papuensis is found to be the sister taxon of a clade containing P. bulbifera, that in itself bears copulatory structures resembling the "remota" group. The close association of both the "papuensis" and "remota" groups is not surprising given the similarly enlarged copulatory organ proper and ventral cirri, contributing to close morphological affinities. Yet, again, the close association of these species breaks up the linear evolutionary scenario of Yamanishi ( Fig. 2a), which is further compromised by the two unrelated origins of the "crassa" forms (P. cf. vestigalis and Pisione wolfi -Pisione hartmannschroederae) within clade-2. This varying morphology in our analyses seen by our recovered clades (Fig. 3) is attributed to several homoplasious characters observed during multiple species examination with detailed microscopy. It appears that the more complex evolutionary scenario found herein may already have been suspected by Yamanishi, who himself excluded several species from his classification of male copulatory structures.

| Distribution patterns
Outside of taxonomical descriptions, no other study attempts to describe biogeographical patterns within Pisione. In an attempt to understand patterns of diversity, Yamanishi (1998) was the first to address the widespread distribution ranges and diversity in Japan.
Aside from the more limited transfer by storms or other migrating sand events, these wide range dispersal events may be due in part to a relatively long planktonic stage and/or perhaps dispersal by anthropogenic means. Pisione larvae drift about the pelagic, capable of collecting food by mucoid or slime nets, and can last up to 10 days in search of suitable substrate (Åkesson, 1961). The patterns of oceanic currents coupled with long planktonic larval stages suggest this as the most plausible means of dispersal within Pisione. While many other dispersal vectors have been proposed for additional groups of interstitial annelids (see Weidhase, Bleidorn, & Simon, 2016), ballast sand or water may also play a role in the dispersal of Pisione individuals.
Our analyses revealed two large clades within Pisione, one that appears to be restricted to the East Atlantic, including the Mediterranean, and the one spanning from the Red Sea, across the Pacific, and into the West Atlantic including the Caribbean. To date, Pisione has always been described as being most commonly distributed within tropical and subtropical areas (Aguado & San Martín, 2004;Salcedo et al., 2015;Yamanishi, 1998). Our geographic analyses confirm this statement, showing maximum diversity across the latitudes −20° to 30° (Fig. 4). We considered that the number of sampling efforts and records might generate false-positive correlations between diversity, and that of latitude and longitude preferences; however, we attempted to address this issue by using a confidence interval for the optimal value (peak of the unimodal curve, 20°) to fit the curve. This method, albeit exploratory, appears to have mitigated our concerns, as the historically high sampling efforts in European waters were not displayed as the highest estimates of diversity, nullifying any false-positive correlation. Furthermore, over half of the described species, including all of the subspecies, are reported from within the West Pacific and Indian Oceans (Salcedo et al., 2015). The effect of area in regard to high diversity at lower latitudes (tropics) has been historically debated in ecology. Given our distribution patterns (Fig. 5) of species richness within Pisione, at first glance it appears to possibly resemble that of the so-called mid-domain effect [MDE; a controversial null model that generates patterns by random overlap of geographic ranges in low latitudes (Colwell & Lees, 2000;Hawkins, Diniz-Filho, & Weis, 2005)].
While addressing this issue was not the focus of this manuscript, our implementation of the bootstrap (implemented when estimating the optimum peak) provides an additional tool to assess the uncertainty of our estimates and our sampling efforts (Kulesa, Krzywinski, Blainey, & Altman, 2015). While we agree that this method lacks a traditional statistical inference, it may represent a more rational way to find and fit models than the traditional linear predictive models that are questionable when it comes to low sampling efforts and area effects.
Unfortunately, obtaining significant optimal values in our longitudinal distribution patterns was not possible.  (Steininger & Rögl, 1984) as shown for corals (Budd, 2000). However, lack of fossils for calibration prevents aging of Pisione clades, and the observed distribution patterns and diversity may very likely just reflect sampling bias in this first phylogeny of Pisione.

| CONCLUSIONS
Our phylogenetic analyses revealed a well-supported monophyletic Presence of a single pair (per segment) of male copulatory organs are found in Pisione and Pisionidens (e.g., Day, 1967;Martín et al., 1999). The number and ornamentation of paired male copulatory segments vary with taxa. Copulatory segments lack dorsal cirri; however, an elongation and/or modification of ventral cirri occurs in some taxa of Pisione. Terminals with only female specimens were coded as present due to separate sexes observed in both the infra-acicular chaetae is the third chaetae, counting from the dorsal-most chaetae on the neuropodia (Yamanishi, 1998