Complete mitochondrial genome of the endangered Mary River turtle (Elusor macrurus) and low mtDNA variation across the species’ range

Additional keywords: control region, D loop, freshwater turtle, MiSeq, next generation sequencing.

The Mary River turtle (Elusor macrurus) represents an ancient lineage of short-necked turtles restricted to a single coastal drainage of eastern Australia (Georges and Thomson 2006). Deep phylogenetic relationships among Australo-Papuan short-neck chelids, including placement of the monotypic genus Elusor, are not yet resolved (cf. Georges et al. 1999; Le et al. 2013; Spinks et al. 2015). A recent treatment shows Elusor as sister lineage to either Myuchelys purvisi or Rheodytes leukops, depending on whether mitochondrial or nuclear DNA datasets are analysed (Spinks et al. 2015). The narrow distribution of E. macrurus, which is limited to the Mary River in south-eastern Queensland, combined with a history of exploitation for the pet trade and recruitment failure due to nest predation, have conspired to make it one of the most threatened species of freshwater turtle in Queensland (Limpus 2012; Micheli-Campbell et al. 2013). The species is listed as ‘Endangered’ at state level under the Nature Conservation Act 1992 and at national level under the EPBC Act 1999, and recognised as endangered globally in the IUCN Red List of Threatened Species 2015-4. Complete mitochondrial genomes (mitogenomes) have proven utility in molecular systematics at both deep and shallow taxonomic scales (e.g. Duchene et al. 2012). Few mitogenomes are available for Australo-Papuan chelid turtles and here we document the first complete mitochondrial genome sequence of E. macrurus and also conduct an assessment of mtDNA variation within the restricted natural range of this highly threatened taxon.

A total of 22 E. macrurus individuals were sampled from eight sites along the main trunk of the Mary River and a tributary separated by estuary, Tinana Creek (Fig. 1; Supplementary Material Table S1). All samples were collected using double-winged fyke nets set overnight in both upstream and downstream directions. A small section of skin (5 mm²) was taken from along the webbing of the hind foot, and preserved immediately in 100% ethanol. Tissue from voucher MRT4 was used for extraction of whole genomic DNA using the DNeasy blood and tissue kit (Qiagen). DNA was sheared to an approximate mean length of 400 bp using the M220 Focused-ultrasonicator (Covaris) and an Illumina MiSeq-compatible sequencing library was prepared. A double-index two-step library preparation was used, with all steps performed in the presence of solid-phase reversible immobilisation (SPRI) beads (iTru protocol: Travis Glenn, pers. comm.). The final full-length library construct was compatible with an Illumina TruSeq library except for unique 8-bp indexes within read 1 and read 2 adapters (Faircloth and Glenn 2012). Sequencing was performed on the Illumina MiSeq platform at Griffith University DNA Sequencing Facility using a 600-cycle MiSeq reagent kit v3, running 2 × 300-bp paired-end reads. The E. macrurus library was sequenced in parallel with four other libraries.

Fig. 1.  Map of the study area. Sample localities for each specimen are denoted by circles. Locality details for Sites 1–8 are given in Supplementary Material Table S1. The catchment boundary for the Mary River is highlighted in bold.

A total of 1.38 × 10⁷ paired-end reads were generated from the E. macrurus shotgun library. Mitogenome assembly was performed using Geneious v8.1.7 (Kearse et al. 2012). Overlapping paired reads were first merged using the FLASH v1.2.9 plugin (minimum overlap 20, maximum overlap 200) (Magoc and Salzberg 2011). Then 5 × 10⁵ merged reads in size range 301–550 bp were sampled. De novo and reference-guided assembly approaches were used for comparison (Schmidt 2015). De novo assembly was implemented in Geneious using default medium-sensitivity settings and allowing contigs with matching ends to circularise. A putative mitogenome was identified from the resulting contigs as the longest contig with circular topology. The putative mitogenome assembly was visually checked for errors derived from heteroplasmy or paralogues using the Geneious genome browser. A low-coverage region occurred towards the start of the putative control region from position 15 553 to 15 688 (mean coverage 3.2×; minimum 1; maximum 5) (Fig. 2). Subsequent mapping of an extra 680 000 merged reads (size range 200–300 bp) to the draft de novo–assembled mitogenome improved coverage in this region to a mean of 10×. A set of primers was designed to amplify across this low-coverage region in two overlapping fragments, spanning a total of 1270 bp from position 15 129 to 16 398 (Fig. 2). Primer pair MRT15129F (5ʹ-TTCGCCTATGCCATCCTACG) and MRT15705R (5ʹ-TGCTAGAGGTAAATAAATTTATGCACG) amplify a 603-bp fragment including partial cyt b, complete trnTh^r, trn^Pro and partial control region. Primer pair MRT15599F (5ʹ-CAACCACACCCTATCCGACA) and MRT16379R (5ʹ-TCGACACTGCACTTGGTGTA) amplify a 800-bp fragment of the control region. Overlap between these two fragments was 133 bp (Fig. 2). Standard PCR conditions of 35 cycles with annealing temperature of 58°C were used. Sequencing of PCR amplicons was performed by Macrogen Inc. (Seoul, Korea). Mitogenome annotation was performed with the Mitos webserver (Bernt et al. 2013) followed by adjustment of gene boundaries in Geneious v8.1.7. Folding of a putative light strand origin of replication (O_L) was examined with the Vienna RNAfold webserver (Gruber et al. 2008), and a search for tandem repeats in the control region was performed with Tandem Repeat Finder (Benson 1999).

Fig. 2.  Distribution and depth of coverage of the Elusor macrurus mitochondrial genome. Gene annotations are listed for protein-coding genes, rRNAs and the control region. Black triangles depict position and orientation of tRNAs. See GenBank accession KC692461 for complete annotation details. Arrow indicates low-coverage region. Black bars represent position of two overlapping PCR amplicons used for mitogenome validation and for assessment of range-wide mtDNA variability. Fragment 1 corresponds to primer pair MRT15129F/MRT15705R, Fragment 2 corresponds to primer pair MRT15599F/MRT16379R.

De novo assembly of 5 × 10⁵ merged raw reads produced a 16 499-bp circular contig from 1053 reads with 23× mean coverage (s.d. 5.4; minimum 1; maximum 41) (Fig. 2). Mapping to two chelid reference mitogenomes (Elseya branderhorsti GenBank accession: NC_026047; Emydura subglobosa GenBank accession: NC_026048) produced the same 16 499-bp contig from either 1057 reads (NC_026047), or 1044 reads (NC_026048), both spanning the entire reference at 23× mean coverage (s.d. 5.4; minimum 1; maximum 41). The mean quality score for 1057 mapped reads was Q36, with 93% at Q30 or higher. The three independent assemblies (one de novo and two map-to- reference) produced contigs 100% identical in length and sequence identity. The top 10 blastn matches for this sequence from the NCBI nucleotide database were turtle complete mitochondrial genomes, with pairwise identity to the query sequence ranging from 78.1% to 88.9%. The mitogenome of Elseya branderhorsti used for reference-guided assembly had 88.4% sequence identity with the new E. macrurus mitogenome. Validation of the E. macrurus mitogenome was achieved by sequencing two overlapping PCR amplicons using primers MRT15129F/MRT15705R and MRT15599F/MRT16379R, which produced an edited 1129-bp contig that aligned with 100% identity to position 15 196–16 324 of the mitogenome consensus sequence.

The complete mitogenome of E. macrurus is available from GenBank under accession number KU736930. It was 16 499 bp in length, with 13 protein-coding genes, 22 tRNAs, 2 rRNAs and a 918-bp control region (Fig. 2; Supplementary Material Fig. S1A). All 37 genes and the control region were arranged in the order observed in other Chelidae (Wang et al. 2012; Zhang and Georges 2014). Base composition was A (34.9%), C (26.6%), G (12.7%), T (25.8%). A portion of tRNA^Asn and tRNA^Cys along with an 18-bp intergenic region between these tRNAs folded into a stem and loop secondary structure (Supplementary Material Fig. S1B). This structure may be homologous with the putative origin of light strand replication (O_L) proposed for Chelodina longicollis by Zhang and Georges (2014). There was no evidence for extensive tandem repeats in the control region as observed in mitogenomes of some turtle taxa (e.g. Wang et al. 2012). The longest tandem repeat detected was (AT)₅. The new mitogenome was compared with publicly available sequences of mitochondrial origin matching E. macrurus. Pairwise identity >99% was found for partial fragments of 12S rRNA (GenBank accession: U40639), 16S rRNA (AF113622), COI (GenBank HQ329617, AF113646, KP876791; KP876792), ND4 (GenBank KC755124; Dryad doi:10.5061/dryad.tf8q1), and cyt b (Dryad doi:10.5061/dryad.tf8q1). All differences between the mitogenome consensus sequence and GenBank accessions were due to the presence of ambiguous base calls in GenBank sequences (16S rRNA from Georges et al. 1999), or to base differences within the first 22 bp at the 5ʹ-end of GenBank sequences (12S rRNA, COI from Georges et al. 1999; Seddon et al. 1997). Sequencing of PCR amplicons for 21 additional E. macrurus samples, using primer sets MRT15129F/MRT15705R and MRT15599F/MRT16379R yielded a full-length edited fragment of 1127–1129 bp, including 246-bp cyt b, complete trnTh^r and trn^Pro, and 743 bp of the control region. Two haplotypes, designated A and B, were identified. Haplotype A (1129 bp; GenBank Accession: KX369542) was found in 18 individuals and is identical to the mitogenome reference sequence KU736930. Haplotype B (1127 bp; GenBank accession: KX369543) was found in 4 individuals and differed by one substitution A > G at Position 16 059 as well as by contraction of the (AT)₅ tandem repeat unit to (AT)₄ at Position 16284–16293. Both of these substitutions were located in the control region. Haplotype A was distributed throughout the study area and Haplotype B was located at Sites 2 and 3 in Fig. 1 (see Supplementary Material Table S1).

We examined a large portion (81%) of the mitochondrial control region, which is the most hypervariable mtDNA region in the Chelidae (Zhang and Georges 2014). Our sample included most of the known geographic range of E. macrurus, including Tinana Creek, which represents a distinct genetic subpopulation in both the Australian lungfish and Mary River cod (Huey et al. 2013; Hughes et al. 2015). Samples from Tinana Creek shared the same control region haplotype (Haplotype A) found throughout the main stem of the Mary River. Limited mtDNA variation means that further data such as polymorphic microsatellites will be required to assess fine-scale genetic subdivision within the Mary River catchment. Overall, control region variation detected across the range of E. macrurus was low relative to most chelid species examined at the within-drainage scale (Souza et al. 2003; Todd et al. 2013, 2014a, 2014b; Georges et al. 2014; Hodges et al. 2014, 2015). One exception is the Bellinger River turtle (Myuchelys georgesi), a critically endangered species with a single control region haplotype known from a sample consisting of 428 bp sequenced in 16 individuals (Georges et al. 2011). Like that of E. macrurus, the range of M. georgesi is restricted to a single small drainage.

Comparison of the new mitogenome sequence with publicly available partial mtDNA sequences of E. macrurus did not yield convincing evidence for additional mtDNA nucleotide variation in this species. All differences detected were ambiguous base calls or variation clustered at the beginning of old sequence accessions, which may be due to editing issues rather than real variation. The distribution of E. macrurus is restricted to a single coastal drainage so we should expect natural levels of mtDNA variability to be limited, as observed in other chelids at the within-drainage scale (e.g. Georges et al. 2011; Todd et al. 2014a). In addition to this, E. macrurus has undergone a dramatic population decline over the last 50 years, attributed to egg harvesting for the pet trade compounded by habitat modification and increased predation (Limpus 2012; Micheli-Campbell et al. 2013). Extreme population bottlenecks of this kind can remove most or all mtDNA variation within a population (Wilson et al. 1985).

Previous studies suggested that absence of an origin of light strand replication located between tRNA^Asn and tRNA^Cys may be characteristic of suborder Pleurodira (Wang et al. 2012). However, evidence for folding of this intergenic region in the new mitogenome of E. macrurus, along with similar evidence for Chelodina longicollis, suggests that this is not the case (Zhang and Georges 2014). Deep phylogenetic relationships among Australo-Papuan chelid genera have proved challenging to resolve (Georges et al. 1999; Georges and Thomson 2006; Le et al. 2013; Spinks et al. 2015). Uncorrected pairwise distance between complete mitochondrial genomes of E. macrurus and those of other Chelidae currently available on the NCBI nucleotide database was >10%, so the new E. macrurus mitogenome will be a useful resource, particularly for resolving the deep mtDNA relationships among Australo-Papuan short-necked Chelidae.