Material and methods
The data set for this study is reported in Table 1. The sampling localities of the studied species in the Mediterranean basin are reported in Fig. 1. Genomic DNA was extracted from entire specimens preserved in absolute ethanol using the Qiamp Tissue kit (Qiagen Inc., USA), following the manufacturer’s protocol.
Fig. 1. Distribution of the analysed Mediterranean taxa, Antrolana lira and Speocirolana bolivari from America are excluded. Abbreviations as in Table 1.
mtDNA sequencingnext section
PCR amplification products were obtained from the 12S and 16S mitochondrial gene portions for all the taxa included in the analysis, with a few exceptions (Table 1). For 12S, we used the primers 12Sbi (5’–AAGAGCGACGGGCGATGTGT–3’) (Simon et al., 1994) and 12SL4 (5’–GTGCCAGCMGCCGCGGTTA–3’) (Schubart et al., 2006). Primers used to amplify the 16S gene portion were 16Sar (5’–CGCCTGTTTATCAAAAACAT–3’) and 16Sbr (5’–CCGGTCTGAACTCAGATCACACGT–3’) (Palumbi et al., 1991). The taxa T. haouzensis, M. delamarei and S. virei were also analysed with a portion of Cytochrome Oxidase I using the primers mtd10 universal primer 5’ –TTGATTTTTTGGTCATCCAGAAGT -3’ (Roehrdanz, 1993) and Florence 5’- CCTAAAAAATGTTGAGGGAA-3’ (Baratti et al., 2005). The PCR profiles for 12S/16S rRNA amplifications were described by Baratti et al. (2004) and for COXI by Baratti et al. (2005). PCR products were run on a 1.5% agarose gel, containing 0.5 g/ml ethidium bromide. The amplification patterns were analysed with the Gel Analysis Program v. 2.0 (Ultraviolet Products Ltd.). PCR products were purified (ExoSAP-IT, Amersham Biosciences), sequenced with a Perkin-Elmer sequencing kit (ABI Big Dye Terminator Cycle Sequencing v. 2.0-ABI PRISM, Applied Biosystems, Foster City, USA) and analysed with an ABI 310 automated sequencer (Applied Biosystems). All the sequences are deposited in GenBank with the accession numbers reported in Table 1.
Table 1. Species used in this analysis. The new data set is underlined; the other is from a previous study (Baratti et al., 2004). Asterisks indicate taxa from other authors (see GenBank accession number).
Sequence analysis and nucleotide diversity
Approximately 480 bp of the 12S rRNA gene and 575 bp of the 16S rRNA gene were sequenced. Unreliable sequences were obtained for the 16S gene portion in six individuals of the taxon Turcolana sp., which was not included in the combined data set tree. For COXI, a portion of 432 nucleotides was amplified. All regions were sequenced in both directions. Electrophenograms were visualised with CHROMAS 1.45 (http://www.technelysium.com.au). The sequences were manually corrected and analysed with ProSeq 2.9 Beta (http://helios.bto.ed.ac.uk/evolgen/filatov/proseq.html) and then aligned using ClustalX 1.81 (Thompson et al., 1997). Multiple alignments were obtained with ClustalX (Jeanmougin et al., 1998) by assigning different gap penalty values (opening: 5, 15, 25; extension: 0.2, 5, 10). The resulting alignments were used to construct tree topologies with the different combinations of gap values. Polymorphic sites, nucleotide statistics and genetic distances (K2P) were analysed with MEGA 2.1 (Kumar et al., 2001).
A chi-square test of homogeneity of base frequencies across taxa was carried out using PAUP 4.0, b. 10 (Swofford, 2001). We executed the likelihood mapping method (Strimmer and von Haeseler, 1997) using TREE-PUZZLE (Schmidt et al., 2002) to test the a priori phylogenetic signal in the mtDNA portions studied. Phylogenetic congruence among 12S and 16S data partitions was also performed using the partition-homogeneity test implemented in PAUP* (Swofford, 2001).
Testing of the evolutionary model that best fits our data was conducted with MODELTEST 3.04 (Posada and Crandall, 1998), based on a likelihood ratio test. Different models of nucleotide substitutions were fitted to each data set and for the combined data set. For 12S, the TRN model was the best one selected (Tamura and Nei, 1993), while for 16S the GTR model was selected (Tavaré, 1986). Both models were corrected for rate heterogeneity among sites with a Gamma (G) distribution (Yang, 1993). For the 12S and 16S data sets, the GTR+G model was selected. We carried out a phylogenetic reconstruction by Maximum Parsimony (MP) (Kluge and Farris, 1969) and Neighbour-Joining (NJ) (Saitou and Nei, 1987) analyses using PAUP and by the Bayesian method using MrBayes 3.0B4 (Huelsenbeck and Ronquist, 2001). MP and NJ analyses were performed separately for each set of DNA sequences and in a combined analysis (total evidence approach) for the taxa sequenced for both genes. Neighbour-Joining trees were constructed with distances computed with the best-fit model obtained with MODELTEST. Parsimony analysis was carried out using the heuristic search algorithm, using 100 random-taxon-replicates for all analyses. The analysis was performed with ACCTRAN optimization and tree bisection TBR branch swapping, considering all characters as unordered and equally weighted, and gaps treated as fifth state. A strict consensus tree was calculated when there was more than one tree. Branch supports were assessed by 1000 non-parametric bootstrap replicates. Non-parametric bootstrapping with heuristic searches of 2000 replicates for MP and NJ was used to assess confidences of branches in MP and NJ. A Bayesian analysis was also performed with MrBayes 3.0B4, with clade support assessed by posterior probability. Four Markov chains, one heated and three cold, were allowed to run for two million generations using random starting trees, trees were sampled every 100 generations, with a burnin amounted to 20%.
Regarding the choice of outgroups for the Mediterranean stygobitic cirolanids, no marine cirolanid isopods were available to us. However, the direct marine ancestors of stygobitic cirolanids are unknown at present. The marine cirolanids described for the Mediterranean sea are 12 species belonging to three genera (Eurydice, Cirolana, Natatolana, Bruce, 1986; Keable, 2006). No samples of these taxa were collected during our research activities and we used as outgroup a marine cirolanid taxon present in the GenBank Database: Cirolana rugicauda Heller, 1861.
Divergence time estimates
A likelihood-ratio test (LRT) was performed to test the molecular clock hypothesis based on 16S gene sequences. The likelihoods with the molecular clock assumption (L0) and without the molecular clock assumption (L1) were calculated using PHYLIP (Phylogeny Inference Package 3.66, Felsenstein, 2005). The molecular clock hypothesis was never rejected between all the Typhlocirolana species and Marocolana, while it was always rejected in comparisons involving Sphaeromides, Antrolana, Speocirolana, Faucheria and Saharolana.
Since the clock-like behaviour of the data was rejected for the complete dataset, a relaxed clock (Uncorrelated Lognormal) was applied, as implemented in BEAST 1.5.3 and accompanying utilities (Drummond and Rambaut, 2007) using a Bayesian MCMC approach. The model of sequence evolution was set according to MODELTEST, and a Yule tree prior was applied. In the absence of suitable fossil calibration points, a strong prior was set for the rate of sequence evolution, applying a mean rate of 0.00325 changes per site per myr. This corresponds to a pair-wise divergence rate of 0.65% per myr, proposed for the 16S gene of crustaceans, in particular of isopods, by Held (2001). Following preliminary runs, the final analysis consisted in two independent runs of 10 million generations each, logging trees and parameters every 1000 generations. Effective sample size (ESS) values, parameter traces and frequency plots were examined with Tracer 1.5 (Rambaut and Drummond, 2003) to confirm stability as well as consistency between the two runs. The initial 10% of generations were removed from the analysis as burnin. Log and trees files were combined with LogCombiner and the posterior distribution of node presence and node height were summarized with TreeAnnotator over the maximum clade credibility tree.