A novel genus and cryptic species harboured within the monotypic freshwater crayfish genus Tenuibranchiurus Riek, 1951 (Decapoda: Parastacidae)

Degree of molecular divergence

Preliminary analyses of both individual and combined gene trees showed a prominent separation between Qld and NSW populations. In light of this, genetic distances between Qld and NSW populations, distances between these two groups and the outgroups, and distances between the outgroups were calculated using both COI and 16S data to compare the degree of separation. These distances were calculated in MEGA5 (Tamura et al., 2011) using the net between group mean distances with 1,000 bootstrap replicates (gamma distribution with shape parameter = 1, Maximum Composite Likelihood (MCL) model; positions containing gaps and missing data were eliminated).

Combined gene tree

Combined gene trees were inferred using both Maximum Likelihood (ML) and Bayesian analyses. Specimens were included in the data set if they were sequenced for at least four of the five genes (see Table S1). The program RAxML v. 7.4.4 through the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2010) was used to infer the ML tree, and MrBayes v. 3.2.0 (Ronquist et al., 2012) for the Bayesian tree. Within the ML analysis, each gene was entered as a separate DNA-partition, the GTR+CAT model used, and bootstrapping automatically halted with the majority rule criterion. For the Bayesian analysis, each gene was entered as a separate partition with the Nst set as 6 for all genes except 16S (where Nst = 2), and Rates set as InvGamma for COI and GAPDH, and Gamma for 16S, H3, and AK, as determined by jModeltest v. 0.0.1 (Posada, 2008). Additional parameters were set as follows; two replicate Markov chain Monte Carlo analyses with four chains in each analysis (one cold, three heated), the statefreq, revmat, shape, and pinvar all unlinked, the ratepr set as variable, and the analysis set to stop when the standard deviations of the partition frequencies <0.0099 (all effective sample size values >100, PSRF+ ≈ 1.000, and the final Ngen was 1,715,000). The same analysis was performed at least twice to verify topological convergence and homogeneity of posterior clade probabilities between runs. The first 25% of samples were discarded as burnin, with the resulting trees visualised using the program Figtree v. 1.4.0 (Rambaut, 2012).

Intra- versus inter-cluster variation

An analysis of molecular variance (AMOVA) was used to calculate variation within and among clusters of sequences, as implemented in Arlequin v. 3.1 (Excoffier, Laval & Schneider, 2005). To determine what the most likely lineages were, the clades identified by the combined gene tree analysis, as well as additional splits evident within the tree that were deemed to plausibly represent lineages, were also tested, as well as groups based on the geographic division of populations (i.e., collection locality; see Table S2). The AMOVA calculates three statistics; ϕ_ST, ϕ_SC, and ϕ_CT, all of which are based on both the haplotype frequency and genetic divergence. ϕ_ST measures variation among all populations, and ϕ_SC measures variation among populations within groups, and ϕ_CT estimates variation among groups. It has been suggested that an F_CT value >0.95 can represent evidence for accurate species groupings (i.e., >95% of the genetic variation can be attributed to differences among groups) (Monaghan et al., 2005). Using the ϕ_CT estimate as a surrogate for F_CT (as this estimate includes genetic divergence as well as haplotype frequency), this can provide an approach to delineate taxa based on population genetic analyses by interpreting the AMOVA results used to calculate intra- versus inter-cluster variation in a way analogous to F-statistics (Wright, 1978). The criterion to determine the appropriate number of lineages using this method is where an increase in the number of suggested lineages does not appreciably increase the ϕ_CT estimate for those lineages.

Barcoding gap

The genetic distances between the hypothesised lineages and between specimens for both COI and 16S were calculated and visualised to determine whether a barcoding gap existed. As the intent of this method was to provide support for, or refutation of, the lineage hypotheses formed through the species discovery approaches, lineages were pre-defined based on those results and genetic distances categorised as representing either intra- or inter-lineage distances. For the purposes of this study, a barcoding gap was defined as a clear separation (or ‘gap’) between the highest intra-lineage and lowest inter-lineage genetic distances measured between the suggested lineages. Although a standard threshold has been suggested by Hebert et al. (2004) for recognising distinct species (10× average intraspecific difference), this approach was not followed as it has been shown that there are vastly different rates of divergence for both different taxa and different genetic markers (Avise, 2009). Rather, a recognisable distinction between the inter- and intra-lineage distances was considered potential evidence for distinct species. Analyses were undertaken for Qld and NSW specimens separately.

Relative divergences between genetic groups were calculated in MEGA5. To determine inter-lineage divergence, the number of base substitutions per site was estimated from the net average between groups of sequences and the diversity between specimens was determined by calculating the number of base substitutions per site between each pair of sequences, both using a MCL model with 1,000 bootstrap replicates. The rate variation among sites was modelled with a gamma distribution with a shape parameter of 1, with positions containing gaps and missing data eliminated. This was performed separately for both COI and 16S, with all unique haplotypes included.

K∕θ Method

The species discovery hypotheses were also tested using the K∕θ method (Birky et al., 2010; Birky & Barraclough, 2009; Birky et al., 2005). Although this method was originally developed for asexually-reproducing organisms and termed the 4X rule, it has been further developed and shown to be effective for the mtDNA region for sexually-reproducing organisms (Birky, 2013). This method provides a simple way of defining species groups based on specimens/populations that form clusters (i.e., clades) that are separated by genetic gaps too deep to be ascribed to random genetic drift within a species and, therefore, must be due to diversifying selection or long-term physical isolation (Apte, Smith & Wallis, 2007).

Using the groups from the species discovery hypotheses, sister clades were identified and statistical support for these was tested. Sequence divergences were estimated within (d) and between each sister clade using uncorrected p-distances calculated in MEGA5. Nucleotide diversity (π) was then calculated using π = dn∕(n − 1), where n is the number of samples per clade. Theta (θ) was then estimated as θ = 2Neμ (where Ne is the effective populations size and μ is mutation rate per base pair per generation) by calculating π∕(1 − 4π∕3) within each clade. If d = 0 (as it did for one clade in this study), then π can alternatively be calculated as 2∕Ln(n − 1), where L is the length of the sequence. K was then calculated for each sister-clade comparison (using MEGA5) as the uncorrected net between group mean distance, with this divided by the highest θ in the comparison (as this is the more conservative approach) to provide K∕θ. Where sister clades were poorly defined in the tree, K was estimated between all potential sister clades in the polytomy, with the clade of the lowest K considered to be the sister clade. Finally, following the method of Birky (2013), if the K∕θ value was greater than 4, then the sister clades were accepted as different lineages.

Combined gene tree

Although not all groupings were statistically supported, both the ML and Bayesian combined gene trees suggested the presence of multiple groups within Qld and NSW, and displayed the same topologies (Fig. 2). Six clades were evident within the Qld populations, with the monophyly of all but two highly supported (as these were represented by single specimens). The first clade included Maryborough and some Tuan State Forest specimens (Lineage 1; BS 90%, Pp 1), and the second contained the remaining Tuan State Forest specimens as well as Bribie Island, Type Locality, and some Beerburrum specimens (Lineage 2; BS 96%, Pp 1). The two groups for which monophyly could not be established were represented by the remaining Beerburrum specimens (Lineage 3) and Hervey Bay (Lineage 4). The final two clades consisted of Tewantin and Lake Weyba specimens (Lineage 5; BS 100%, Pp 1) and Gold Coast specimens (Lineage 6; BS 100%, Pp 1). There was also some geographic structuring evident within each of the clades.

The two monophyletic clades evident within the NSW populations were strongly supported, and form Lineage 7 (Lennox Head) and Lineage 8 (Lake Hiawatha, Broadwater National Park 1 & 2) (Fig. 2). Although there was some structuring evident within Lineage 8, the branching patterns were very shallow and were therefore not explored as potential distinct lineages.

Intra- versus inter-cluster variation

A total of eight lineage arrangements was deemed plausible based on apparent genetic groupings and collection localities, and were tested using AMOVAs (Table 4). The process of assigning the potential lineages is outlined in Table S2, where a hierarchical approach was taken to split the tree into major genetic groups, minor genetic groups, and geographic localities. As there was no logical reason for combining the NSW lineages for the AMOVA analysis based on either the phylogenetic or geographic information, the NSW populations were considered to consist of the LH lineage and the LakeH/BNP lineage. Further testing, however, was considered appropriate to determine the lineages present within Qld. Figure 3 shows an increase in the ϕ_CT estimate, with a plateau reached at six lineages for both the COI and 16S estimates. These six Qld lineages represent the most parsimonious arrangement of the specimens into lineages.

Species discovery hypothesis

As the combined gene tree was inferred using only specimens that were successfully sequenced for at least four of the five genes, not all collection localities were represented on the tree (i.e., TSFN, KNP, Moo, Eu). Of these localities, only TSFN was represented in the AMOVA analysis, as the remaining localities were represented by a single sequence and therefore could not be included in the AMOVA. In order to assign these populations to a lineage for further testing, the individual gene trees and haplotype networks were examined and the localities were designated through the closest phylogenetic connection (data not shown). Both of the species discovery approaches suggested the presence of eight lineages (six in Qld and two in NSW; Table 5), and formed the lineages to be validated.

Barcoding gap

The COI data showed some overlap of the intra- and inter-lineage estimates within Qld, resulting in no usable barcoding gap for lineage separation (Fig. 4A). Where the overlap occurred, the low inter-lineage estimates were attributable to the Lineage 1 vs. Lineage 2 comparison, and the high intra-lineage estimates were seen between specimens within Lineage 1. However, many estimates between these two lineages fell in the higher range of the inter-lineage estimates as well as the low range.

The 16S data for Qld populations showed a clearer relationship between lineages (Fig. 4C). Although there was a very small overlap between the intra- and inter-lineage distances (occurring between two specimens from Lineage 1), this represented an overlap of less than 0.01%. When the existence of this overlap was disregarded, there was a small gap at 2.8–3.0%. However, despite there not being a distinguishable gap due to the overlap, identification of the majority of lineages through the comparison of intra- and inter-lineage distances was clear and distinguishable.

When the estimates within and between Lineage 1 and 2 specimens were removed from both the COI and 16S data (with the comparison between these two lineages and all other lineages remaining), a clear barcoding gap was seen in both data sets (Fig. 4B, Fig. 4D). For COI, the gap occurred between 1.7–4.7%, and between 0.9–3.5% for 16S. This shows that all other Qld groups (i.e., Lineage 3 through 6) represent clear lineages based on the barcoding approach using both COI and 16S data.

For NSW populations, there was a clear barcoding gap between the two lineages (i.e., Lineage 7 and 8), occurring between 1.5 and 6.6% for the COI data and 0.7–3.0% for the 16S data (Fig. 5).

K∕θ method

The sister clades within Qld and NSW were tested using the K∕θ method for a delimitation of eight lineages (six from Qld, two from NSW) using both COI and 16S data (Table 6). In some instances, sister clades that were defined by the lowest K-distance (as they were ambiguous based on the combined gene tree) differed between the COI and 16S datasets. In these cases, only the relevant K∕θ comparison for the applicable gene was calculated.

Lineage assignment

Although there was some ambiguity in the barcoding analysis of the Qld COI data regarding the separation of Lineage 1 and 2, the 16S data showed support for the species discovery lineage hypothesis. Because of the deeper phylogenetic inferences provided by 16S in addition to the fact that there were many genetic distances within and between Lineage 1 and 2 falling within the expected distributions, the lineage hypothesis for Qld populations was considered supported by this analysis (Table 7). The two NSW lineages were clearly separate based on both the COI and 16S data and were therefore also supported (Table 7). In the K∕θ analysis, all lineages were supported by both genes with the exception of the split between Lineage 1 and 2 (both genes), and Lineage 1 and 3 (16S) (Table 7).