Molecular phylogeny of Kazimierzus Plisko, 2006 (Clitellata, Kazimierzidae) from the Western and Northern Cape Province inferred from mitochondrial DNA sequences

Species identification of earthworms using morphology can be challenging and inconclusive as homoplasy in many characters is high. The use of molecular DNA technology, such as the use of conserved regions in mtDNA and nuclear DNA has unravelled the phylogenetic background of several earthworm species. The current study utilised the cytochrome c oxidase subunit I (COI) mitochondrial marker to reconstruct the phylogeny of Kazimierzus Plisko, 2006 species from the Western and Northern Cape provinces of South Africa. Phylogenetic reconstructions were implemented using Bayesian Inference, as well as Maximum Likelihood. Both tree building methods adhered to the monophyly of the majority of the taxa. Results showed that species fell into two clades and validated eleven currently known Kazimierzus species (K. circulatus (Plisko, 1998), K. franciscus (Pickford, 1975), K. guntheri (Pickford, 1975), K. hamerae (Plisko, 1998), K. kleinoodi Nxele & Plisko, 2017, K. nietvoorbiji Nxele & Plisko, 2017, K. nieuwoudtvillensis Nxele & Plisko, 2017, K. occidualis (Plisko, 1998), K. pearsonianus (Pickford, 1975), K. phumlani Nxele & Plisko, 2017, K. sophieae (Plisko, 2002)). Cryptic diversity is evident in K. occidualis with genetic diverAfrican Invertebrates 61(2): 83–92 (2020) doi: 10.3897/AfrInvertebr.61.53380 http://africaninvertebrates.pensoft.net Copyright Thembeka Clara Nxele et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Thembeka Clara Nxele et al. / African Invertebrates 61(2): 83–92 (2020) 84 gence greater than 12% amongst populations. Kazimierzus franciscus and K. ljungstroemi (Pickford, 1975) have a low genetic variability suggesting close relatedness or probably conspecificity. A group of specimens from Clanwilliam are morphologically identical to K. sophieae, but are genetically distinct and may belong to undescribed species. This study demonstrates the importance of integrative taxonomy in earthworms in order to present reliable taxonomic and biogeographic data.


Introduction
Earthworms constitute a large component of soil invertebrates and are regarded as soil engineers (Jouquet et al. 2006). They alter soil properties and enhance nutrient cycling (Lavelle 1988) which determines plant community composition. Regardless of their importance, African earthworm fauna have not received much attention and their evolutionary relationships are subject to debate. Species identification is possible by investigating their anatomy (Chang et al. 2007, Csuzdi 2010, Csuzdi and Zicsi 2003, Plisko and Zicsi 1991, but the structural simplicity of their body plan and existence of cryptic species may hinder taxonomic classification (Domínguez et al. 2015). Traditional morphology-based identification also requires substantial taxonomic expertise in this group, because it involves observation of minute morphological characters . As such, the use of molecular DNA technology in species delineation is an important complement to morphological classification of earthworms.
The use of DNA sequences has increased in the recent past, because it is less subjective than morphological characters, allows for the analysis of several characters (Scotland et al. 2003) and is applicable at all developmental stages (Decaëns et al. 2013;Chang and James 2011). The molecular studies of earthworms that have used the mitochondrial cytochrome c oxidase subunit I (COI) gene in integrative taxonomy have shown good results (Blakemore 2013a, b;Blakemore et al. 2010;Bantaowong et al. 2011;Huang et al. 2007;King et al. 2008;Rougerie et al. 2009;Chang et al. 2009;James et al. 2010). This marker is able to distinguish the intra-and interspecific genetic variation and groups conspecifics together, because COI sequences are variable enough to differentiate between taxa, but are less variable in conspecifics Stoeckle and Hebert 2008;Valentini et al. 2008).
According to Nxele (2012), knowledge of the taxonomic diversity of the indigenous megadriles in South Africa is incomplete. Therefore, integration of taxonomic methods is vital to improve the knowledge and understanding of the megadrile fauna. Furthermore, Plisko (2013) stressed that a molecular study on indigenous South African megadriles is essential in order to reveal the evolutionary relationships amongst them. Against this background, the phylogenetic relationships in Kazimierzus, a genus occurring in the Western and Northern Cape provinces of South Africa, was investigated.

Sampling
In order to obtain earthworms of Kazimierzus, qualitative sampling was carried out in 2011 and 2015 during the rainy season (July-September) in the Western and Northern Cape, South Africa. Besides the focus on type localities, potential sites other than the type localities were also sampled. Earthworms were collected by digging and hand sorting. Collected specimens were anaesthetised in 20% ethanol solution, fixed in 4% formalin solution and preserved in 75% ethanol for taxonomic purposes. A subsample was preserved in absolute ethanol for molecular analysis. All specimens were examined using a Wild Heerbrugg stereomicroscope and were identified according to the descriptions in Plisko (1996Plisko ( , 1998Plisko ( , 2002, Pickford (1975), Michaelsen (1913) and Nxele et al. (2017). All specimens are at the KwaZulu-Natal Museum, Pietermaritzburg, South Africa.

Genomic DNA extraction, amplification and sequencing
Tissue from the posterior section of the earthworm was used. All DNA extractions were performed using the ZR Genomic DNA TM Tissue MicroPrep kit, following the manufacturer's standard protocol. The concentration of DNA in each sample was estimated using the NanoDrop 2000 (Thermo Scientific). A fragment of the mitochondrial cytochrome c oxidase subunit I (COI) gene was amplified using LCO1490 (5' GGT-CAACAAATCATAAAGATATTGG 3') and HCO2198 (5' TAAACTTCAGGGT-GACCAAAAAATCA 3') primers (Folmer et al. 1994). Polymerase chain reactions (PCR) were performed in a final volume of 25 µl using the BIO RAD T100 Thermal Cycler and contained: 2 µl of DNA template (approximately 35 ng/ul), 12.5 µl One Taq Quick-Load 2X Master Mix with standard buffer, 0.5 µl of 10 uM forward and reverse primers and sterile water. The thermocycler conditions were as follows: 95 °C for 2.30 min for initial denaturation, followed by 35 cycles at 95 °C for 30 sec denaturation, 50 °C to 52 °C for 45 sec annealing and 72 °C for 75 sec extension. A final extension step at 72 °C for 10 min completed the reactions.
Sequencing of the 675 bp fragment of the COI mtDNA was conducted at Inqaba Biotechnical Industries (Pty) Ltd.

Sequence alignment and phylogenetic analysis
Identity of sequences was verified by the Basic Local Alignment Search Tool (BLAST) in the National Centre for Biotechnology (NCBI). Sequences of Amynthas minimus (Horst, 1893) and Amynthas corticis (Kinberg, 1867) were included as outgroup taxa. The sequences for outgroup taxa were obtained from GenBank (Accession nos: AB542509.1, AB542469.1) and current sequences are added as supplementary data.
All specimens are at the KwaZulu-Natal Museum. The sequences were aligned using CLUSTAL X 2.1 (Larkin et al. 2007). These alignments were then manually edited using BIOEDIT 3.3.19.0 (Hall 1998). Unreliable nucleotides (low signal strength), as well as primers sequences, were trimmed off at both the 5' and 3' ends. The programme JMODELTEST v.0.1.1 (Darriba et al. 2012) was used to select the best-fit evolutionary model using the AKAIKE Information Criterion (AIC; Akaike 1973). Phylogenetic analyses were based on two approaches, Bayesian Inference (BI) was performed using MRBAYES 3.2 (Huelsenbeck and Ronquist 2001) and Maximum Likelihood (ML) analysis was performed using GARLI (Zwickl 2011). In each case, the best-fit evolutionary model selected by JMODELTEST was specified.
Clade support was evaluated by 1000 bootstrap replicates for the ML analysis and posterior probability values for the Bayesian analysis. For Bayesian analyses, all MR-BAYES analyses were run for 5000000 generations with a sampling frequency of 1000. The deviation of split frequencies was less than 0.01 at the conclusion of all analyses which confirmed that the MCMC chains had converged. The programme TRACER v1.5 (Drummond and Rambaut 2007) was used to check that the Effective Sampling Size > 200 and that posterior distribution for all parameters was unimodal. Consensus trees were generated using PHYLIP 3.69 (Felsenstein 2005) and viewed in FIG  TREE v1.3.1 (Rambaut 2009). Uncorrected p genetic distances were obtained for the sequenced specimens using MEGA 6 (Tamura et al. 2013). Each species is represented by one specimen, except where the species appeared in more than one clade on the phylogenetic tree. Sequences have been deposited in GenBank (Table 1).

Phylogenetic analyses
The sequences were 675 bp. Variable sites were 431 bp and conserved sites 233 bp showing great differentiation amongst taxa.
Most species pair comparisons showed a genetic distance above 13%, except for K. ljungstroemi and K. franciscus which have one percent genetic distance between them ( Table 2).
The MI and BI trees were congruent, therefore, support values were annotated on to the branches of the most likely trees generated for each of the datasets analysed (ML run with no bootstrap, rooted using outgroup species (Amynthas minimus and Amynthas corticis).
Bootstrap values and posterior probabilities below 50% and 0.50, respectively, were not shown on the trees. Two clades, A and B, are distinct; clade A separates further to clades C and D whilst clade B separates to clades E and F (Fig. 1). Both tree building methods support monophyly of the majority of taxa (Fig. 1). Kazimierzus sophieae is paraphyletic, found in two distinct clades; one shared with K. sp and the second clade  with Kazimierzus hamerae (clade D). The two species in clade C, K. circulatus and K. nieuwoudtvillensis, are sister taxa. Kazimierzus phumlani and K. occidualis are sister taxa comprising clade E. Clade F has six species including a strongly-supported sub-clade, composed of K. franciscus and K. ljungstroemi. The other species (K. guntheri, pearsonianus, kleinoodi and nietvoorbiji) form the other sub-clade.

Discussion
Morphological examination revealed eleven currently known species (circulatus, franciscus, guntheri, hamerae, kleinoodi, nietvoorbiji, nieuwoudtvillensis, occidualis, pearsonianus, phumlani and sophieae). However, the current phylogenetic analysis resulted in additional lineages that suggest cryptic diversity (Fig. 1) and a new species, K. sp. Most lineages had strong support but weak branch support was noted in deeper nodes with some clades having bootstrap values less than 50%, hence not annotated on the tree. Genetic distances also support the relationships observed in the phylogenetic tree. Genetic distances between species were all greater than 13% and generally much greater than that, consistent with the recommendations of Hebert et al. (2003), Decaëns et al. (2013) and Huang et al. (2007) for separation of species by genetic distance using the CO1 barcode (Table 2).
Kazimierzus hamerae is similar in appearance to K. sophieae, but molecular data confirmed that they are separate species. The phylogenetic tree (Fig. 1), however, highlights, unexpectedly, that K. sophieae specimens that were collected from two localities, one in Clanwilliam and the other in van Rhyn's Pass, fall into two clades. The two clades may be because K. sophieae is polyphyletic or there is no morphological divergence of these taxa, but are different species genetically (25.5% divergence). The type locality of K. sophieae is Nieuwoudtville, which is close to van Rhyn's Pass. It is likely that the sp. may possibly be an undescribed species; however, a comparison of sequences of K. sophieae specimens and those from the type locality, Nieuwoudtville is necessary. It seems that morphological evolution is not rapid in this complex of taxa and there is presence of cryptic diversity. Clade D has poor ML support (52%) but BI support is high (0.99). The low ML support may suggest that the present data are not sufficient to resolve this polytomy fully. Cryptic diversity was observed in K. occidualis. Finding cryptic diversity is common in earthworms, Novo et al. (2010) reporting five cryptic species within the complex of Hormogaster elisae.  reported cryptic species in Lumbricus terrestris and Pérez-Losada et al. (2009) found cryptic diversity within Aporrectodea caliginosa species complex. The observations support the view that diagnosis, based on morphology, only underestimates taxonomic diversity.
The specimens of Kazimierzus franciscus and Kazimierzus ljungstroemi were collected in two neighbouring forests, Duiwelbos and Koloniesbos in Marloth Nature Reserve, Swellendam. Duiwelbos is the type locality of K. franciscus, whilst the type locality for K. ljungstroemi is Great Winterhoek, Tulbagh District. Therefore, it is possible that all specimens belong to K. franciscus, but it would be of benefit to collect specimens from Winterhoek and compare them with the ones collected from Marloth Nature Reserve.
Although analysis of other conserved genomic regions in both mtDNA and nuclear DNA in the future would benefit the study of Kazimierzus species, the phylogenetic analysis of COI recovered several well-supported phylogenetic relationships, some of which are congruent with existing classification.