Geographically structured genomic diversity of non-human primate-infecting Treponema pallidum subsp. pertenue.

Many non-human primate species in sub-Saharan Africa are infected with Treponema pallidum subsp. pertenue, the bacterium causing yaws in humans. In humans, yaws is often characterized by lesions of the extremities and face, while T. pallidum subsp. pallidum causes venereal syphilis and is typically characterized by primary lesions on the genital, anal or oral mucosae. It remains unclear whether other T. pallidum subspecies found in humans also occur in non-human primates and how the genomic diversity of non-human primate T. pallidum subsp. pertenue lineages is distributed across hosts and space. We observed orofacial and genital lesions in sooty mangabeys (Cercocebus atys) in Taï National Park, Côte d'Ivoire and collected swabs and biopsies from symptomatic animals. We also collected non-human primate bones from 8 species in Taï National Park and 16 species from 11 other sites across sub-Saharan Africa. Samples were screened for T. pallidum DNA using polymerase chain reactions (PCRs) and we used in-solution hybridization capture to sequence T. pallidum genomes. We generated three nearly complete T. pallidum genomes from biopsies and swabs and detected treponemal DNA in bones of six non-human primate species in five countries, allowing us to reconstruct three partial genomes. Phylogenomic analyses revealed that both orofacial and genital lesions in sooty mangabeys from Taï National Park were caused by T. pallidum subsp. pertenue. We showed that T. pallidum subsp. pertenue has infected non-human primates in Taï National Park for at least 28 years and has been present in two non-human primate species that had not been described as T. pallidum subsp. pertenue hosts in this ecosystem, western chimpanzees (Pan troglodytes verus) and western red colobus (Piliocolobus badius), complementing clinical evidence that started accumulating in Taï National Park in 2014. More broadly, simian T. pallidum subsp. pertenue strains did not form monophyletic clades based on host species or the symptoms caused, but rather clustered based on geography. Geographical clustering of T. pallidum subsp. pertenue genomes might be compatible with cross-species transmission of T. pallidum subsp. pertenue within ecosystems or environmental exposure, leading to the acquisition of closely related strains. Finally, we found no evidence for mutations that confer antimicrobial resistance.

Data Summary

The authors confirm that all supporting analyses and external sequence data have been provided within the article or through supplementary data files. All new sequences generated in this study have been deposited with the National Center for Biotechnology Information (NCBI) under BioProject PRJNA588802.

We were able to reconstruct three new subsp. genomes from swabs/biopsies and three partial genomes from bone samples, adding to the limited number of subsp. genomes existing today. Phylogenomic analyses revealed that simian subsp. strains did not form monophyletic clades based on the host species or the types of clinical symptoms manifested, but rather clustered by geographical origin. This is compatible with cross-species transmission of subsp. within ecosystems where this disease occurs. Additionally, this study joins a growing body of evidence that human and wildlife bones are useful for tracing the history of treponemal diseases. Using bones, we confirmed that subsp. infected previously undescribed non-human primate hosts and has been present in Taï National Park for at least 28 years.

Introduction

Spirochete bacteria belonging to the species (TP) have affected humankind since at least the late 15th century [1] and cause a large global disease burden in humans [2, 3]. Three pathogenic subspecies are currently recognized that are morphologically similar, but genetically and epidemiologically distinguishable [4]. Clinically, pathogenic TP subspecies cause three distinct disease syndromes; yaws (subsp. pertenue; TPE), venereal syphilis (subsp. pallidum; TPA) and bejel (subsp. endemicum; TEN) [5, 6]. Though usually treatable with antibiotics, these treponematoses remain major public health threats across the globe [2-4]. For example, in 2012, nearly 8 million new cases of venereal syphilis were reported globally [7], and in 2018, an estimated 80 000 new yaws cases occurred in the 15 countries where the disease remains endemic [8]. Exact estimates of the number of bejel cases from the Sahel region and Arabian Peninsula where the disease is endemic are not available [9]. Efforts to reduce the prevalence of these treponematoses are underway, particularly for yaws, which an ongoing campaign aims to eradicate globally by 2030 [10]. Eradication efforts face a number of challenges, such as issues concerning the availability of diagnostics, the distribution of treatment options and resistance to available antibiotics.

One of the key questions for the potential success of eradication efforts and for understanding TPE evolution is the degree to which other animals are infected with these pathogens and whether cross-species transmission occurs. To date, only TPE has been shown to infect non-human primates (NHPs), as suggested by the finding of TPE DNA in lesion samples from a number of different NHP species [11-13]. The TPE strains infecting humans and NHPs are extremely similar. For example, the first genome sequenced from a TPE infecting a wild NHP – the strain Fribourg-Blanc obtained from a Guinea baboon (Papio papio: RefSeq ID: NC_021179.1) – exhibits a 99.97 % sequence identity with human-infecting strains [13]. There is no clear evidence for phylogenetic separation of NHP-infecting and human-infecting TPE strains, as strains do not form well-supported reciprocally monophyletic groups [12, 13]. Although this pattern might be suggestive of past transmission events between NHPs and humans, whether such events really happened remains uncertain due to a paucity of genomic data from both humans and wildlife. However, in an experimental setting, the Fribourg-Blanc strain induced classical yaws symptoms in humans [14] and human-infecting TPE strains were reported to elicit yaws-like symptoms in NHPs [15]. This suggests that molecular compatibility barriers to cross-species transmission of TPE are low, though other barriers to spillover might exist [16]. In fact, recent multilocus sequence typing (MLST) of genes unique to TPE showed that TPE strains infecting wild NHPs cluster according to geography and not by host species in Tanzania, which is compatible with local cross-species transmission of TPE in natural systems [17].

Here, we aimed to further explore TPE genomic diversity across space, hosts and clinical manifestations. For this, we first focused our analyses on a single group of sooty mangabeys (Cercocebus atys) in Taï National Park (TNP), Côte d’Ivoire, whose affected individuals showed orofacial or genital lesions. We then drew on recent findings that TP can be detected in NHP bones [18] to screen NHP bones from TNP and examine whether closely related strains of TPE circulate in this ecosystem, as seen in Tanzania [17]. In addition, we also analysed non-symptomatic bones from 11 species sampled across 8 sub-Saharan countries to further explore the range of hosts that this pathogen infects and be in a position to detect potential broad-scale spatial patterns. Finally, we examined the TPE genomes we generated for evidence of mutations associated with antibiotic resistance, as ongoing eradication efforts are based on antibiotic therapy.

Methods

Study sites and samples

In January 2014, sooty mangabeys from a habituated social group in TNP were observed with orofacial lesions and lesions of their distal extremities; samples were collected from two symptomatic animals and through whole-genome sequencing, TPE was determined to be the cause of the infection [12]. Over the next 2 years, other individuals started showing genital ulcerations and necrotizing dermatitis on the inner parts of the thighs and abdomen, often with visible yellow crusts. Orofacial and genital lesions were still observed on other animals in the group during the study period (Fig. 1). Three more individuals with visible lesions were chemically immobilized using a combination of xylazine (1 mg kg−1) and ketamine (10 mg kg−1) administered by blowpipe. Biopsy and swab samples were collected from orofacial and genital lesions (Table 1). Samples were preserved in RNAlater Thermo Fisher Scientific, Waltham, MA, USA) and shipped to the Robert Koch Institute, Berlin, Germany for molecular analysis.

Fig. 1.

Lesions due to TPE infection in sooty mangabeys. (a) Necrotizing dermatitis of inner parts of the thighs and ventral abdomen with yellowish crusts. (b) Necrotic orofacial lesions. (c) Genital lesions in females. (d) Genital lesions in males.

Table 1.

Lesion and sample types observed at TNP and PCR screening results for each animal

Animal ID	Species	Lesions observed	Sample type	PolA PCR	cfpA PCR
5847	C. atys	Anogenital lesions	Lesion biopsy	Negative	Negative
		Anogenital lesions	Lesion swab	Positive	Positive
		Facial lesions	Lesion biopsy	Negative	Negative
		Body ventral lesions	Lesion biopsy	Negative	Negative
2116	C. atys	Genital lesions	Lesion swab	Positive	Positive
2117	C. atys	Genital lesions	Lesion swab	Positive	Positive
1864	C. atys	Orofacial lesions	Lesion biopsy	Positive	Positive

Animal ID 1864 was previously reported in Knauf et al. [12] and in this study we prepared a new library and resequenced it to improve genome coverage.

To obtain insights into NHP treponematoses in TNP over the last three decades and to obtain an understanding of the host range of this pathogen, we screened non-symptomatic bones collected opportunistically at TNP (n=67; Table S1, available in the online version of this article). Previous studies have shown that informative TP sequences can be obtained from such skeletal remains [18]. We also screened NHP bones collected from 11 additional field sites in sub-Saharan Africa (n=83; Table S1) through the support of the Pan African Programme: The Cultured Chimpanzee (PanAf; http://panafrican.eva.mpg.de/). NHP bones were assigned to particular species using molecular methods targeting the 16S ribosomal RNA gene as described previously [19], coupled with morphological assignment by experts in the field.

DNA extraction

DNA was extracted from skin biopsies (n=4) and swabs (n=3) using the DNeasy Blood and Tissue kit (Qiagen, Hilden, Germany) following the manufacturer’s instructions (Table 1). DNA was extracted from bones using a silica-based method. Briefly, non-lesioned bones were drilled using a fine drill bit at slow speed to produce ~150 mg of bone powder. The drilling was performed in a designated sealed glove box, both to prevent any contamination of the bones and to prevent exposure of researchers to pathogens that might be present in the bones [e.g. in TNP, cases of sylvatic anthrax ( biovar anthracis) are frequent, and this pathogen can be cultured from bones] [20]. The box was UV-sterilized and surfaces were bleached following the drilling of each bone and extraction. Drill bits, the metal pieces that actually come into contact with the bones, were changed for each bone to prevent cross-contamination. DNA extraction from bone powder was performed following a protocol that was modified from one that was used previously [21, 22] and is described in detail elsewhere [23]. Extracted DNA from all samples was quantified with a Qubit fluorometer using a double-stranded DNA high-sensitivity assay kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s instructions. DNA was subsequently stored at −20 °C.

Screening for

To screen for TP DNA in swabs and biopsy samples, we performed an end-point polymerase chain reaction (PCR) assay targeting the 67 bp of the polA gene fragment, using primers previously developed for screening human clinical specimens (Table S2) [24]. PCR reactions were performed in 25 µl reactions; up to 200 ng of DNA was amplified using 1.25 U of high-fidelity Platinum Taq polymerase (Thermo Fisher Scientific, Waltham, MA, USA), 10× PCR buffer (Thermo Fisher Scientific, Waltham, MA, USA), 200 µM dNUTPs, 4 mM MgCl2 and 200 nM of both forward and reverse primers. The thermal cycling profile was as follows; denaturation at 95 °C for 5 min, followed by 40 cycles of 95 °C for 15 s, 60 °C for 30 s and 72 °C for 1 min, with a final elongation step at 72 °C for 10 min. Known TPE-positive DNA extracts and negative controls were included.

The 67 bp amplified product is too short for direct Sanger sequencing, so to confirm the results of this initial screening positive samples were further tested with a semi-nested assay targeting the cytoplasmic filament protein gene (cfpA) [25]. These primer pairs yield a 352 bp outer product in the primary PCR and a 189 bp inner fragment in the second round PCR (Table S2). In both primary and second round PCR assays, the reactions were performed as follows: ~200 ng of DNA was amplified in a 25 µl reaction using 1.25 U of high-fidelity Platinum Taq polymerase, 10× PCR buffer, 200 µM dNUTPs, 4 mM MgCl2 and 200 nM of both forward and reverse primers. For the second round PCR, 2 µl of a 1 : 20 dilution of the first round PCR product was used as input template. The thermal cycling profiles for both rounds of the nested PCRs were the same as in the initial screening assay described above.

PCR products were visualized on a 1.5 % agarose gel stained with GelRed (Thermo Fisher Scientific, Waltham, MA, USA). Positive bands were purified using the PureLink Quick Gel Extraction kit (Thermo Fisher Scientific, Waltham, MA, USA) following the manufacturer’s protocol. Purified products were stored at −20 °C until they were sequenced using the BigDye Terminator v3.1 Cycle Sequencing kit (Thermo Fisher Scientific, Waltham, MA, USA) and sequences were compared to publicly available sequences in EMBL through blast [26]. All samples that tested positive in the confirmatory cfpA PCR assay were selected for whole-genome in-solution hybridization capture and high-throughput sequencing.

Potential DNA degradation in bone samples precluded the use of the cfpA PCR assay (demonstrated in [18]). To select the most promising samples, we estimated copy numbers in the bones using a real-time qPCR also targeting the 67 bp fragment of the polA gene [24]. Samples were tested in duplicate. Briefly, 5 µl of total DNA was amplified in a 25 µl qPCR reaction containing 10× PCR buffer, 200 µM dNUTPs, 4 mM MgCl2, 300 nM of both forward and reverse primers, 100 nM of a specific probe and 0.5U of high-fidelity Platinum Taq polymerase. The thermal cycling profile was set as follows DNA denaturation at 95 °C for 10 min followed by 45 cycles at 95 °C for 15 s and 60 °C for 34 s. All bone samples that had detectable TP DNA in duplicate reactions were selected for in-solution hybridization capture and high-throughput sequencing.

Library preparation, genome-wide capture and high-throughput sequencing

For biopsy and swab sample extracts, we sheared 1000 ng of DNA per sample to 400 bp fragments using the Covaris S2 ultrasonicator (intensity: 4; duty cycle: 10%; cycles per burst: 200; treatment time: 55 s; and temperature: 4–5 °C). Bone samples were not sheared due to the potentially fragmented nature of DNA in these older specimens. Two library preparation methods were applied (Table 2). Single-indexed Illumina libraries had previously been built and sequenced from samples 1864 and 22_52 using the Accel NGS ds DNA Library Preparation kit (Swift Biosciences, Ann Arbor, MI, USA) [12, 18]. In an effort to improve genome coverage for these samples, we generated new dual-indexed libraries using the NEBNext Ultra II DNA Library Prep kit (New England Biolabs, Ipswich, MA, USA). The different library preparation methods used reflect the adoption of a new method in the laboratory because of higher library conversion rates with the NEBNext Ultra II DNA Library Prep kit; we have no reason to think that these different methods will systematically affect the genomes generated. For all other samples, DNA was converted into dual-indexed Illumina libraries using the NEBNext Ultra II DNA Library Prep kit (New England Biolabs, Ipswich, MA, USA; samples: Tai_105, Boe_092, 2116, 2117, 5847). All generated libraries were quantified using the KAPA library quantification kit (KAPA Biosystems, Wilmington, MA, USA) following the manufacturer’s instructions. In all library preparations, chicken DNA was included as a control, as no subspecies of TP is known to infect birds.

Table 2.

Summary table for library preparation, capture and sequencing methods used for each sample analysed in this study

Sample ID	NHPs species	Sample type	Library preparation method		Site/country
			NEBNext	Accel NGS
Boe_092	C. atys	Bone	Yes	No	Boe/GB
2117	C. atys	Swab – genital lesion	Yes	No	TNP/CI
5847	C. atys	Swab – genital lesion	Yes	No	TNP/CI
2116	C. atys	Swab – genital lesion	No	Yes	TNP/CI
1864	C. atys	Biopsy – face lesion	Yes	Yes	TNP/CI
Tai_105	P. badius	Bone	Yes	No	TNP/CI
22_52	P. badius	Bone	Yes	Yes	TNP/CI
Control	Chicken DNA	Tissue	Yes	No	na

Samples 1864 and 22_52 had been sequenced previously, but in this study we prepared new libraries and resequenced them to improve genome coverage [12, 18].

na, not applicable.

Libraries were enriched for sequences using in-solution hybridization capture with biotinylated RNA baits following the manufacturer’s protocol (myBaits, Arbor Biosciences, Ann Arbor, MI, USA). The baits spanned the simian-derived Fribourg-Blanc reference genome (RefSeq ID: NC_021179.1) with a two-fold tiling. In-solution hybridization capture was performed for two rounds of 48 h each. After each round of capture, a post-capture amplification step was performed using the KAPA HiFi HotStart library amplification kit (KAPA Biosystems) with P5 and P7 Illumina primers to generate ~200 ng of enriched DNA per sample. The post-capture amplification thermal profile was as follows: initial hot start at 98 °C for 2 min followed by 12 to 16 cycles at 98 °C for 20 s, 65 °C for 30 s and 72 °C for 45 s, with a final elongation step at 72 °C for 5 min. Enriched libraries were quantified using the KAPA library quantification kit (KAPA Biosystems, Wilmington, MA, USA). Prior to sequencing, libraries were diluted to 4 nM and pooled for sequencing on an Illumina MiSeq (Illumina, San Diego, CA, USA) with 300 bp paired end reads (V3 chemistry; Table 2) at the Robert Koch Institute (Berlin, Germany).

Bioinformatics analysis

The paired-end reads generated here, along with published SRAs from prior TP sequencing efforts (Table S3), were trimmed using Trimmomatic v0.38, removing the leading and trailing reads below a quality score of 30, clipping any part of the read where the average base quality across 4 bp was less than 30, and removing reads less than 30 bp in length [27]. The surviving read pairs were merged using Clip and Merge version 1.7.8 with default settings [28]. Merged reads and surviving single-end reads were combined and mapped to TPE Fribourg-Blanc (RefSeq ID: NC_021179.1) using BWA-MEM [29] with a minimum seed length of 29. Mapped reads were sorted using Picard’s SortSam and subsequently deduplicated with Picard’s MarkDuplicates (https://broadinstitute.github.io/picard/index.html). Alignments with a mapping quality score smaller than 30 and a mapping length lower than 30 were also removed using SAMtools [30]. We called two consensus sequences for each genome using Geneious v11.1.5 [31]: (i) the first consensus required a minimum of 10 unique reads per position (i.e. 10× coverage) with at least 95 % identity for a base to be unambiguously called; (ii) the second consensus required 5× coverage and at least 95 % identity. Positions that did not meet the coverage requirement were subsequently treated as missing data. In the following, we mostly report on the results obtained with the 5×, 95 % dataset, which includes all the (partial) genomes generated in this study. To confirm these findings, we repeated the analyses using the 10×, 95 % dataset that included only the higher-quality genomes available, thereby excluding two partial genomes of lower quality that we assembled from bone material (Tai 105 and Boe 092).

Each set of whole-genome consensus sequences was aligned using the multiple sequence alignment programme MAFFT [32]. We then removed all previously described paralogous and putative recombinant genes (Table S4) and selected conserved blocks using the Gblocks tool [33] in SeaView v4 [34]. Phylogenetic inference was performed on the resulting alignments of informative positions after stripping all identical sites and ambiguities in the final datasets: 5×, 95 % (1737 positions: 37 sequences including the two TPE genomes from bones) and 10×, 95 % (1739 positions: 35 sequences). To account for the effects of potential long branch attraction, we ran analyses on versions of these alignments that either included or excluded non-TPE sequences.

Final alignments were uploaded to the online ATGC PhyML-SMS tool (http://www.atgc-montpellier.fr/phyml-sms/) for construction of a maximum-likelihood (ML) phylogeny using smart model selection [35] with the Bayesian information criterion and subtree pruning and regrafting (SPR) for tree improvement, but otherwise using default settings. Branch robustness was estimated using the Shimodaira–Hasegawa approximate likelihood ratio test (SH-like aLRT) [36]. We also used RAxML Next Generation [37], which allowed for an accounting of the number of invariant sites that had been stripped from the alignment to correct for acquisition bias, using 200 bootstrap replicates and the model selected by PhyML-SMS [38, 39]. Both ML analyses retrieved very similar trees; here we report the ML trees from this analysis that accounted for the number of invariant sites. The ML trees were then rooted using TempEst (version 1.5.1), which estimated the best-fitting root of these phylogenies using the heuristic residual mean squared function, which minimizes the variance of root-to-tip distances [40]. Evolutionary pairwise distance between TP strains were extracted from the X5, 95 % ml phylogeny using the Patristic program [41].

To explore the robustness of our phylogenetic analysis, we also ran Bayesian phylogenetic analyses with beast (version 1.10.4) using four different models. These models all used the general time-reversible (GTR) nucleotide substitution model identified in the PhyML-SMS analysis, but combined either a strict clock model or an uncorrelated lognormal relaxed clock model with demographic models assuming a coalescent process and a constant population size or a birth–death process. These models covered a plausible range of clock and tree priors, similar to those explored in another recent publication on syphilis genomics [42]. We incorporated the number of invariant sites into these models to account for sampling and ascertainment bias. For all models, we examined the output of multiple runs for convergence and appropriate sampling of the posterior using Tracer (version 1.7.1) [43] before merging runs using Log Combiner (version 1.10.4) [44]. The maximum clade credibility (MCC) tree was then identified from the posterior set of trees and annotated with Tree Annotator (version 1.10.4: distributed with beast). The resulting ML and MCC tree files were edited using iTOL (https://itol.embl.de/) [45].

For those bone samples for which the aforementioned phylogenetic pipeline was unable to resolve the position of their partial genomes in the phylogeny, we performed phylogenetic read placement using the evolutionary placement algorithm tool EPA-ng to determine the position of individual reads on the TP phylogeny [46]. We performed read placement on a red colobus bone sample (22-52; this particular sample exhibited a suspiciously long branch in initial ML and Bayesian analyses, which we considered to be potentially indicative of an assembly problem, e.g. the inclusion of reads from TPE and an environmental spirochete) and two previously sequenced chimpanzee bone samples (11786 and 15028) that were collected in TNP, for which the presence of TP was confirmed but no subspecies assignment had been performed [18]. Briefly, we selected filtered merged and single-end reads that mapped to TPE and were 50 bp or longer. The surviving reads were then aligned to the whole-genome alignment using the parsimony-based phylogeny-aware read alignment program (PaPaRa) [47]. The resulting PaPaRa alignment was then split into the query reads and the original alignment using the split function available in the EPA-ng toolkit [47]. To estimate the best fitting evolutionary model for the phylogenetic placement of the query reads, both the reference tree and reference alignment were evaluated using the RaxML-ng toolkit [48]. The best fitting model was then used to place query reads on the reference tree with the EPA-ng tool. The resulting EPA-ng jplace tree files were visualized as heat trees depicting the percentage of reads placed on each branch using the gappa toolkit [49]. For tree visualization, every query read was treated as a point mass concentrated on the highest-weight placement and the multiplicity of each query read was set to 1.

Results

Screening PCR and whole-genome capture

All symptomatic animals sampled by biopsy or swab tested positive for in at least one of the sample types collected (Table 1). Sequences generated from the respective assays were all identical (a representative sequence was uploaded to Zenodo: doi.org/10.5281/zenodo.3540499). In addition, based on the polA gene qPCR assay, we detected TP DNA in NHP bones from the following species/subspecies and field sites: 14 western chimpanzees (Pan troglodytes verus), 5 western red colobus (Piliocolobus badius), 2 Diana monkeys (Cercopithecus diana), and 1 sooty mangabey from TNP in Côte d’Ivoire; 5 eastern chimpanzees (Pan troglodytes schweinfurthii) from Bili-Uere in the Democratic Republic of the Congo; 2 western chimpanzees from Loango National Park, Gabon; 1 sooty mangabey from Boe in Guinea-Bissau; 1 western chimpanzee from East Nimba, Liberia; 1 western chimpanzee and one lesser spot-nosed monkey (Cercopithecus petaurista) from the nationwide survey in Liberia; 2 western chimpanzees from Sapo National Park, Liberia; 1 eastern black-and-white colobus (Colobus guereza) from Gashaka Gumti National Park, Nigeria; 1 blue monkey (Cercopithecus mitis) and 1 eastern black-and-white colobus from Budongo in Uganda (Table S1).

From the C. atys study group in TNP, we were able to sequence three new TPE genomes from biopsy and swab samples with 5× coverage of 47.5–97.5 % of the genome (Table 3). We reanalysed a previously sequenced sample from a C. atys with facial lesions (1864) [12], and improved the 5× genome coverage from 82.4 to 95.5 %. Further, we recovered partial TPE genomes from two P. badius bones (Tai 105 and 22_52) from TNP, Côte d’Ivoire, as well as one C. atys (Boe 092) bone from Boe, Guinea Bissau, with 76.1, 62.4 and 69.6%, respectively at 1× genome coverage and 39.3, 26.9 and 21.7 % at 5× genome coverage, respectively. In the control library (chicken DNA), only a single TP read survived quality control. All genomes were distinct, including those obtained from individuals belonging to a single social group of C. atys. None of the genomes generated in this study had the A2058G and A2059G mutations in the 23S ribosomal RNA gene, which have been demonstrated to confer antimicrobial resistance to macrolide antibiotics [50].

Table 3.

Mapping results for non-human primate subsp. strains from TNP and BNP determined in this study and those sequenced previously

Sample ID	NHPs species	Sample type	Deduplicated reads mapped	Positions covered [1×]	Positions covered [5×]	% genome coverage [1×]	% genome coverage [5×]
Boe 092	C. atys	Bone	36 407	793 392	247 484	69.6	21.7
2117	C. atys	Swab – genital lesion	44 143	1 009 366	541 728	88.5	47.5
5847	C. atys	Swab – genital lesion	457 063	1 117 566	1 111 969	98.0	97.5
2116	C. atys	Swab – genital lesion	131 550	1 119 023	1 066 350	98.1	93.5
1864	C. atys	Biopsy – face lesion	202 416	1 116 210	1 089 159	97.9	95.5
Tai 105	P. badius	Bone	42 719	867 318	448 209	76.1	39.3
22_52	P. badius	Bone	43 373	711 608	306 789	62.4	26.9
Control	Chicken	Tissue	1	52	0	0.005	0.0

Phylogenetic analyses

Phylogenetic analysis of genomes generated as part of this study, all other TPE and TEN genomes, and a representative selection of TPA genomes from GenBank (Table S3) yielded tree topologies largely consistent in both the RAxML next-generation and beast-based approaches and the two minimum coverage thresholds used. The ML and MCC tree topologies resolved into distinct reciprocally monophyletic groups representing the TP subspecies (TPA, TPE and TEN). The TPE clade included both human- and all NHP-infecting strains, while TPA and TEN clades consisted only of human strains (Fig. 2). These analyses clearly showed that all NHP-derived strains were TPE, but the relationships within the TPE clade were poorly resolved, potentially due to the long branches separating these clades. Long branches have the potential to complicate phylogenetic analyses [51]; thus we focused our examination of relationships of strains with the TPE clade on the ingroup phylogenetic analyses presented below.

Fig. 2.

Maximum clade credibility tree of strains (relaxed clock model assuming a birth–death process) and map showing the sampling locations for all TP genomes that appear in the tree. All simian-infecting strains are shown in bold with tip labels showing the host species, location of the sampled lesion and sample ID. Genomes generated in this study are shown in red with a minimum coverage of 10× to call a base, and a threshold of 95 % identity for a base to be called. Branches supported by SH-like aLRT values <0.90 in the maximum-likelihood tree and posterior probabilities <0.95 in the maximum clade credibility tree are indicated in grey. The scale shows nucleotide substitutions per variable site.

Our ingroup phylogenetic analysis revealed that neither the ML nor the MCC tree topologies showed monophyly of TPE strains infecting NHPs based on their clinical manifestations, regardless of the minimum coverage thresholds used. In other words, TNP strains causing different clinical manifestations did not form statistically supported, reciprocally monophyletic groups based on clinical manifestations sampled (Figs 2 and 3, Table 4).

Table 4.

Summary table of TPE ingroup phylogenetic analyses. The support for monophyletic groups of simian and human strains based on the different probabilistic methods, molecular clocks and tree priors is shown

Dataset	Probabilistic method	Clock	Tree prior	Based on clinical manifestations	Senegal–The Gambia– Guinea-Bissau	Tai National Park	Human-infecting strains
5×, 95 %	ML	na	na	No	Yes	Yes (weak)	No
	BMCMC	Strict	Coalescent	No	Yes	No	No
	BMCMC	Strict	Speciation	No	Yes	No	No
	BMCMC	Relaxed	Coalescent	No	Yes	No	Yes (weak)
	BMCMC	Relaxed	Speciation	No	Yes	No	Yes (weak)
10×, 95 %	ML	na	na	No	Yes	Yes (weak)	Yes (weak)
	BMCMC	Strict	Coalescent	No	Yes	Yes	Yes
	BMCMC	Strict	Speciation	No	Yes	Yes	Yes
	BMCMC	Relaxed	Coalescent	No	Yes	Yes	Yes
	BMCMC	Relaxed	Speciation	No	Yes	Yes	Yes

na, not applicable.

Rather, it appeared that simian strains formed reciprocally monophyletic groups based on geography, as strains from the same location generally seemed to cluster together. We observed the most consistent support for a clade of NHP-infecting TPE strains from Senegal, The Gambia and Guinea-Bissau, with all analyses and models showing this geographical clustering, in most cases with strong statistical support (Figs 3 and S1, Table 4). This clade consistently included a sooty mangabey bone-derived TPE (Boe_092) collected in Guinea-Bissau and included the Fribourg-Blanc strain isolated from a baboon in Guinea and green monkeys (Chlorocebus sabaeus) from neighbouring Senegal and The Gambia (Fig. 3, Table 4). We found less consistent support for a clade including NHP-infecting TPE strains from TNP, specifically TPE-infecting western red colobus and sooty mangabeys (Figs 3 and S1, Table 4). For the 5×, 95 % dataset, statistical support for this TNP clade was relatively weak, predominantly due to inconsistent placement of the low-coverage 1864 genome; support for a monophyletic TNP clade was much stronger in the 10×, 95 % dataset (Fig. S1, Table 4). In addition, we observed that human-derived TPE strains formed a monophyletic group with relatively strong statistical support in the 10×, 95 % dataset (Fig. S1, Table 4).

Fig. 3.

Maximum clade credibility tree from the ingroup analysis of TPE strains (relaxed clock model assuming a birth–death (speciation) process). All simian-infecting strains are shown in bold with tip labels showing the host species, location of the sampled lesion and sample ID. Genomes generated in this study are shown in red with a minimum coverage of 10× to call a base, and a threshold of 95 % identity for a base to be called. Branches supported by SH-like aLRT values <0.90 in the maximum-likelihood tree and posterior probabilities <0.95 in the maximum clade credibility tree are indicated in grey. The sampling locations that appear on this tree are those shown on the map in Fig. 2. The scale shows nucleotide substitutions per variable site.

The analysis of bone-derived TPE strains further supported a geographical structuring of NHP TPE genomic diversity. Read placement of bone-derived TPE reads from one western red colobus specimen (22_52) and two western chimpanzees specimens (11 786 and 15 028) from TNP showed that the majority of reads fell within the TNP clade (Fig. 4).

Fig. 4.

Phylogenetic read placement of bone samples. Heat tree visualization of phylogenetic placement of TPE mapped reads from bone samples on to the TP MCC (10× coverage and 95 % threshold) reference tree using the evolutionary placement algorithm (EPA-ng). The approximate percentage of reads placed on to a particular branch of the cladogram is shown as a linearly scaled colour density. Genomes generated in this study are shown in red. (a) Sample 11786 (P. troglodytes verus: total number of reads, 517), (b) sample 15028 (P. troglodytes verus: total number of reads, 3581) and (c) sample 22_52 (P. badius: total number of reads, 19 389).

Within these geographical clusters, for TPE-infecting particular species at each field site, there appeared to be differences in average genetic distances between strains. The sooty mangabey-infecting strains at TNP had an average patristic distance of 4.23333E-05 %, compared to 0.000152, 0.000002 and 0.000002 % for strains infecting Papio anubis in Tanzania and Chlorocebus sabaeus monkeys in Senegal and The Gambia, respectively (Table S5), although sample sizes were small, precluding a formal statistical analysis.

Discussion

In agreement with previous studies [11, 13, 17], our phylogenetic analyses demonstrated that the NHP-infecting TP strains all belonged to the TPE clade. We confirmed the presence of TPE in both orofacial and genital lesions, adding to a growing body of evidence that TPE causes a diversity of symptoms in NHPs [11, 25, 52–55]. We found no evidence that TPE strains causing these different symptoms formed separate monophyletic groups, suggesting that TPE pathology may not be determined by bacterial properties alone, but also by host factors (e.g. individual immune status) or the route of exposure.

While TPE genomes did not cluster by symptoms, we found that simian TPE strains from different ecosystems in sub-Saharan Africa formed monophyletic groups largely based on their geographical origin. These results support and extend the findings of Chuma et al. [17], who used multilocus sequence typing (MLST) data to reveal that TPE sequences from NHPs in Tanzania did not form monophyletic clades based on host species, but rather clustered based on geography [17]. Here, considering a larger geographical scale, we find that TPE strains from TNP cluster separately from strains infecting NHPs in Guinea, Guinea-Bissau, Senegal and The Gambia (all countries in close proximity to one another). This geographical signal was supported by phylogenetic read placement of TPE reads derived from the bones of two NHP species collected in TNP, which fell predominantly in the TNP clade. This observed incongruence between the host phylogeny and the phylogeny of their TPE infections could be indicative of cross-species transmission or infection from some common unknown source within a habitat [56, 57].

We observed high TPE diversity in some ecosystems, which could lend support to the hypothesis of exposure of NHPs to a diversity of strains, potentially from other species sharing an ecosystem. We observed the highest level of diversity within one region in the TPE strains sampled from a single sooty mangabey group in the TNP ecosystem. Indeed, the average number of single-nucleotide polymorphisms (SNPs) separating any two TNP sooty mangabey-infecting TPE genomes was 97 (5×, 95 % identity threshold). While one pair differed at only 3 positions, potentially suggesting an epidemiological link, the most divergent pair had 264 SNPs. Considering that the TPE mutation rate in humans has been estimated to be ~1.21×10−7 per nucleotide site per year (or lower [58]), equalling fewer than 0.138 mutations per year across the 1.14 MB TPE genome, the 264 SNPs between the 2 most divergent TNP strains would likely have required hundreds to thousands of years to accumulate (Table S6). High levels of divergence precluding a direct epidemiological link were also observed for two TPE strains sampled from the Lake Manyara National Park ecosystem in Tanzania [17]. In other ecosystems, TPE genomes from the same species were very similar, compatible with an epidemiological link (green monkeys in Bijilo Forest Park, The Gambia and Niokolo National Park, Senegal [12]). The factors driving the high diversity of TPE in some ecosystems and not others are unknown, but between-species transmission in diverse primate ecosystems, such as TNP, could play a role.

Interspecies interactions that could facilitate transmission via direct contact between NHP species inhabiting TNP are well documented; these include a strong predator–prey relationship between chimpanzees and red colobus, as well as direct contact between monkeys that spend large amounts of time in mixed-species associations (e.g. grooming, fighting, play, mating [59, 60]). Another transmission mode that has long been suggested for TPE is vectorial transmission; under experimental conditions viable TP spirochetes were transmitted by flies between different host species causing clinical disease [61, 62]. Knauf et al. amplified TP DNA from flies in ecosystems where TP infections in NHPs are common [63] and Gogarten et al. showed that flies carrying TP DNA formed high-density persistent associations with NHP social groups in TNP [64]. Primate-associated flies were observed to move between groups of different species, suggesting that they could be involved in transmitting the yaws pathogen between species, even when NHP are not found in mixed-species associations [64]. Further work is needed to confirm whether flies or other arthropods actually transmit TPE in the wild [65] and to understand the routes and rates of transmission of TPE between NHP species.

Wide spatio-temporal differences in sampling of human- and simian-derived TPE strains precludes a robust assessment of whether zoonotic transmission occurs. The sampling of NHP TPE genetic diversity has improved over the last several years, with 14 of the 29 TPE genomes included in the current analyses originating from NHPs in sub-Saharan Africa. Unfortunately, in countries where NHP-infecting TPE strains have been sequenced, no genomic data from human infections are available. Given the geographical signal observed for simian isolates [17], future studies may benefit from investigating human yaws infections in these regions. Such data will help determine whether zoonotic inter-species transmission of TPE between humans and NHPs actually occurs. Our analysis found some support for the monophyly of human-infecting strains; if larger datasets support this finding, it will suggest that zoonotic spillover does not occur frequently (or at all). If zoonotic spillover of TPE from NHPs to humans does occur, it is noteworthy that all NHP-derived TPE strains characterized to date [17], including those generated in this study, have no mutations in the 23S ribosomal RNA gene, which is known to confer antimicrobial resistance to macrolide antibiotics [50]. Therefore, the available antibiotic treatments are expected to be effective against NHP-derived strains. However, we caution against overinterpreting this finding, as sample sizes are small, meaning that low-frequency mutations would likely remain undetected.

This study joins a growing body of evidence that human and wildlife bones are a useful resource for generating sufficient TPE reads to inform phylogenetic analyses that can both extend the known host range of TPE and push back the date of emergence for TPE in particular populations [18, 66]. Despite a lack of clinical evidence from western chimpanzees and western red colobus in TNP, the phylogenetic placement of bone-derived TPE reads into the TNP clade confirms that sooty mangabeys are not the only species affected by TPE in this ecosystem. The oldest bone from which we were able to phylogenetically assign reads to the TPE clade was collected in 1992, complementing clinical evidence from TNP that only started accumulating from 2014 onward [12]. Our results suggest that, together with archaeological collections [67-69], natural history specimens can provide important genomic information regarding TPE and thereby further our understanding of its ecology and evolution.

Conclusions

We found that genomically diverse TPE circulate in NHP, even in a single social group [17]. On a larger spatial scale, whole-genome sequences and bone-derived TPE sequences allowed phylogenetic analyses that revealed that the genomic diversity of TPE strains derived from NHPs is geographically structured. This pattern is compatible with cross-species transmission of TPE within ecosystems, although how often and by what means this transmission occurs remains an important area of future research.

Click here for additional data file.

Click here for additional data file.