Molecular approaches uncover cryptic diversity in intertidal Ligia isopods (Crustacea, Isopoda, Ligiidae) across the southern Africa coastline

Field sampling, preservation, and identification

We hand-collected Ligia individuals from 18 localities around the coastline of southern Africa between 2014–2017. Detailed locality information is provided in Table 1. All samples were field-preserved and stored in 70% ethanol until molecular analyses were carried out. In the laboratory, specimens were identified to species by visual inspection of key characteristics (e.g., appendix masculina of the second pleopod of males) and comparing these traits to those reported for the various recognized Ligia species in southern Africa (Barnard, 1932; Ferrara & Taiti, 1979). Field collections were carried out under Scientific Collection Permit RES2017/53 issued by the South African Department of Environmental Affairs.

Molecular laboratory methods

We extracted total genomic DNA from several pleopods for 2–10 Ligia individuals per location using the Quick g-DNA MiniPrep Kit (Zymo Research, Irvine, CA, USA), following standard protocol instructions. For each individual, we PCR-amplified a 658-bp fragment of the Cytochrome Oxidase I (COI) mitochondrial gene using the LCO-1490 and HCO-2198 primers and previously published conditions (Folmer et al., 1994). We also PCR-amplified a 661-bp region of the sodium-potassium ATPase alpha subunit (NaK) gene using the NaKFb and NaKR2 primers and standard conditions (Tsang et al., 2008). Positive PCR amplifications were determined by visualizing PCR products on 1% agarose gels stained using SYBR Safe (Invitrogen, Carlsbad, CA, USA). Positive amplicons were sequenced at the University of Arizona Genetics Core, with sequences and assembled and edited (i.e., primer removal) using Geneious R8.0.5.

Sequence alignments, phylogenetic analyses, and estimation of molecular divergence

The mitochondrial COI and nuclear NaK sequence datasets were aligned independently using the MAFFT server (Katoh & Standley, 2013) under standard settings for nucleotide sequences. Visual inspection of the resulting alignment produced no evidence suggestive of pseudogenes (e.g., stop codons, high rates of amino acid substitutions) or indels. Due to the limited phylogenetic signal within the NaK dataset, we did not concatenate the two datasets and carried out phylogenetic searches only on the COI resulting alignment. Relationships within the NaK dataset were estimated using haplotype network reconstructions.

We carried out preliminary phylogenetic analyses incorporating sequences produced in this study, as well as all publicly available COI sequences for other Ligia species and two Ligidium species (Ligidium germanicum and Ligidium hypnorum; COI accession numbers DQ182795 and DQ182781, respectively) so as to identify the most appropriate outgroup for our phylogenetic analyses. These preliminary analyses recovered the monophyly of southern Africa Ligia species; however, relationships amongst the ingroup were poorly resolved. As a result, we carried out our phylogenetic excluding outgroups and opting instead for rooting using a midpoint-root approach.

Prior to phylogenetic searchers, we determined the most appropriate model of nucleotide evolution for our COI dataset using the Modeltest script (Posada & Crandall, 2001) as implemented by FindModel (http://hiv.lanl.gov/content/sequence/findmodel/findmodel.html). Model selection was made by comparing the likelihood scores of all 28 models available based on an initial Weighbor tree. We then carried out phylogenetic searches under both Maximum Likelihood and Bayesian inference approaches. Maximum Likelihood phylogenetic searches were carried out in RAxML v8.1.2 (Stamatakis, 2014; Stamatakis, Hoover & Rougemont, 2008) and consisted of 1,000 thorough bootstrap replicates, followed by a thorough ML search under the GTR +Γ model. We produced a majority-rule consensus tree of all bootstrap replicates using the Sumtrees command of DendroPy v4.1.0 (Sukumaran & Holder, 2010).

We carried out Bayesian phylogenetic searches in MrBayes v3.2.5 (Ronquist & Huelsenbeck, 2003) and Phycas v2.2.0 (Lewis, Holder & Swofford, 2015). Searches in MrBayes consisted of two simultaneous searches of four chains, each sampled every 5,000th tree, while Phycas searches consisted of a single search, sampled every 50th tree. All Bayesian searches were carried out under the GTR +Γ model. For each Bayesian analysis, we estimated node support values by discarding all samples prior to stationarity (10–25% of sampled trees) and calculating a majority-rule consensus tree using the Sumtrees command of DendroPy v4.1.0 (Sukumaran & Holder, 2010).

Lastly, we used MEGA v7.0.7 (Kumar, Stecher & Tamura, 2016) to estimate COI Kimura 2-Parameter distances (K2P) within and amongst sampled localities and major lineages observed in the above phylogenetic reconstructions.

Haplotype network reconstructions

We used the ancestral parsimony algorithm proposed by Templeton, Crandall & Sing (1992) as implemented in PopART v1.7 (Leigh & Bryant, 2015) to visualize relationships between all COI haplotypes recovered in this study. We estimated branch connections using the TCS network option (Clement, Posada & Crandall, 2000) of PopArt with networks considered separate if connections between them exceeded 33 steps (i.e., a 95% connection limit). We repeated this approach to visualize the relationships amongst NaK alleles.

Molecular species delimitation analyses

We implemented molecular species delimitation analyses (hereafter MSDAs) using ABGD (Puillandre et al., 2012). This approach identifies putative species by producing pairwise differences amongst haplotypic data and identifying a putative inter-species gap in the distance distribution. Individual haplotypes are assigned to groups if their genetic distances are less than the proposed gap. Analyses were carried out on the COI data for each Ligia species separately using the ABGD server (http://wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html) under standard settings, with the exception of the relative gap width (set to 1.0) and the distance setting (“Kimura (K80) TS/TV” = 2.0).

Phylogenetic results

We observed a basal split between two well supported clusters of highly divergent clades in all analyses: a “Western” cluster (reds and greens in all figures; Bootstrap support (BS): 100; Maximum Posterior Probability (MPP): 100%; Fig. 2) with a geographic distribution from Namibia to the Cape Agulhas region, and an “Eastern” cluster (blues, yellows, and purples in all figures; Bootstrap support (BS): 100; Maximum Posterior Probability (MPP): 100%) that was distributed from Knysna, on the south coast of South Africa (hereafter SA), to the KwaZulu-Natal region of SA. Each of these clusters was composed of two or more highly divergent clades (clades A–F; amongst clade COI K2P divergences 3.1–17.2%, Table 2).

COI haplotypes assigned to the “Western” cluster were further divided into two highly divergent clades (amongst clade COI K2P divergences: 8.5–10.7%, Table 2). Clade A (reds in all figures; BS: 93%; MPP: 100%) included all Ligia individuals sampled in Namibia (A1), as well as from two locations in SA: Jacob’s Bay (A2) and Ganzekraal (A3) and corresponds to the species morphologically identified as L. glabrata. Within this clade, we observed three lineages that correspond with the sampled localities and that were moderately divergent from each other (COI K2P: 5.1–5.6%; Table 3). The relationships between these lineages were not well supported; however, our analyses suggest a sister-taxon relationship between the lineage found in Ligia from Luderitz, Namibia (A1) and that found in Jacob’s Bay (A2) (BS, MPP < 60). The second clade part of the “Western” cluster, Clade B (greens in all figures; BS: 100; MPP: 100), comprised all Ligia individuals collected from localities between the Cape of Good Hope and Cape Agulhas (B1–B5) and morphologically corresponds to the species L. dilatata. Clade B, contrary to the Clade A, does not appear to be composed of any further divergent lineages and within-clade divergences within it were low (COI K2P: 0.0–1.2%; Table 2).

The “Eastern” cluster, which contained all Ligia collected from Knysna to the KwaZulu-Natal region of South Africa, was composed of three highly divergent and well supported monophyletic clades (C–E; COI K2P 3.1–12.2%; Table 2) which morphologically correspond to the established species L. natalensis. Within this cluster, clades D, E, and F (blues and purples in all figures) are placed in a well-supported clade (BS: 82; MPP: 100) with Clade C (yellows in all figures) in turn sister to this group. Relationships between D, E, and F are not well resolved. Clade C containing all Ligia individuals collected in the Port Elizabeth area (C1–2), was highly supported across analyses (BS: 86; MPP: 98–100), and exhibited low within-clade divergences (COI K2P 0.0–1.1%, Table 2). This clade was highly divergent from all other clades in the “Eastern” cluster (COI K2P 9.4–12.2%; Table 2) and appears genetically distinct enough to be considered a separate and previously unrecognized species within the natalensis group. Clade D (BS: 86; MPP: 99–100) includes most COI haplotypes obtained from Ligia individuals collected in Knysna (D1) and the Port Edward area (D2–4). Within clade divergence for Clade D was low (COI K2P 0.0–1.9%, Table 2). Clade E contained all COI haplotypes recovered from individuals from the Kenton-on-Sea area (E1–2), those from Kidd’s Beach (E3), as well as one each from Knysna (D1) and Salmon Bay (D2). Although this clade was not strongly supported by any phylogenetic analyses (BS, MPP < 50), we denote it as a separate clade given the very low levels of divergence between all haplotypes in it (COI K2P average ∼1.3%; Table 2), moderate amongst clade divergence when compared to haplotypes from clades D and F (COI K2P 3.5%–6.4%; Table 2), and the results of haplotype network reconstructions (Fig. 3). Lastly, the well supported Clade F (cyan in all figures, BS: 92; MPP: 99–100) contained all individuals collected at the East London Harbor (F1) and exhibited low levels of within clade divergence (COI K2P 0.0–0.6%; Table 2).

COI haplotype network reconstructions

The results of our COI haplotype network reconstructions (Fig. 3) largely match patterns produced by our phylogenetic analyses, as we recovered four separate networks (i.e., connections of <95%) largely corresponding to clades observed in phylogenetic reconstructions.

Network I (Fig. 3A) contained four haplotypes recovered from Ligia individuals from Luderitz (A1), Jacob’s Bay (A2) and Ganzekraal (A3). In Luderitz, we recovered a single haplotype that diverged by 32–33 steps from haplotypes recovered from Ganzekraal, which in turn diverged by 33–34 steps from the single haplotype recovered in Jacob’s Bay. This network closely parallels the patterns observed for Clade A in our phylogenetic reconstructions and contains all individuals morphologically identified as L. glabrata.

Network II (Fig. 3B) contained all 15 haplotypes recovered from Ligia collected between the Cape of Good Hope and Cape Agulhas (B1–B5), closely matches Clade B, and include all individuals identified as L. dilatata. Divergences in this network were low, with most connections between haplotypes being only 1–2 steps and the maximum connection between haplotypes being 10 steps. Despite such short connections, the network suggests some isolation between localities, as no sharing of haplotypes is apparent. Furthermore, haplotypes recovered within a single location appear to be much more similar (1–2 steps) than to those found at other locations in the region (∼4 steps).

Network III (Fig. 3C) consisted of six haplotypes recovered from Ligia collected in localities near Port Elizabeth (C1–2) and corresponds with Clade C from our phylogenetic findings. As observed in Network II, connections between haplotypes are very short, as most haplotypes are connected by <3 steps and the maximum span between haplotypes is seven steps. Morphologically, all members of this network were identified as L.  natalensis.

Lastly, Network IV (Fig. 3D) contained 34 haplotypes divided into three sub-networks separated by <18 steps. These sub-networks appear to correspond with clades D–F from our phylogenetic results and include all other L. natalensis individuals. One sub-network (blues in Fig. 3) contained all but two haplotypes from localities around Knysna (D1) and Port Edward (D2–4). Another sub-network (purples in Fig. 3) contained all the haplotypes collected in the localities of Kenton-on-Sea (E1–2) and Kidd’s Beach (E3). Intermediate to these two subnetworks is a small subnetwork of four haplotypes recovered from individuals collected in East London (F1; cyan in Fig. 3). In general, haplotypes collected from the same locality are much more similar to each other (<6 steps) than those from other localities (>10 steps), with two exceptions. A COI haplotype recovered from a Ligia individual collected in Salmon Bay (D2) was much more similar to those found in Kidd’s Beach (E3; 4–9 steps) than others from its own location (>29 steps). This haplotype was not observed in any other Ligia individual from any other locality. The other exception was a COI haplotype collected from an individual collected in Knysna (D1) that was shared with individuals from the Kenton-on-Sea area (E1–2). These patterns are congruent with the amongst-locality divergences where these lineages were found (Tables 4 and 5).

NaK haplotype network reconstructions

NaK haplotype network reconstructions (Fig. 4) were congruent with the above results; however, they produced much simpler patterns. We uncovered four NaK alleles separated by 1–10 steps: one allele shared by all surveyed individuals within Clade A (L. glabrata), one shared by all individuals within Clade B (L. dilatata), and two alleles from individuals from within the “Eastern” Cluster (L. natalensis). These latter two alleles were much more similar to each other (1 step) than to the other two recovered alleles (7–9 steps). The alleles found in the other clades were also highly divergent from those found in other clades, with the Clade A allele separated by 5–10 steps from other alleles and the one found in Clade B being separated by 5–9 steps. These patterns are largely concordant with those produced by mitochondrial analyses when considering the low levels of variation of this marker.

Molecular species delimitation analyses

MSDAs suggested the presence of cryptic species for all three nominal Ligia species in the region. ABGD identified three putative species for L. glabrata individuals regardless of the prior maximal distance (i.e., barcode gap) assumed, with groups corresponding with the locations in which these individuals were found (Group 1: A1, Group 2: A2, Group 3: A3). For L. dilatata two groups were identified (Group 1: B1, Group 2: B2–5); however, this split was based on a prior maximal distance of only 0.001668. Another two putative species were identified for L. natalensis individuals, with those from C1 and C2 placed in Group 1, and those from D1–4, E1–3, and F1 placed in Group 2. These results were based on prior maximal distances of 0.02154 and above. Exclusion of the Port Elizabeth individuals (C1–2) from ABGD analyses did not produce any additional groupings within the L. natalensis individuals.