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 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.

Figure 1. 

Sampling localities of Stenotaenia in the Italian region. Greek letters refer to the species tentatively recognized after the analyses (see text).

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).

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 and Giribet 2007, Edgecombe and Giribet 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.

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 and Hebert 2007).

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).

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.

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 p-distances) 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.

Figure 2. 

Frequency distribution of COI pairwise distances. A K2P distances. B p-distances.

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.

Figure 3. 

Maximum likelihood phylogeny. ML tree obtained from concatenated COI and 28S sequences, by the GTR+I+G model, and manually rooted. The following support values are indicated at the nodes (only for those present in the topology obtained from the concatenated sequences): ML bootstrap for the analysis of concatenated genes (upper left); Bayesian posterior probabilities (upper right, in italics); ML bootstrap for the analysis of COI sequences (lower left); ML bootstrap for the analysis of 28S sequences (lower right). 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.

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.

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.

Figure 4. 

Dated phylogeny. 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.

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 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.

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.