The gut microbiome of the Sunda pangolin (Manis javanica) reveals its adaptation to specialized myrmecophagy

Sample collection

Fecal samples were collected from nine Sunda pangolins kept in captivity at the Pangolin Research Base for Artificial Rescue and Conservation Breeding of South China Normal University (PRB-SCNU). The sample details are shown in Table 1. Their artificial diet consists of mealworms (58.94%), silkworms (20%), dry ants (15%), yeast powder (2%), vitamin complex (0.06%), crushed stone or coarse sand (4%), and vinegar (0.85 mL/100 g) (Zou, 2016). We performed comparative analyses of myrmecophages, carnivores, herbivores, folivores, insectivores, and omnivores using 47 samples that were originally sequenced by Muegge et al. (2011) and Delsuc et al. (2014) (Table S1). We first used a sterile swab to obtain a small amount of fecal material and transferred it to a vial containing a lysis and stabilization buffer that preserved the DNA for transport at ambient temperatures (Almonacid et al., 2017). The samples were then frozen in liquid nitrogen and stored at −80 °C until DNA extraction.

Data acquisition

Bacterial genomic DNA was extracted from the feces samples using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany). The concentration of each sample was tested with a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA). The hypervariable V3 and V4 regions of the 16S rRNA gene were amplified with target-specific primers (338F and 806R) that specifically amplify bacterial and archaeal 16S genes (Caporaso et al., 2012; Allali et al., 2017). Three replicate PCRs were performed for each DNA sample. After amplification, the PCR products from the same sample were mixed, purified using the Agencourt AMPure XP system (Beckman Coulter Genomics, Brea, CA, USA), and examined by 2% agarose gel electrophoresis. DNA was extracted from the gels using the AxyPrep DNA Gel Extraction Kit (Axygen, Tewksbury, MA, USA) and examined by 2% agarose gel electrophoresis. Illumina sequencing adapters and dual-index barcodes [index 1(i7) and index 2(i5)] (Illumina, San Diego, CA, USA) were added to each amplicon in a subsequent limited cycle PCR. The final libraries were repurified using the Agencourt AMPure XP system, quantified with a Quant-iT™ PicoGreen^® kit (Invitrogen/Life Technologies, Green Islands, NY, USA), and normalized before pooling. The DNA library pool was denatured with NaOH and diluted with hybridization buffer. The libraries were validated using an Agilent 2100 biological analyzer (Agilent Technologies, Palo Alto, CA, USA) and quantified with a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA). Finally, the libraries were multiplexed and loaded on an Illumina MiSeq instrument according to the manufacturer’s instructions (Illumina). Sequencing was performed using a 2 ×250 paired-end (PE) configuration and MiSeq Control Software embedded in the MiSeq instrument did the image analysis and base calling. Our data were submitted to the NCBI sequence read archive under accession number PRJNA558047.

Sequence processing and analysis

Clean reads were obtained after paired-end reads assembly and quality control from raw data. Initially, the PE reads were merged using FLASH (Magoc & Salzberg, 2011) and quality filtering was performed using Trimmomatic (ver. 0.39) (Bolger, Lohse & Usadel, 2014). Potential chimeras were moved using the UCHIME algorithm (Rognes et al., 2016). Then, the clean reads were aligned and clustered into operational taxonomic units (OTUs) with a sequence similarity threshold of 97% using Usearch (ver. 10) (Edgar et al., 2011). Chimeras and singletons were removed in the clustering process. The Ribosomal Database Project (RDP) (ver. 2.2) classifier was used to assign taxonomic categories to all OTUs at a confidence threshold of 0.7 (Wang et al., 2007). The RDP classifier uses the Silva 119 database, which has taxonomic categories predicted to the species level (Quast et al., 2013).

Sequences were rarefied before calculating the alpha and beta diversity statistics. We calculated the OTU abundance and Ace, Chao1, Simpson, and Shannon indices using the Mothur software package (ver. 1.30.1) to measure alpha diversity based on the number of different data categories and their respective abundances in the dataset (Schloss, Gevers & Westcott, 2011). All statistical analyses were conducted in R 3.5.1 (R Core Team, 2015). We calculated Good’s coverage (Good, 1952) and plotted rarefaction and rank abundance curves using Mothur to determine whether the sequencing effort was sufficient to describe bacterial communities. We compared the beta diversity to identify differences in the microbial communities between each sample. Principal coordinates analyses (PCoA) and non-metric multidimensional scaling (NMDS) were performed using the Bray–Curtis distance (Bray & Curtis, 1957). To explore whether pangolins underwent convergent evolution with other groups, we performed PCoA on the UniFrac distance (Jiang et al., 2013) and constructed a neighbor-net network based on the UniFrac distance (Huson & Bryant, 2006). Permutational MANOVA (PERMANOVA) was performed to identify significant differences between these groups (“Adonis” in the R package “vegan”).

Analysis of the sequencing data

We generated 16S rRNA gene sequences from nine fecal samples from nine Sunda pangolins. All samples were sequenced. Initially, we generated 735,783 PE reads. After all quality filtering steps, we had 704,251 clean reads with an average read length of 412.10 bp and an average of 78,250 sequences per sample (Table S2). Figure 1 shows the distribution of the lengths of the clean sequences. Using a minimum identity of 97% as the threshold for any sequence pair, we identified 536 OTUs (Table S3). The rarefaction curve for the OTUs flattened gradually, almost plateauing, indicating that the number of OTUs for each sample was sufficient and reasonable (Fig. 2A). The rank abundance curves indicated the evenness and abundance of species in fecal samples horizontally and vertically (Fig. 2B). Good’s coverage was nearly 99.9% for all samples, showing that most bacterial species in the samples were detected (Table 2). We generated a gene clustering heatmap using the 50 most abundant OTUs (Fig. 3).

Composition of Sunda pangolin gut microbiota

The detected bacteria were classified into 14 phyla, 24 classes, 48 orders, 97 families, and 271 genera. The gut microbiomes were dominated by bacteria in the phylum Firmicutes (73.71%), followed by Proteobacteria (18.42%), Actinobacteria (3.44%), and Bacteroidetes (0.51%) (Fig. 4A). It was rich in the families Clostridiaceae (31.97%), Peptostreptococcaceae (21.68%), Enterobacteriaceae (15.34%), Lachnospiraceae (7.58%), Streptococcaceae (4.96%), Planococcaceae (1.47%), Ruminococcaceae (1.43%), Lactobacillaceae (1.39%), and Nocardiaceae (1.2%), followed by Erysipelotrichaceae, Enterococcaceae, Bacillaceae, Xanthomonadaceae, Rhodobacteraceae, Corynebacteriaceae, and Brucellaceae. The principal genera were Clostridium_sensu_stricto_1 (29.42%), Escherichia-Shigella (13.57%), Terrisporobacter (10.39%), Lactococcus (4.92%), Epulopiscium (4.41%), Clostridium_sensu_stricto_13 (1.78%), and Morganella (1.73%), followed by Lactobacillus, Sporosarcina, Cellulosilyticum, Staphylococcus, Bacteroides, and Bacillus (Fig. 4B). We also identified several genera that are likely involved in chitin digestion. Lactococcus (Vaaje-Kolstad et al., 2009) and Bacteroides (Borrelli et al., 2017) were abundant (Fig. 5), followed by Bacillus (Gooday, 1990) and Staphylococcus (Wadström, 1971).

Table 2 shows the alpha-diversity indices (including Chao1, Ace, Shannon, and Simpson). These indices reflected the number of OTUs and the diversity in the nine samples was high. Two-dimensional PCoA and NMDS based on the Bray–Curtis distance suggested that seven of the nine Sunda pangolin gut microbiomes cluster, which means that the results are reliable and have statistical significance (Fig. 6).

Comparison of the Sunda pangolin with mammals with different diets

To explore whether pangolins underwent convergent evolution with other groups, we compared them to mammals with various diets, including myrmecophages, carnivores, herbivores, folivores, insectivores, and omnivores. Using PCoA and the neighbor-net network of weight UniFrac distance, these data confirmed that Sunda pangolins form a group (PERMANOVA: pangolins vs. other diets, p<0.01), which was close to myrmecophages (aardvarks, anteaters, aardwolves, and armadillos), and distant from carnivores and herbivores (Figs. 7A, 8A). In addition, beta diversity analysis of 56 samples by PCoA of unweighted UniFrac distance revealed that the gut microbiomes of Sunda pangolins formed a distinct cluster clearly separated from other diets (PERMANOVA: pangolins vs. other diets, p<0.001) (Fig. 7B). Moreover, the network reconstructed from the unweighted UniFrac distances also showed that Sunda pangolins differed from other diets (Fig. 8B). These findings all demonstrate that Sunda pangolin gut microbiomes are close to those of myrmecophages and carnivores, and distant from herbivores.