Genomic analysis reveals a cryptic pangolin species

Abstract

Eight extant species of pangolins are currently recognized. Recent studies found that two mitochondrial haplotypes identified in confiscations in Hong Kong could not be assigned to any known pangolin species, implying the existence of a species. Here, we report that two additional mitochondrial haplotypes identified in independent confiscations from Yunnan align with the putative species haplotypes supporting the existence of this mysterious species/population. To verify the new species scenario we performed a comprehensive analysis of scale characteristics and 138 whole genomes representing all recognized pangolin species and the cryptic new species, 98 of which were generated here. Our morphometric results clearly attributed this cryptic species to Asian pangolins (Manis sp.) and the genomic data provide robust and compelling evidence that it is a pangolin species distinct from those recognized previously, which separated from the Philippine pangolin and Malayan pangolin over 5 Mya. Our study provides a solid genomic basis for its formal recognition as the ninth pangolin species or the fifth Asian one, supporting a new taxonomic classification of pangolins. The effects of glacial climate changes and recent anthropogenic activities driven by illegal trade are inferred to have caused its population decline with the genomic signatures showing low genetic diversity, a high level of inbreeding, and high genetic load. Our finding greatly expands current knowledge of pangolin diversity and evolution and has vital implications for conservation efforts to prevent the extinction of this enigmatic and endangered species from the wild.

Keywords: pangolins; population genomics; species delimitation.

Fig. 1.

Thirty-three scales from multiple confiscations in Hong Kong and Yunnan representing the published H04 and H08 haplotypes and the novel YN01 and YN02 haplotypes of the cryptic new species. (A) Five scales attached to skin fragments. (B) Three scales attached to claws, including two from the hindfoot (H04) and one from the forefoot (H08). (C) A scale from the dorsum of the tail whose origin was deduced from the fact that it is smaller than scales from other body parts, has a V-shaped posterior margin, and lacks a ridge bulge. (D) Thirteen scales from the dorsal side, whose origins were deduced from their large size, their fan-like, rhombic, or inverted trapezoidal shape, and their lack of ridge bulges. (E) Three scales from the ventral side, whose origins were deduced from their small size and the presence of a ridge in the middle of the scale. (F) Two scales from the head, whose origin were deduced from their small size and fan-like or almost round shape. (G) One scale from the lateral rows of the tail, whose origin was deduced from its small size and reduced triangular shape. (H) Five scales from unknown body parts. All body part assignments (C–G) are based on previously reported characteristics (–23).

Fig. 2.

Phylogenomic analyses, genetic difference evaluations, species delimitation, and identification of genomic island of differentiation between the cryptic new species and other Asian pangolin species. (A) Gene tree (ML) inferred from autosomal and X chromosome SNPs. They both obtained a well-supported and consistent phylogenetic tree topology. (B) Gene tree (ML) inferred from mitogenome sequences. (A and B) The number of individuals for each species is indicated (SNP/mitogenome). (C) Species tree (SNAPP) inferred from the autosomal SNPs. (D) Gene tree (ML) and species tree (MP-EST and ASTRAL) inferred from concatenated single-copy coding sequences. Interspecific divergence times (nodes from a to j) are shown on the branches, including median values and their 95% HPDs (purple bars). The red boxes on the nodes (A–D) indicate full bootstrap values (ML, MP-EST, and ASTRAL BS = 100%; SNAPP BP = 1.0) at the nodes representing relationships among pangolin species. (E) Pairwise comparisons of nucleotide site differences between the cryptic new species and the four recognized Asian pangolin species across 50-kb genomic windows. (F) Species delimitation analyses based on autosomal SNPs using the BFD* method. (G–L) Genomic regions of elevated differentiation identified by performing Dxy and Fst analyses between the cryptic new species and the Chinese pangolin (G), Malayan pangolin (H), Indian pangolin (I); and between the Indian pangolin and the Chinese pangolin (J), Malayan pangolin (K); and between the Chinese pangolin and the Malayan pangolin (L). Data points above the horizontal dotted line (corresponding to the 5% right tail of the empirical Fst distribution) and on the right side of the vertical dotted line (the 5% right tail of the empirical Dxy distribution) were identified as genomic islands of differentiation for each species pair. Reproduction-related genes with significant signals are indicated by arrows and red circles.

Fig. 3.

Genetic diversity, inbreeding levels, and genetic load evaluation. (A) Whole-genome genetic diversity based on the average heterozygosity (He) among all individuals for each species. (B) Inbreeding levels calculated based on the average inbreeding coefficient (F_ROH) among all individuals for each species. (C) Long-ROHs (>1 Mb; ROH_1Mb) detected in the new pangolin species and three other Asian pangolin species. Evaluations of genetic load based on LOF mutations (D), missense mutations (E), deleterious GS mutations (F), and synonymous mutations (G) among all individuals for each species. Error bars indicate SD values based on 100 resamples.

Fig. 4.

Demographic history reconstruction and gene flow analyses. (A) Results of a PSMC analysis of Asian pangolin species, including the new species (shown in the upper half of the figure), and African pangolin species (lower half) based on a mutation rate per generation of 1.47 × 10⁻⁸ (43) and an average generation time of 1 y (44). Note that 100 bootstrap replicates were performed to assess variance in the demographic trajectories. Glacial periods, including the PG, LGP, and last interglacial period are indicated by background shading. Surface temperatures (temp.) over the last 1 Mya based on δ ¹⁸O data (45) are indicated by dotted lines. (B–D) The gene flow analysis results for Asian pangolins, including the new species, using the D statistic (B), f4-ratio test (C), and the f-branch method (D). The depth of the color for each block in the matrix represents the degree of migration ratio. The darker color means the higher migration ratio. The gray shading corresponds to tests that are not applicable to the provided phylogeny.