Specimen collection and identification

Centipede specimens were collected from both natural and urban habitats through the course of biodiversity surveys in Thailand and adjacent countries under grants to SP since 2010. Permission to enter remote and natural reserve areas was granted by The Plant Genetic Conservation Project under the Initiative of Her Royal Highness Princess Maha Chakri Sirindhorn (grant proposals 2010–2015). During field surveys, some individuals were photographed with either a Nikon D700 or D90 camera equipped with a Nikon AF-S DX Micro-Nikkor 40 mm f/2.8G lens and two Nikon SB 600 Speedlight flash sources to record their living colouration pattern. Collected samples were relaxed with 50% ethanol concentration mixture for 5–20 minutes and then transferred to a higher ethanol concentration for setting their posture for photography. Species identification was made based on previous literature on scolopendrid taxonomy, i.e. studies by Attems [29, 36], Chao [37], Lewis [22, 38, 39] Shelley et al. [24, 35] Schileyko [40–42], Koch [34], Kronmüller [30, 43] and Verhoeff [44]. Terminology applied to taxonomic characters used the standardized nomenclature for centipede morphology by Bonato et al. (2010) [45]. All taxonomic characters were observed under an Olympus stereo microscope connected with a montage imaging system operated under the Cell’D program. All examined materials are housed at Chulalongkorn University, Museum of Zoology, Bangkok, Thailand. Molecular samples are kept in absolute ethanol at -40 degrees Celsius to prevent DNA degradation.

DNA extraction, amplification and sequencing

Sixty Scolopendra samples were dissected to separate tissue from locomotory legs for DNA extraction. The commercial Geneaid DNA extraction and Nucleospin purification kits were used for harvesting genomic DNA. Genomic DNA yields were determined using DNA quantification loading dye blue, loaded in 1X TBE-Agarose gel and run under a 135 V electrical gradient impulse for 15 minutes. Standardized conditions of PCR amplification were edited based on previous molecular works on centipede phylogeny as follow; Edgecombe and Giribet, 2008 [13], Joshi and Karanth, 2011; 2012 [46, 47], Murienne et al. 2010; 2011 [48, 49], Vahtera et al. 2012, 2013 [50, 51] and Siriwut et al. 2015 [52]. Three partial gene fragments were used in this study; the barcode region of cytochrome c oxidase subunit I (COI), 16S ribosomal DNA (16S rRNA), and 28S nuclear ribosomal DNA (28S rRNA).

The PCR mixture consisted of the following: 0.6–1 μl of DNA template, 2.5 μl of forward and reverse primers, 25 μl of Ultra-Pure Taq PCR Master Mix combined with ruby and emerald loading dye and 18 μl of ddH₂O. All gene fragments were amplified using the selected primers for each region as follow: COI was assembled using forward universal primers for COI for animal phylogeny LCO1490 [53] and the modified reverse primer for centipede phylogeny HCOoutout [49, 54], 16S rRNA used primers 16Sa and 16Sb [51, 55], and 28S rRNA was amplified by 28SF4 and 28SR5 [56].

All PCR mixtures were activated by an Eppendorf Mastercycler Pro S machine with manual and gradient functions. The COI, 16S and 28S gene amplifications were performed under standard conditions of PCR reactions cycled at 94°C for 5 min of an initial step, followed by 35 cycles of 94°C for 30 s in a denaturation step. The selected temperatures of the annealing step were 42.5–44.1°C for COI, 45–48°C for 16S, and 55–58.1°C for 28S, for 30 s, 72°C for 15 s in an extension step, and then a final extension step at 72°C for 10 min. PCR cycler was installed at a holding temperature at 4°C as the final step.

The PCR products were inspected under 1% (w/v) agarose gel electrophoresis in 0.5x TBE buffer. The fluorescence of PCR bands was enhanced with SYBR Safe illuminant and observed under a UV light source. The gene target products were purified using a QIAquick purification kit (QIAGEN Inc.). The purified PCR products were directly cycle-sequenced using the original amplification primers with an Applied Biosystems automatic sequencer (ABI 3730XL) at Macrogen and Bioneer Inc. (Korea). Sequences were aligned with libraries in GenBank using the BLASTn algorithm to verify the correct group of organisms from product sequences.

Phylogenetic reconstruction

DNA sequences were assembled in Sequence Navigator [57]. Double strand sequence comparisons were made by a shadow pair-wise alignment function analysis to detect missing sites and gaps in nucleotide sequences, and correlated with chromatograms for each sequence sample. Seven Scolopendra species were used in this analysis: S. dawydoffi Kronmüller, 2012 [30], S. dehaani Brandt, 1840 [58], S. japonica Koch, 1878 [59], S. morsitans Linnaeus, 1758 [60], S. pinguis Pocock, 1891 [61], S. subspinipes Leach, 1816 [23] and Scolopendra sp. Moreover, sequences of two scolopendromorph taxa from GenBank, Cormocephalus monteithi Koch, 1983 [62] and Cryptops doriae Pocock, 1891 [61], were chosen as outgroups to root the trees. DNA alignment was carried in MEGA 6 [63] using MUSCLE [64] with the default parameter set. File format preparation, i.e. FASTA, PHYLIP and NEXUS, for further phylogenetic analysis was implemented in MEGA 6 and Mesquite 3.03v [65]. The heterogeneity of nucleotide substitution model fit was calculated by JModelTest v.1.7 [66] based on the PhyML likelihood algorithm of heuristic search [67] and MEGA 6. Eleven nucleotide substitution schemes [68] and 88 candidate models were set at the beginning of an analysis for each gene fragment in JModelTest. Gaps and missing data were discarded. For the best fit DNA substitution model, the three fragments were analyzed independently. Kakusan 4 [69], implemented in MOLPHY, was used to assemble the final concatenated file.

In this study, Maximum likelihood (ML) and Bayesian inference (BI) were applied to construct phylogenetic trees from the combination of the partitioned DNA dataset. For ML analysis, the concatenated files were analyzed with Treefinder [70] and RAxML 8.0.0v [71]. A single search was conducted to find the starting tree in Treefinder. Fast likelihood-based analyses were performed with 1,000 bootstrap pseudo-replicates. Bayesian inference was conducted in MrBayes, ver. 3.2.5. [72, 73] using four metropolis-coupled, Markov chain Monte Carlo runs [74]. The program was ordered for random sampling of starting trees before more exhaustive analyses. The number of pseudo-replicates was set at 10,000,000 generations, with simultaneous tree sampling at every 500 random replicates. Sixteen nucleotide substitution schemes and the invisible gamma parameter were applied. Seventy percent of harvested trees were removed as burn-in. The analyses were terminated after the standard deviation of proportional frequency reached below 0.01. The consensus tree implemented from 50% majority rules was obtained at the final stage, and draft tree topology files were reconstructed by FigTree [75]. Node support values have been depicted on the trees in those instances where bootstrap values exceed 70% (ML) and posterior probabilities exceed 0.95 (BI) [74, 76]. A Kimura 2-Parameter model [77] was used to calculate corrected distance of all gene fragments in MEGA 6. Genetic distance was compared for both interspecific and intraspecific variation within and between populations. Finally, species justification and validation of the selected Scolopendra species based on genetic affinities are discussed in relation to morphological identification in previous literature.

Geometric-morphometric analysis

One-hundred seventeen Scolopendra specimens were used in this analysis (see S1 Table). Because of morphological changes through the course of ontogeny, a minimum size for sampled individuals was set at 40 mm, following a suggested standard for scolopendrid taxonomy [78, 79]. Three morphological features were examined in this analysis: the cephalic plate, the forcipular coxosternite, and the tergite of the ultimate leg-bearing segment (tergite 21). All samples were photographed in the same orientations and magnifications under a light stereo-microscope. Each feature was analyzed independently by using landmark geometric methods [80, 81]. The landmark points were digitized from a set of stable, conserved parts of each feature, the position detail of each landmark point being as follows:

1. Cephalic plate
   1. Landmark 1:. anterior end of median sulcus of cephalic plate
   2. Landmark 2:. interior basal part of first antennal article (right side)
   3. Landmark 3:. anterior end of anterior ocellus (right side)
   4. Landmark 4:. posterior end of posterior ocellus (right side)
   5. Landmark 5:. intersection point between pleurite of forcipular segment and trochanteroprefemur (right side)
   6. Landmark 6:. intersection point between pleurite of forcipular segment and Tergite 1 (right side)
   7. Landmark 7:. intersection point between pleurite of forcipular segment and Tergite 1 (left side)
   8. Landmark 8:. intersection point between pleurite of forcipular segment and trochanteroprefemur (left side)
   9. Landmark 9:. posterior end of posterior ocellus (left side)
   10. Landmark 10:. anterior end of anterior ocellus (left side)
   11. Landmark 11:. interior basal part of first antennal article (left side)
2. Forcipular coxosternite
   1. Landmark 1:. median diastema
   2. Landmark 2:. inner end of oblique suture (left tooth-plate)
   3. Landmark 3:. outer end of oblique suture (left tooth-plate)
   4. Landmark 4:. left upper corner of forcipular coxosternite
   5. Landmark 5:. coxosternal condyle (left side)
   6. Landmark 6:. inner end of coxosternite collar (left side)
   7. Landmark 7:. left junction between presternite and sternite of first leg-bearing segment
   8. Landmark 8:. right junction between presternite and sternite of first leg-bearing segment
   9. Landmark 9:. inner end of coxosternite collar (right side)
   10. Landmark 10:. coxosternal condyle (right side)
   11. Landmark 11:. right upper corner of coxosternite
   12. Landmark 12:. outer end of oblique suture (right tooth-plate)
   13. Landmark 13:. inner end of oblique suture (right tooth-plate)
3. Tergite of ultimate leg-bearing segment
   1. Landmark 1:. anterior interior margin (right side)
   2. Landmark 2:. anterior exterior margin (right side)
   3. Landmark 3:. posterior exterior margin (right side)
   4. Landmark 4:. posterior exterior margin (right side)
   5. Landmark 5:. distal point of postero-median margin
   6. Landmark 6:. posterior exterior margin (left side)
   7. Landmark 7:. posterior exterior margin (left side)
   8. Landmark 8:. lower exterior margin (left side)
   9. Landmark 9:. upper interior margin (left side)

All landmark points were marked manually by WS with tpSDig2 [82]. The standard image of each constant character of all samples was randomly chosen with tpsUtil [83] in order to avoid personal bias. MorphoJ 1.06b [84] was used for testing shape variation. Procrustes superimposition was calculated to minimize effects such as sample size, orientation and depth [85, 86]. The covariance metric was generated as two-dimensional axes for each feature. Multivariate regression was performed using the Procrustes superimposed data to define allometry and statistically test for correlation between centroid origin and shape variation [87]. A category of sampled specimens was classified based on morphological identification. Canonical variates analysis (CVA) implemented the relative determination of two or more classified groups under Mahalanobis and Procrustes distance values [88, 89]. A permutation test for pairwise distance was set at 10,000 permutation rounds for calculation of Mahalanobis distance in both between- and among-classified groups. From the CVA results, the shape variation detected from landmark positions was linked serially as Wire-frame outlines to visualize shape reformation between negative and positive canonical variates groups on a three-dimensional axis. Confidence ellipses were calculated to indicate the centroid origin of each defined sample group in a three-dimensional CVA graph. Comparative CVA plots were generated separately in each dimension for CV1-CV2 and CV2-CV3. All graphs were exported in Encapsulated Postscript Vector Graphics format (EPS) for processing in Adobe Illustrator.

Morphological identification

This study is based on 176 centipede specimens collected from 134 localities in mainland Southeast Asia and one voucher specimen from the Japanese archipelago. All specimens were observed by light microscopy, with taxonomy based on traditional external morphological characters of scolopendrids. These are as follow: the number of antennal articles, as well as the number of those that are sparsely hirsute (“glabrous”); number of teeth on the forcipular coxosternal tooth-plates; the first tergite to possess complete paramedian sutures; the first tergite with complete margination; the extent of paramedian sutures on the sternites (complete or confined anteriorly to a variable extent); the number of spines on the coxopleuron (specifically, the number of apical spines and the presence/numbers of subapical and dorsal spines); the prefemoral spine arrangement on the ultimate legs; the presence or absence of tarsal spurs on legs 19 and 20; and, presence or absence of a gonopod (“genital appendage”) on the first genital segment of the male. The taxonomic results show six nominal species that can be identified as named species and one putative new species in the sampling area as follow: Scolopendra dawydoffi, S. dehaani, S. japonica, S. morsitans, S. subspinipes and Scolopendra sp. Diagnostic character combinations of all assigned species are summarized in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Diagnostic description of all examined species based on external morphology and common colouration schemes of voucher specimens in this analysis, with references to recent taxonomic descriptions with additional information.

https://doi.org/10.1371/journal.pone.0135355.t001

Sequence annotation

Sixty nucleotide sequences from partial gene targets for cytochrome c oxidase subunit 1, 16S rDNA and 28S rDNA were obtained (Table 2). All raw nucleotide sequences were blasted with other available scolopendromorph sequences in GenBank as a check for contamination. The compatibility values of all sequence reached up to 80% of available scolopendromorph sequences, suggesting that outgroup contamination is not affecting the genomic DNA. The final aligned sequences obtained from sequence editing and the alignment program, consisted of 814 bp for COI, 446 bp for 16S, and 638 bp for 28S. Sequence annotation (Table 3) of each gene is as follows: COI sequences consist of 334 parsimony-informative sites, and 403 and 411 variable and conservative sites, respectively; 16S sequences comprise 197 parsimony-informative sites, and 271 and 175 variable and conservative sites, respectively; 28S sequences include 110 parsimony-informative sites, and 206 and 432 variable and conservative sites, respectively. Corrected genetic distances were calculated under the Kimura-2-parameter model for DNA sequence alignment. Interspecific variation in each partial sampling of genes is 15–24.2% for COI, 10.6–22.4% for 16S, and 0.8–10.8% for 28S. A summary of inter-intra specific variation and best fit scores for the nucleotide substitution model is given in Tables 4 and 5.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. List of voucher specimens of seven Scolopendra species and selected outgroups used in phylogenetic analyses.

Each sample includes the collecting locality, GPS co-ordinates, CUMZ registration numbers, and GenBank accession number for three selected genes (COI, 16S and 28S).

https://doi.org/10.1371/journal.pone.0135355.t002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Characteristics of nucleotide sequence for three amplified genes and best fit models of heterogeneous nucleotide substitution for each gene calculated from jModel Test under AIC and BIC criteria.

https://doi.org/10.1371/journal.pone.0135355.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Corrected distance of interspecific variation in seven Scolopendra species from COI and 16S partial gene analyses under calculation model of K-2 parameter.

https://doi.org/10.1371/journal.pone.0135355.t004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Corrected distance of intraspecific variation in six Scolopendra species from COI, 16S and 28S partial gene analysis under calculation model of K-2 parameter.

https://doi.org/10.1371/journal.pone.0135355.t005

Phylogenetic analysis

The best fit models of heterogeneous nucleotide substitution under the two optimality criteria (Maximum likelihood and Bayesian inference) for each partial gene analysis are JC, GTR+G+I and GTR+G (Table 3). The output trees for the concatenated analyses have congruent topologies in both analyses (Fig 1; see S1 Fig for node support). They depict the expected monophyly of the subfamily Scolopendrinae Kraepelin, 1903 [90], with all species of Scolopendra and the outgroup OTU of Cormocephalus monteithi nesting together with strong posterior probability support in BI (Fig 1A). Inside the Scolopendrinae, however, the monophyly of the genus Scolopendra is contradicted by the interpolation of C. monteithi within it (Fig 1D). In the case of examined Scolopendra taxa, seven genetically delimited taxa can be discriminated that are congruent with their morphological identification as species, and phylogeographic structure is resolved within each of these species. The tree separates members of Scolopendra into three main clades. One clade includes S. dehaani, a putatively new species that we refer to as Scolopendra sp., and S. subspinipes (labeled as clade B in Fig 1), and another unites S. morsitans and S. pinguis (clade C in Fig 1). The third clade (clade E in Fig 1) groups S. japonica and S. dawydoffi. From the tree topology, short internal branch lengths are characteristic of populations within S. dehaani (Fig 1, clade F therein) while the remaining species showed greater amounts of genetic diversity between their populations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Phylogenetic tree of Scolopendra mainland Southeast Asia.

Relationships among Scolopendra and two outgroups indicated similarly both in Maximum likelihood (ML) and Bayesian inference (BI) of the concatenated COI, 16S and COI partial gene analyses. Significant support values in ML and BI are indicated by three colouration circles; black circle = support both in ML and BI (above 70% bootstrap in ML and 0.95 posterior probability in BI), white circle = support only in ML, grey circle = support only in BI. The gradient colouration bars on the tree represent the genetic affinities of populations relative to morphological identification in each species.

https://doi.org/10.1371/journal.pone.0135355.g001

In S. dehaani (Fig 2), the phylogenetic tree indicates three major groups that have clear relationships to geographical zones as follow: the Chao Phraya Basin population (CPB; Clade A), the Mekong River Basin population (MRB; Clade C), and the Lower Tenasserim Range population (LTR; Clade B). The CPB clade unites populations from the northern, western (upper part of the Tenasserim Range) and central parts of Thailand. The northern population is separated from the others with strong node support in ML and BI analysis (73/1) while the relationship between MRB and LTR remains unresolved, with weak support values in ML and BI (64/0.86). The MRB clade consists of all populations from northeastern and eastern parts of Thailand, Laos and some western Cambodian populations. Inside the MRB clade, node support values indicate that the western Cambodian population and the lower northeastern Thailand population are closely related (98/1; ML and BI) while relationships in the other populations are undefined. The LTR clade includes southern populations starting from the Isthmus of Kra through the coastal part of Thailand to the northern part of the Malay Peninsula. In this clade, the tree topology indicates two main populations, one northern and other southern. The split between these two clades received strong support both in ML and BI (98/1). In addition, genetic relationships are congruent with colouration among regional populations of S. dehaani in the Southeast Asian mainland. During field surveys, five colouration patterns have been recorded (Fig 2A–2E). Colour morphs 2–4 are found throughout most regions whereas morphs 1 and 5 (Fig 2A and 2E) are specific to only some populations. Colour morph 1 (Fig 2A) has been reported only from some specimens of the CPB population. This colour morph has a dark brown body and dark violet legs, with most specimens coming from the north western part of Thailand, close to the Thailand-Myanmar border. Also geographically restricted is colour morph 5 (Fig 2E), which has been found only in the LTR population. In this, the body is bright reddish with a dark band on the posterior part of the tergites while all legs are usually reddish or yellowish.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2.

Phylogenetic relationship of S. dehaani population (left) based on genetic structure among its populations relative to regional distribution in mainland Southeast Asia (middle); colours indicate the major populations. Five patterns of live colour morphs in S. dehaani were found (right); A. Dark colour morph; B. Light brownish colour morph; C. Reddish-brown body color morph, yellowish legs with reddish on distal part; D. Dichromatic pattern; cephalic plate, tergite 1, 20 and 21 reddish, tergites 2–19 brownish with yellowish legs; E. Reddish colour morph with dark band on anterior and posterior parts of tergites.

https://doi.org/10.1371/journal.pone.0135355.g002

Scolopendra sp. and S. subspinipes are the closest relatives of S. dehaani (Fig 1; clade B), although the precise sister group relationships are equivocal because the putative clade that unites Scolopendra sp. and S. dehaani to the exclusion of S. subspinipes has low support values in both ML and BI (39/0.63). The union of the two sampled populations of Scolopendra sp. from northern and southern Laos supports the validity of this species, as is indicated by diagnostic morphological characters, notably the incomplete paramedian sutures on the tergites. In case of S. subspinipes, the tree topology indicated the separation of S. subspinipes from S. dehaani as well as Scolopendra sp. according to both support values (56/0.97) and corrected distances (12.9–13.7% and 17–19.9% in COI and 16S, respectively).

S. morsitans was nested in the same clade with S. pinguis and an outgroup (C. monteithi) (Fig 1; clade C). All S. morsitans populations were nested together as monophyletic. However, the genetic affinity divided S. morsitans into two minor populations (Fig 3), one of which is found in the northern part of Thailand (clade A in Fig 3) whereas the other is located in northeastern Thailand and some populations from Cambodia. This separation was supported with strong nodal support in both ML and BI (100/1) whether between or within the two populations. Furthermore, observation of colour patterns indicated two colour morphs for S. morsitans in this region (Fig 3). The most common pattern, called colour morph 1, covers most populations in this region. The cephalic plate, tergite 21 and ultimate legs of this morph are red-yellowish, while colour morph 2 has the cephalic plate, tergite 21 and ultimate leg blackish. This latter morph is only found in the north-western population.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Phylogenetic tree of S. morsitans.

Colour gradient indicates population structure; blue gradient indicates the northern population, yellow gradient the eastern population. Scolopendra morsitans exhibited two colour morphs: colour morph 1—antenna, cephalic plate, tergites 1, 20 and 21 and ultimate legs orange; colour morph 2—antenna, cephalic plate, tergites 1, 20 and 21 and ultimate legs blackish.

https://doi.org/10.1371/journal.pone.0135355.g003

The tree shows marked genetic diversity between populations of S. pinguis, as indicated both by branch lengths and the corrected within-species distances (18% and 11% in COI and 16S, respectively). Monophyly of the species is supported (Fig 4), as is consistent with its diagnostic morphological characters, and groupings within the species received strong support in both ML and BI analysis. The genetic structure of S. pinguis showed that Thai populations exhibited a genetic connection with all examined populations from Laos. Moreover, the tree topology depicts a relationship between western and northwestern Thailand, which showed closer affinities with each other than with other populations. However, one population (Fig 4; P1) from the western part of Thailand did not group with the adjacent population (Fig 4; P4) and the former is resolved basally relative to all other S. pinguis populations. This species also exhibited colour variation among its populations, as for the previously discussed species. Four colouration patterns (Fig 4) have been recorded both in the Thai and Laotian faunas. The major colouration pattern is seen on the cephalic plate, which differs between yellowish and blackish colour morphs. Minor variability has also been detected on the legs, these exhibiting two colour forms, either being monochrome or dichromatic. The tree topology indicated that all populations with the yellowish cephalic plate were grouped together but one of the two blackish populations is nested within the yellowish population. The blackish populations are divided into two lineages, one of which is sister lineage to a clade that unites all other sampled individuals. With regards to their geographical distribution, the blackish cephalic plates are specific to the north western area of Thailand whereas the other forms are distributed both in Thailand and Laos.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Phylogenetic tree of S. pinguis based on genetic structure of its populations.

Colour gradient bar indicates colour morphs of sampled individuals that divide into two patterns; colour morp 1—blackish population (monochromatic); colour morph 2- yellowish—black population (dichromatic). Four live colour morph pictures depict the variability of colouration on legs of the two colour morphs in S. pinguis; colour morph 1A and 2A—animal with dark blue legs, colour morph 2A and 2B—animal with yellowish legs with blue stripes on distal part.

https://doi.org/10.1371/journal.pone.0135355.g004

Scolopendra dawydoffi and S. japonica unite as a clade in which both species are reciprocally monophyletic (Fig 5). The taxonomic validity of these two species has been corroborated by strong node support (96/0.99). The corrected distances indicated low genetic diversity among S. dawydoffi populations as depicted by short branch lengths in the phylogenetic tree (Fig 5; clade A). In case of S. japonica (Fig 5; clade B), individuals from the sampled populations unite as well-supported clades, the corrected distance within the species being 12.4% and 5.3% in COI and 16S, respectively. A Japanese specimen was resolved as sister taxon to a group composed of the S. japonica populations from Laos. With regards to colouration, there is no evidence in our collections for S. dawydoffi exhibiting colour variation. In contrast, S. japonica has two morphs that can be distinguished by colour of the legs: in colour morph 1, the legs are reddish, whereas colour morph 2 has green-yellowish legs. The Laos fauna includes both colour morphs whereas the specimen from Japan shows similarity to the Laotian population classified as colour morph 1. However, the phylogenetic tree indicated that the colouration does not precisely correlate with genetic affinities in this species.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Phylogenetic tree of S. dawydoffi and S. japonica.

Clade A, S. dawydoffi, clade B, S. japonica. In S. japonica, colour gradients indicate the colour morph of sampled individuals; colour morph 1—greenish body with reddish antenna, cephalic plate and legs; colour morph 2—greenish body with blue antenna, yellowish cephalic plate and legs.

https://doi.org/10.1371/journal.pone.0135355.g005

Geometric morphometrics

Six species of Scolopendra in the Southeast Asian mainland (excluding S. subspinipes) have been identified based on morphological taxonomy. Specimen numbers used in this analysis are as follow: five specimens of S. dawydoffi, 84 specimens of S. dehaani, eight specimens of S. japonica, 10 specimens of S. morsitans, 12 specimens of S. pinguis, and two specimens of Scolopendra sp. The Procrustes ANOVA statistical test found no measurement error in all analyses (p<0.0001; S2 Table). CVA plots were performed based on group classified datasets with 10,000 replicates of permutation testing for pairwise distance. Eigenvalue and variable percentages of selected characters are summarized in S3 Table. The three dimensions of CVA plots have been used for species discrimination from shape variation (Fig 6). The discriminant results of Scolopendra species from CVA analysis with statistical testing in three selected features are described as follows.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Diagram of landmark locations on three constant characters and the CV plot of individual scores on each CV axis from canonical variates analysis (CVA).

A. Cephalic plate; B. Forcipular coxosternite; C. Tergite 21. The CV plots represent the discrimination of classified individuals scored from CV axis comparison, showing comparisons of CV1 and CV2 axes (middle column) and CV2 and CV3 axes (right column)

https://doi.org/10.1371/journal.pone.0135355.g006

Cephalic plate (Fig 6A): The CV1 axis captured 46.468% of shape variation while the CV2 and 3 axes exhibited 26.269% and 15.575%, respectively. The CVA plot clearly indicated the differentiation of centroid origin in the CV1-CV2 axis. From this CVA plot, all individuals of the same species were grouped together. The clusters of individuals of both S. morsitans and S. japonica separated from other species, whereas S. dawydoffi, S. dehaani and Scolopendra sp. were grouped closely together. In contrast, a poorer discriminant resolution was found in the CV2-CV3 axes, with all species apart from Scolopendra sp. overlapping with the range of variation in S. dehaani. From the statistical result (S4 Table), the p-values from permutation tests of Mahalanobis distances among groups supported the discrimination of four Scolopendra species, these being S. dehaani, S. japonica, S. morsitans and S. pinguis (p<0.0001). However, no significant difference among groups was indicated by p-values of Procrustes distances.

Forcipular coxosternite (Fig 6B): Forcipular coxosternal shape showed appreciable variability in most Scolopendra species except in Scolopendra sp., although it need be noted that the sample size of the latter species is small. The percentages of variability contributed by the CV1, CV2 and CV3 axes are 54.298%, 30.728% and 7.995%, respectively. The CV1-CV2 pairwise comparison plot showed that individuals of each species show a marked clustering. Variation in the shape range of the forcipular coxosternite in S. japonica and S. dawydoffi overlapped and extended into the range of S. dehaani, whereas S. morsitans and S. pinguis are well delineated. In the CV2-CV3 plot, all individuals of S. pinguis were pooled separately with each other while the remaining species were grouped closely. All individuals of S. morsitans and S. dawydoffi were interpolated inside the range of variation of S. dehaani. From permutation tests (S5 Table) of Mahalanobis distances among groups, p-values indicated the distinctness of four Scolopendra species, i.e., S. dehaani, S. japonica, S. morsitans and S. pinguis, with significant support (p<0.0001). The Procrustes distances significantly differed only in three examined species, S. dehaani, S. morsitans and S. pinguis.

Tergite 21 (Fig 6C): The shape variation percentages of the three CV axes are as follow: CV1 explained 64.358%, CV2 28.756%, and CV3 5.558%. The CVA analysis exhibited the discriminant centroid origin of CV1-CV2 comparison plot in all examined species. Scolopendra dehaani, S. pinguis and Scolopendra sp. depicted wholly or mostly unique ranges of variation. The three remaining species exhibited shape variability against each other. However, the CV2-CV3 comparison plot allowed S. pinguis and Scolopendra sp. to be distinguished from other species. In this comparison, most individuals of S. dawydoffi, S. japonica and S. morsitans were gathered inside the variation of S. dehaani. The Mahalanobis and Procrustes distances among groups (S6 Table) indicated four distinct species, i.e., S. dawydoffi, S. dehaani, S. japonica, S. morsitans and S. pinguis, with statistical support (p<0.0001).

Based on the CVA analyses for each assigned feature, the Eigenvalue and percentage of variance were generated in S3 Table. The variation of the landmarks in the three selected morphological features has been recorded and described as follows:

Cephalic plate: In the CV1 axis, variable reformation was found on Landmarks 1, 2, 3, 4, 5, 6, 7 and 8 (see column A in Fig 7). The position of Landmark 1 shifts posteriorly as CV scores trend positively. This position describes a deeper median sulcus on the anterior part of the cephalic plate. Landmarks 2, 3 and 4 (as well as Landmarks 9, 10 and 11) are displaced medially as CV scores trend from negative to positive. This variation relates to the shape of the cephalic plate, the size of basal antennal article, and relative length of the basal article and the ocelli. The movement of Landmarks 5 and 6 (also 7 and 8) record an increase in width of the posterolateral part of the cephalic plate from positive to negative. In the CV2 axis (column B in Fig 7), variable sites were found to relate to Landmarks 1, 4, 5, 6, 7, 8, 9 and 10. The most extreme change was recognized both in CV negative and positive groups, this being located between Landmarks 5 and 6 (and 7 and 8), affecting length of the cephalic plate. In the CV3 axis (column C in Fig 7), only six landmarks showed slight variability, these being Landmarks 1, 2, 3, 6, 7 and 8 in both CV positive and negative score datasets. This variation again impacted on the length of the cephalic plate.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Wireframe diagram from continuous linkage of all landmark positions in three features derived from CVA scores on three axes.

CV1, CV2 and CV3 arranged vertically, respectively. In the wireframes, dotted lines represent shape changes relative to CV score moving from both negative and positive directions, solid lines represent the shape consensus in negative and positive groups.

https://doi.org/10.1371/journal.pone.0135355.g007

Forcipular coxosternite: In the CV1 axis, shape variation was detected in all landmark points (column A). Movement in Landmarks 3–12 describes the shape of the coxosternite in both CV negative and positive groups. The forcipular coxosternite in the CV positive group is relatively shorter than in the negative. Landmarks 3 and 4 (also 11 and 12) shift laterally from negative to positive, describing a broadening of the anterolateral corner of the coxosternite. Landmarks 6, 7, 8 and 9 moves anteriorly from negative to positive. In the case of the CV2 axis (column B in Fig 7), the outline between Landmarks 12–13 and 2–3 showed variation. Specifically, landmarks describing the course of the anterior parts of the forcipular coxosternite shift from more curved to straighter towards the positive group. Moreover, the posterior part of the forcipular coxosternite also showed high variability at Landmarks 7 and 8, which shift anteriorly from negative to positive groups. Simultaneously, the coxosternal condyles are situated more posteriorly from negative to positive groups. From this variability, the outline of the forcipular coxosternite in the CV positive group is more trapezoidal and straighter / more transverse across both its anterior and posterior parts whereas the CV negative group exhibited a more curved diamond shape on its anterior and posterior parts. From the CV3 axis (column C in Fig 7), the shape variation of the forcipular coxosternite in both the CV negative and positive groups showed only slight change, this being captured by Landmarks 2, 3, 4, 5, 6, 7, 8 and 10.

Tergite 21: In the CV1 axis, shape variation exhibited substantial changes in both CV negative and positive groups (Column A in Fig 7). Tergite 21 depicts relatively elongate versus broad rectangular shapes in the CV negative and positive groups, respectively. According to Landmark 5 in particular, the posterior margin of tergite 21 shifts posteriorly from positive to negative, and combined with movements of other landmarks corresponds to the change in length: width ratio of the tergite. Landmarks 1, 2, 8 and 9 shift laterally from negative to positive groups, whereas Landmarks 2, 3, 7 and 8 shift so that the tergite transforms from longer than wide to the reverse. In CV 2 axis (column B in Fig 7), a stable landmark was identified in Landmark 5, which showed consistency in both CV extremes. However, the visualized shape outline of tergite 21 also reflected the distinctness between CV populations, similar to the case for the CV1 axis. The CV positive group shows a displacement of Landmarks 3 and 4 (also 6 and 7) so that the posterior margin of the tergite changes from a broad V-shape to having the short extent delimiting the margination being transverse, and then flexing to its V-shaped extent. In the CV3 axis (column C in Fig 7), the shape of tergite 21 showed less variation than did the CV1 and CV2 axes. Landmarks 1, 3, 5, 7, 8 and 9 all exhibited variability on this axis. In the CV negative group, tergite 21 has the lateral margin (between Landmarks 2 and 3 / 7 and 8) more strongly diverging posteriorly than in the CV positive group.

From the results of discriminant function analysis among classification categories (Table 6), five assigned categories of morphometric samples were confirmed by the percentage of correct classification in two features, the forcipular coxosternite and tergite 21 (>80% in all categories). However, the percentage of correct classification based on variation in the cephalic plate was lower (<50%) in S. dawydoffi and Scolopendra sp. The cross-validation of discriminant function analysis showed a low percentage of correct classification in all sampled characters and taxa except for S. dehaani, which received 89% correct classification when defined by the forcipular coxosternite and tergite 21.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Results of CV discriminant function in three selected characters; the total number and percentage of correction of leave-one-out cross validation tests in CV discriminant function are in parentheses.

https://doi.org/10.1371/journal.pone.0135355.t006

The diversity of Scolopendra in mainland Southeast Asia

The field survey in this study identified the occurrence of S. japonica in mainland SE Asia, its character data conforming to a recent taxonomic review [30]. The distribution range of all species has been refined and it is now possible to make inferences on several species usually reported as regional widespread species [22]. In this study, seven Scolopendra species have been found in both natural and anthropogenic areas. Previous records of Scolopendra in SE Asia [42, 90, 91] indicated that there are three additional species that can potentially be found in this area: S. gracillima Attems, 1898 [92], S. calcarata Porat, 1876 [93], and S. hardwickei Newport, 1844 [94]. However, very few specimens were reported in the relevant studies and the species listed above were treated as likely introductions [42]. In mainland SE Asia, S. dehaani has been found to be the dominant, widespread species throughout the sampling territory, whereas S. subspinipes (S. subspinipes subspinipes of most previous studies) is rarely found, as conforms to previous work [91, 95, 96]. In contrast to those species, some species seem to be endemic and have scattered distributions and sparse populations in nature, i.e., S. dawydoffi, S. japonica, Scolopendra sp. and S. pinguis. From phylogenetic analysis, the regional populations of some examined species such as S. dehaani, S. morsitans and S. pinguis suggested the genetic affinities among geographically neighboring members, infraspecific structure that might be affected by the geographical richness of the region. Several areas have been promoted as corridors for dispersal or land bridges [1, 97, 98], and the mechanisms that sculpt the gene pool among regional populations may derive from these boundaries, as has been proposed for other invertebrates and vertebrates [99–101].

Taxonomic validity of some Scolopendra members in SE Asia

Currently, the species diversity of Scolopendra in Southeast Asia comprises 13 species that are distributed in the mainland and insular faunas [21]. Among them, morphological examination is adequate for species delimitation in some species, such as S. morsitans and S. dehaani. However, there are other species that show high morphological variability, indeed more than previously estimated, such as S. dawydoffi, S. japonica and S. pinguis. As noted above, their variability might overlap with other related species in this region such as S. multidens (in the case of S. dawydoffi) and S. gracillima (in the case of S. pinguis) [29]. Scolopendra dehaani was originally established as a full species [56] but because of its morphological similarity with S. subspinipes, subsequent taxonomists treated it as a subspecies of S. subspinipes, until it was only recently revalidated as a separate species [30]. Another two former subspecies of S. subspinipes, S. dawydoffi and S. japonica (previously known as S. subspinipes cingulatoides Attems, 1938, and S. subspinipes japonica Koch, 1878, respectively), have likewise most recently been elevated to the ranks of separate species [30]. Our molecular examination confirmed the validation of these two species, not the least because they are resolved as more closely related to each other than either is to S. subspinipes (Figs 1 and 3), and was supported by new additional characters for species discrimination. For these reasons, we consider that three former subspecies of S. subspinipes sensu Lewis, 2010 are valid species, as likewise determined based on external morphology alone [30].

Phylogeography of Scolopendra in mainland SE Asia

Phylogeography has been introduced for centipede systematic studies in the past decade [19, 20, 27, 109]. Geological events in the past that potentially bear on the distribution range of these animals include the drift and collision of former fragments of Gondwana during the Jurassic [46], micro-refugia in the last glacial period during the Pleistocene [26, 27], and the debated hypothesis of biotic shuttle for some insular centipedes [48]. Phylogenetic results in this study depict genetic relationships of some Scolopendra populations that are congruent with geographical barriers in mainland Southeast Asia. These findings may relate to the two sub-regional faunas of the Indo-Burmese biodiversity hotspot, Indochina and Malesia [1, 111–113]. The genetic structure of S. dehaani populations suggested a separation into three lineages. The CPB population occupies the northern, central and some western elements of Thailand while the MRB population contained the entire eastern element of the Indochina sub-region, including the northeast and east of Thailand, Laos and Cambodia. The separation between these two major populations is located on the western margin of Korat Plateau along the Phetchabun, Dong Paya Yen and Sankambeng mountain ranges. Uplift of the Korat Plateau has been estimated to date to the middle-late Triassic [114, 115]. This plateau is delimited by the mountain block between northeastern Thailand and Cambodia, which is likely to have initiated speciation and divergences between these two faunas, as exemplified by amphibians [16] and reptiles [116]. However, our analysis groups together the entire population of S. dehaani from northern Laos to the eastern coast of Thailand, demonstrating less sensitivity to vicariance than in some other animal groups.

In the case of the CPB population, the northwest and central Thai samples were united as a core group. The one sample from the upper northern region exhibited genetic difference from the rest, being resolved as the basal clade of the CPB population. The genetic distance may reflect the genetic transfer among neighbouring populations along a geographical gradient. The parallel mountain range along the western and northern parts of Thailand may not limit the dispersal of this population only to the northwestern part of Thailand. Collision of the Indian subcontinent with the Laurasian plate during the Eocene (55 Mya) [98, 117, 118] is likely to have contributed to dispersal of S. dehaani, which is known from India [119]. Other elements of the SE Asian and Indian biotas, both extinct [120] and extant [47], show similar patterns of diversification. Moreover, a hypothesis that refers to some Indian centipede taxa dispersing into SE Asia is consistent with molecular dating and genetic composition among populations [52].

Similar to the previous population, the LTR population is divided into two main groups. The first covered the population along the Isthmus of Kra, which is assumed to be a transitional zone for animals in this region [111, 121, 122]. The second lineage is the lower Tenasserim Range population, which extends from the central part of southern Thailand to the Malay Peninsula. These two minor populations were separated by the mountain valley between the Phuket Range in the west coast and the Titiwangsa Range in the lower east coast of the Thai-Malay Peninsula. Sea level fluctuation in the South China Sea and its role in habitat change/loss has been discussed as a factor impacting on population-level patterns in this region [112, 123–127], and also contributed to exposure of the Sunda Shelf along the coastal areas of Indochina and the Malay Archipelago [101, 112, 128–131]. These geological events seem to be driving mechanisms for the migration of both flora [132] and fauna [16, 111, 116, 133, 134] in this area. Our result showed that one sample from the eastern coast of Thailand exhibited genetic affinity with the entire southern population. This result is congruent with data from taxonomic studies on other animal groups that report relationships between populations from the same two regions, including butterflies [135] as well as centipedes [42, 136]. These repeated patterns are suggestive of a cause that is linked to the geographical history of Sundaland. However, passive dispersal from anthropogenic activity is also a possibility because of the habitat preference of this species, S. dehaani often being collected in disturbed sites with human modification.

Contrary to the widespread S. dehaani, the rare species S. pinguis exhibited high genetic diversity among its populations. The updated distribution from this study showed a limit to dispersal along the northwestern and northeastern mountain ranges of mainland SE-Asia. Genetic distances indicated that population structure relates to geography in that western, central and eastern populations are differentiated along the northern region. Genetic variability likewise reported in other morphologically conservative groups of vertebrates and invertebrates has often been found to signal cryptic species, as exemplified by amphibians [137], molluscs [3], diplopods [12], and even in some scolopendrid centipedes [47]. Accordingly, some of the distinct populations of S. pinguis in the west should be monitored to test the hypothesis that cryptic species may be involved.

In case of S. japonica, the phylogenetic result showed two regional populations. This genetic relationship provides evidence for shared animal diversity between the eastern coast of the Palearctic (here sampled by Japan) and Indo-Malay ecozones, as reported in other invertebrates, e.g. dragon flies [138], scorpions [139], and spiders [140, 141]. Populations within S. dawydoffi depict close relationships to each other (i.e., short branch lengths). However, the population from the east coast of Thailand was resolved as sister group to the remaining conspecific populations. The continuous geography of this area may facilitate genetic transfer among these populations, as is likewise the case for populations of S. morsitans from Northeastern and Eastern Thailand and Cambodia. The ecological richness of this area likely contributes to it being a migration route for some organisms throughout the Indochina sub-region [142]. However, the precise distribution range of S. dawydoffi is unknown (records are scattered in Laos and Vietnam) and additional material will be needed to reveal the fine detail of its genetic structure.

Colouration patterns of Scolopendra

The colour variation on animal bodies has been promoted as an interesting topic for evolutionary studies of various kinds of organisms, such as nudibranchs [143], land snails [144], and butterflies [145], as well as in centipedes [26]. Recently, molecular phylogeny has been used to explain the relationship between two different colour morphs distributed sympatrically with each other [26]. In our survey, colour variation in populations was recorded for all examined species, Scolopendra being the only centipede genus in the biota that shows highly diverse patterns of colour morphs [34, 35]. Several species exhibited extreme difference as bright or contrasted colour morphs between populations. The number of patterns differs in various species [95], for example S. laeta Haase, 1887 from Australia exhibited five colour morphs [34], two colour morphs of S. cingulata were found in the Mediterranean [26], and four colour morphs have been documented in S. morsitans amazonica Bücherl, 1946 [146] (= S. morsitans) from Africa [147]. Moreover, the significance of colour variation in several species has been debated in the literature with regards to its taxonomic value. In previous studies, authors usually assigned the different colour morphs of centipedes to new species [148, 149] or treated them as subspecies or varieties of species [59]. Recently, several taxonomic revisions argued against this approach. The number of fresh specimens for examination in several species is limited, and when colours have faded in preserved museum specimens the taxonomic utility of colour is compromised. For these reasons, several nominal species in the past, which had been identified by colouration pattern, have fallen into synonymy with other widespread species [30, 35]. However, our study revealed the relationship of these morphological changes to the genetic structure of some species, in several cases showing that colour morphs correspond to monophyletic groups.

In this study, the species with most diverse colour morphs is S. dehaani. This centipede exhibited five colouration patterns. Colour morphs 2, 3 and 4 are standard colour morphs which are usually found in all populations of S. dehaani whereas colour morphs 1 and 5 are specific to only two regional populations, CPB and TLH. Environmental factors such as humidity and other habitat characteristics were assumed to affect this variability. Comparable variability in colouration has been reported and discussed in systematic studies on tropical organisms, e.g., snakes [150], amphibians [151], and sea cucumbers [152]. These results suggest a correlation between colouration and such factors such as habitat type [153] and predation [154].

The distinct colouration patterns in S. pinguis populations differ from other species in our analysis. In general, juvenile Scolopendra show different colour patterns compared to adult stages, e.g., in S. morsitans and S. dehaani. In the case of S. pinguis, however, colouration patterns are consistent through an observed series of post-embryonic developmental stages, such that paedomorphosis may account for this variation. Ontogenetic data such as these have been suggested as necessary to clarify some taxonomic problems in scolopendromorphs [155, 156].

The remaining examined species, S. morsitans and S. japonica, exhibited specific colour morphs in different populations. Interestingly, the black and reddish patterns found in S. morsitans were also reported in African and Australian populations [34, 157]. The pattern of colouration change is specific to particular morphological characters such as the cephalic plate and tergite of the ultimate leg-bearing segment [34, 95], as often noted in taxonomic descriptions of Scolopendra. According to the colour variation of these samples, we can separate the colouration pattern among SE Asia Scolopendra into either monochrome or dichrome (Table 1), following similar descriptions in some previous Scolopendra studies [37]. The colour polymorphism of Scolopendra populations might be useful as a model for molecular ecology and evolutionary studies in SE Asian centipedes.

Taxonomic implications of shape variation

Shape variation tests in morphology have been widely used in fossilized material because of a frequent shortage of morphological characters for species identification [158–160]. These techniques have also been applied to systematic study of varied organisms for species discrimination [3, 161, 162]. In this study, the three surveyed features show results that appear to be useful as taxonomic signal for species-level taxonomy, as exemplified by plots from canonical variates analysis (Fig 6) and CV discriminant functions (Table 6). Each of the examined characters had been documented in some scolopendrid literature over the past century [29, 61], but this study showed the first indications of taxonomic value of these features based on their discrete shape variation among Scolopendra species. Shape variation in tergite 21 seems to be the most variable of the features explored in our analysis and, indeed, this character has been used as diagnostic for some scolopendrid species, e.g., Cormocephalus cupipes Pocock, 1891 [102] and Otostigmus caudatus Brölemann, 1902 [163]. To reduce confounding effects and improve accuracy, further shape analysis techniques such as combined outline and semi-landmark methods that have been used in scutigeromorph centipedes [31, 32], as well as meristic measurements [33], are likely to be promising topics for further study.