Introduction

Millipedes (Diplopoda) are a large, understudied group. Twelve thousand species belonging to sixteen orders have been described and it is estimated that as many as 80, 000 species exist. Little is known about the ecology, life history, and phylogeny of the great majority of species (Sierwald and Bond 2007). Millipedes are detritivores, feeding on decaying organic matter in the leaf litter and as such play an essential role in nutrient cycling. They are also model organisms for microendemic studies with many species confined to small geographic areas near populations of closely related taxa. Such distributional patterns are probably due to their inability to travel long distances (Sierwald and Bond 2007). Millipedes are the fourth largest group of cave-obligate (troglobiotic) invertebrates in North America, after the insects, arachnids and crustaceans (Culver et al. 2000).

Troglobiotic millipedes inhabit caves in temperate and tropical areas (Mauries 2003; Culver and Pipan 2009). Troglobiotic millipedes belong to a variety of orders including Chordeumatida, Callipodida, Glomerida, Polydesmida, Julida and Spirostrepsida (Mauries 2003; Culver and Pipan 2009). Most have evolved troglomorphic characters associated with cave-adapted organisms including reduced ocelli, elongated appendages, lack of pigmentation, and decalcification of the integument (Mauries 2003; Culver and Pipan 2009). Like other cave organisms, troglobiotic millipedes are generally dependent on resources that have been carried into caves from the surface. Troglobiotic millipedes may feed on leaf litter, scrape bacteria off of rocks or feed on bat guano. Some species are amphibious, and filter particulate matter from the water. Troglobiotic millipedes can be abundant in caves with dozens of individuals present within a small area (Culver and Pipan 2009).

More than 85 troglobiotic millipede species have been described in North America, most of which belong to the family Chordeumatidae. Some genera are exclusively troglobiotic: e.g. Zygonopus Ryder, 1881 from West Virginia, Mexiterpes Causey, 1963from Mexico and Scoterpes Cope, 1872 from southeastern North America, whereas other genera, such as Pseudotremia Cope, 1869, contain surface-dwelling species as well (Shear 2008). Troglobiotic millipede species may either have broad ranges or be limited to a single cave or several caves within a given karst system. Range size may be influenced by regional geological differences. For example, Scoterpes species tend to have small ranges where the limestone strata are highly folded and broad ranges where limestone strata are flat and extend over long distances (Shear 2010).

Genetic research on troglobiotic millipedes is limited. Laing et al. (1976) used isozymes to study two populations of Scoterpes inKentucky and found high genetic diversity between populations but low genetic diversity within each population. These populations were later described as belonging to different species (Shear 2010). Woo et al. (2007) sequenced the entire mitochondrial genome of a Korean cave millipede.Despite the lack of genetic studies on cave millipedes, there are a number of genetic studies on other cave organisms. These studies have shown that both vicariance and dispersal play a role in the speciation of cave organisms. Such research has also revealed an abundance of cryptic species among cave taxa (reviewed in Porter 2007; Juan et al. 2010).

Troglobiotic millipedes belonging to the genus Tetracion Hoffman, 1956(Callipodida, Abacionidae) inhabit caves of the southern Cumberland Plateau in Tennessee and Alabama. Three species have been described: Tetracion jonesi Hoffman, 1956, Tetracion antraeum Hoffman, 1956 and Tetracion tennesseensis Causey, 1959. They are relatively large (up to 8 cm in length) and exhibit troglomorphic characters including lack of pigmentation, reduced ocelli and elongated appendages (Peck 1989; Figure 1). Peck (1989) noted that they may be the most abundant scavenger species in Alabama cave communities. As in other callipodid species, Tetracion jonesi secretes p-cresol as a defensive compound (Peck 1989; Shear et al. 2010). It appears this defensive mechanism is effective as these millipedes are not part of the diet of Eurycea lucifuga Rafinesque, 1822 or Plethodon glutinosus (Green, 1818), two common cave salamanders (Peck 1974; Peck and Richardson 1976).

The type species for the genus, Tetracion jonesi, was described from Bat Cave, near Grant, Marshall County, Alabama (Hoffman 1956). In the same publication Hoffman also described the subspecies Tetracion jonesi antraeum, with a type locality of Barclay Cave, Madison County, Alabama. The subspecies were distinguished on the basis of differences in the structure of the male gonopods (Hoffman 1956). Tetracion tennesseensis was subsequently described from Warren County in central Tennessee; this species is smaller and inhabits caves at the northernmost range of the genus (Causey 1959). (Causey (1959, 1960) noted new records for Tetracion jonesii antraeum, extending its range into Jackson County, Alabama and Franklin County, Tennessee (Figure 2). Shear (1969) claimed that the morphological differences between Tetracion jonesii jonesi and Tetracion jonesii antraeum were insignificant and that subspecific status was not justified. Tetracion jonesii antraeum was subsequently elevated to species level status since its range did not overlap with that of Tetracion jonesii jonesi and therefore was completely isolated from this species (Shelley 1996). Hoffman (1999) reported that Tetracion is found in Georgia; we believe this is an error.

Because classification within Tetracion is based solely on morphological characters, we wanted to test if morphological species definitions corresponded to patterns of genetic variation. We also wanted to understand the population structure of these cave-obligate millipedes and determine if gene flow was occurring between populations. To do this, we sequenced part of the mitochondrial cytochrome oxidase 1 (CO1) gene from across the range of Tetracion tennesseensis and in the northern part of the range of Tetracion jonesi. We expected to find high genetic divergence between species and significant population structure within species. This is the first molecular phylogeographic study of troglobiotic millipedes.

Materials and methods Collecting

Tetracion specimens were collected from eleven caves on the southern Cumberland Plateau between 1996 and 2009. Specimens were collected from Franklin, Grundy, Warren, White and Van Buren Counties in Tennessee and Jackson County in Alabama. Due to the uncertain status of Tetracion antraeum, we refer to our specimens from Jackson County, Alabama and Franklin County, Tennessee as Tetracion jonesi. Our sampling covered the full range of Tetracion tennesseensis and the northern portion of the range of Tetracion jonesi (Figure 1). After the arrival of White Nose Syndrome to northern Tennessee in early 2010, we canceled all sampling in Alabama to avoid the possibility of spreading the fungus south. As a result we did not sample the southern portion of the range of this genus. Further sampling across the range of Tetracion in Alabama is required to determine if there is genetic evidence for the distinctiveness of Tetracion antraeum. We refer to all caves by their Tennessee Cave Survey or Alabama Cave Survey names. Due to the sensitive nature of cave habitats, cave coordinates are not published in this manuscript.

As Tetracion and Abacion Rafinesque, 1820are the only genera of the tribe Abacionini in the family Abacionidae (Shelley 1979; Hoffman 1999), we used an individual of the surface-dwelling species Abacion magnum (Loomis, 1943) as our outgroup. This individual was collected from Howard County, Maryland by J. Shultz. All millipedes were stored in 95% ethanol at -80°C.

DNA extractions, PCR, and DNA sequencing

We followed the manufacturer’s protocol from the DNeasy Blood and Tissue Kit (Qiagen; P/N: 69506) for DNA extractions. Millipede tissue was taken from the legs, antennae, or a section of the trunk. We used polymerase chain reactions (PCR) to amplify the mitochondrial CO1 gene. Amplitaq Gold PCR Master Mix (Applied Biosystems; P/N: 4318739) was added to all reactions and several different primer combinations were necessary for successful amplifications. The most successful amplifications for Tetracion jonesi were obtained using the primer combinations HCO1-tet (5’-GATATAGAATAGGATCTCCTCCAGC-3’) and LCO1-milli (5’-TCCACAAACCACAAAGACATTGG-3’), and for Tetracion tennesseensis HCO1-Tetten (5’-TCCTCCAGCGAGCAGGATCAAAGA-3’) and LCO1-Tetten (5’- ATTTTTGGAGCTTGAGCTGCCATG-3’). We cycled all reactions once for 5 min at 95°C and 35 times for 15 s at 95°C, 15 s at 50°C, and 1 min at 72°C. Occasionally, we lowered the annealing temperature from 50°C to 45°C to improve amplification. Successful PCR reactions were purified following the manufacturer’s protocol for the QIAquick PCR Purification Kit (Qiagen; P/N: 28106) and both strands were sequenced on an Applied Biosystems 3730 sequencer. All sequences have been submitted to GenBank (Accession #JN656558-JN656611).

Genetic analyses

Sequences were aligned and edited using Sequencher (v. 4.9; Gene Codes Corp., Ann Arbor, MI) and the number of indels, transitions and transversions were counted by eye. We used TCS (v. 1.21; Clement et al. 2000) to analyze haplotypes within species. This analysis required a standardized sequence length and therefore all sequences were trimmed to the length of our shortest haplotype. Mean pairwise distances between haplotypes and between populations were determined using MEGA (v. 4.0.2; Tamura et al. 2007). We also tested for population structure between all populations where we sampled five or more individuals using F statistics in Arlequin (v. 3.1; Excoffier et al. 2005).

To test for monophyly of Tetracion tennesseensis and Tetracion jonesi, we constructed a phylogenetic tree.We used MrBayes (v. 3.1.2; Ronquist and Huelsenbeck 2003) to conduct Bayesian phylogenetic analyses on a matrix of all haplotypes. We partitioned the data by codon position, and for each partition we used a General Time-Reversible (GTR) model with six substitution rates, estimated nucleotide frequencies, and a gamma distribution of rates. These model parameters were linked between partitions with the exception of the gamma parameter, which was unlinked for the 3^rd codon position partition due to the large number of changes in the 3^rd position relative to the 1^st and 2^nd codon positions. We calculated clade credibility values from 8000 trees by sampling every 1000^th tree from two runs of 5, 000, 000 trees after discarding the first 1001 sampled trees of each run. We also conducted a branch-and-bound parsimony bootstrap analysis (1000 replicates) in PAUP* (v. 4.0; Swofford 2001).

Results Sequencing results

We obtained sequences from 54 individuals: 17 from Tetracion jonesi, 36 from Tetracion tennesseensis, and 1 from Abacion magnum (Table 1). Average sequence length was 568 bp, with a maximum length of 571 bp and a minimum length of 541 bp; 30, 090 bp were sequenced in total. No indels, gaps, or stop codons were present and seven ambiguous bases were located. Within Tetracion, there were 55 variable sites of which 40 were transitions, 12 were transversions, and 3 were sites where both a transition and transversion had occurred. Of those variable sites, 41 were fixed differences between Tetracion jonesi and Tetracion tennesseensis including 27 transitions, 11 transversions, and 3 sites where both a transition and transversion had occurred. There was one fixed amino acid difference between Tetracion jonesi and Tetracion tennesseensis and there were five fixed amino acid differences between Abacion and Tetracion.

Intraspecific variation

Six haplotypes were present in Tetracion tennesseensis (T1 to T6; Figure 3). Genetic variation within Tetracion tennesseensis was low as haplotypes differed by no more than eight nucleotides. Clustering was apparent in the Tetracion tennesseensis network; haplotypes T1-T4 formed one cluster whereas haplotypes T5-T6 formed another cluster which differed from the first cluster by four nucleotides (Figure 3). Four haplotypes were limited to a single cave (T2, T3, T4, T6) and two were found in multiple caves (T1 and T5; Figure 4).The most common haplotype (T1) was shared among 14 out of 36 individuals. The least common haplotype (T2) was present in only one individual. Five of seven cave populations were fixed for a single haplotype (Table 1; Figure 4). Woodlee and Coppinger Caves each had populations containing two different haplotypes. In both cases the two haplotypes present in a single cave differed by a single nucleotide (Table 1; Figure 3). Consistent with the observation that most caves had distinct haplotypes, FST values between caves were high and significant (FST > 0.90, P < 0.05; Table 2) for twelve of fifteen comparisons. The only exceptions were the comparisons between Coppinger, Crystal and Woodlee Caves, which shared haplotype T1 (Table 1) and were non-significant (Table 2). Although caves generally did not share haplotypes, mean CO1 divergence between caves was low (<1.4% for all comparisons; Table 2).

Four haplotypes were present within Tetracion jonesi (J1-J4; Figure 3). Genetic variation within Tetracion jonesi was lower than within Tetracion tennesseensis as haplotypes differed by no more than four nucleotides. The most common haplotype (J2) was present in 7 out of 17 individuals; the least common haplotype (J3) was found in only one individual. Haplotypes were not shared between caves except J2 which was present in two caves (Jess Elliot and Little Crow Creek) and all cave populations were fixed for a single haplotype except for the Grapevine population which contained two haplotypes (J1 and J3; Figure 4), that differed by a single nucleotide (Figure 3). FST values between caves were high and significant for all comparisons (FST > 0.80, P < 0.01; Table 2). Mean genetic distance between caves ranged from 0.2-0.7% (Table 2).

Interspecific variation

Within Tetracion, interspecific genetic variation was much greater than intraspecific genetic variation. Kimura 2-parameter (K2P, Kimura 1980) distances between Tetracion tennesseensis and Tetracion jonesi haplotypes ranged from 8.2-9.2% with a mean of 8.8%. Our phylogenetic tree reflected these results as we found strong support for the monophyly of Tetracion tennesseensis and Tetracion jonesi (Figure 5). Our Bayesian and parsimony analyses supported the same topology (Figure 5). Mean K2P distance between Tetracion and Abacion was 20.9%.

Discussion

We found low intraspecific genetic variation in Tetracion, both in the number of haplotypes observed and in the amount of divergence between haplotypes. Only ten haplotypes were observed and populations were generally fixed for a single haplotype. As our sample size was small (mean = 4.8 individuals/cave), it is possible that further sampling would reveal greater diversity within populations. The maximum intraspecific divergence between haplotypes was 1.4%. Genetic variation was lower within Tetracion jonesi populations but that may be due to sampling bias as we did not survey the entire range of this species. Such low intraspecific genetic variation contrasts with other troglobiotic organisms from the Cumberland Plateau. For example, Snowman et al.(2010) found high genetic variation (up to 4%) between populations of the cave spider Nesticus barri Gertsch, 1984, and Dixon and Zigler (2011) found high genetic variation in several cave species (including a fly, a beetle, and an isopod) on a small scale. High levels of genetic variation were also observed in the cave crayfish genera Orconectes Cope, 1872 and Cambarus Erichson, 1846 (Buhay and Crandall 2005; Buhay et al. 2007), and in cave salamanders of the genus Gyrinophilus Cope, 1869 (Niemiller et al. 2008).

Although the differences between Tetracion populations were not great, populations were generally fixed for a single haplotype. Across eleven populations, only three haplotypes were shared between caves, and in each case the haplotypes were shared by geographically proximate caves. In the three caves with two haplotypes, the haplotypes in those caves differed by a single nucleotide, which is consistent with in situ evolution through mutation, as opposed to migration from a genetically distinct population. These patterns are similar to those observed in other terrestrial troglobionts from the Cumberland Plateau (Snowman et al. 2010; Dixon and Zigler 2011), but contrast with those of the cave crayfish and cave salamanders, where haplotypes were commonly shared among populations and across large distances (Buhay and Crandall 2005; Buhay et al. 2007; Niemiller et al. 2008). Aquatic subterranean habitats may be better connected than terrestrial subterranean habitats, or the longer lifespans of crayfish and salamanders may permit greater levels of migration and gene flow between populations.

Our molecular results were consistent with previous morphology-based species definitions. Interspecific genetic variation in Tetracion was high (8.8%). Several studies have used the CO1 gene to study millipede phylogeny and population structure, and the divergence between Tetracion jonesi and Tetracion tennesseensis was greater than that observed for most interspecific comparisons within Appalachioria Marek and Bond, 1996 (Swafford and Bond 2010), Parafontaria (Verhoeff, 1936) (Sota and Tanabe 2010), and Pseudotremia (Dixon and Zigler 2011). Although we do not have a molecular clock for Tetracion, if we use the often cited insect mitochondrial clock of 2.3% divergence/my (Brower 1994), or Papadopoulou et al.’s (2010) beetle CO1 clock of 3.5% divergence/my, we can estimate that Tetracion jonesi and Tetracion tennesseensis diverged several million years ago. Interestingly, the boundary between the ranges of Tetracion jonesi and Tetracion tennesseensis in northeast Franklin County coincides with species boundaries in two other genera of troglobiotic millipedes. Pseudotremia barri Lewis, 2005 and Scoterpes ventus Shear, 1972 are found north of the break point whereas Pseudotremia minos Shear, 1972 and Scoterpes stewartpecki Shear, 2010 are found south of it (Lewis 2005; Shear 2010). CO1 divergence between Pseudotremia barri and Pseudotremia minos, which may or may not be sister taxa, across that boundary is 3.7%, less than half that observed in Tetracion (Dixon and Zigler 2011).

Tetracion is an interesting contrast to other genera of troglobiotic millipedes known from the Cumberland Plateau and adjacent areas. Tetracion is not speciose, and each of its species has a large range spanning several counties and dozens of caves. Scoterpes, with fourteen troglobiotic species (Shear 2010), and Pseudotremia, with more than fifty species, most of which are troglobiotic (Lewis 2005; Shear 2008), are much more speciose. Many Scoterpes and Pseudotremia species are known from one or a few caves, and it has been suggested that Scoterpes species with large ranges may be ‘superspecies’ containing multiple distinct lineages worthy of recognition as species (Shear 2010). It is not clear why Tetracion exhibits this distinct pattern. Further research on the phylogeography and population structure of other troglobiotic millipedes, as well as sampling the full range Tetracion, may help unravel this mystery.