Phylogeny and Evolutionary Patterns in the Dwarf Crayfish Subfamily (Decapoda: Cambarellinae)

The Dwarf crayfish or Cambarellinae, is a morphologically singular subfamily of decapod crustaceans that contains only one genus, Cambarellus. Its intriguing distribution, along the river basins of the Gulf Coast of United States (Gulf Group) and into Central México (Mexican Group), has until now lacked of satisfactory explanation. This study provides a comprehensive sampling of most of the extant species of Cambarellus and sheds light on its evolutionary history, systematics and biogeography. We tested the impact of Gulf Group versus Mexican Group geography on rates of cladogenesis using a maximum likelihood framework, testing different models of birth/extinction of lineages. We propose a comprehensive phylogenetic hypothesis for the subfamily based on mitochondrial and nuclear loci (3,833 bp) using Bayesian and Maximum Likelihood methods. The phylogenetic structure found two phylogenetic groups associated to the two main geographic components (Gulf Group and Mexican Group) and is partially consistent with the historical structure of river basins. The previous hypothesis, which divided the genus into three subgenera based on genitalia morphology was only partially supported (P = 0.047), resulting in a paraphyletic subgenus Pandicambarus. We found at least two cases in which phylogenetic structure failed to recover monophyly of recognized species while detecting several cases of cryptic diversity, corresponding to lineages not assigned to any described species. Cladogenetic patterns in the entire subfamily are better explained by an allopatric model of speciation. Diversification analyses showed similar cladogenesis patterns between both groups and did not significantly differ from the constant rate models. While cladogenesis in the Gulf Group is coincident in time with changes in the sea levels, in the Mexican Group, cladogenesis is congruent with the formation of the Trans-Mexican Volcanic Belt. Our results show how similar allopatric divergence in freshwater organisms can be promoted through diverse vicariant factors.


Introduction
The freshwater crayfish subfamily Cambarellinae is comprised of the unique genus Cambarellus, with 17 recognized species and a disjunctive distribution across the freshwater streams of the Gulf Cost of the United States and North and Central México (Fig. 1) [1]. The subfamily is unique because of the exceptionally small body size of its species. They typically reach only 4 cm compared to most crayfish averaging a maximum body size of .5 cm; hence, the reference to the genus as the ''Dwarf'' crayfishes. Their distribution goes from the Swanee River in northern Florida, eastward through the southern Mississippi River watershed to southern Illinois and continues southwest to the Nueces River in Texas [2,3]. In México, Cambarellus has a discontinuous distribution with three distant and isolated populations from the northern states of Chihuahua, Coahuila and Nuevo León and then along the Trans-Mexican Volcanic Belt (TMVB) [4,5]. The genus contains species largely inhabiting lakes and lentic habitats. The evolutionary history of such a broad and disjunct distribution of species is unclear and our goal with this study is to shed some light on the biogeography of the Cambarellinae.
A series of apomorphic morphological characters define the subfamily and, therefore, monophyly has been accepted since its proposal. These include, as for other crayfish groups, genital morphology, which is particularly important, but also a small body size, specific branchial formula, movable and enlarged annulus ventralis (female genitalia) and the absence of the cephalic process in the first pair of pleopods (male genitalia) [1,2,6]. The morphological unity of these characters that define the subfamily contrasts with the wide morphological variation in other characters described for populations of several species [2,5,7]. This diversity within and among species makes designation and identification difficult, especially for widely distributed species [2,5].
Despite the intriguing geographic distribution and species diversity in the Cambarellinae, the only phylogenetic hypothesis for species relationships in the group is based on phenotypic information and genital morphology [2]. With this hypothesis (Fig. 2) three subgenera were proposed; Pandicambarus (containing seven species), the monotypic Dirigicambarus, both comprised of species occurring north of the Rio Grande (the Gulf Group), and Cambarellus, containing species south of the Rio Grande (the Mexican Group) [2]. However, no apomorphic characters have  Table 1. Sample locations are colored to represent different clades recovered by phylogenetic analyses (see Fig. 3). Open circles correspond to the only locality records for the three species not included in the analyses as they were not found during sampling, or did not amplify during PCR reactions. Gray background refers to elevation (500-6000 m). doi:10.1371/journal.pone.0048233.g001 Evolutionary Patterns in Cambarellinae Table 1. Sampling localities and Genbank accession numbers from individuals of Cambarellus used in this study. + Species id from this study Subgenus (Fitzpatrick, 1983) GeneBank accession numbers  been proposed to support these subgeneric classifications and no formal phylogenetic hypothesis has been evaluated using either molecular or morphological characters. Therefore, we propose to estimate a robust phylogenetic hypothesis for the group using an extensive molecular data set. We then use this phylogenetic framework to evaluate a coherent taxonomy for the group and to test biogeographic hypotheses regarding the origin and spread of the dwarf crayfish.
We also examine diversification patterns in the subfamily through the estimated phylogenetic history of the species within the subfamily. Phylogenetic diversity patterns are impacted by geographic features and geologic history due to their effects on allopatric speciation [8]. Given the contrasting geographical features (Fig. 1) coupled with their distinct geological histories occupied by the different groups in Cambarellinae, we will use reconstructed molecular phylogenies to serve as models of lineages through time (LTT), that will allow us to test the tempo and pattern of change across lineages [9,10,11]. In the present study, we used our molecular dataset on the subfamily Cambarellinae to infer the timing and mode of lineage accumulation (patterns of speciation minus extinction) which allows us to determine whether there have been contrasting patterns in rates of diversification between the two geographical components of this group; namely, those defined as the Gulf and Mexican Groups, as a result of contrasting biogeographic histories. Finally, we identify a geological timescale consistent with biogeographic factors and cladogenetic events in this group.

Sampling and Sequencing
No specific permits were required for the described field studies, as none of the studied species were included in any endangered list, at national or international levels at the time of sampling (comprising the years 2005 and 2006). Including field and museum localities, 59 geographic locations covering 14 of the 17 species were collected throughout the distributional range of the subfamily Cambarellinae ( Fig. 1). Taxonomic identification was carried out using existing keys [12]. The two main ranges for the subfamily were covered, along the Neartic and the Transition zone of North America, from the Mississippi River basin to the TMVB in central México. Most of the species could be sampled, but those tissues from species with very restricted distribution ranges and/or being collected in a reduced number of times in wild were obtained from museum specimens (National Museum of Natural History, Smithsonian Institution) ( Table 1). Detailed data about samples included are summarized in the Table S1.
The central goal of this work is to estimate a robust phylogenetic hypothesis for relationships among the species within the subfamily to test taxonomic hypotheses, biogeographic hypotheses, and speciation hypotheses. As phylogenies are most accurately estimated using broad taxonomic sampling as well as extensive character sampling, we attempted to sample all species within the subfamily (but are missing three of them) and collected sequence data from five different gene regions (three mitochondrial and two nuclear). We sequenced the mitochondrial genes 16S rDNA (16S), 12S rDNA (12S) and Cytochrome Oxidase subunit I (COI). These genes have good phylogenetic signal in crustaceans [13] and are considered optimal choices to characterize the genetic variation in crustacean groups. Nuclear genes sequenced were 28S rDNA large ribosomal unit (28S) and Histone 3 (H3) gene, which also have some variation among species and are particularly good at discerning deeper nodes [13].
PCR amplifications using gene specific primers ( Table 2) were carried out in 25 mL reactions containing: 16PCR buffer, 0.5 mM of each primer, 0.2 mM of each dNTP, 1.5 mM MgCl 2 , 1 U Taq Figure 2. Morphologic hypothesis tested. Phylogenetic hypothesis based on morphologic analysis of the monotypic subfamily Cambarellinae (genus Cambarellus), indicated are the subgenera previously proposed, mainly based on genital morphology (Fitzpatrick, 1983). doi:10.1371/journal.pone.0048233.g002 Evolutionary Patterns in Cambarellinae polymerase (Biotools), and about 10-50 ng of template DNA. The cycling profile for PCR amplifications was 3 min at 94uC (1 cycle), 30 s at 94uC, 30 s at the primer-specific melting temperature and 60 s at 72uC (30 cycles), followed by a final extension of 4 min at 72uC. PCR products were visualized in 1.0% agarose gels (16TBE) and stained with SYBR-Safe (Invitrogen). Fragments were sequenced on an ABI 3730XL DNA Analyzer. Sequences of the different gene fragments were aligned using MUSCLE [14]. In the case of the COI gene, recommendations to detect the occurrence of possible nmtDNA were carried out for each sequence. These included the identification of stop codons, repeated sequencing of samples, nonsynonymous substitution and unusual levels of genetic divergence in samples from the same population [15,16].
In order to identify the most appropriate evolutionary model of nucleotide substitution (Table 2), we considered the Akaike corrected information criterion (AICc) [18], and the Bayesian Information Criterion (BIC) [19] as estimated using the program jModeltest [20]. A phylogenetic tree was constructed under ML using PHYML 3.0 [21] and AICc-selected parameters for the concatenated matrix. The tree search was started with an initial BIONJ tree estimation followed by a Subtree Pruning and Regrafting (SPR) topological moves algorithm. We assessed confidence in branches using 1000 nonparametric bootstrap [22] replicates under the best-fit evolutionary model.
Bayesian inference of phylogeny was implemented in MrBayes v. 3.1.2 [23], following the BIC-selected parameters and applying a Monte Carlo Markov Chain (MCMC) search procedure for 10 million generations. Sequences were partitioned by codon position for COI and by gene for the rest of fragments, using the parameters found by BIC as priors and unlinking the run parameters. Convergence between the different run parameters in paired simultaneous runs (4 chains by run), trees were sampled every 100 generations and run length was adjusted considering an adequate sampling based on average standard deviation of split frequencies being ,0.01 [24]. We examined the results and determined the burn-in period as the set of trees saved prior to log likelihood stabilization and convergence as estimated using Tracer 1.4.1 [25], eventually the first 10% trees. Tracer was also used to check for convergence between chain runs and optimal values of run parameters. Confidence in nodes was assessed from the posterior probabilities along the MCMC run. Highly supported nodes are termed herein as those with a value of 95% or more in posterior probabilities and bootstrap values.
We tested our resulting topology against the phylogenetic hypotheses put forth by Fitzpatrick [2]; namely, the three subgenera are monophyletic and show the following relationships ((Dirigicambarus, Pandicambarus),Cambarellus). Topology constrained ML scores were estimated for each hypothesis in PAUP*. Congruence with alternative hypotheses was evaluated in a ML framework applying the Shimodaira-Hasegawa (SH; [26]) test and the Approximate Unbiased (AU) test [27] with 50,000 RELL bootstrap replicates as implemented in TreeFinder [28]. We also tested these hypotheses using a Bayesian approach by identifying the alternative hypothesis within the set of Bayesian tree topologies and testing for significant differences. To do so, we filtered the post-burnin Bayesian topologies included in the set of trees with the constraint topology in PAUP* [17].

Divergence Dating
In order to propose an accurate time frame for phylogenetic divergence processes, we estimated mean node ages and their 95% highest posterior densities (HPDs) using Bayesian relaxed molecular clock methods [29] as implemented in BEAST ver. 1.6.1 [30]. In this method, tests of evolutionary hypotheses are not conditioned on a single tree topology, which allows for simultaneous evaluation of topology and divergence times while incorporating uncertainty in both. A uniform Yule tree prior was specified, as appropriate for hierarchical rather than reticulate relationships, and a subsampling of one representative of every lineage was included to avoid over-representation of certain individual lineages with more sampling. We applied the optimal model of data partitioning and DNA substitution identified by BIC for each gene (COI, 16S, 12S, 28S and H3) and for codon positions for COI. An uncorrelated relaxed lognormal molecular clock was applied to model rate variation across branches, and pertinence of a relaxed estimation was checked after verifying that the distribution of the coefficient of variation was .1. The dating analysis was performed with the total matrix, but calibration of the molecular clock was done using COI and 16S mutation rates only, as information on rates of mutation of these two fragments is widely described in multiple groups and for which there is extensive fossil calibrated divergence time data in crustaceans 16S-1472 AGATAGAAACCAACCTGG 46 [76] 12S 12sf GAAACCAGGATTAGATACCC 53 [77] 12sr TTTCCCGCGAGCGACGGGCG 53 [77] 28S 28s-rD1a CCCSCGTAATTTAAGCATATTA 52 [78,79] 28s-rD3b CCYTGAACGGTTTCACGTACT 52 [78,79] 28s-rD3a AGTACGTGAAACCGTTCAGG 52 [78,79] 28s-rD4b CCTTGGTCCGTGTTTCAAGAC 52 [78,79] 28sA GACCCGTCTTGAAGCACG 52 [78,79] 28S B TCGGAAGGAACCAGCTAC 52 [78,79] H3 H3 AF ATGGCTCGTACCAAGCAGACVGC 57 [80] H3 AR ATATCCTTRGGCATRATRGTGAC 57 [80] doi:10.1371/journal.pone.0048233.t002 Evolutionary Patterns in Cambarellinae [31,32]. As a representation of these substitution rates, we considered the range to include extreme values reported, which extends between 0.23-1.1% per million years (PMY) for 16S [33,34] and 0.7-1.3% PMY for COI [35,36,37]. These sets were introduced as uniform prior distributions, as no evidence justifies a specific distribution of rates in our data, avoiding the introduction any additional bias to the rate values assumed. Considering the geographic distribution of the genus, a geological calibration was also included as identified with the uplifting of the TMVB, which began around 12 MYA [38]. This age was set as a maximum for MRCA of the Mexican species. Additionally, fossil calibration was included in one point as the minimum age to account from the oldest fossil from the genus Procambarus [a Procambarus primaevus, 52.6-53.4 MYA, [39]]. Monophyly was not enforced for any node. Analyses were run for 20 million generations with a sampling frequency of 2000 generations. Tracer was used to determine the appropriate burn-in by monitoring run parameters by ensuring all effective sample sizes (ESS) were larger than 200 and independent runs converged. Two million generations were discarded before recording parameters and four independent runs were performed to ensure values were converging on similar estimates.

Diversification Patterns
The two main components of the subfamily occupy two regions highly contrasting in topography and biogeographic history. Thus, a second objective in this study was to describe the patterns of cladogenesis involved in the evolutionary history of Cambarellinae and to test the hypothesis that the different biogeographic histories from the two different geographic ranges of the subfamily (i.e., the Mexican and Gulf Groups), could lead to contrasting cladogenetic patterns evidenced by possible diversification shifts. Shifts in birth and death rates can leave distinctive signatures in phylogenies, resulting in departures from linearity in semi-log LTT plots [9,11]. We compared diversification rates from the reconstructed phylogeny of the entire subfamily and of the two main clades (Mexican Group vs. Gulf Group) to different null models of diversification by using the Birth-Death Likelihood method (BDL). This temporal method was used to test different hypothesis of cladogenesis rate shifts [40]. BDL uses maximum likelihood estimates of speciation rate parameters and a likelihood score per tree, and test different rate-variable models against null models of rate-constancy under the Akaike Information Criterion (AIC) [18]. To provide an indication of the diversification rates in each case, we generated a logarithm LTT plot using the LASER package version 2.2 [41]. The LTT plot was generated from the Maximum Clade Credibility tree from BEAST, after pruning the terminals not included in each clade tested using TreeEdit v1.0a10 [42] and rooting the basal age to the one observed from the dating analysis. To test for significant departures from the null hypothesis of rateconstancy, observed DAIC RC from our data was compared to those from the different rate diversification models using BDL as implemented in the LASER package version 2.2 [41]. The test statistic for diversification rate-constancy is calculated as: where AICR C is the Akaike Information Criterion score for the best fitting constant-rate diversification model, and AICRv is the AIC for the best fitting variable-rate diversification model. Thus, a positive value for DAIC RC indicates that a rate-variable model best approximates the data. We tested five different models, of which two are rateconstant and three are rate-variable: 1) the constant-rate birth model (Yule) [the Yule process; [43]] with one parameter l and m set to zero; 2) the constant-rate birth-death model with two parameters l and m (BD); 3) a pure birth rate-variable model (yule2rate) where the speciation rate l1 shifts to rate l2 at time ts, with three parameters (l1, l2, ts); density-dependent speciation models with two variants, 4) exponential (DDX) and 5) logistic (DDL). Significance of the change in AIC scores was tested by generating a distribution of scores. This was done through simulation of 9000 trees using yuleSim in LASER, for the entire Cambarellinae subfamily and each geographic group, reflecting our sampling size in each case and having the same speciation rate as under the pure-Birth model.

Phylogeny
We sequenced three mitochondrial (16S (519 bps), 12S (365 bps) and COI (1527 bps)) and two nuclear (28S 1100 bps and H3 322 bps) gene fragments resulting in 3833 characters (2411 mitochondrial and 1422 nuclear) and giving a series of substitution models (Table 3). These new data have been deposited in GenBank (Table 1). COI-like sequences were found in seven cases, identified by the occurrence of one or several stop-codons along the sequence and an unusual sequence divergence, which affected position in the tree and divergence regarding the other sequences coming from the same population. These sequences were removed from data sets and not considered for any analysis. As previously reported [15], when working with COI sequences in crayfish these sequences have to be specially checked to ensure they are mitochondrial.
The most variable fragment was 12S, followed by COI and 16S (variable sites: COI = 530/1527, 16S = 199/519 12S = 143/365; besides this, COI showed the highest proportion of parsimony informative (PI) sites: COI = 419, 16S = 121, 12S = 80) ( Table 3). As expected, nuclear fragments were the most conservative (for the The topologies recovered by mitochondrial and nuclear analyses based on ML and BI methods were similar (Figure 3), although some discrepancies can be found in some terminal taxa arrangements and between genera-outgroup relationships, principally concerning the relative positions of Cambaridae genera representatives. Both topologies show Cambarellus as a monophyletic clade (Figure 3). Within Cambarellus we found two divergent clades which correspond to the two distinct geographic ranges of the genus based on a highly supported node by ML and BI analyses (more than 95% of nodal support values). The first lineage included the species from the Mexican Group, coincident with the TMVB in México. The second lineage included the Gulf Group, containing the species distributed in USA. Only results from the combined analyses of mitochondrial and nuclear information are shown, as nuclear evidence did not have enough phylogenetic signal to distinguish relationships within each geographic group (Mexican and Gulf Groups). As shown in different studies, mitochondrial and nuclear information could resolve different portions of the phylogeny (i.e., shallow vs. deep levels of tree, [44,45] ) and that was one of the major reasons for combining these data types in this study. The hypothesis explaining this is that long-branch attraction might be more common among deeper nodes, and that slow-evolving nuclear DNA might help to resolve such issues [46,47].
Topology tests rejected the null hypothesis of an equally good explanation for all the constrained and the unconstrained topologies. The topology obtained in this study showed a significantly better Likelihood score (L = 227483.1) than the monophyletic grouping of Pandicambarus subgenus. Our phylogenetic estimate resulted in a monophyletic subgenus Cambarellus and Dirigicambarus, but Dirigicambarus was nested within the paraphyletic Pandicambarus (Fig. 3). We tested the monophyly of the Pandicambarus by forcing this alternative topology and we can reject this hypothesis by the results of SH and AU tests (likelihood values for the alternative hypothesis/p values for SH and AU -27565.1/ 0.043, 0.047). Except for the division within Pandicambarus, Fitzpatrick's notion of relationships among the subgenera is supported by our resulting topology, except for the nonmonophyletic Pandicambarus as Pandicambarus and Dirigicambarus are nested together as a sister clade that is then sister to Cambarellus as proposed by Fitzpatrick. Bayesian inference also failed to support the monophyly of Pandicambarus failing to find a monophyletic Pandicambarus in 9900 trees resulting from the MCMC search.
Species were generally well recovered as monophyletic groups for most of those included in the Gulf Group, but a different situation is depicted for the Mexican Group ( Figure 3). The clades highly supported by phylogenetic analyses have a geographic concordance, supporting the hypothesis that geographic events could have been important factors influencing cladogenesis in the genus, especially those regarding geographic features of the TMVB. Phylogenetic structuring between all Mexican taxa did not support the monophyly of some of the species currently recognized, as the highly supported clades showed representatives of multiple named species, suggesting that some of the named species did not form monophyletic assemblages.
Low 16S divergences can be observed between taxa. Divergences obtained between those contained in the Gulf Group were higher than those from the Mexican Group. The mean sequence divergence considering the likelihood model within the former was D HKY = 4.13%, and that within the latter was D HKY = 1.18% (Table 4).
The Mexican Group is composed of several clades highly supported by ML and BI analyses (95-100% support, termed with roman numerals in Figure 1), which also show geographic concordance. Some geographic overlapping between clades was observed, mainly along the Lerma Basin. The Clade I included populations from the Cuitzeo and Middle-Lerma basins, morphologically assigned to C. montezumae. C. zempoalensis from type locale was placed inside this clade as well. Cambarellus patzcuarensis from the basins of Pátzcuaro and Zirahuén were contained in Clade II and sister clade to Clade I. The third and more divergent clade (Clade III) consisted of a population from La Mintzita, geographically close to the Cuitzeo basin.
Clade IV consisted of populations from the basin of Chapala and its tributaries (Duero River), as well as its neighboring basins, Cotija and Zapotlán. This group included two species, C. chapalanus and C. prolixus, both found in Lake Chapala and associated with different habitat conditions. Also included here were populations from up-stream tributaries of the Santiago River, which originates as an outflow of the Chapala Lake. Clade V contained populations from the river Ameca basin. Clade VI contained the population from Zacapu Lagoon. The Clade VII included two populations from the eastern-limits of the distribution of the genus in the TMVB, the populations of Xochimilco (type locality for C. montezumae) from the Valley of México basin and the crater lake Quechulac. The Clade VIII was composed of two populations from the northern margin of the Middle-Lerma basin and the Clade IX by populations from the basins of the Santiago and Magdalena rivers, in the west part of TMVB.
Gulf Group relationships depict a phylogenetic structuring corresponding to geographic ranges. C. diminutus corresponds to the most divergent lineage, while two clades were recovered with high ML and BI support corresponding to a west-east pattern. The first clade contained most of the species from the Central and East Gulf Coast (CEG), except C. diminutus, and included four recognized species. Populations of C. shufeldtii from the Mississippi river basin form a monophyletic group, while C. blacki, C. lesliei, and C. schmitti are grouped together in a sister clade to the latter, geographically covering the eastern extreme distribution range of the genus in the Gulf Group from the Mobile Bay, Alabama to the Swuanee River, Florida. A similar grouping is observed in the second clade of the Gulf Group, containing populations from the West Gulf Coast (WG), mainly in the south-west part of Texas, where C. puer was recovered as a sister lineage to the clade grouping C. texanus and C. ninae.

Diversification Patterns and Dating
Log-likelihood scores with the molecular clock enforced and not enforced were 213.893 and 213.767, respectively. As the LRT rejected the null hypothesis of a global molecular clock (x2 252, P = 0.001), the sequences analyzed did not evolve at a homogenous rate along all branches and we proceeded to use a relaxed molecular clock (Fig. 4) as a result.
Ages from the dating analysis were recovered with consistency through repetitions (Figure 4). The crown age for the tree was 53 Myr (95% highest posterior density [HPD] interval for node heights/ages: 52.6-53.7 Myr), which corresponds to the separation of the genus Procambarus from the rest of the groups. We estimated an approximate age of 31.0 Myr (27.4-34.9 Myr 95% HPD) for the TMRCA of clade containing the Cambarellinae. MRCA for the terminals included in the two lineages of the Gulf Group is approximately 16.7 Myr (13.9-19.7 Myr 95% HPD). MRCA of the Mexican Group was dated around 11.1 (9.8-11.9 Myr 95% HPD). We propose some major biogeographic events inferred from the phylogenetic structure, which depicted different vicariant and dispersion events along the evolutionary history of Cambarellinae (depicted in Figure 4.).
The LTT plots track the temporal accumulation of lineages in a clade and indicate that the subfamily Cambarellinae did not significantly deviate from a constant model of diversification during its evolutionary history, as evidenced in the LTT analyses for the entire subfamily (including both, Gulf and Mexican Groups, see Fig. 5). LTTs rate-constancy models received better AIC scores, and they were not significantly different from the best rate-variable model for all analyses ( Table 5). The pure birth speciation rate model was identified as having the lowest AIC value amongst the other models tested for the subfamily together and the two groups separately. Although the Mexican Group showed the highest diversification rate (under pureBirth model r = 0.174), it is still a low value as compared to recognized shifts in diversification in other animal groups ranging from 0.4 to 0.8 speciation events per million years [48,49].
Quick inspection of the LTT plots shows some differences between the cladogenesis of the entire subfamily and that of the Gulf and Mexican Groups alone ( Figure 5). However, according to the BDL analysis, the diversification rate-constancy statistic DAICRc was found to be similar between them, being 20.

Phylogenetic Relationships
Our results are consistent with the monophyly of the Cambarellinae subfamily, previously proposed from morphology and a set of apomorphic characters [2]. The combination of mitochondrial and nuclear markers provide sufficient information to resolve the relationships between highly supported clades, namely the Gulf (Pandicambarus/Dirigicambarus) and Mexican (Cambarellus) Groups and included clades (Figure 3). Less resolution is observed at the deeper nodes of the Mexican Group, where several clades were not supported by all analyses. It is possible, as commonly argued for polytomies, that such patterns could be related to an acceleration of speciation rates in a short period of time [50]. Species sampling in this study is not complete, as three species are still to be added to the phylogenetic analysis. These correspond to C. alvarezi, C. areolatus and C. chihuahuae from North of Mexico and have almost no collection records. Populations from the aforementioned species are currently under serious threat or possibly extinct, as we did not find any specimens in our attempts to collect them. Their rarity is possibly due to extreme habitat alteration or drought, a situation reported as critical for freshwater fauna in some of the localities from where they have been recorded [51,52]. Their future inclusion, if possible (mostly through museum collections or captive populations), could provide valuable insight into the phylogenetic relationships within the subfamily, especially between the Mexican and Gulf Groups defined here.
Several differences can be found between the phylogenetic relationships emerging from this work and the previous hypothesis [2]. First, relationships between species in the Gulf Group are not congruent with several assumptions made from morphology, especially regarding the phylogenetic meaning of genitalia variation. Although species are generally well recovered as monophyletic, their relationships are not congruent. As evidenced by topology tests carried out in this study, sister relationships proposed by genital morphology between the two subgenera from the Gulf Group (Pandicambarus and Dirigicambarus) is not supported. Instead, Dirigicambarus (composed by C. shufeldtii) is recovered as a sister taxon of a clade containing C. lesliei and C. schmitti. This would leave the subgenus Pandicambarus as paraphyletic, ultimately questioning also its phylogenetic validity. Maintaining of the subgenus Dirigicambarus for C. shufeldtii could be also questioned, as no phylogenetic evidence supports it, pointing out that genital distinctiveness in this species could be the result of drift events or selective processes along its history. Besides its proposition as a member of a separate subgenus, C. shufeldtii has been recognized as a derived rather than a plesiomorphic representative [2], an assumption supported in this study. Therefore, we recommend that the subgenus Dirigicambarus be disregarded and that the genus Cambarellus should contain only two subgenera, namely Cambarellus and Pandicambarus that correspond to the Mexican and Gulf clades, respectively (resulting in Cambarellus shufeldtii being considered a member of the subgenus Pandicambarus). Our phylogenetic results support the hypothesis of C. diminutus as having plesiomorphic character states for the Gulf Group. Its unique morphological traits (outlined in [2]) are in agreement with this hypothesis.

Taxonomic Implications
Numerous species concepts have been proposed that emphasize different features for delimiting species. Sometimes, this has led to contrasting conclusions regarding species limits and the number of species in many groups. A 'unified species concept' was advocated that emphasizes the common element found in many species concepts, which is that species are separately evolving lineages [53]. This unified concept also allows the use of diverse lines of evidence to test species boundaries [e.g., monophyly at one or multiple DNA loci, morphological diagnosability, ecological distinctiveness, etc. [53,54] and is the species concept we follow in this study.
There were two cases in which the inferred topology did not recover species' monophyly in the Gulf Group. The first one shown by one individual morphologically assigned to C. shufeldtii (Locality 5, Colorado Basin), which grouped with individuals of C. Bootstrap support from ML (above) and Posterior Probabilities from Bayesian Inference (bellow) are indicated on each node. ***Stands for 95 or more, **for 85-94 and *for 75-84 support values from ML analyses. Drawings correspond to male genital morphology, which is the base for traditional taxonomy of subgenus and species in the group. Individual 5-1 was morphologically identified as C. shufeldtii, but is considered here as C. ninae based on the phylogenetic position in tree. doi:10.1371/journal.pone.0048233.g003 Evolutionary Patterns in Cambarellinae Table 4. Uncorrected (below diagonal) and ML 16S rDNA distances (above diagonal) between phylogenetic groups of Cambarellus.  Evolutionary Patterns in Cambarellinae ninae and the other by one individual morphologically assigned to C. puer (Locality 12, San Bernard Basin), grouped with individuals of C. texanus. The most plausible explanation for this could be the finding of introgression of C. shufeltii, supported by the overlapping ranges of these species in east Texas. As a common consequence, introgression between species with smaller ranges could be favored when they share similar regions with widely distributed species like C. shufeldtii and C. puer (e.g., [55]). The aforementioned hypothesis needs to be supported with faster-evolving nuclear markers, which allow the differentiation between species, and could be approached in the near future. For the Mexican Group, the phylogenetic structure shows a geographic correspondence. This observation supports the hypothesis that cladogenesis in the group has been influenced by geological history. This geographic correspondence could explain why instead of recovering species, cladogenetic structure recovered different hydrological units as monophyletic. This is the case for the widely distributed C. montezumae, which is not recovered as monophyletic, as several populations morphologically assigned to this species were located in different clades in the Mexican Group. In fact, several populations morphologically assigned to C. montezumae form a paraphyletic group, as C. zempoalensis is recovered inside this group. Another example concerns C. prolixus, included inside the wider genetic variation of C. chapalanus. However, the striking morphological distinctiveness of C. prolixus suggests a very recent processes of divergence between this species and C. chapalanus which may be missed by the genetic markers used here [see [56] for discussion on the relative importance of genetic markers versus selected morphological differences in species studies]. Based on an unified species criterion, we found support for all described species in the Mexican Group from TMVB, which match to the terminal clades in tree ( Figure 3) plus C. prolixus, which possesses contrasting morphological and ecological features. These clades correspond to six described species: 1) C. zempoalensis, corresponding to the population from Zempoala. Temporarily, we consider this species as valid but this needs to be confirmed with an analysis including the 'lermensis form' (in the terms of Villalobos' proposal) [57]. This is because when considering the range of this clade, it probably includes the aforementioned form, from the upper Lerma Basin. As such, C. montezumae lermensis would be raised to species rank and C. zempoalensis would stand as a junior synonym; all populations found in clade I would temporarily correspond to C. zempoalensis, until confirmation of the above mentioned issue regarding its synonymy with C. montezumae lermensis; 2) C. patzcuarensis, for those populations from the Patzcuaro basin; 3) C. chapalanus, from the basin of Chapala and adjacent basins; 4) C. prolixus, from certain habitat conditions at Chapala Lake; 5) C. montezumae, from the Valley of Mexico and adjacent basins and 6) C. occidentalis, from the lower part of the Río de Santiago basin, at the western extreme of the distribution in México. In addition, we found several monophyletic clades, and in congruence to the same criterion, we propose they correspond to no recognized species, those from the terminal clades in tree ( Figure 3): 1) clade III, for the population from La Mintzita spring; 2) clade V, for populations from Ameca basin; 3) clade VI, for populations from the Zacapu Lagoon and 4) clade VIII, for certain populations from the northern side of the Middle Lerma Basin (populations of La Laja basin and Vegil, see Table  S1).

Rates of Cladogenesis and Contrasting Cladogenetic Forces
Unlike the cladogenetic structure, the rate at which cladogenesis took place in Cambarellus does not seem to be affected by geologic events. Even when most of the cladogenetic events in the Mexican Group are probably the result of vicariance corresponding to geological features as the formation of the TMVB, a geologic region that has been proposed to affect cladogenesis in different freshwater groups [58,59], this study has found no effect of geologic events on speeding or reducing cladogenesis rates. Although different in nature, and affected by contrasting geographic ranges, vicariant events in both groups lead to similar cladogenetic trajectories, demonstrating the impact of climatic and geologic forces on allopatric speciation.
All these lines of geological evidence indicate that the historical geographic range of the hypothesized ancestral species of Cambarellinae in México and the Southeast of the United States have changed dramatically over time. Additionally, some other effects could have played a roll in speciation in both groups. Although the continental ice sheets during the Pleistocene glacial periods in North America never extended into the study area, these glaciations had some profound indirect effects in freshwater faunas in México and are hypothesized to have permitted dispersal by stream captures, local inland or estuarine flooding, and interconnecting drainages due to lowered sea levels during the late Neogene [60].

Biogeography
Our results support that MRCA for the Cambarellinae existed in the Eocene, ,40.4 MYA (35.2-45.7 MYA). A singular biogeographic event inferred from this study comes from the separation of the two major clades, which could be related to the Eocene-Oligocene boundary, a transition documented to strongly affect terrestrial, marine and freshwater dwellers, as evidenced by significant extinctions and taxonomic turnovers in a wide range of groups [61,62]. In this case, the formation of the Rio Grande Rift could have vicariant effects on the ancestors of both groups.
We postulate that historical vicariant events are related to change in geographical barriers and climate in both groups while dispersal events of some species are responsible for occupying the current wider range, with their current absence related to extinction periods. While these biogeographic events are present in both groups, contrasting vicariance and dispersal impacts on distributions are not unusual for freshwater crayfishes [63]. Here we explain some possible alternatives, inferred from congruence in timing of cladogenetic events. Species cladogenesis in the Gulf and Mexican Groups are best explained by an allopatric mechanism of speciation because no overlap is observed between sister taxa in Cambarellus [8]. Estimation of divergence times provides a temporal scenario of these events, allowing for a relationship of earth history with hypothesized vicariant mechanisms proposed to promote allopatric speciation. We postulate that divergence patterns in these groups are contrasting in several ways. First, date estimates agree with a more ancient diversification in the Gulf Group than in the Mexican Group, even when possible events related to species diversification from Northern Mexico could not be inferred. It is possible that the latter predate the ones observed for the Mexican Group. Second, while the Gulf Group cladogenesis could be more related to climatic oscillations, orogenic impacts could have been more important for diversification of the extant species in the Mexican Group. The absence of Cambarellus from the Río Grande basin could be explained by a generalized extinction of its populations related to the high desiccation rate since the Tertiary [64]. The ultimate evidence of that would be the presence of C. chihuahuae from the Guzman basin in the Southern part of the Río Grande Rift. This high extinction rate could explain the current disjunct geographic pattern between the Mexican and Gulf Groups. It seems reasonable to consider that as a consequence of the extinction rate along the former contact zone between the groups. It would not be surprising to find relict populations from both groups if further sampling efforts in this region could take place, which could modify their known range and find regions containing both lineages.

Proposed Vicariant Events Promoting Speciation
The first diversification event in the Gulf Group was dated to Early Miocene, ,16.7 MYA (13.9-19.7 MYA), and corresponded to the separation of the C. diminutus lineage. It is possible that extinction events could explain the observation of a unique well differentiated branch leading to C. diminutus, although a wider genetic variation not yet sampled from this lineage could be possible, which would be consistent with the wide morphological variation previously observed [2]. Orogenic activity dating to this period corresponds at the SE of United States with the formation of the Edwards Plateau, and the Miocene increased activity along the Balcones Fault. The next cladogenesis recorded for the Gulf Group is congruent with Late Miocene times, approximately 8.7 MYA (7.3-10.3 MYA), originating the Western (WG) and Eastern Gulf Coast (EG) species groups. These speciation events are consistent with sea levels along the gulf coast driven by climatic oscillations since the Middle Miocene. These were characterized by a dramatic rise in sea level between 80 and 100 m above the present day sea level [65,66]. As a consequence, a marine incursion took place along the coast of the Gulf of Mexico, which could be important in the split of West and Central-East distribution ranges and keep them separated long enough to induce a strong speciation event.
In the Late Miocene there was a sharp drop of 80-100 m below present sea levels, extending Gulf of Mexico tributaries further south. This southward extension of Gulf Coastal rivers created connections between tributaries that were isolated during periods of higher sea level. The Late Miocene drop in sea level correlates with the estimated age of the first speciation events among the extant species of Cambarellus from the Gulf Group. Later, the Pliocene (2.5-5.5 MYA) was characterized by a 50-80 m rise above current sea level, but this incursion lasted for only a short time, approximately one million years [66]. Sea levels dropped in the Late Pliocene, and during the Pleistocene there were at least three major fluctuations in sea level, none rising higher than 10-20 m above the current level [66,67]. All these events could affect the most recent speciation events and possible inter-basin connection between the Gulf species of Cambarellus could allow for dispersal of some of the today widely distributed taxa along the coastal drainages.
As a lentic-habitat dweller is the widespread way of life for the subfamily, it is reasonable to think that this could be the same situation for the ancestors of the different groups. In this case, formation of Paleolakes during the Middle Miocene (,10.8 MYA) along the Northern Central Plateau of México and South-East United States could be important features driving early cladogenesis in the Mexican Group.
The pattern of distribution observed in Cambarellus agrees with those proposed for other freshwater organisms, like the Plateau Track and western Mountain Track [64]. To explain similar distributions in fish genera such as Ictalurus, Moxostoma, and Micropterus, these patterns suggest former hydrographic exchanges across the present arid plateau. Based on faunal composition and the finding of sister taxa between those regions like Tampichthys/ Codoma and Algansea/Agosia sister pairs [68], possible connections between drainages of the South Western Gulf Slope (Nueces, Colorado and Guadalupe rivers) and those from the northern Río Grande tributaries have been suggested [69]. These connections could explain the presence of the Northern Central Plateau species (C. alvarezi, C. areolatus and C. chihuahuae), especially joined to lake habitats. Extensive lakes associated with the past Río Grande inflow have been documented to cover much of north-western Chihuahua and southern New Mexico in Pleistocene times, like the Lake Cabeza de Vaca [70].
The reduction in volume of lacustrine habitats in the Central Plateau by climatic events may have resulted in a high rate of late Cenozoic extinction [64], and the patchy distribution pattern of Cambarellus in this region. This high rate of desiccation, now increased by human activities [51], could have eroded diversity in this region, driving to extinction most of the Cambarellus populations in the Northern Central Plateau of México and could also explain the current absence of Cambarellus from the rivers south of the West Gulf Coast drainages. Partial extirpation from a formerly continuous range due to increasing dry rate during the Tertiary, has also been seen in different fish groups with a similar range, like Goodeidae and Cyprinodontidae [71]. Additionally and continuing southward, former connections between Northern and Western Central Plateau rivers could explain the presence of C. occidentalis in the Lower Santiago basin and from there a connection to the rest of the TMVB could be inferred.
Along the TMVB, the Lerma-Santiago river system is the main drainage. Previous connections between the Lerma River and northeastern and western drainages have been suggested for Goodeidae and Cyprinid fish [71]. Similar to what has been postulated for freshwater fish groups like the families Atherinidae, Goodeidae Cyprinidea, diversification in Cambarellus along the TMVB could have been related to an ancient and successive fragmentation of the Lerma-Santiago drainage across extensive lacustrine systems from the Miocene to Pleistocene [71,72].
Separation of the main clades of Cambarellus in the TMVB is dated along the late Miocene and Pliocene (10.8-4.6 MYA), a period of high geological activity in México [73]. Formation of the TMVB advanced in a West-East direction [74], and this could influence the separation of clades from the main groups of Cambarellus. This formation could have begun before the presence of Cambarellus in TMVB, given its absence on the Pacific basins south to the Zapotlán basin, at its western margin. Major diversification of the genus took place in an interval of time of less than 9 MYA.
This study found evidence consistent with a long and complex evolutionary history of the Cambarellinae. The group's distribution has been modified extensively by geologic and climatic factors. Although there appear to be contrasting causes for cladogenesis between the two groups, they have similar diversification rates. In addition, these results showed that genital and morphological changes widely used in the subfamily in particular and in crayfish in general, should be compared with other kinds of evidence in order to make more robust use of morphological differences for evolutionary inferences.

Supporting Information
Table S1 Voucher numbers and localities of the individuals analyzed. Proposed new taxa are indicated as Cambarellus sp., and correspond to the terminal clades indicated with roman numerals in phylogeny (see Figure 3). (DOCX)