Total DNA was extracted from tissue culture cells of M. javanica (Thailand, tissue code: 1999.508) and muscle samples of P. tricuspis (Gabon, tissue code: 2005.GLC14) and S. temminckii (Tanzania, tissue code: SUN106) using the QIAGEN DNeasy Tissue Kit (Qiagen, Hilden, Germany).

Twenty-five overlapping fragments of the mitochondrial genome were amplified by PCR using primers published in previous studies [5–8], as well as four new primers specially designed for this study: DLU400 or DLU400M1/12SL41 or 12SL200; 12SU1230/12SL2226M1LA; 12SU829/16SL518; 16SU365/16SL1056; 16SU946/N1L64; Uleu or U16S1421/LMet or LMet2; N1U840/N2L492; N2U354M2 or IleU/AsnL; TrpU/C1L705; UTyr or C1U246/C1L1017; C1U897M1/C2L15M1; SerU or SerUM1/A8L1; C2U603M4/C3L168; A6U654M1 or A6U654M2/GlyLM1; C3MANIU (5′- CAA TAT ATC AAT GAT GAC GTG A-3′)/LARGMANI (5′- GTT GAY TTG TTT GTG ATG CTC A-3′); C3U780 or C3U780M3/N4L366 or N4L366M2; UArg or UArgM1/N4L918M1 or N4L1071R; N4U681/Leu2L or Leu2LM1; Ser2U/N5L652 or N5L652M1; N5U501/N5L1214; N5U1146M5/N5L154M5; N6UT298 (5′- GGC TCA ATC AAA CTA TAC TTC CT-3′)/CBLT298 (5′-AGT TTC ATC ATG CTG AGA TG-3′); N6RU102 or N6RU102M1/CBL402 or CBL402M1; GluMA/ProMA; CBU162 or U844/L482; and UThr13/L482. Amplifications were done in 20 μl using 3 μl of Buffer 10× with MgCl₂, 2 μl of dNTP (6.6 mM), 0.12 μl of Taq DNA polymerase (2.5 U, Qiagen, Hilden, Germany) and 0.75 μl of the two primers at 10 μM. The standard PCR conditions were as follows: 4 min at 94 °C; 5 cycles of denaturation/annealing/extension with 45 s at 94 °C, 1 min at 60 °C and 1 min at 72 °C, followed by 30 cycles of 30 s at 94 °C, 45 s at 55 °C, and 1 min at 72 °C, followed by 10 min at 72 °C.

Both strands of PCR products were sequenced using Sanger sequencing on an ABI 3730 automated sequencer at the Centre national de séquençage (Genoscope) in Évry (France). The sequences were edited and assembled using Sequencher 5.1 (Gene Codes, Corp., Ann Arbor, MI, USA).

The analyses based on complete mitochondrial genomes included five species of the order Pholidota: M. javanica, M. pentadactyla, P. tetradactyla, P. tricuspis, and S. temminckii. Four genera, representing different orders of Laurasiatheria, were used to root the Pholidota tree: Boselaphus (Cetartiodactyla), Canis (Carnivora), Ceratotherium (Perissodactyla), and Rhinolophus (Chiroptera). GenBank accession numbers of the DNA sequences are provided on the tree of Fig. 2. Complete mitochondrial genomes were aligned using MUSCLE [9] and then further adjusted by eye with Se-Al v2.0a11 [10]. All ambiguous regions, i.e., involving ambiguity in the position of gaps, were excluded from the analyses to avoid erroneous hypotheses of primary homology. For this reason, the control region was not included in the final alignment. The length of the reduced alignment is 14,926 nt (available upon request to AH). The phylogenetic tree was constructed using PAUP* version 4 [11] based on the Maximum Likelihood (ML) method and the GTR + I + G model selected by the Akaike information criterion under jModelTest 2 [12]. Bootstrap percentages (BP) were computed using 1000 replicates.

BLASTN searches were performed on NCBI [13] using the mtDNA genome of P. tricuspis as a query to extract all available mitochondrial sequences of Manidae. The sequences were filtered by mitochondrial markers (COI, Cytb, 12S, 16S, and CR) and aligned as explained above. Uncorrected nucleotide distances (D) were calculated in PAUP* for each datasets and the Neighbour-Joining (NJ) method was applied to reconstruct clusters of similar sequences, with BP computed after 1000 replicates.

Protein-coding, ribosomal and transfer RNA genes were identified by comparison with the two published genomes of Pholidota, M. pentadactyla [3] and P. tetradactyla [4]. The annotated sequences of M. javanica, P. tricuspis, and S. temminckii were deposited in GenBank under accession numbers KP306514-KP306516. The three new mitochondrial genomes of Manidae represent circular double stranded DNA molecules that range between 16,559 and 16,576 bp in length, and that contains the 37 genes described in Table 1. Among these genes, only ND6 and eight tRNAs are encoded by the L-strand, whereas all the other genes are encoded by the H-strand. The gene arrangement is the same as that found in other mammalian genomes. All protein-coding genes of the mtDNA have a methionine start codon (ATR), except ND5 (ATT). Most protein-coding genes appear to be terminated by TAR, although this stop codon is incomplete in the COIII, ND3 and ND4 genes. In M. javanica, the COI, Cytb, and ND6 genes are terminated by AGA or AGG.

The control region is highly conserved among the five species of Manidae, both in structure and length (sequence length varies between 1156 and 1169 bp), with three tandem RS2 repeats of 78/79 bp in the L domain, which are located just after a partial RS2 motif, i.e., “YATGTATAATCGTGCAT” (pos. 15,453 in M. javanica).

The overall base composition of the mt genome showed marked differences between the two Asian species (Manis) and the three African species (Phataginus and Smutsia), with higher percentages of C (29.1 vs. 23.2–24.0%) and G nucleotides (14.7 vs. 13.3–13.8%), and lower percentages of A (33.1 vs. 33.2–33.9%) and T nucleotides (23.1 vs. 28.7–29.9%) (Fig. 1).

Maximum Likelihood analyses based on an alignment of 14,926 bp resulted in a robust tree (Fig. 2) with maximal BP (100%) for interspecific relationships within the family Manidae. The results are in full agreement with the morphological study of Gaudin et al. [1]: the genera Manis and Phataginus are found to be monophyletic and the family Manidae is divided into two geographic groups corresponding to the Asian subfamily Maninae (Manis) and the African subfamily Smutsiinae (Smutsia and Phataginus). However, the terminal branch lengths of the two species of Manis are unexpectedly short in the ML tree (Fig. 2). Nucleotide distances indicate that the mitochondrial genome of M. javanica is 99% identical to that of M. pentadactyla published by Qin et al. [3] (NC_016008), suggesting species misidentification for one of them. To solve this issue, the mitochondrial genomes were compared to all homologous fragments available for Pholidota in GenBank. Five mtDNA alignments were analysed using the NJ clustering method, i.e., COI (38 taxa and 658 nt), Cytb (35 taxa and 402 nt), 12S (38 taxa and 395 nt), 16S (25 taxa and 521 nt), and CR (56 taxa and 1264 nt). In Fig. 3, the COI tree shows the existence of six robust clusters (BP = 100) corresponding to the following species: M. javanica, M. pentadactyla, P. tricuspis, P. tetradactyla, S. gigantea and S. temminckii. On the one hand, the NC_016008 genome of M. pentadactyla fell into the M. javanica cluster (0.2 < D < 3.5 %), which is composed of representatives from China, Malaysia and Thailand. On the other hand, the NC_016008 genome appears to be highly divergent from COI sequences of M. pentadactyla from Taiwan and China (provinces of GuangDong, Hainan, and Hunan) (13.2 < D < 13.7 %). Similar results were found for three other markers (Cytb, 12S and CR) using mitochondrial sequences produced by independent teams (Appendices 1, 2 and 4). The analyses suggest, therefore, that the NC_016008 genome published in Qin et al. [3] belongs to M. javanica rather than M. pentadactyla. However, the specimen was collected from Guangxi, a Chinese province where M. javanica is not currently recorded (Fig. 4). Accordingly, two hypotheses can be proposed to explain the presence of M. javanica in Guangxi: (1) the specimen was not native to Guangxi, but was imported in China from a Southeast Asian country (Vietnam, Lao PDR, Myanmar, Cambodia, Malaysia, Thailand or Indonesia) to satisfy the persistent demand for meat and scales used in traditional medicines [14]; (2) alternatively, the species M. javanica may also occur in sympatry with M. pentadactyla in the Guangxi province, implying that the geographic range currently provided by the IUCN [2] for M. javanica is incomplete. To resolve this and to ensure that optimal conservation strategies are in place, the possible occurrence of M. javanica in Guangxi should be reassessed urgently.

A similar problem was found for the mitochondrial genome of P. tetradactyla published by Arnason et al. [4]. Indeed, the NC_004027 genome is highly divergent from mitochondrial sequences of P. tetradactyla produced in different labs (COI: 13%; Cytb: 15%; 12S and 16S: 6%), whereas it shares 99–100% of identity with COI, Cytb, 12S, and 16S sequences of P. tricuspis obtained by independent teams (Fig. 3 and Appendices 1–3). Here again, the most likely hypothesis involves species misidentification of the specimen sequenced by Arnason et al. [4].

Within P. tricuspis, the analyses of mitochondrial data revealed high nucleotide divergence between the Gabonese pangolin sequenced for this study and the specimens collected in Ghana, Nigeria and Cameroon: 5.9–6.7% for COI, 8.0–9.5% for Cytb, 1.5–2.6% for 12S, 2.4–3.3% for 16S (Fig. 3 and Appendices 1–3). Given that mitochondrial markers are available for all species of pangolin currently recognized in Africa, and that the holotype of P. tricuspis was described from Guinea in West Africa [15], these comparisons suggest the existence of a new species of tree pangolin in Gabon. This finding could be consistent with recent studies supporting the existence of a strong biogeographic barrier (named charnière climatique) between Cameroon and Gabon for both animals and plants endemic to tropical rainforests [16,17]. Further studies are needed to determine the morphology, geographic distribution and conservation status of this new taxon.