Rickettsia felis DNA recovered from a child who lived in southern Africa 2000 years ago

Abstract

The Stone Age record of South Africa provides some of the earliest evidence for the biological and cultural origins of Homo sapiens. While there is extensive genomic evidence for the selection of polymorphisms in response to pathogen-pressure in sub-Saharan Africa, e.g., the sickle cell trait which provides protection against malaria, there is inadequate direct human genomic evidence for ancient human-pathogen infection in the region. Here, we analysed shotgun metagenome libraries derived from the sequencing of a Later Stone Age hunter-gatherer child who lived near Ballito Bay, South Africa, c. 2000 years ago. This resulted in the identification of ancient DNA sequence reads homologous to Rickettsia felis, the causative agent of typhus-like flea-borne rickettsioses, and the reconstruction of an ancient R. felis genome.

Fig. 1. The provenience of the Later Stone Age hunter-gatherer skeletal remains recovered from a mound formed by a shell midden overlooking the beach in Ballito Bay, South Africa.

a The provenience of the Later Stone Age hunter-gatherer skeletal remains recovered from a mound formed by a shell midden overlooking the beach in Ballito Bay, KwaZulu-Natal Province, South Africa (map modified from GISGeography ‘https://gisgeography.com/south-africa-map/’). b The cranial remains of the BBayA male child indicating aDNA sample sources, i.e., DNA was extracted and sequenced from bone samples acquired from the left petrous bone (LPB), right petrous bone (RPB) and the upper left premolar (ULPM) (image adapted from Pfeiffer et al., 2019). c The C14 date (1980 ± 20 cal. BP) obtained for the remains of the child.

Fig. 2. Genome reconstruction of the ancient BBayA R. felis genome.

a Genome reconstruction of the ancient BBayA R. felis genome and mapping of the ancient genome to the genome of R. felis LSU-Lb. The R. felis BBayA genome consists of 1,512,774 bases and 69 contigs, with a N50 of 42,410 bases and the longest contig comprising 121,989 bases. It exhibits a GC value of 32.5% and a coding density of 84.58% for 15161 predicted proteins. A comparison was performed with the chromosome and plasmids (pRF and pLBaR). Rings (from outer to inner ring) show ORFs, SNPs, GC skew, GC content and coverage. Following the initial mapping of our datasets against all available (i.e., 126) NCBI Rickettsia reference genomes (Supplementary Data 3), the genome coverage analysis of BBayA R. felis was performed using the reads mapping to the R. felis LSU-Lb genome (the closest phylogenetic homologue to the ancient BBayA R. felis strain), with average coverage estimated at 58.01-fold. We estimated that 94.73% of the LSU-Lb genome was covered by the BBayA R. felis reads. The RHS-like toxin plasmid (pLbaR-38) occurs only in R. felis LSU-Lb, with R. felis URRWxCal2 and the ancient BBayA being devoid of it. The coverage for the LSU-Lb pRF plasmid within the BBayA read dataset is 132.26 (mean) and 106.19 (trimmed mean), and 0/0 for pLBaR. b DNA damage pattern analysis for the BBayA R. felis reads using mapDamage. G-to-A and C-to-T misincorporations are plotted in blue and red, respectively, and grey lines indicate all possible misincorporations. c DNA fragment read-length distributions of the BBayA R. felis reads, exhibiting a mean read-length of 82.30 base-pairs (bp). d Distribution of edit distance of high quality BBayA R. felis reads mapping to R. felis LSU-Lb and R. felis URRWxCal2.

Fig. 3. Phylogenomic reconstruction of the Rickettsia genus including the ancient R. felis strain derived from BBayA.

a Maximum likelihood (ML) phylogenomic tree including 127 Rickettsia genomes, b pruned ML phylogenomic tree including 31 genomes from the TG and TRG Rickettsia groups and the five R. felis genomes and c inset of the R. felis clade. Purple dots in leaf nodes represent confident ultrafast bootstrap support values of >95%. Multiple sequence alignment of the codons of 138 concatenated core genes for the 127 Rickettsia genomes was used for the phylogenomic reconstruction using the maximum likelihood method (ML). The ML reconstruction was obtained with the IQtree software using the best substitution model (GTR + F + R10) selected using the ModelFinder option and utilising all the sites from the codon alignment of the 138 selected core genes (103,280 nucleotide sites), from which 53.4% of the sites were constant across the 127 genomes (55,120 nucleotide sites). The parsimony informative sites were 33,343 nucleotide sites which present 17,683 distinct nucleotide composition patterns. The legend displays the branch length. The black star (⋆) indicates the position of the Rickettsia felis BBayA genome. The trees indicate the four Rickettsia groups. i.e., the spotted fever group (SFG), typhus group (TG), transitional group (TRG) and ancestral group (AG) (see Fig. S4 and Supplementary Data 6 for additional maximum-likelihood trees constructed using the same codon alignment). The final resulting tree was analysed in the iTol tree visualisation tool, with phylogenetic reconstructions visualised and managed using the iTOL web server.