Pangolin genomes and the evolution of mammalian scales and immunity

Pangolin genome, divergence, and heterozygosity

A female Malayan pangolin derived from a wild specimen from Malaysia and a female Chinese pangolin from Taiwan were sequenced using whole-genome shotgun sequencing strategies to a coverage ∼145× and ∼59× based on the estimated genome size of 2.5 gigabases (Gbp) and 2.7 Gbp, respectively (Supplemental Information 1.0). We identified 23,446 and 20,298 protein-coding genes in the Malayan and Chinese pangolin genomes, respectively, using the MAKER annotation pipeline based on several sources of evidence from ab initio gene prediction, transcriptomic data, and protein evidence from the Canis familiaris reference genome (Broad CanFam3.1/canFam3) and transcriptome (Lindblad-Toh et al. 2005; Cantarel et al. 2008). The assemblies of both Malayan and Chinese pangolins capture at least 99.7% and 94.8% of the expressed genes assembled from RNA-seq data generated from eight different organs, respectively (Supplemental Information 2.1). The lower overall mapping rate of the Chinese pangolin is partially explained by the divergence of Malayan and Chinese pangolin, although it could also be impacted by the different sequencing depths. Furthermore, Core Eukaryotic Genes Mapping Approaches (CEGMA) and genome comparison analyses indicated that these genomes are good candidates for genome annotation (Supplemental Information 2.2, 2.3, 3.0).

Coalescent dating analysis suggests that pangolins diverged from their closest relatives, the Carnivora, ∼56.8–67.1 million years ago (MYA) (Supplemental Information 4.0) and that Malayan and Chinese pangolin species diverged from each other ∼4–17.3 MYA. Malayan pangolin has a more than three times higher heterozygosity rate than Chinese pangolin (1.55 × 10⁻³ and 0.4 × 10⁻³, respectively), which likely reflects the demographic impacts of a founder population and smaller effective population in the latter species (Supplemental Information 5.0). Interestingly, analysis of the trajectory of historical effective population size of the two pangolins using the pairwise sequential Markovian coalescent model showed an inverted pattern during the Middle Pleistocene, reflecting the influence of climate and sea level on both pangolin species (Supplemental Fig. S6.1).

Pangolin-specific phenotypes

To gain insight into possible genomic patterns linked with the unique traits of pangolins, we first searched for pseudogenized (presumably loss-of-function) genes. Gene loss through pseudogenization is increasingly recognized as an important factor in lineage and adaptive phenotypic diversification (Wang et al. 2006; Albalat and Canestro 2016). Given that pangolins are edentulous mammals, we screened 107 tooth development-related genes and found three pseudogenized candidate genes in pangolins. Most notable of these was ENAM, the largest protein of the enamel matrix and essential for normal tooth development (Meredith et al. 2009). ENAM mutations can lead to enamel defects (Hart et al. 2003). We identified several frameshift indels and premature stop codons in the ENAM genes of both the Malayan and Chinese pangolins (confirmed by Sanger sequencing; sample size = 8) (Fig. 1A; Supplemental Information 7.0; Meredith et al. 2009). ENAM pseudogenization has also been reported in the African tree pangolin (M. tricuspis) (Meredith et al. 2009). We also identified two other pseudogenized genes that have been linked with normal tooth development, amelogenin (AMELX) and ameloblastin (AMBN). These genes shared common mutations in both Asian species and all four African species (M. tricuspis, M. tetradactyla, M. gigantea, and M. temminckii) that we examined, suggesting that pseudogenization occurred early in the evolutionary history of the Pholidota lineage (Supplemental Figs. S7.1, S7.2). Loss of function of these genes through pseudogenization has also been reported in other edentulous vertebrates such as toothless baleen whales, birds, and turtles (Meredith et al. 2014).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Case studies of pseudogenized genes in pangolins. (A) Three tooth development-related genes were pseudogenized and may be related to the lack of teeth in pangolins. There was a frameshift deletion in positions c1311–c1324, a single base pair deletion in c1476, and an insertion of AGAT at position c1621, resulting in another premature stop codon in the AMBN gene. The AMELX gene contains a large deletion at position c438–c497. (B) Two genes were pseudogenized and may be related to the poor vision of pangolins. Blue = insertion, green = deletion, pink = stop codon.

Pangolins are mostly nocturnal (the only exception being the long-tailed pangolin) and are thought to have poor vision (Soewu and Ayodele 2009). We screened 217 vision-related genes and identified several candidate genes that were pseudogenized (Fig. 1B). One of these, BFSP2, encodes a lens-specific intermediate filament-like protein that is a unique cytoskeletal element in the ocular lens of vertebrates (Sandilands et al. 2003). Interestingly, one frameshift insertion and three premature stop codons were consistently identified in the BFSP2 genes of both pangolin species (confirmed by Sanger sequencing; sample size = 8), suggesting possible loss of optical clarity, which may lead to progressive cataracts (Alizadeh et al. 2004). Another candidate gene, guanylate cyclase activator 1C (GUCA1C), expresses proteins that stimulate photoreceptors GC1 and GC2 at low concentrations of free calcium but inhibit GCs when the concentrations of free calcium ions are elevated, which is important for regulating the recovery of the dark state of rod photoreceptors following light exposure (Haeseleer et al. 1999). We found identical frameshift mutations and premature stop codons in GUCA1C in both pangolin species (confirmed by Sanger sequencing; sample size = 8), suggesting reduced stimulation of GC1 and GC2 and reduced rates of phototransduction in pangolins (Stephen et al. 2006). The GUCA1C gene was also pseudogenized in all African species with identical mutations (Supplemental Fig. S7.3). Interestingly, the pseudogenization of BFSP2 and GUCA1C have been reported in mice with reduced vision, supporting the hypothesis that the two genes are likely associated with the poor vision of pangolins (Imanishi et al. 2002; Sandilands et al. 2004; Song et al. 2009).

As an evolutionary consequence of being covered by scales, it is plausible that the pangolin immune system evolved differently than in other mammals. For example, genes such as IFNE, a unique interferon exclusively expressed in skin epithelial cells and inner mucosa-protected tisssues (e.g., lung, intestines, and reproductive tissues), establish a first line of defense against pathogens in other placental mammals (Day et al. 2008; Ponten et al. 2008; Xi et al. 2012; Fung et al. 2013; Demers et al. 2014; Uhlen et al. 2015). Interferons (IFNs) are a cluster of highly conserved gene families that encode for cytokines expressed by host cells for communication between cells, leading to the activation of the immune system in the presence of pathogens (De Andrea et al. 2002; Fensterl and Sen 2009). Strikingly, the single copy intronless IFNE gene is pseudogenized in both pangolin species (confirmed by Sanger sequencing; sample size = 8) but is intact in 71 other mammalian species (Fig. 2A). We found an insertion from positions 195 to position 219 and a point mutation at position c235C > T causing a premature stop codon. Other frameshift deletions at positions 264 and 540–546 would cause the loss of function. Although an alternate splicing pattern might avoid the premature stop codon and preserve the protein, high-throughput RNA sequencing of lung, cerebellum, cerebrum, and skin transcriptomes failed to detect IFNE expression (fragments per kilobase million [FPKM] = 0.0), although it is expressed in these organs in other mammals (Uhlen et al. 2005; Demers et al. 2014). We were also unable to identify any sequencing reads that mapped to the remainder of the gene after the premature stop codon, which is consistent with the hypothesis of a loss of function even if the stop codon were rescued. Interestingly, many putative functional domains or key signatures are deleted in the pangolin IFNE protein due to the premature stop codon (Fig. 2B). This suggests that the IFNE protein is unlikely to function properly if expressed, indicating that resistance to infection may be impacted in pangolins. We also examined other IFN families (IFNB, IFNK, IFNA, IFNG, and IFNL) and found that all families have intact gene copies in both pangolin species. It should be noted that the IFNA family consists of a cluster of 13 functional genes in human. Although we found that several intact genes from the IFN family in both pangolin species were retained, the number of IFNA genes is relatively low compared to human, suggesting diminished functions of this gene family in pangolins.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Multiple sequence alignment of all mammalian IFNE genes. Seventy-three mammalian species have available IFNE sequences, used for alignment. (A) Nucleotide sequence alignment of IFNE genes across different mammalian species. Blue = insertion, green = deletion, pink = stop codon. Protein sequence alignment of IFNE genes across different mammalian species is also shown. We identified an insertion (blue) in IFNE starting from nucleotide position 195 but not in other mammalian species, indicating that this insertion is specific to pangolins and possibly a marker to differentiate pangolins and other mammalian species. A premature stop codon or frameshift (orange) in the IFNE gene was consistently detected in both pangolin genomes. These mutations were validated by Sanger sequencing in eight Malayan pangolins. Protein sequences highlighted in orange are the frameshift mutations and premature stop codon. (B) Comparison between pangolin and human (reference) IFNE protein sequences. Pangolins have a short putative protein sequence because of a premature stop codon. Predicted functional domains and signatures in the human IFNE gene are represented by colored boxes. Yellow and pink boxes represent the predicted binding residues to IFNAR2 and IFNAR1, respectively, and which are bounded by interferons, which we obtained from a previosuly published paper (Fung et al. 2013). IFabd (SM00076) is a conserved functional domain in known interferons. The main conserved structural feature of interferons is a disulfide bond. INTERFERONAB (PR00266) is a three-element fingerprint that provides a signature for alpha, beta, and omega interferons. The elements 1 and 3 contain Cys residues involved in disulfide bond formation.

To test whether IFNE was also pseudogenized in the African pangolins, we sequenced the gene in all African species using the primers that we designed based on the genome sequence of Malayan pangolin. Our data showed that the IFNE gene was also pseudogenized in all African species and shared ancestral mutations with the Asian pangolin species (Supplemental Fig. S7.4). Therefore, we suggest that the pseudogenization of IFNE in pangolins most likely occurred before the divergence of the African and Asian pangolins, or at least 19–26.9 MYA, based on prior estimates (Bininda-Emonds et al. 2007; Fritz et al. 2009; Meredith et al. 2011).

Comparative analyses among pangolins and mammals

To further explore genetic differences between pangolins and their closest relatives, we compared the overlap of gene families among the two pangolin species and dog, cat, and giant panda. Our analysis identified 8325 ancestral gene families common to all five species (Fig. 3A). Interestingly, many gene families were unique to either Malayan pangolin (4958) or Chinese pangolin (3465), confirming the high level of divergence within pangolins, likely resulting from evolutionary adaptations to different ecological environments. We also found 1152 pangolin-specific gene families, of which 61% are expressed in at least one of the eight organs assessed, further suggesting that at least these genes are not assembly or annotation artifacts (Fig. 3B). Functional enrichment analyses suggest that the pangolin-specific genes are significantly overrepresented in signal transduction (496 genes; GO:0007165; FDR P-value = 2.53 × 10⁻⁹), neurological system processes (262 genes; GO:0050877; FDR P-value = 2.51 × 10⁻²⁰), cytoskeleton organization (99 genes; GO:0007010; FDR P-value = 1.53 × 10⁻³), and cell junction organization (28 genes; GO:0034330; FDR P-value = 1.32 × 10⁻²) (Fig. 3C). The large number of pangolin-specific genes significantly enriched in signal transduction and neurological system processes suggests that pangolins evolved complex and unique signal transduction and neurological system networks. The large number of cytoskeleton organization-associated genes may be linked with the assembly, formation, and maintenance of the cytoskeletal elements involved in cell shape, structural integrity, and motility of the pangolin scales.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Comparative analysis between pangolins and mammals. (A) Venn diagram showing the unique and shared gene families among pangolins and their closest relatives (cat, dog, and giant panda). (B) Heat map showing the expression level of pangolin-specific genes across different pangolin organs, represented by FPKM values. Any FPKM values >5 were set to 5 in the heat map for visualization purposes. Only genes expressed (FPKM ≥ 0.3) in at least one organ were shown. (C) GO enrichment analysis of 1152 pangolin-specific genes. Significantly enriched GO terms are shown for the categories of cellular compartment (blue), molecular function (yellow), and biological process (light orange). (D) Phylogenetic tree and gene family expansion and contraction. Expanded gene families are indicated in blue, whereas contracted gene families are indicated in red. The proportion of expanded and contracted gene families is also shown in pie charts. (E) Phylogenetic tree showing significant contraction of the interferon gene family in pangolins.

Gene family analysis identified 147 families that are significantly expanded and 18 that are contracted in the pangolin lineage, indicating that gene expansion and contraction events played a major role in the functional diversification of pangolins (Fig. 3D; Supplemental Information 8.0). Notably, pangolins have a highly reduced number of interferon genes, which have a major role in responding to infections, inflammation, and healing of skin (Xi et al. 2012; Fung et al. 2013). Of the 10 genes in the interferon family, only three were found in Malayan pangolin and two in Chinese pangolin (Fig. 3E). The heat shock gene family also contracted significantly, possibly contributing to the sensitivity of pangolins to stress (Hua et al. 2015). However, we also document gene expansion in several important families, including ribosomal genes (17 families), olfactory receptor (OR) genes (six families), cathepsin genes (one family), and septin genes (one family), which were recently shown to deter microbial pathogens, preventing them from invading other cells (Supplemental Information 8.0; Mascarelli 2011). The significantly expanded OR gene families suggest that pangolins have an enhanced sense of smell, possibly helping locate prey and counterbalancing poor vision (Supplemental Fig. S8.2). Strikingly, functional enrichment analysis of the 147 expanded gene families reveal these genes are significantly overrepresented in neurological system processes (including neuromuscular processes and sensory perception) (207 genes; GO:0050877; FDR P-value = 6.43 × 10⁻¹⁴), which is consistent with our finding in the functional analysis of pangolin-specific genes (Supplemental Fig. S8.3). The gene families associated with symbiosis and host-parasite interactions were also overrepresented (82 genes; GO:0044403; FDR P-value = 2.67 × 10⁻¹³).

Positive selection analysis

To identify signatures of natural selection, we used a set of 8250 protein-coding orthologs shared among the dog, cat, panda, cow, horse, mouse, human, and megabat genomes (Supplemental Information 9.0). Branch site tests identified evidence of positive selection along the pangolin lineage in 427 genes that had significant signals (P < 0.05 adjusted) (Supplemental Table S9.1).

We hypothesize that the reduced IFN-mediated immunity from the loss of IFNE and the contraction of the interferon gene family in pangolins imposed strong selective pressure on immunity-related genes. We identified a large proportion of genes under selection in the pangolin lineage that involve a wide range of immunity-related pathways, including hematopoietic cell lineage, cytosolic DNA-sensing pathway, complement and coagulation cascades, cytokine–cytokine receptor interaction, and the phagosome pathway (Supplemental Table S9.2). For instance, in the hematopoietic cell lineage, the colony stimulating factor 3 receptor (CSF3R) gene encodes a transmembrane receptor for a cytokine (granulocyte colony stimulating factor 3) (Fig. 4A). CSF3R proteins are present on precursor cells in the bone marrow, playing important roles in the proliferation and differentiation into mature neutrophilic granulocytes and macrophages which are essential for combating infection (Liongue and Ward 2014). As predicted by Protein Variation Effect Analyzer (PROVEAN), CSF3R protein contains a critical pangolin-specific mutation at position 247(G→H) in the functionally relevant Fibronection Type III domain that likely affects its protein function in pangolins (Fig. 4C). Another gene, integrin subunit alpha M (ITGAM), mediates leukocyte activation and migration, phagocytosis, and neutrophil apoptosis. It has been reported that ITGAM can migrate neutrophils across intestinal epithelium to maintain intestinal homeostasis and eliminate pathogens that have translocated across the single layer of mucosal epithelial cells (Parkos et al. 1991). We identified two pangolin-specific amino acid changes in the functional Integrin alpha-2 domain of ITGAM which likely affect protein function (Fig. 4C). In the cytosolic DNA-sensing pathway, transmembrane protein 173 (TMEM173), also known as Stimulator of Interferon Genes, plays an important role in innate immunity by sensing cytosolic foreign DNA and inducing type I interferon production when cells are infected with pathogens (Fig. 4B; Ran et al. 2014). We identified two critical amino acid changes in pangolins at positions 216(D→G) and 217(P→L) in the functionally relevant region, suggesting that they may have a functional impact on TMEM173 (Fig. 4C).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Evolution of immunity-related pathways and scale-related genes in the pangolin ancestor. Genes under positive selection in the pangolin lineage in hematopoietic cell lineage (A) and cytosolic DNA-sensing pathway (B) are highlighted in red. (C) Several critical pangolin-specific mutations were identified in the functionally relevant Fibronection Type III signatures (INTERPRO ID = IPR003961) of CSF3R, the Integrin alpha-2 signatures (IPR013649) of ITGAM, and the Stimulator of Interferon Genes Protein region (INTERPRO ID = IPR029158) of TMEM173. For the CSF3R, ITGAM, and TMEM173, P-values of the detection of positive selection were 1.58 × 10⁻², 3.75 × 10⁻², and 2.92 × 10⁻⁵, respectively. (D) Several critical pangolin-specific amino acid changes were detected in hair-/scale-related keratin proteins, KRT36 and KRT75, which are located in functionally relevant regions that may affect protein functions. For the KRT36 and KRT75, P-values of the detection of positive selection were 2.24 × 10⁻² and 1.33 × 10⁻⁵, respectively. We also examined whether African species (M. tricuspis, M. tetradactyla, and M. temminckii) have the critical pangolin-specific amino acid changes that we observed in the Asian pangolins. Circles at the bottom indicate our preliminary results: Red circle = all African species that we examined have the identical amino acid change, yellow circle = not all African species that we examined have the identical amino acid change (but they all have amino changes/deletion that likely affect protein function as predicted by PROVEAN), brown circle = all African species that we examined have the same amino acids like the human reference sequence, white circle = data not available. The alignment results are shown in Supplemental Figure S7.5.

We further postulate that genetic changes during the evolution of hair-derived pangolin scales likely involved genes related to hair formation in general, and particularly keratins, which are essential components of scales and hair (Meyer et al. 2013). Among the positively selected candidate genes associated with hair formation are KRT36, KRT75, KRT82, and KRTAP3-1. For instance, type II keratin 75 (KRT75), a hair follicle–specific keratin, has an essential role in hair and nail integrity (Chen et al. 2008). Type I keratin 36 (KRT36), a hair or “hard” sulfur-rich keratin, is mainly responsible for the extraordinarily high degree of filamentous cross-linking by specialized keratin-associated proteins (Moll et al. 2008). Many critical pangolin-specific amino acid changes were identified in keratin genes in functionally relevant domains (Fig. 4D). These critical amino acid changes are also present in both KRT36 and KRT75 genes in the distantly related African pangolin species that we examined, suggesting that they potentially contributed to the development of pangolin scales (Fig. 4D; Supplemental Fig. S7.5). We also identified genes related to energy storage and metabolism (12 genes) and mitochondrial metabolism (11 genes), perhaps an adaptive response to reduce their metabolic rate given their large body size but low energy diet of ants and termites (da Fonseca et al. 2008). Other candidate genes are associated with the development of the robust muscular system of pangolins, which assists its fossorial lifestyle (10 genes), and the nervous system (10 genes) (Supplemental Table S9.2).

Furthermore, we explored whether there was evidence that genes under selection were associated with any response to diseases, particularly those related with infection or skin/mucosal organs. We identified genes associated with inflammation (18), cancer or viral infection (25), bacterial infection (10), pneumonia and gastrointestinal disease (21), and skin diseases (15) (Supplemental Table S9.3). Interestingly, one gene, lactotransferrin (LTF), encodes a multifunctional immune protein found at mucosal surfaces, providing the first line of defense to the host against inflammation and infection (Ward et al. 2002). It has been shown that LTF has antimicrobial activities and induces both systematic and mucosal immune responses, for example, against lung and gut-related systemic infections (Debbabi et al. 1998). Two pangolin-specific amino acid changes were detected in the LTF protein sequence, potentially having an impact on LTF function (Supplemental Fig. S9.1). Together, evidence of selection on these genes provides possible links with the unique phenotypes and physiological requirements of pangolins. Importantly, the immunity- or inflammation/infection-related genes provide avenues to further explore the evolution of mammalian immunity.

Pangolins have a noncanonical repeat layout

Transposable elements accounted for approximately 28.9%–30.0% of the pangolin genomes, lower than in other Carnivora such as the giant panda (35.3%), tiger (36.7%), cat (37.2%), and dog (36.2%), suggesting divergence in genome architecture between pangolins and other placental mammals (Supplemental Fig. S10.1a). Pangolins also have relatively low proportions of SINE repeats (2.6%–2.8% vs. 8.2%–11%), mainly due to a relatively smaller number of tRNA SINEs (Supplemental Fig. S10.1a,b). A search using a de novo constructed repeat library supports the comparatively low proportion of this repeat family (Supplemental Table S10.2). Furthermore, both pangolin species have relatively low proportions of intronic repeats (SINEs, LINEs, and LTRs) (Supplemental Fig. S10.1c).

Discussion and conclusion

Our analyses have revealed that, similar to observations in other lineages (Wang et al. 2006; Meredith et al. 2014), gene loss through pseudogenization has had a significant impact on the evolution and diversification of pangolin genomes and their biology. We detected and validated pseudogenes in pathways related to multiple morphological or physiological functions including dentition, vision, and immunity. The pseudogenization of the IFNE gene in pangolins, otherwise functional in all other mammal genomes surveyed so far, is especially noteworthy. In mice, a deficiency of IFNE in epithelial cells of the female reproductive tract (FRT) can increase susceptibility to microbial infection (Fung et al. 2013). Therefore, the pseudogenization of this interferon gene in pangolins suggests that innate immunity may be compromised, resulting in an increased susceptibility to infection, particularly in the skin and mucosa-protected organs. Since pangolins may have intrinsically weak mucosal immunity, Burkholderia spp. might easily infect the lung by penetrating lung epithelial cells and the brain by penetrating nasal mucosa and migrating to the brain through olfactory nerves, which has previously been shown with other Burkholderia (Owen et al. 2009; Sim et al. 2009; Taylor et al. 2010; Dando et al. 2014; St John et al. 2014). Furthermore, captive pangolins are prone to frequently fatal gastrointestinal disease, pneumonia, and skin maladies (Clark et al. 2008; Hua et al. 2015). Pangolins are notoriously difficult to maintain in captivity, and the stress of captivity and/or poor husbandry might render pangolins even more vulnerable to infections or diseases by suppressing immune responses. Whether or not these afflict free-ranging pangolins where pathogen exposure is rampant remains unknown.

Our study raises the fundamental evolutionary question of how pangolin adaptations occurred given their seemingly increased risk to infection if their skin and mucosal surfaces are constantly being exposed to pathogens. Our analyses suggest several evolutionary changes in genes that may respond to or interact with pathogens, including immunity-related genes/pathways, the significant expansion of the septin gene family, and the enrichment of the genes in symbiosis. These changes may enhance the immunity system of pangolins, although we do not know the relative efficiency of these adaptations compared to interferon-mediated immunity. We propose that pangolin scales may be an important morphological innovation to compensate for the decrease in immunity normally provided by the skin. These hard and overlapping scales may act as defensive armor to protect pangolins against injuries (or stress) which would make pangolins even more vulnerable to infection or the invasion of pathogens. Moreover, pangolins curl into near impregnable balls, covering their scaleless abdomen using their well-developed neuromuscular system, during sleep or when threatened, which would also support the hypothesis that these adaptations serve to protect pangolins from skin injuries. It is noteworthy that besides the antimicrobial effects, plasmacytoid dendritic cell (pDC)-produced type 1 IFNA/B helps heal wounds through re-epithelization of injured skin (Gregorio et al. 2010). It is unknown whether keratinocyte-produced type-I IFNE is also involved in skin wound healing (e.g., recruiting pDCs to the wound). Pangolins have a marked contraction of interferon genes that likely facilitates wound healing. Future studies are needed to test the relationship between IFNE, skin healing, and scale development.

Our analyses showed that pangolin-specific genes were significantly overrepresented in the GO categories governing cytoskeleton organization, neurological system process, and signal transduction. Further, certain genes involved in the muscular and nervous systems showed evidence of positive selection. These could be linked to the unique adaptations and behaviors of pangolins, particularly those related with their highly sophisticated musculoskeletal system (Supplemental Information 12).

The Chinese pangolin had a significantly larger predicted genome size but significantly fewer genes than the Malayan pangolin, despite the lower sequencing coverage in the former species. We found segmental duplications were identified in a higher proportion in the Chinese pangolin compared to the Malayan pangolin, suggesting that the segmental duplications partly contribute to the larger genome size of the Chinese pangolin (Supplemental Information 3.4). It should be noted that a previous comparative chromosome map between the Malayan (2n = 38) and Chinese (2n = 36–42) pangolins revealed many chromosomal rearrangements between the two species (Nie et al. 2009). Besides that, we did not observe significant recent repeat family activity unique to Chinese pangolins, suggesting that the genome expansion is unlikely directly accounted by repeats (Supplemental Fig. S10.2). However, we cannot rule out the possibility that the recent segmental duplications and repeats were underestimated (e.g., some might collapse together due to high sequence similarity during assembly since we sequenced the pangolin genomes using the short-read Illumina technology), which would complicate the analyses.

Our data also showed that the number of gene families unique to each of the pangolin species is 3–4 times greater than the number of shared gene families. In this analysis, we faced the common problem of almost all draft genomes—the potential difference between the real and available gene set. We can expect a number of misannotated genes because of fragmented assemblies, and this can partly inflate the number of genes for each species. However, it should not affect subsequent analyses because we had a refined pangolin-specific gene set by orthology clustering as described in the Methods section. It is also worth noting here that the discrepancy in the unique gene family numbers indeed can be explained by biological reasons, such as a high rate of molecular evolution. It was previously shown that the chromosome complements of Malayan pangolin and Chinese pangolin differ by seven Robertsonian rearrangements (three centric fissions and four centric fusions), which could also affect the number of unique gene families for each species after their divergence from the common ancestor (Nie et al. 2009).

In conclusion, this study provides new insights into the biology, evolution, and diversification of pangolins. The annotated reference genomes will be invaluable for future studies addressing issues of species conservation and gene function and may provide a new natural animal model for understanding mammalian immunity.