Molecular-based estimate of species number, phylogenetic relationships and divergence times for the genus Stenotaenia (Chilopoda, Geophilomorpha) in the Italian region

Abstract Stenotaenia is one of the largest and most widespread genera of geophilid centipedes in the Western Palearctic, with a very uniform morphology and about fifteen species provisionally recognized. For a better understanding of Stenotaenia species-level taxonomy, we have explored the possibility of using molecular data. As a preliminary assay, we sampled twelve populations, mainly from the Italian region, and analyzed partial sequences of the two genes COI and 28S. We employed a DNA-barcoding approach, complemented by a phylogenetic analysis coupled with divergence time estimation. Assuming a barcoding gap of 10–16% K2P pairwise distances, we found evidence for the presence of at least six Stenotaenia species in the Italian region, which started diverging about 50 million years ago, only partially matching with previously recognized species. We found that small-sized oligopodous species belong to a single clade that originated about 33 million years ago, and obtained some preliminary evidence of the related genus Tuoba being nested within Stenotaenia.


Introduction
Stenotaenia Koch, 1847 is one of the largest and most widespread genera of geophilid centipedes occurring in the Western Palearctic. Species of Stenotaenia have been recorded mainly from the Italian region, through the Balkan peninsula and the Aegean islands, to Anatolia (Bonato and Minelli 2008). However, there are also sparse records reaching the British isles, to the west (Barber 2009), the Alborz mountains in Iran (R. Zarei, unpublished data), to the east, and the Atlas chain plus Crete, Cyprus and Israel, to the south (Chipman et al. 2013, Simaiakis et al. 2013. Like in most other centipedes, taxonomic recognition and delimitation of species in Stenotaenia have so far been based on scanty morphological evidence. However, all Stenotaenia species are remarkably similar in body anatomy and several fine morphological details, including most of the characters that are traditionally considered diagnostic at the species level in other geophilid genera (e.g., details of the labrum and the maxillary complex, shape and denticulation of the forcipules, arrangement of the sternal pore areas along the trunk, structure and shape of the legs of the ultimate pair and the associated metasternite, and arrangement of the coxal pores). On the contrary, Stenotaenia exhibits high variability in the number of trunk segments and adult body size, so that species-level current taxonomy is based almost exclusively on these two characters. Extreme morphologies are represented by S. romana (Silvestri, 1895), which is reported being less than 17 mm long when fully grown, with some specimens having only 43 leg-bearing segments, and S. sturanyi (Attems, 1903), reaching 77 mm in length, with specimens having up to 115 leg-bearing segments (Bonato and Minelli 2008). According to a recent critical reassessment of all available morphological information, no more than fifteen species have been recognized in Stenotaenia and more than half of these should be considered valid only provisionally, because supporting evidence is inadequate (Bonato and Minelli 2008).
For a better understanding of Stenotaenia species-level taxonomy, we have explored the possibility of using molecular data. As a preliminary assay, we analyzed two genes in a sample of populations from a significant part of the geographic range of the genus. Our aim was to estimate how many species could be recognized on the basis of DNA sequences, especially in comparison with the taxonomic scheme currently in use. We adopted a DNA-barcoding approach, complemented by a phylogenetic analysis coupled with divergence time estimation.

Sampling and DNA extraction
The target of our study was the genus Stenotaenia according to the taxonomic concept and circumscription currently in use (Bonato and Minelli 2008).
Our sampling focused on the western part of the known range of the genus. This area is centred on the Italian region s.l., which extends from the Alps and Istria, through the entire Italian peninsula, to the Italian islands (Fig. 1). So far, this area has been investigated much more intensively than the remaining range of the genus (Bonato and Minelli 2008) and it represents the only part of the range for which a consistent taxonomic scheme for Stenotaenia has been developed. We sampled ten specimens from ten different populations within this area and two specimens from localities at the eastern borders of the known range of the genus, i.e. from Cyprus and from Iran ( Table 1). Most of the specimens were obtained directly by recent sampling in the field, because most of Greek letters refer to the species tentatively recognized after the analyses (see text). Table 1. Sampled specimens of Stenotaenia, arranged west to east and then north to south. The preliminary identification is just a tentative one, based only on the few morphological characters hitherto proposed as diagnostic at the species-level, including number of legs and geographical provenance (Bonato and Minelli 2008 the material already available in collections has proved unsuitable for DNA extraction, amplification and sequencing; in fact, the DNA of materials stocked in ethanol 70% or in other solutions used for morphological studies (i.e., lactophenol) easily degrades. After a preliminary species identification based on morphology (Table 1) following Bonato and Minelli (2008), the specimens were fixed in absolute ethanol to preserve DNA from degradation. Total genomic DNA was extracted from a dissected intermediate portion of the trunk of each specimen, using the DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany).

DNA amplification and sequencing
We sequenced a portion of the cytochrome c oxidase subunit I (COI) using the primer pair LCOI490/HCO2198 (Folmer et al. 1994), which amplifies an approximately 700 bp long fragment, and a portion of 28S rRNA (28S) using the primer pair 28SD1F/28SrD4b (Boyer andGiribet 2007, Edgecombe andGiribet 2006), which amplifies a 1050 bp long fragment. For COI amplification we employed the following PCR conditions: an initial denaturation step at 95 °C for 5 min, followed by a variable number of cycles (27-35) including denaturation at 94 °C for 1 min, annealing (ranging from 40 to 47 °C) for 1 min and extension at 72 °C for 1.5 min, then a final step at 72 °C for 7 min. For 28S amplification we applied the same thermal cycling profile except for the annealing step, ranging from 42 to 50 °C for 1 min. Each PCR product was screened for the potential successful amplification by electrophoresis on 1% agarose gel in 1X TAE and purified using MinElute PCR Purification Kit (Qiagen). Then it was directly sequenced on both strands with the same primer sets as used for amplification, using an ABI 3130 XL automatic capillary sequencer (Applied Biosystems, Branchburg, USA; service provided by BMR Genomics, Padova, Italy).
We obtained COI sequences for all 12 specimens of Stenotaenia (between 642 bp and 647 bp long) and 28S sequences from 11 Stenotaenia specimens (between 952 and 1005 bp long) ( Table 2).
Finch TV 1.4.0 (Geospiza, PerkinElmer) was used to check each chromatogram for nucleotide signal intensity and whole sequence signal strength. Sequences were edited manually to obtain a more accurate reading. For each specimen, forward and reverse sequences of the gene were aligned with default parameters with Clustal W2 (Larkin et al. 2007) and combined in a single sequence. After that, datasets of COI and 28S of Stenotaenia specimens were aligned with the same procedure.

Analyses of sequence variation
For delimiting species in our sample through a DNA-barcoding approach, we calculated the pairwise distances of the COI sequences between all 12 Stenotaenia specimens in two alternative ways, namely by K2P distances (which is a standard for DNA-barcoding; e.g., Hebert et al. 2003, Ratnasingham and Hebert 2007, Chevasco et al. 2014) and p-distances (which is also used in arthropods; e.g., Montagna et al. 2013), treating gaps with partial deletion and estimating standard errors by 500 bootstrap pseudoreplicates. The software MEGA 6 (Tamura et al. 2013) was employed. We analyzed the frequency distribution of both kinds of pairwise distances in order to recognize the putative barcoding gap. This is the interval which is expected to separate within-species distances from between-species distances in the distribution of pairwise distances (Hebert et al. 2003, Ratnasingham andHebert 2007).

Phylogenetic analysis
For the phylogenetic analysis of the Stenotaenia sequences, we chose as outgroups six species of Geophilidae (Table 2) for which sequences of COI and 28S were already available. These species were selected as representative of the genera reputed to be the most strictly related to Stenotaenia according to previous phylogenetic analyses ( Murienne et al. 2010, Bonato et al. 2014) and, as far as possible, including the type species of these genera. Because compositional biases may interfere with the phylogenetic signal of the dataset (Rota-Stabelli and Telford 2008), we checked whether the outgroup sequences were ingroup-like in the GC-content and the GC-skew.
The full set of ingroup and outgroup sequences were aligned using Clustal W2 (with default parameters) and only the positions shared by all sequences were considered for the analyses. The single genes and the concatenated sequences were analyzed by a maximum likelihood (ML) approach with PHYML 3.1 (Guidon et al. 2010) and by a Bayesian approach using BEAST v1.7.2 (Drummond et al. 2012). We estimated the best-fitting replacement model according to the Akaike information criterion (Posada and Buckley 2004) implemented in jModelTest v2.1.1 (Darriba et al. 2012). The ML analyses were bootstrapped using 100 pseudoreplicates. Newick output trees were visualized and manually rooted (assuming, as far as possible, the monophyly of the ingroup) with FigTree v1.4.0 (Rambaut 2009).

Molecular dating
To estimate divergence dates between Stenotaenia sequences, we used the Bayesian method implemented in BEAST v1.7.2 with XML input files prepared using BEAUti v1.7.2 (Drummond et al. 2012). Since no fossil record exists for any of the taxa represented in our dataset, we relied on both a root prior and a replacement rate derived from previous date estimates. As root prior we used an estimate of 200 million years (Ma), with a permissive SD=50, for the basal split in the phylogeny obtained, between the most basal outgroup and the remaining taxa, according to the dated phylogeny of Murienne et al. (2010). As replacement rate we assumed a value of 0.0016 substitutions/site/Ma, with SD = 0.0010, which was previously estimated analyzing a 18S+28S dataset from a wide range of arthropods, including myriapods (Rota-Stabelli et al. 2013). The choice of this rate was motivated by a lack of more taxon-specific rate for centipede 28S. We modelled the molecular replacement with the best-fitting model selected by the Akaike information criterion, and the rate distribution using a random clock model. All other priors were those of BEAST at default settings.

Pattern of sequence diversity
Considering the pairwise distances of COI between the sampled specimens of Stenotaenia (Table 3), the minimum distance (about 0.5% for both the K2P and the pdistances) was recorded between two specimens from the Western Alps (Valdieri, Vernante), which had been previously recognized as belonging to a single species based on morphology ( Table 1). The maximum distance (27.4% for the K2P distances, 22.6% for the p-distances) was estimated between a specimen from Cyprus and a specimen from the Eastern Alps (Barbarano), previously recognized as belonging to distinct species based on morphology (Table 1). Considering the frequency distributions of the pairwise distances of both models (Fig. 2), two major intervals of zero frequency might be recognized: a lower one (0.5-6.5% for K2P, 0.5-6.2% for p-distances) and a higher one (10.4-16.1% for K2P, 9.5-14.4% for p-distances). Assuming a barcoding gap corresponding to the lower observed gap, we would obtain 11 species of Stenotaenia from our sample of 12 specimens, of which nine species in the Italian region (between Alps, Istria and Sicily), with only two specimens from the Western Alps resolved as conspecific. On the contrary, assuming a barcoding gap corresponding to the higher observed gap, we would obtain eight species of Stenotaenia, of which six in the Italian region.

Phylogenetic relationships
After the alignment of ingroup and outgroup sequences, the average nucleotide composition of COI turned to be A = 0.289, C = 0.248, G = 0.163, T = 0.299, and that of 28S resulted as A = 0.219, C = 0.283, G = 0.317, T = 0.181. A bias against G-C in the composition of COI has been observed in other Chilopoda as well (Spelda et al. 2011, Oeyen et al. 2014). The GC-contents (Table 2) of the sequences of the outgroups turned out being within the range of variation of Stenotaenia sequences (38.2-46.8% for COI, 57.5-62.1% for 28S), with minor exceptions for A. glacialis and G. alpinus (slightly higher for COI, slightly lower for 28S). Also the GC-skew values ( Table 2) of most of the outgroup sequences turned out within the range of variation of the Stenotaenia sequences (between -0.246 and -0.136 for COI, between 0.040 and 0.069 for 28S), but with some exceptions (slightly lower for COI in A. glacialis, higher for 28S in G. alpinus and G. electricus).
For the ML phylogenetic analysis, the Generalized Time Reversible model with proportion of invariable sites and a Gamma distribution (GTR+I+G) with four discrete categories was selected as the best-fit model for nucleotide substitution with the Akaike information criterion. We applied this model for the datasets of the COI sequences and the 28S sequences, when analyzed separately and when concatenated.
The ML tree, obtained from the concatenated sequences of the 11 Stenotaenia specimens from which we got workable sequences for both genes, is shown in Fig.  3. Very similar topologies for the ingroup were found in the ML trees obtained from the single genes separately and in the Bayesian tree obtained from the concatenated sequences (trees not shown, but node supports shown on the ML tree, Fig.  3). Three groups of Stenotaenia specimens were found to be strongly supported (at least when analyzing the two genes together) and corresponded to groups obtained by the DNA-barcoding analysis assuming the higher barcoding gap (10.4-16.1% for K2P, 9.5-14.4% for p-distances). These are: a group including three specimens between the Eastern Prealps and the northern part of the Italian peninsula (Barbarano, Frasassi, Isola Fossara), a second group with the two specimens from the southern part of the Italian peninsula and Sicily (Frosinone, Randazzo), and a last group with both specimens from the Western Alps (Vernante, Valdieri). In addition, we found strong support for a clade including the first group listed above and a specimen from Istria (Lovran), and for another clade including the other two groups already listed plus another specimen from the Eastern Prealps (Giavera) and a specimen from Iran.
However, the monophyly of the genus Stenotaenia, as currently circumscribed, was not recovered in our analyses: the specimen representative of the genus Tuoba was found nested within the clade encompassing all Stenotaenia specimens in both ML and Bayesian analyses, and also in the ML analysis of 28S, even though in different positions (Fig. 3); in the ML analysis of COI, instead, the specimen representative of the genus Clinopodes was nested within Stenotaenia.

Divergence time
By applying a molecular clock on a Bayesian analysis of the 28S sequences, we obtained a tree topology (Fig. 4) largely consistent with ML and Bayesian trees obtained from the concatenated dataset (Fig. 3). Differences are limited to the unstable position of Tuoba and two nodes with weak support.
The common ancestor of all Stenotaenia representatives (including Tuoba if it comprises a monophyletic group together with Stenotaenia) was dated at about 64 Ma. The divergence between the two major groups of Stenotaenia present in the Italian region (one represented by specimens from the northern part of the Italian peninsula, besides from the Eastern Prealps and Istria; another represented by specimens from the southern part of Italian peninsula and Sicily, besides from the Western Alps and the Eastern Prealps) was dated at about 51 Ma. Within the first group, the lineage of the specimen from Istria diverged from the remaining lineages around 33 Ma. The estimated age of the subsequent divergences are all younger than 10 Ma. In particular, for the groups obtained by the DNA-barcoding analysis of COI sequences assuming a higher barcoding gap and well supported in the phylogenetic analyses (see above), the . Bootstrap values < 50% and posterior probabilities < 0.50 are not shown. Circles indicate ingroup nodes that are highly supported in the tree based on concatenated sequences. Terminal node groupings indicated by Greek letters refer to the species tentatively recognized (see text and Fig. 1). The specimen from Volpago (species δ) is absent because its 28S sequence was not obtained. common ancestor of each group was dated younger than 1.6 Ma, whereas all divergences between these groups and other specimens of Stenotaenia were all dated older than 3.4 Ma.

Putative species of Stenotaenia
The results of the DNA-barcoding analysis of COI sequences, of the phylogenetic reconstructions on the basis of the two genes, and of the divergence date estimation based on the 28S sequences, all agree in suggesting that at least six species of Stenotaenia are recognizable in our sample from the Italian region ( Fig. 1): α (three populations between the Eastern Prealps and the northern part of the Italian peninsula), β (two populations from the southern part of Italian peninsula and Sicily), γ (two populations from Western Alps), δ and ε (each represented by a single population in the Prealps), ζ (one population from Istria). Additionally, the two populations from Cyprus and Iran probably belong to different species. Estimates of divergence time, calculated using 28S sequences and two priors (age of the root and substitution rate) in the package BEAST v1.7.2 (see text). 95% High Posterior Density intervals are represented by coloured bars for the most robust nodes, emphasized by a circle.
Greek letters refer to the species tentatively recognised (see Fig. 1). The tree has the same topology of the concatenated ML tree of Fig. 3, but for the position of Tuoba sydneyensis and the relationships within the group formed by species β, ε, γ and the specimen from Iran. The specimen from Volpago (species δ) is absent because its 28S sequence was not obtained. Time scale is different in the two intervals 0-100 and 100-200 Ma.

Discussion
The analysis of the pairwise distances of the COI sequences, following our DNA-barcoding approach, suggested two alternative putative barcoding gaps, both consistent with the results of our phylogenetic and dating analyses. However, the two gaps entail different scenarios for the cladogenesis in Stenotaenia, including a different estimate of the minimum number of species inhabiting the Italian region.
Assuming a low barcoding gap (1-6% of K2P distances) would be in agreement with a threshold around 1% between intra-and inter-specific distances, as more frequently found and applied in a large assortment of animal taxa, including vertebrates and some groups of arthropods (Ratnasingham andHebert 2007, Chevasco et al. 2014). However, in this case we should conclude that (i) almost every single specimen in our small sample would actually represent a distinct species, (ii) some of these species separated as recently as a few hundred thousand years ago, and (iii) at least nine species of Stenotaenia actually inhabit the Italian region. In this scenario, a remarkable number of species would be expected for the entire geographic range of the genus, in face of the morphological uniformity between species and departing dramatically from the taxonomic scheme in use (Bonato and Minelli 2008).
As an alternative, assuming a higher barcoding gap (10-16% of K2P distances) would be in agreement with gaps estimated in other centipedes, in the very few studies so far published: a gap in p-distances between 9% and 14% for the scolopendromorph genus Scolopendra (Oeyen et al. 2014), a gap in K2P distances between 7% and 10% for the lithobiomorph genus Eupolybothrus (Stoev et al. 2010(Stoev et al. , 2013 and average intraspecific vs. interspecific K2P distances of 7% vs. 18%, for a centipede sample from the Bavarian fauna (Spelda et al. 2011). This is also in agreement with the divergence time analysis, as under the higher barcoding gap scenario all interspecific nodes would be older than one Ma, which is a more reasonable figure (Coyne and Orr 2004). In this scenario, we should acknowledge a lower number of species in our sample. On the basis of this more parsimonious assumption, the minimum number of six species we obtain for the Italian region is nonetheless higher than the four species hitherto recognized on the basis of morphology, a number which besides is derived by a far more extensive geographic sampling (Bonato and Minelli 2008).
Despite evidence being very preliminary, species delimitation suggested by our study may be tentatively compared and matched with these four morphospecies.
The putative species α might correspond to Stenotaenia romana (Silvestri, 1895), even if it includes also a specimen initially identified as S. sorrentina. S. romana was customarily distinguished from all other species in the Italian region by its remarkably minute body size (total length not surpassing 17 mm) and a distinctly lower range of variation in the number of trunk segments (43-49 pairs of legs, Bonato and Minelli 2008). S. romana inhabits mainly the Tyrrhenian side of the Italian peninsula, from Liguria to Campania, including some minor islands (e.g., Ischia and Elba), and also Sardinia; apparently disjunct populations are present in the Euganei and Berici hills (unpublished data).
Species β confidently corresponds to S. sorrentina (Attems, 1903), which is usually circumscribed by having intermediate body size (at least 20 mm when fully grown) and an intermediate number of trunk segments (usually 53-67 leg pairs, Bonato and Minelli 2008). This morphospecies has been frequently recorded in the Italian peninsula from Maritime and Ligurian Alps to Gargano and Calabria, and more rarely and doubtfully also from the Prealps. It is also reported from several islands, including Sardinia and minor Tyrrhenian islands (e.g., Elba, Ponza, Capri and Ischia). Our sample from Sicily (Table 1) confirms the presence of this species in this island, from where only old and dubious records were available so far (Bonato and Minelli 2008).
The species ζ may correspond to S. palpiger (Attems, 1903), which was known so far only for the holotype (total length 15 mm and 49 leg pairs), collected from Istria, about 30 km from the locality of our sample. The species name palpiger has been most often ignored and only recently resurrected as a potentially valid species (Bonato and Minelli 2008). However, the identity and distribution of S. palpiger remain to be clarified, especially with respect to other nominal species of Stenotaenia reported from the Dinaric region, including S. antecribellata (Verhoeff, 1898) and S. cribelliger (Verhoeff, 1898) (Bonato and Minelli 2008).
The remaining three putative species (γ from Western Alps, δ and ε from Eastern Prealps) cannot be assigned confidently to known species. Different number of legs in their specimens had initially suggested different tentative identification (see Table  1). However, according to their origins, it is impossible to tell which of them, if any, actually belongs to S. linearis, which is the fourth species traditionally reported for the Italian region. The type locality of S. linearis is in Baviera (near Regensburg) and based on morphology the species has been reported also from the Alps, as far south as the Maritime Alps (Bonato and Minelli 2008).
While the specimen from Cyprus is confidently recognizable as a different species, S. naxia (Verhoff, 1901) (Simaiakis et al. 2013), which is also recorded from the Aegean area (Bonato and Minelli 2008), nothing can be said for the specimen from Iran. The number of leg-bearing segments and the morphology of the ano-genital region are unknown because the posterior part of the trunk is missing. This makes it very difficult to associate this specimen to any known species, also because no previous records for the genus are known for this geographic area.
Our preliminary evaluation of the phylogenetic relationships and the divergence dates between these species have implications on the evolutionary history of the group which invite some cautious comments.
All specimens in our sample that are representative of species with smaller body size and lower number of legs (S. romana and S. palpiger) apparently represent a single derived clade within the genus, suggesting that these features are synapomorphic for these two species. Considering the entire clade, as far as known, the body length is less than 20 mm long at full growth and the number of leg pairs ranges from 43 (in S. romana; Bonato and Minelli 2008) to 57 (in the S. romana specimen from Isola Fossara).
The genus Tuoba Chamberlin, 1920, which is currently recognized as distinct from Stenotaenia in different respects (morphological features, ecological niche, pattern of dis-tribution; Bonato 2011) turned out to be strictly related to the latter, suggesting also the possibly that Tuoba is a derived clade within Stenotaenia. The phylogenetic position of T. sydneyensis among the Stenotaenia representatives in our sample is unstable and weakly supported, and could also be an artefact of our analysis, possibly produced by the attraction among sequences sharing similar GC-content and GC-skew values (Table 2).
This study is a very preliminary attempt to address the diversity of Stenotaenia with genetic data. It was limited by a relatively low efficacy in DNA extraction, amplification and sequencing from the available specimens (actually up to 15 specimens from 14 localities had been originally processed). This prevented us from analyzing a larger sample of specimens, and did not permit taking into account variation within populations. Nevertheless, our results could serve as a basis for forthcoming investigations on the genus Stenotaenia, not only in the Italian region, but also in the entire distribution range of the genus.