Genome sequence analysis of Malayan pangolin (Manis javanica) forensic samples reveals the presence of Paraburkholderia fungorum sequences

Ethics statement

Veterinary officers conducted all procedures involving dissecting animals under the oversight of experts from the Department of Wildlife and National Parks (DWNP) Peninsular Malaysia (PERHILITAN, Malaysia), following internationally recognized guidelines and approved by the University of Malaya Institutional Animal Care and Use Committee (reference number DRTU/11/10/2013/RH (R)).

Biological samples

In 2012, more than 40 Malayan pangolins (Manis javanica) were seized in an anti-smuggling operation by PERHILITAN in Malaysia. One of the female pangolins (“UM3”) was alive and pregnant, and weighed 2.73 kg at the time of confiscation. It was euthanized by the professional and experienced veterinarians at PERHILITAN using Dolethal® (Pentobarbitone sodium 200 mg/L; dosage: 1 mL/kg body weight) for animal welfare reasons. Only after the veterinary surgeon had pronounced the animal dead by checking for vital signs such as breathing and pulse, which ceased within 5 min of administering the drug, several organs including the cerebellum, cerebrum, lungs, thymus, liver, blood, heart, spleen and skin of the female and the fetus’s cerebrum, cerebellum, intestine, kidney, cord blood, liver, lung, and gastrocnemius muscle were harvested and used in this study. The samples were immediately flash-frozen using liquid nitrogen and subsequently stored at −80 °C. For DNA extraction, minimal thawing was performed to obtain the tissue samples. Blood was drawn from seven additional live adults from the same seizure. Before conducting blood sampling, an anaesthetic mixture of ketamine-xylazine 1:1 (dosage: 0.5–1 mL per pangolin) or Zoletil 100 (3–4 mg/kg) was injected intramuscularly (IM) into the adult pangolins (276T, 2T9, 12T, 2T2, UM1, UM2, and UM3) to minimize their intervention. Once the muscles had relaxed, blood sampling was carried out via the coccygeal vein located at the tail, using a 5 mL syringe and 21G needle.

Genomic DNA extraction and library preparation

The fetal gastrocnemius muscle tissue was used for genomic DNA extraction using a Qiagen Genomic Tips 20/G kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol.

Genomic DNA libraries were prepared with a fragment length of approximately 300 bp and sequenced using Illumina HiSeq 2000 following the manufacturer’s sequencing protocol.

Discovery of bacterial sequences in the pangolin genome

The genome of the female pangolin was sequenced using the Illumina HiSeq 2000 platform (Choo et al., 2016). The tissue-specific genome assemblies were generated by CLC Assembly Cell using sequencing data from pangolin brain (cerebrum and cerebellum), liver, and lung samples. The sequencing reads were searched against a bacterial nucleotide sequence database using BLASTN (Mount, 2007). We screened the bacterial identity using two criteria: 90% sequence identity and 90% sequence coverage, and 97% sequence identity and 97% sequence coverage.

Average nucleotide identity (ANI) and average amino acid identity (AAI) analyses

The average nucleotide identity (ANI) values between bacterial species were calculated using previously described methods (Goris et al., 2007; Rodriguez-R & Konstantinidis, 2014; Ang et al., 2016; Tan et al., 2016). We used two-way BLAST and only used the forward and reverse-matched orthologs in the calculations. For robustness, the BLAST hits were filtered for at least 50% identity at the nucleotide and amino acid level, and a sequence coverage of at least 70%.

The protein sequences of 18 genomes belonging to the Paraburkholderia and Burkholderia genera were annotated using Rapid Annotation Search Tool (RAST) (Aziz et al., 2008). The RAST-predicted protein sequences for each assembly were retrieved, and average amino acid identity (AAI) values were calculated using the AAI calculator (Goris et al., 2007; Rodriguez-R & Konstantinidis, 2014).

PCR assays

To further validate the presence of bacterial sequences in the pangolin, the frozen tissue samples of pangolin UM3 and her fetus were examined. Genomic DNA was extracted from nine adult tissues (cerebrum, cerebellum, liver, lungs, heart, spleen, thymus, skin, and blood) and eight fetal tissues (cerebrum, cerebellum, intestine, kidney, cord blood, liver, lung, and gastrocnemius muscle), and these were screened using polymerase chain reaction (PCR) assays. Three different target genomic regions that showed top hits to the bacteria identified from the previous BLASTN results were selected to design and synthesize novel PCR primers. We used bacterial universal 16S primers and three in-house designed primer sets (Table S1), targeting bacterial 16S rRNA, Burkholderia-specific transposase genomic region, OI25_7129 hypothetical protein genomic region, and the P. fungorum-specific DNA polymerase genomic region, respectively.

All PCR assays were performed using a total reaction volume of 50 µL containing 160 ng purified organ gDNA, 0.3 mol of each primer, deoxynucleotide triphosphates (dNTP, 400 µM each), 1.0 U Taq DNA polymerase and a supplied buffer. The PCR was performed as follows: one cycle (94 °C for 2 min) for initial denaturation; 35 cycles (98 °C for 10 s; 68 °C for 3 min) for annealing and DNA amplification. After completion of PCR, we visualized the PCR products on a 1% TAE agarose gel at 100 V for 70 min. The PCR products were purified by GeneJET PCR Purification Kit and directly sequenced with the same primers using BigDye^© Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA, USA) for validation.

Tissue preparation and histological staining

We examined the histology of the adult pangolin’s cerebellum and lungs. Each of the thawed organs was excised into two sets of smaller tissue pieces and fixed in 10% buffered formalin at 12 °C for a week, followed by embedding in paraffin wax to produce paraffin blocks. For histology, the tissue blocks were sectioned on a rotary microtome (Leica RM2235; Leica Biosystems, Wetzlar, Germany) using a 3 μm blade. To prevent cross contamination, the blades were cleaned with 99% ethanol between sections. Subsequently, the slices were dewaxed using graded alcohol. Tissue slices were separately counterstained using hematoxilin/eosin (HE; Sigma, St. Louis, MO, USA) for tissue abnormality such as inflammation, and Brown-Hopps Gram stains for bacterial presence, as described by Brown & Hopps (1973). Slices were examined under a light microscope with a Leica DF300 camera.

Assembly of P. fungorum genomes

To further analyse the genomes of P. fungorum, we assembled genomes using three different strategies: (i) mapping reads from UM3’s cerebrum and cerebellum whole-genome data onto P. fungorum reference genome ATCC BAA-463 (accession number: CP010024–CP010027) and assembling them into contigs; (ii) sequencing the DNA extracted from UM3’s fetal muscle and mapping these reads onto the P. fungorum reference genome and assembling them into contigs, and (iii) mapping reads from UM3’s cerebrum, cerebellum, and fetal muscle whole-genome data onto the P. fungorum reference genome and assembling them into contigs.

The core-genome single nucleotide polymorphisms (SNPs) were employed in constructing a sturdy phylogenetic tree to determine the taxonomic classification of these genomes. Core-genome SNPs are found in the core genome of a species or a group of closely related strains (Vignal et al., 2002) and are commonly used in phylogenetic studies to infer evolutionary relationships among bacterial populations. The use of core-genome SNPs has been shown to be a powerful tool for phylogenetic analyses, as they are less likely to be subject to homoplasy (convergence or reversal of nucleotide changes) compared to other markers. Seventeen Burkholderia and Paraburkholderia whole-genome sequences were retrieved from NCBI (http://www.ncbi.nlm.nih.gov) (Pruitt, Tatusova & Maglott, 2007; Table S2), with the sequences used being similar to those used by Tan et al. (2020). The newly generated and reference sequences were uploaded to the PanSeq server to identify the core-genome SNPs in common genomic regions (Laing et al., 2010), and the extracted core SNPs were subsequently merged into a continuous sequence for each genome. Recent studies used conserved sequence indels (CSI) to study the relationships between species of Burkholderia and Paraburkholderia (Sawana, Adeolu & Gupta, 2014). Thus, CSI were used to further verify the taxonomic classification in this study. Conserved sequence indels were identified using the protein sequence of the assembled genome (iii) and 17 closely related species that were detected by ProteinOrtho (Lechner et al., 2011). The resultant 27 CSI were aligned using ClustalW (Thompson, Gibson & Higgins, 2003). The phylogenic trees using sequences of core genome SNPs and CSI were reconstructed using MEGA-X (Kumar et al., 2018). Neighbour-joining trees were inferred using the Kimura’s two parameter model and nodal support was estimated using 1,000 non-parametric bootstrap replicates.

Presence of bacterial sequences in UM3

In our previous pangolin genome sequencing project, we sequenced the genomes of pangolin cerebellum, cerebrum and liver using the Illumina HiSeq2000 platform (Choo et al., 2016). During an in-silico bacterial sequence screening of the contig sequences of the tissue-specific assemblies using BLAST, we found many exogenous DNA sequences. The bacterial sequences were found in the assemblies of the cerebrum and cerebellum, but not in the liver assembly (Table S3). Specifically, in the assembled cerebral genome, there were 6,730 contigs mapped to bacterial genomes, where 6,635 of them (98.58%) had best matches with P. fungorum. Similarly, in the cerebellum-specific genome, 3,533 contigs mapped to bacterial genomes, among which 3,452 (97.7%) were from P. fungorum. These results indicate that the cerebrum and cerebellum tissues were predominantly colonised or infected by P. fungorum even though they should be sterile.

PCR screening and sanger sequencing across different tissues of UM3

Our results showed clear positive bands in the lung, cerebrum, cerebellum, and blood of UM3 (Fig. 1). A weak positive band was also visible in the liver whereas no clear bands were observed in the other tissue types (Fig. 1).

Histological examinations

To further confirm the presence of P. fungorum in UM3, the lung and cerebellum tissues were dissected and stained using Brown-Hopps Gram stains. Our staining revealed the presence of gram-negative and rod-shaped bacteria with a size of approximately 6–7 microns, supporting the notion that the lungs and cerebellum were invaded by P. fungorum (Figs. 2A, 2B, 2D, and 2E).

Histopathology screening was also performed using hematoxylin/eosin (H&E) staining to confirm P. fungorum infection in the lung and brain (cerebrum and cerebellum). The histological presentation of Paraburkholderia infection observed from lung tissues in other mammals is an abscess composed of cellular debris, numerous degenerate neutrophils, and macrophages that contain abundant intracytoplasmic basophilic material composed of rod-shaped bacteria (Glaros et al., 2015). However, our investigations showed no significant pathological signs in the dissected organs (Figs. 2C and 2F).

Presence of P. fungorum in other adult pangolins

To examine whether the presence of P. fungorum in the UM3 pangolin organs was an isolated case, we randomly selected and screened the blood of four live adult pangolins (26T, 2T9, 12T, and 2T2) that were seized in the same batch as UM3. We also tested the blood of two live adult pangolins (UM1 and UM2) seized in a separate operation. UM3 was included as a positive control. All samples were screened for the presence of P. fungorum using the same primer sets as were used for UM3 and her fetus (Table S1). Of the seven blood samples, four (2T9, 12T, 2T2, and UM3) showed positive PCR bands for P. fungorum, and the bacterial identity was confirmed by Sanger sequencing and phylogenetic analyses (see Choo et al., 2020; Fig. 4 and Table S4), again indicating that the pangolins’ blood was infected by P. fungorum (Fig. 3).

Paraburkholderia fungorum may have the capability to infect the fetus

Since UM3 was pregnant, we wondered whether its fetus was also infected by P. fungorum. To examine this, we harvested and screened tissues from its fetus (including cord blood, lungs, intestine, kidney, liver, and brain) by performing PCR with the primer sets. We found that the fetal gastrocnemius muscle showed clear positive PCR bands (Fig. 4), and the bacterial identity was confirmed by Sanger sequencing and phylogenetic analyses (see Choo et al., 2020; Fig. 4 and Table S4). No significant bands were observed in other tissues (cerebrum, cerebellum, kidney, lung, cord blood, intestine and liver; Fig. 4). Furthermore, the fetal muscle genome was also sequenced using the Illumina HiSeq 2000 platform with a 20X sequencing coverage (after removing the pangolin sequences) and we found a substantial amount of P. fungorum DNA sequences.

Assembly of P. fungorum genomes

The newly assembled genome sequence of P. fungorum (accession numbers: CP028829–CP028832) can be downloaded from the GenBank database. Our whole-genome results confirm that the sequences are most similar to the P. fungorum reference sequence ATCC BAA-463 (Figs. 5 and 6) which further validates the assignment of our sequences to P. fungorum.

To further validate this species assignment, we compared the P. fungorum genomes to other Burkholderia and Paraburkholderia species using ANI and AAI values. These comparisons indicated that the identified P. fungorum was closely related to the reference P. fungorum ATCC BAA-463, with an ANI value of 98.49% (Fig. 7). Other species had ANI values below the threshold of 97% used to define a species (Goris et al., 2007). The identified P. fungorum found in UM3’s brain and fetal muscle had almost identical AAI and ANI values, indicating that they were from the same source. The ANI, AAI, and core-genome SNP-based phylogenetic analyses provided consistent evidence for the presence of P. fungorum in UM3.