Phylogenetic analysis of the Australian trans-Bass Strait millipede genus Pogonosternum (Carl, 1912) (Diplopoda, Polydesmida, Paradoxosomatidae) indicates multiple glacial refugia in southeastern Australia

Abstract This study documents the first detailed phylogenetic analysis of an Australian paradoxosomatid millipede genus. Two mitochondrial genes (partial COI and 16S) as well as partial nuclear 28S rDNA were amplified and sequenced for 41 individuals of the southeastern Australian genus Pogonosternum Jeekel, 1965. The analysis indicates that five species groups of Pogonosternum occur across New South Wales, Victoria and Tasmania: Pogonosternum nigrovirgatum (Carl, 1912), Pogonosternum adrianae Jeekel, 1982, Pogonosternum laetificum Jeekel, 1982 and two undescribed species. Pogonosternum coniferum (Jeekel, 1965) specimens cluster within Pogonosternum nigrovirgatum. Most of these five species groups exhibit a pattern of high intraspecific genetic variability and highly localized haplotypes, suggesting that they were confined to multiple Pleistocene refugia on the southeastern Australian mainland. The phylogenetic data also show that northwestern Tasmania was colonized by Pogonosternum nigrovirgatum, probably from central Victoria, and northeastern Tasmania by an as yet undescribed species from eastern Victoria.

Thus, Pogonosternum occurs on both sides of Bass Strait, which separates mainland Australia from Tasmania. The paradoxosomatid genus Somethus Chamberlin, 1920 also has a trans-Bass Strait distribution (Jeekel 2006), as do the paradoxosomatid species Dicranogonus pix Jeekel, 1982 and Notodesmus scotius Chamberlin, 1920(Mesibov 2014. Many soil invertebrates, including millipedes, have limited active dispersal capabilities. Phylogenetic studies of southeastern Australian soil invertebrates can give important insights into the impact of glacial periods during the Pleistocene (Byrne 2008, Endo et al. 2014, Garrick et al. 2004, Schultz et al. 2009, Sunnucks et al. 2006 and assist in identifying biogeographic barriers (Chapple et al. 2011). Unfortunately, phylogenetic studies of Australian millipedes are rare and restricted to a few taxa from a small number of localities (Adams and Humphreys 1993, Nistelberger et al. 2014, Wojcieszek and Simmons 2012. For the australiosomatine species Orocladosoma kosciuskovagum (Brölemann, 1913) from the Australian Alps a hypothesis of multiple glacial refugia has been proposed (Endo et al. 2014) to explain the results of such studies. Similarly, the australiosomatine genus Somethus in South Australia was found to have high morphological and genetic variability within species was discovered: it seems probable that isolation in multiple glacial refugia during the Pleistocene was the evolutionary driving force for this variability (Decker 2016).
The present study documents a molecular phylogenetic analysis of the antichiropodine genus Pogonosternum, using specimens from across the genus range, and with molecular evidence indicating past isolation in multiple Pleistocene refugia. Finally, the identity and origin of Tasmanian Pogonosternum populations are clarified.

Specimen collecting and preservation
Pogonosternum specimens were collected by hand in Victoria and New South Wales in August 2014 by the author, Karin Voigtländer and Robert Mesibov, and by Mesibov in Tasmania in May 2014 and May 2015 (Fig. 1). Most sites were searched for 1-5 hours with the aim of finding 1-3 adult males. At only a few localities were Pogonosternum found to be abundant. Specimens were killed and stored in 95% ethanol, with a change of ethanol after 1-2 months. Full details of locality, date, collector, collection number and coordinates (WGS84 decimal degrees) are provided in Suppl. material 1.  Table 1 and Suppl. material 1 for further details). P. adrianae (light green), P. laetificum (green), P. nigrovirgatum s. l./coniferum (yellow), P. sp. A (red), P. sp. B (blue).

Illustrations
Maps were created with ArcGIS 10. The final phylogenetic trees were edited using Adobe Illustrator CS4.

Molecular analysis
DNA was extracted from 2-4 legs from each of 41 Pogonosternum specimens and from the three paradoxosomatid species Archicladosoma magnum Jeekel, 1984, Somethus scopiferus Jeekel, 2002 andS. castaneus (Attems, 1944), which were chosen as outgroups (Table 1). Total genomic DNA was extracted using the Qiagen DNAeasy Blood&Tissue kit following the standard protocol except that tissue was incubated for 48h.
Glom primer cocktail pairs (Decker 2016, Macek et al. 2014) were used to sequence a 618 bp fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene. Primer pairs 28S D1a (Fw) and 28S D3b (Rv) (Dell'Ampio et al. 2009) were used to amplify 1225 bp of the D2 fragment and adjacent areas of D1 and D3 on the nuclear 28S ribosomal RNA gene.
All fragments were sequenced in both directions by the BiK-F Laboratory Centre, Frankfurt, Germany. All obtained sequences were checked via BLAST searches of GenBank; no contamination was discovered. The sequences were aligned by hand in ClustalX ver. 1.83 (Chenna et al. 2003) and uploaded to GenBank (Table 1).
Some homologisation problems in the 16S rRNA sequences arose mainly because of the highly variable expansion loops. As a result, selected alignment positions (272-297) were excluded from the 16S rRNA dataset for all further analyses using MEGA6.
The final alignments consisted of 618 bp of COI mtDNA, 540 bp of 16S rRNA and 1206 bp of 28S rRNA. The combined datasets after these exclusions consisted of 1158 bp for COI+16S. Individual partial alignments can be obtained from the author upon request. The alignment of the combined dataset can be found in the Suppl. material 2 as a FASTA file.
COI and 16S sequences were combined as a single dataset and incongruence assessed between the mtDNA intergenic spacer sequences with the incongruence length difference (ILD) test (Farris et al. 1994) implemented as the partition homogeneity test  in PAUP* version 4.0b10 using a full heuristic search, 10 random taxon addition replicates, tree-bisection-reconnection (TBR) branch swapping, and with MaxTrees set to 100 (Swofford 2002). The best-fit model of nucleotide substitution for the individual COI and 16S dataset was determined by MrModelTest 2 (Nylander 2004). The bestfit model of nucleotide substitution selected using MrModelTest 2 was the General Time Reversible model with gamma distribution and proportion of invariant sites (Nei and Kumar 2000) for the individual COI and 16S dataset. The trees constructed from individual genes did not show significant conflicts in topology (nodes different among trees with support > 70% in ML) and no significant incongruence among the three genes was revealed by the ILD test (P > 0.83 in all of the pairwise comparisons), so the sequences were concatenated into a dataset containing 1158 characters for phylogenetic analysis.
The combined dataset of COI and 16S was analysed under maximum likelihood (ML) using MEGA6 (Tamura et al. 2011) and Bayesian inference (BI) using MrBayes version 3.2 (Ronquist et al. 2012). For ML analysis, three independent runs were performed with nodal support estimated from 1000 bootstrap (BP) pseudoreplicates using the best-fit model for the concatenated dataset. For Bayesian analysis, two independent runs were carried out with four differentially heated Metropolis-coupled Monte Carlo Markov chains for 10 000 000 generations started from a random tree and chains were sampled every 100 generations.
Multiple runs of ML and BI converged in trees with the same topology and similar likelihood score so that only the result of the first run is presented. The topology resulting from ML and BI analyses was largely congruent except for the arrangements of several terminal nodes with low support. Thus, results from the ML and BI analyses are shown together based on the ML tree with bootstrap (BP) and posterior probabilities (PP) of the major lineages shown on the corresponding branches with BP values > 70 (Fig. 2).
An appropriate DNA substitution model was determined for 28S under the Bayesian Information Criterion (BIC) in Modeltest implemented in MEGA 6 (Tamura et al. 2011). The lowest Bayesian Information Criterion score (BIC) was obtained for 28S rRNA (BIC 3875.11) with the Tamura 3-parameter model (Tamura 1992).
A phylogenetic hypothesis was inferred for COI+16S and 28S by using the maximum likelihood method conducted in MEGA6 (Tamura et al. 2011). The phylogenetic tree with the highest log likelihood (COI+16S: -7237.4280; 28S: -1831.9238) is shown (Figs 2, 3). Initial trees for the heuristic search were obtained by applying the neighbor-joining method to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach (Tamura et al. 2004). A discrete Gamma distribution was used to model evolutionary rate differences among sites (five categories (+G, parameter = COI+16S: 0.2338)). The bootstrap consensus tree inferred from 1000 replicates (Felsenstein 1985) is here used as the best estimate of the phylogeny of the analyzed taxa (Figs 2, 3).
Mean uncorrected pairwise distances between terminals (transformed into percentages) were determined using MEGA6 (Tamura et al. 2011) and can be found in Suppl. material 3.

Phylogenetic and distance analysis
The monophyly of the genus Pogonosternum is strongly supported (ML BP = 97; BI PP = 1.0) in the mitochondrial tree and shows five clades within Pogonosternum, resembling five species groups (Fig. 2).
Pogonosternum nigrovirgatum sensu lato with a trans-Bass Strait distribution formed a well-supported (ML BP = 89; BI PP = 1.0) sister clade to the new species P. sp. A (ML BP = 98; BI PP = 1.0) that also has a trans-Bass Strait distribution. Pogonosternum sp. A also occurs in New South Wales (Car 2010) and in northeast Tasmania (Mesibov & Churchill 2003). Pogonosternum nigrovirgatum s. l. occurs on mainland Australia (Otway Ranges to eastern Victoria) and in northwest Tasmania. Pogonosternum coniferum clusters with another form with intermediate gonopods (referred to as P. cf. nigrovirgatum in Fig. 2) between P. nigrovirgatum sensu stricto and P. coniferum.
Both P. nigrovirgatum s. l. and P. sp. A show high intraspecific distances ranging from 1.8 to 6.8% within P. nigrovirgatum s. l. and 1.1 to 5.9% within P. sp. A.
Within the P. nigrovirgatum s. l. species-group, the greatest genetic distances were observed between populations in the Strzelecki Ranges (S60, S63; ML BP = 100; BI PP = 1.0) and more western populations, with values ranging from 5.0 to 6.8%. Specimens from the Otway Ranges (S77, S78, S81) all formed a well-supported cluster (ML BP = 86; BI PP = 1.0). The Tasmanian specimen (X2) was distinct from both the Strzelecki Ranges (5.4-6.0%) and central and western Victorian specimens (3.7-3.8%). In the case of Pogonosternum sp. A the largest distances (4.2-5.8%) were between the Eastern Gippsland populations (S42, S47; ML BP = 100; BI PP = 1.0) and all other specimens. The status of the northeast Tasmanian specimen is not well resolved; it is closest to a population from Kosciuszko National Park (S31, 3.0%), the two forming a poorly supported sister clade with a specimen from Gippsland (S52; ML BP = 55; BI PP = 0.6).
All species show considerable intraspecific genetic distances and high phylogeographic structure, especially P. laetificum, and, except in the case of P. adrianae, no haplotypes are shared between different populations. Additional one to three sequenced specimens from eight sampling sites (S14, S15, S22, S58, S59, S78, S83, S87) always showed the same haplotype in Pogonosternum (data not published).
Owing to the general lack of variability within the nuclear 28S rRNA dataset, the phylogenetic relationships among species were largely unresolved. Distances for 28S rRNA within Pogonosternum are very low, with a maximum of three base pair differences noted for P. sp. B (Fig. 3). Only the two condensed sister clades of P. nigrovirgatum + P. sp. A and P. adrianae + P. laetificum, as well as P. sp. B are shown.

Morphology
In a separate paper (Decker, in preparation), the morphology of the Pogonosternum species groups is described in detail and new species are described, based on the specimens used here and from ca 130 additional localities. Here I note briefly that several common morphological features were observed in the gonopods of P. nigrovirgatum s. l., P. laetificum, and P. sp. A: some specimens also showed intermediate states of those features (Fig. 2). It was found, however, when additional material was examined from each population that the morphology of each population was locally stable. It was only in rare cases in the Otway Ranges and NW Tasmania populations that two gonopod morphs occurred in one place.
Surprisingly, gonopod morphology did not appear to agree well with the phylogenetic tree (Fig. 2). Various gonopod forms were distributed with no apparent phylogeographical correlation. Only the species P. adrianae and P. sp. B showed stability in both gonopods and some other non-gonopodal characters over their distribution area, even when material from other museum collections was included (Decker, in preparation).

Phylogenetic analysis
The mitochondrial tree (Fig. 2) shows five main clades, suggesting five species. Pogonosternum coniferum clustered within P. nigrovirgatum, and its taxonomic status needs re-examination (Decker, in preparation).
The 28S tree shows little or only little resolution at the species level (Fig. 3), but was useful in identifying sister clades. This result contrasts with that from a study of the paradoxosomatid genus Somethus in South Australia, in which the 28S gene was used successfully for species identification (Decker 2016). Future studies on other Australian Paradoxosomatidae will reveal if 28S is useful as a diagnostic nuclear gene at the species level.

Morphological variability
With the exception of P. adrianae and P. sp. B, Pogonosternum species show significant variability in gonopod form, with local morphs occurring throughout each species' distribution area.
Interestingly, P. adrianae is morphologically distinct (in size, spiracles, male tibiotarsal brushes and gonopods, female coxal process) from P. laetificum despite their close genetic distance.
Gonopod variability was also documented for some species of Somethus in South Australia (Decker 2016) and Stygiochiropus Humphreys & Shear, 1993 from Western Australia (Humphreys and Shear 1993). Another good example of variability is seen in the trans-Bass Strait (eastern Victoria, NE Tasmania) paradoxosomatid millipede, Dicranogonus pix: while this species shows only slight variability in gonopods there is marked variation in the development of their paranota. Individuals with no paranota are separated from those with keels by a gap between the Kent and Furneaux Groups of islands (Mesibov 2014).
This study has shown that in the area of southern and southeast Australia, there are at least two genera, Pogonosternum and Somethus (Decker 2016), which both show variability in morphology and genetics. Poor sampling and too few specimens could lead to incorrect conclusions and unnecessary multiple species descriptions.

Multiple glacial refugia in southeastern Australia
The results indicate that there is high intraspecific genetic divergence, with high genetic distances and haplotype diversity in the mitochondrial genes between populations of Pogonosternum, even those adjacent to each other. The P. laetificum clade, which has been sampled extensively in the Central Highlands, shows particularly high intraspe-cific genetic differences (mean genetic distance of 3.9%), apparently without corresponding geographic patterning, or morphological variation (Decker, in preparation).
The phylogenetic patterns with high intraspecific divergence, high genetic distances, and haplotype diversity with unique local haplotypes, resulting in long branches, shown by Pogonosternum, indicate multiple Pleistocene refugia according to Byrne (2008). These refugia provided suitably moist habitats in which species could persist during the dry, cold climate cycles of the Pleistocene period in southern Australia, while glaciation was limited to the alpine areas of the Great Dividing Range and Tasmania (Barrows et al. 2002). Moderate to high genetic diversity prior to these cycles can be assumed for poorly dispersing millipedes, through isolation by distance, and it is likely that populations were isolated within refugia, leading to further genetic diversification. In contrast, contractions to one or few major refugia during cold, arid periods would result in a low genetic diversity, few divergent lineages and low haplotype diversity, with few haplotypes in areas of postglacial recolonisation (Byrne 2008).
The phylogenetic patterns shown by Pogonosternum suggest that in Victoria and New South Wales there were large areas with multiple local refugia during the Pleistocene. No region in the study area on mainland Australia showed results which indicate rapid postglacial resettlement of Pogonosternum.
Evidence for multiple glacial refugia was also identified in the spirostreptidan millipede Atelomastix bamfordi Edward & Harvey, 2010 in Western Australia (Nistelberger et al. 2014) and for some species of Somethus in South Australia (Decker 2016). Similar phylogeographic patterns seem to occur in other soil invertebrates with limited dispersal capacities in southern Australia, for example flatworms (Sunnucks et al. 2006) and springtails (Garrick et al. 2004). Endo et al. (2014) have suggested, however, that glacial periods have had less of an impact on the distribution and genetic diversity of invertebrate groups (Coleoptera, Orthoptera, Collembola, Diplopoda) in the Australian Alps than they have in alpine systems in the Northern Hemisphere.
However, further studies on genetic and morphological variability on a finer geographical scale could lead to a better understanding of the pattern and impact of isolation in multiple glacial refugia during the Pleistocene, also as an evolutionary driving force for morphological variability in some species.

Gippsland phylogeography
There is a notable high genetic distance gap within P. nigrovirgatum sensu lato between specimens from the Strzelecki Ranges (S60, S63), West Gippsland, and those sampled in the central and western regions in Victoria, but some specimens of adjacent populations from the latter (S64, S65) were morphologically indistinguishable from specimens from the Strzelecki Ranges. A similar genetic gap was observed in P. sp. A for the populations in Eastern Gippsland east of Orbost (S42, S47) and all other populations. These two cases indicate that these areas may have been isolated for long periods from neighboring regions, possibly before the Pleistocene, perhaps during a marine incursion in the Gippsland Basin and other parts of southeast Australia close to the Miocene-Pliocene boundary (Dickinson et al. 2002).

Trans-Bass Strait distribution
The genus Pogonosternum shows a trans-Bass Strait distribution and most likely originated in mainland southeast Australia, since the highest species diversity is found on the mainland and the two Tasmanian branches occupy only very subordinate positions on the tree (Fig. 2). Tasmanian populations of this genus are restricted to the northeast and northwest corners of the Tasmanian mainland and neighboring islands, and presumably dispersed from Victoria when it was largely connected with Tasmania during the Pleistocene (Lambeck and Chappell 2001). Mitochondrial data suggest that the sequenced population of P. nigrovirgatum s. l. in northwest Tasmania was most likely derived from one in central Victoria or the Otway Ranges. While the results for P. sp. A from northeast Tasmania do not show a close relationship to coastal Victorian populations, analysis of 16S (data not included here) including sequences from two other localities in the western part of East Gippsland showed the Tasmanian specimen clustering with the latter. This indicates that the settlement of Tasmania by this species started in the Gippsland region. A remarkably similar distribution to that of P. sp. A across Bass Strait is also known for the paradoxosomatid millipedes Dicranogonus pix and Notodesmus scotius (Mesibov 2014).
Further studies using more sampling localities in Tasmania and its islands could indicate points of origin in Victoria and the timing of millipede settlement of Tasmania.