In quest of contact: phylogeography of helmeted terrapins (Pelomedusa galeata, P. subrufa sensu stricto)

Sampling and laboratory procedures

Fieldwork and sampling in South Africa were permitted by the Limpopo Provincial Government (permit ZA/Lp/80202), Ezemvelo KwaZulu-Natal Wildlife (permit OP 139/2017), CapeNature (permit AAA007-00212-0056), the Department of Environmental Affairs, Eastern Cape (permit CRO117/13C & CRO 118/13CR), and Biodiversity Northern Cape Province (permit 245/2015). Fieldwork and sampling in Namibia were permitted by the Ministry of Environment and Tourism (permit 1910/2014). Terrapins were hand-collected or captured using baited traps, and blood and saliva samples were taken as approved by the Ethics Committee of the University of the Western Cape under ethical clearance number ScRiRC2008/39. Terrapins were released after sampling at the capture sites.

Using alcohol-preserved blood or saliva samples and wet laboratory approaches for fresh material as described in Fritz et al. (2014), we generated sequences of the partial 12S rRNA gene (360 bp) for 36 Pelomedusa galeata and seven P. subrufa sensu stricto from South Africa and Namibia. Another sequenced mtDNA fragment of these samples comprised the partial ND4 gene plus adjacent DNA coding for tRNAs (816 bp), and a third mtDNA fragment corresponded to the partial cytochrome b (cyt b) gene (674 bp). In addition, we used homologous mtDNA sequences of helmeted terrapins from other investigations (Vargas-Ramírez et al., 2010; Wong, Fong & Papenfuss, 2010; Fritz et al., 2011; Fritz et al., 2014; Fritz et al., 2015; Petzold et al., 2014; Nagy et al., 2015). Including previously published material, we studied mtDNA sequences of 116 P. galeata and 41 P. subrufa sensu stricto and georeferenced their collection sites (Table S1).

In addition, we sequenced two nuclear loci of samples representing all mitochondrial lineages and subclades and almost all sampling sites (Table S1). One locus, the intron 1 of the RNA fingerprint protein 35 gene (R35) has been previously shown to be species-diagnostic for P. galeata and P. subrufa (Vargas-Ramírez et al., 2010). The other locus, including coding and non-coding parts of the ornithine decarboxylase gene (ODC), is relatively variable in chelonians (Fritz et al., 2012; Praschag et al., 2017) and therefore looked promising to be also species-specific. Laboratory procedures for the nuclear loci followed Praschag et al. (2017) except that we applied newly designed internal primers for sequencing the R35 gene (forward: GCAAGGAAAAATGTTTG, reverse: ACGCTGACTCCATGCACA). The resulting R35 sequences were 1,101 bp long. The ODC sequences comprised a hardly readable simple-sequence-repeat (SSR) region, which could not be sequenced for all samples. We excluded this region from further analyses, yielding 739 bp length. Including some previously published data (Vargas-Ramírez et al., 2010), R35 sequences were available for 37 P. galeata and 16 P. subrufa. For the ODC gene, sequences were available for 31 P. galeata and 10 P. subrufa.

Phylogenetic analyses and uncorrected p distances of mtDNA

We concatenated individual mtDNA fragments for phylogenetic analyses and merged this dataset with previously published sequences, resulting in an alignment of 1,850 bp length that included 233 Pelomedusa sequences (also including sites outside South Africa and Namibia, and other species). Pelusios sinuatus served as the outgroup. For Pelomedusa galeata, P. subrufa and the outgroup, accession numbers and collection sites are given in Table S1; for other species, see Petzold et al. (2014), Fritz et al. (2015), and Nagy et al. (2015). We assessed the best partitioning scheme using PARTITIONFINDER (Lanfear et al., 2012) and the Bayesian Information Criterion (BIC). Accordingly, we partitioned the dataset using each codon position of the protein-coding genes, the 12S gene and the lumped DNA coding for tRNAs as a distinct partition. We inferred phylogenetic relationships under Maximum Likelihood using RAxML 7.2.8 (Stamatakis, 2006) and the GTR + G substitution model across all partitions. We performed five independent ML searches using different starting conditions and the fast bootstrap algorithm to explore the robustness of the results by comparing the best trees. Then, we calculated 1,000 non-parametric thorough bootstrap replicates and plotted the values against the best tree. In addition, we calculated uncorrected p distances for each mtDNA fragment using MEGA 7.0.21 (Kumar, Stecher & Tamura, 2016) and the pairwise deletion option.

Parsimony networks

We calculated for concatenated mitochondrial sequences a parsimony network using POPART (http://popart.otago.ac.nz). Since the underlying algorithm is sensitive to missing data, we excluded all individuals with lacking genes. In addition, we removed all individual missing sites and homologous data, resulting in an alignment of 1,602 bp length, which contained 87 sequences of P. galeata and 18 sequences of P. subrufa.

For network construction of nuclear data, we phased heterozygous R35 and ODC sequences using the PHASE algorithm in DNASP 5.10 (Librado & Rozas, 2009). For R35, we built two networks because the sequences of Vargas-Ramírez et al. (2010) were approximately 300 bp shorter than ours. One network comprised only our 86 phased sequences of 1,101 bp length, whereas the second also contained the previously published data (in total 106 phased sequences). It was based on an alignment trimmed to 699 bp to match the sequence lengths of Vargas-Ramírez et al. (2010).

Climatic niche models

To assess whether historical climate fluctuations influenced the distribution of Pelomedusa galeata and P. subrufa, we computed climatic niche models for present conditions as well as for the Last Glacial Maximum (LGM) and the mid-Holocene using the machine-learning algorithm MAXENT 3.3.3k (Phillips, Anderson & Schapire, 2006; Phillips & Dudík, 2008). To remove spatial autocorrelation, we filtered the genetically verified occurrences for each species, as well as for subclades Ia-Ic and lineage II of P. galeata, retaining only one record per sampling site, and supplemented the dataset with unambiguously assignable records from VertNet (http://vertnet.org/). The resulting datasets contained 31 localities for P. subrufa (23 own and eight VertNet records) and 51 for P. galeata (47 own and four VertNet records). Of the latter, 14 localities corresponded to subclade Ia, 10 to Ib, 17 to Ic, and 10 to lineage II.

We obtained eight uncorrelated bioclimatic predictors (R₂ < 0.75) with a spatial resolution of 2.5 arc minutes (∼5 km at the equator) from WorldClim (http://www.worldclim.org/) for current and past climatic conditions (mid-Holocene, ∼6,000 BP; LGM, ∼21,000 BP). For the mid-Holocene and LGM, we used datasets for three different general circulation model scenarios, namely the Community Climate System Model (CCSM4), the Model for Interdisciplinary Research on Climate (MIROC-ESM), and the Max-Planck-Institute Earth System Model P (MPI-ESM-P) that were statistically downscaled to a spatial resolution of 2.5 arc min. The following predictor variables were selected and clipped to a rectangular study extent: Bio 3 = isothermality (Bio2/Bio7)(*100), Bio 5 = maximum temperature of warmest month, Bio 7 = temperature annual range, Bio 8 = mean temperature of the wettest quarter, Bio 9 = mean temperature of the driest quarter, Bio 15 = precipitation of the wettest quarter, Bio 17 = precipitation of the driest month, and Bio 18 = precipitation of the warmest quarter.

We trained models separately for P. subrufa, P. galeata, and for subclades Ia-Ic and lineage II of P. galeata using circular buffers of 200 km surrounding the respective records and projected the results onto the full study extent for current and past conditions. In MAXENT, we selected linear, quadratic and hinge features to reduce model complexity and applied a bootstrapping approach with 100 replicates randomly splitting the records into 80% used for training and 20% for model evaluation. We performed a maximum number of 5,000 iterations and used the area under the curve AUC (Swets, 1988) for model evaluation. We used the average projection across 100 replicates for further processing, wherein we applied the “minimum training presence logistic threshold” as presence-absence threshold.

Phylogenetic analyses of mtDNA

The general branching pattern of our phylogenetic tree (Fig. 1) was in agreement with previous publications (Vargas-Ramírez et al., 2010; Wong, Fong & Papenfuss, 2010; Fritz et al., 2014; Fritz et al., 2015; Petzold et al., 2014) in that a well-supported major clade included all species and candidate species from the northern part of the distribution range of the genus Pelomedusa. This northern clade was sister to a weakly supported southern clade comprising Pelomedusa galeata and P. subrufa sensu stricto. Whilst the monophyly of P. subrufa was well supported, the monophyly of P. galeata received only weak bootstrap support of 61.

Pelomedusa galeata showed clear genetic structuring, with two well-supported clades (clade I and clade II) corresponding to deeply divergent mitochondrial lineages. One of these mitochondrial lineages (clade I) comprised three weakly supported subclades; Ia, Ib, and Ic. Subclade Ia was sister to a weakly supported, more inclusive clade with subclades Ib and Ic. Localities for subclade Ia were from the interior of South Africa, at high elevations above the Great Escarpment, where summer-rainfall prevails (Fig. 2; records in the provinces of the Free State, Gauteng, Limpopo, North West, Northern Cape). Samples in subclade Ib approximated geographically to the subtropical (low-elevation) summer-rainfall region along the northeast coast of South Africa (Eastern Cape, KwaZulu-Natal). The samples of subclade Ic were from the south coast and adjacent inland regions, mostly below the Great Escarpment, with all-year (aseasonal) rain (Eastern Cape, Western Cape). The other mitochondrial lineage (clade II) was from the western, winter-rainfall region of South Africa, represented by several samples from the southwestern Cape and the historical type specimen of Pelomedusa galeata devilliersi Hewitt (1935) from a site in the arid northwest of South Africa, close to the Namibian border. Even though P. subrufa also showed sequence variation, there was no obvious geographic pattern (Fig. 1 and Fig. S1). Considering the wide distribution of P. subrufa (Fig. 2), this is unexpected and contrasts with the phylogeographic structure found in P. galeata. Our new records for P. subrufa from the Lapalala Wilderness Area in the Waterberg region (Limpopo) currently represent the westernmost known occurrences for this species in South Africa. Before it was only known from the western border region of the Kruger Park.

Uncorrected p distances of mtDNA

Mean uncorrected sequence divergences between Pelomedusa galeata and P. subrufa amounted to 5.6% for the 12S fragment and 10.3% for the cyt b fragment, two genes that have previously been used for species delimitation in Pelomedusa (Petzold et al., 2014). Within all P. galeata (lineages I and II together) and within P. subrufa, the sequence divergence for 12S was 0.7% for each. For cyt b, the divergences were 2.0% for P. galeata and 1.3% for P. subrufa (Table 1). Between lineages I and II of P. galeata, the sequence divergence was 2.2% for 12S and 7.1% for cyt b. Additional values, and divergences for the mtDNA fragment comprised of the ND4 and tRNA genes, are summarized in Table 1.

Parsimony networks

The concatenated mitochondrial sequences were grouped in five distinct clusters (Fig. 3). Among the four clusters corresponding to the clades and subclades of Pelomedusa galeata, a maximum of 100 mutation steps occurred; among the haplotypes of P. subrufa, a maximum of 38 steps. Clade II of P. galeata was separated from the most similar haplotype of clade I (subclade Ia) by a minimum of 73 steps; subclade Ia differed from subclade Ib by a minimum of 13 steps, and from subclade Ic by 14 steps. Subclades Ib and Ic diverged by a minimum of 11 steps from one another. The haplotype clusters of P. galeata and P. subrufa were connected in a loop, with subclades Ia and Ib differing from P. subrufa by a minimum of 119 mutation steps. Clade II of P. galeata was separated by a minimum of 126 mutations from P. subrufa.

Even though much less variation occurred in the nuclear data, no haplotype sharing was observed between P. galeata and P. subrufa (Fig. 4). In the ODC network, haplotypes of P. galeata differed by a minimum of two mutations from haplotypes of P. subrufa. In the R35 network comprising the long sequences (1,101 bp), the two species differed by a minimum of seven mutations; in the R35 network with the short sequences (699 bp), the minimum was two steps. For the samples corresponding to the mitochondrial lineages within P. galeata, haplotype sharing was observed for the ODC gene between subclades Ia and clade II, and between subclades Ia, Ib, and Ic. However, unique haplotypes occurred in each subclade. Both R35 networks consisted of three haplotype clusters, one corresponding to P. subrufa, and the two others to lineages I and II of P. galeata, respectively. No haplotypes were shared between lineage I and lineage II. The number of mutations separating the two lineages of P. galeata resembled (1,101 bp network) or clearly exceeded (699 bp network) the divergence between P. subrufa and the two haplotype clusters of P. galeata. Haplotype sharing was observed for subclades Ia, Ib, and Ic of P. galeata, but to a lesser extent in the network based on the longer sequences. Unique haplotypes occurred for each subclade. The occurrence of shared haplotypes of the two lineages or the individual subclades of P. galeata showed no correlation with geography (Fig. S3), rather suggestive of ancestral polymorphism than of gene flow.

Climatic niche models

The models for current climatic conditions revealed distinct areas of climatically suitable space for the two Pelomedusa species (Fig. 5), with some overlap in north-eastern South Africa. The discrimination capability, expressed as AUC_test scores, of both models was high (P. subrufa = 0.79, P. galeata = 0.73), indicating a good discrimination of suitable and unsuitable areas. Across all 100 replicates, contributions of variables beyond 20% showed both models to be strongly affected by temperature and precipitation during the driest month and quarter (Table 2: Bio 9, Bio 17). However, the model for the widespread P. subrufa was also influenced by temperature during the wettest quarter (Bio 8 > 10%), whilst P. galeata was revealed to be strongly dependent on temperature annual range (Bio 7 >20%) and precipitation during the warmest quarter (Bio 18 > 10%), two variables less important for P. subrufa. Potentially suitable areas for the subclades of P. galeata were highly distinct (Fig. 5) and described by different predictors. For subclade Ia, inhabiting a highland summer-rainfall area, temperature annual range and mean temperature of the wettest and driest quarter contributed the most (Bio 7 and Bio 8 > 20%, Bio 9∼20%). However, the model was also impacted by precipitation during the driest month (Bio 17 > 10%). For subclade Ib, inhabiting a subtropical low elevation summer-rainfall area, the model was exclusively shaped by temperature-related variables (Bio 7∼40%; Bio 3 and Bio 9 >20%), whilst precipitation-related predictors contributed less than 5%. For subclade Ic, inhabiting the southern coast and adjacent inland regions characterized by aseasonal rainfall, variable contribution was highest for precipitation during the driest month (Bio 17 = 36%), precipitation seasonality, and temperature during the wettest quarter (Bio 15, Bio 8 > 10%). For subclade II from the western winter-rainfall area of South Africa, precipitation of the wettest quarter was the variable with the highest importance (Bio 16∼77%), whilst no other predictor exceeded 10%. All models for the subclades of P. galeata received excellent AUC_test scores (Ia = 0.81, Ib = 0.75, Ic = 0.78, II = 0.93) and suggest mutually exclusive climatic niches.

Projections onto past climatic conditions found similar patterns across all three scenarios (Figs. 6 and 7, Figs. S4–S7) with high discrimination abilities (MIROC-ESM: P. subrufa = 0.78, P. galeata = 0.73, Ia = 0.84, Ib = 0.77, Ic = 0.78, II = 0.95; CCSM4: P. subrufa = 0.73, P. galeata = 0.80, Ia = 0.86, Ib = 0.77, Ic = 0.80, II = 0.91; MPI-ESM-P: P. subrufa = 0.71, P. galeata = 0.80, Ia = 0.85, Ib = 0.81, Ic = 0.80, II = 0.91).

Projections onto climatic conditions of the LGM revealed for P. subrufa a suitable area resembling the present situation, whilst the space for P. galeata shifted considerably, excluding a large area in the center of the extant distribution range. When the models for subclades Ia-Ic and lineage II of P. galeata were inspected individually, the models for subclade Ic and lineage II matched well with extant conditions. However, the models for subclades Ia and Ib shifted to the northeast, outside the current niche, but remained geographically mutually exclusive (Figs. 6 and 7). All projections onto LGM conditions yielded high AUC_test scores (MIROC-ESM: P. subrufa = 0.80, P. galeata = 0.71, Ia = 0.85, Ib = 0.77, Ic = 0.76, II = 0.91; CCSM4: P. subrufa = 0.77, P. galeata = 0.73, Ia = 0.83, Ib = 0.77, Ic = 0.78, II = 0.94; MPI-ESM-P: P. subrufa = 0.79, P. galeata = 0.72, Ia = 0.85, Ib = 0.77, Ic = 0.76, II = 0.94).