Mycobacterium leprae genomes from a British medieval leprosy hospital: towards understanding an ancient epidemic

Archaeology and osteology

Excavations at St Mary Magdalen, Winchester, are amongst the most extensive of any British leprosy hospital to date [11, 12]. The hospital is one of the earliest such institutions known in Britain and, as such, provides insights into both early institutional care and the nature and status of such communities at a relatively early date. Such hospitals may well have been a model for succeeding charitable institutions of social care. Thirty-eight burials have been excavated from the site, of which 33 (87%) show osteological signs of leprosy, a much higher percentage than observed in other British material. Of these, five M. leprae strains have been SNP genotyped (Sk2, Sk7, Sk8, and Sk14, and Sk19) [10] and 2 whole genomes (Sk2 and Sk8) sequenced [7].

The current study focuses on typing a further burial, Sk27, and generating whole genome sequences of the M. leprae from both Sk14 and Sk27. Both were excavated from the earlier and more northern of two cemeteries at the site (together with skeletons Sk7, Sk8, Sk9 and Sk19 [10]). Based upon comparative and diagnostic material this cemetery was thought to have been in use primarily in the decades immediately following the Norman Conquest, between 1070 AD and 1100 AD. This interpretation was confirmed by ¹⁴C dating, with dates of between 940 cal AD and 1160 cal AD [11]. Sk14 and Sk27 were ¹⁴C dated to between 955–1033 cal AD (WK 28629) and 1020–1162 cal AD (SUERC-39676), respectively.

The remains of Sk14 are those of an adolescent, probably male, displaying clear pathological features of leprosy: circumferential wasting of the foot phalanges, loss of bone of the distal ends of the distal foot phalanges, porosity of the shafts of the metatarsals and porous, woven and compact bone on the shafts of the fibulae, thinning of the palate and nasal bones and rounding of the margins of the nasal aperture [13]. In contrast, Sk27 was a middle-aged adult, with only minimal evidence of characteristic leprosy lesions. Indeed, Sk27 was originally identified as a possible control, non-lepromatous skeleton and only on closer examination was porosity and resorption of bone noted on the dorsal surface of the distal ends of the distal foot phalanges, particularly those for the first metatarsal. The grave of Sk27 was also notable for containing a scallop shell (Figure 1), indicating that the remains are likely those of a pilgrim who had completed the Way of St James to the Cathedral of Santiago de Compostela in Galicia, Spain [11].

Confirmation that skeleton Sk27 had been infected with M. leprae

Samples of bone from Sk27 were scraped from the rhino-maxillary area under clean and controlled conditions using disposable gloves, scalpels and tubes. A control sample was taken from around the vomer (the thin bone that separates the nostrils) of a skeleton, Sk12, which showed no osteological evidence of M. leprae. The vomer is a site that often has a high M. leprae burden in lepromatous leprosy, and so is a good site from which to isolate M. leprae DNA. The presence of M. leprae DNA in the Sk14 extract was previously confirmed by real-time PCR using primers for the multicopy RLEP and the single copy locus, ML1795 (encoding the 18 kDa antigen), both of which gave single bands of the expected size (111 bp and 114 bp, respectively) [13]. The presence of M. leprae DNA in extracts from SK27 was similarly confirmed by PCR. PCR for the multicopy M. tuberculosis loci IS1081 gave no products (79 bp expected) for either sample, indicating an absence of tuberculosis. The control sample Sk12 gave no products with primers for loci of either M. leprae or M. tuberculosis, indicating that contamination of the samples from extraneous sources or cross-contamination between the samples was negligible.

Whole genome sequencing of M. lepraegenomes from skeletons SK14 and Sk27

The control sample, Sk12, was characterized by very low coverage of the M. leprae genome, punctuated by short regions with high coverage and high diversity. This is consistent with a sample containing little or no M. leprae sequence, but with some carry-over of DNA from environmental organisms. This again confirms that there was little or no cross-contamination between samples or from modern M. leprae DNA during the sample preparation. Many of the reads that were retrieved from Sk12 had similarity to environmental organisms that contain genomic regions that are likely to be conserved across mycobacteriaceae and beyond. Such regions apparently have sufficient similarity to the microarray probes to hybridize, particularly in the absence of competing M. leprae DNA, but are sufficiently dissimilar to generate the high diversity observed. Regions with high coverage included approximately 50 genes such as rpoB, tuf, fusA and certain smaller regions within genes that often represent conserved motifs, for example within ABC transporters. Base calls from regions with high coverage in the control sample were not considered reliable and were discounted from the data.

Analysis of the ancient M. lepraegenomes

We have previously used PCR and Sanger sequencing of a subset of SNPs [8] to genotype the M. leprae from Sk14 as type 2F [10]. Using identical methods, we also identified Sk27 as a type 2F strain (Table 2).

Independent analysis of the sequence reads at the University of Tübingen and the University of Surrey (Table 1), generated differing values for coverage and read depth. However, both centres reconstructed very similar genomes with only a single SNP discrepancy (not included in subsequent analysis) between the comparable datasets. Comparison of the Sk14 and Sk27 M. leprae genomes with other modern and ancient genomes identified the presence of very few novel polymorphs (Additional file 1: Table S1). Sk14 contained 3 novel indels, all of which were in pseudogenes, and 3 novel SNPs, only one of which was non-synonymous, encoding a change of Gly₁₆₃ to valine in hypothetical protein ML1637. Sk27 had a single novel synonymous SNP and no identifiable novel indels. Phylogenetic analysis generated trees comparable to those of Schuenemann et al.[7] and placed Sk14 and Sk27 in the highly conserved type 2F branch of M. leprae (Figure 2 and Additional file 2: Figure S1). It is notable that the strain most closely related to Sk14 is another isolate from Winchester, Sk8, described in our previous studies [7, 10].

Bone retrieval and methods

The skeletons were excavated between 2009 and 2011 by a team from the University of Winchester. The burials were recovered as articulated remains from single graves and associated grave fills were subject to 100% sampling. Pathological conditions were recorded and photographed in detail and diagnosed with reference to appropriate sources [11, 12]. All necessary permits were obtained for the described field studies, including a licence (−0070) to exhume and retain human remains, provided by the Ministry of Justice, 102 Petty France, London, SW1H 9AJ.

Bone sampling and DNA extraction

Samples of bone were taken from the rhino-maxillary region of skeletons Sk14 and Sk27 to maximize the chances of recovering M. leprae DNA. Measures were taken to minimize the opportunities for cross-contamination between cases and, use of modern DNA as a positive control was avoided completely [10]. The vomer region of a skeleton, Sk12, with no signs of leprosy was sampled as a negative control and was treated identically to the other samples throughout. The bones were ground to a fine powder and DNA extracted using the NucliSens™ extraction kit (bioMériux Limited, Boxtel, The Netherlands) as described previously [13] in laboratories physically separate from rooms in which leprosy DNA has been previously amplified.

PCR amplification of specific loci

Real-time PCR assays were carried out on a Mx3005P qPCR system (Agilent Technologies, Wokingham, UK). Methods for the PCR amplification of the multi-copy RLEP locus, the single copy locus coding for the 18 kDa antigen (ML1795), and for the M. tuberculosis specific multi-copy loci IS1081 locus are all described by Taylor, 2013 [10].

Whole genome amplification, enrichment and analysis

DNA was extracted, enriched and sequenced either at the University of Surrey alone, for Sk27, or independently at both the University of Surrey and University of Tübingen for Sk14. Both protocols are fully described by Schuenemann et al. [7]. Both involved removing uracil residues with USER enzyme (New England Biolabs, Hitchin, UK), repairing the DNA, ligating linkers to the ends, and PCR amplifying before enriching for M. leprae sequences using microarrays (once at the University of Tübingen and twice at the University of Surrey). Samples were pair-end sequenced with double indices on HiSeq2000 or MiSeq machines.

At the University of Surrey, reads were quality controlled and aligned to the M leprae TN genome [9] using Bowtie2 [24]. Duplicate reads were removed with MarkDuplicate script [25] and SNPs assigned. The criteria for inclusion in the SNP/indel list were that polymorphs had to be outside of annotated repeat regions, to have a read depth of more than 5 reads in the sample, but less than 5 reads in the control sample Sk12, have a QUAL score of > =50 (not required for longer indels) and appear in more than 80% of reads. Indels had an additional requirement that all qualifying reads had to span the indel.

At the University of Tübingen the reads were first quality controlled using FastQC, followed by adapter clipping, merging of corresponding paired-end reads and, finally, quality trimming of the resulting reads. The reads were mapped against the M. leprae TN genome as a reference using the Burrows-Wheeler Aligner (BWA). Mapped reads were subject to duplicate removal using Samtools’ rmdup method. The Genome Analysis Toolkit (GATK) was then used to generate a mapping assembly and to call SNPs for each strain. The TN reference base was called if the quality score was at least 30 and the position was covered by at least 5 reads. A variant position (SNP) was called if the position was covered by at least 5 reads and the fraction of mapped reads containing the SNP was at least 90%.

To be included in phylogenetic analysis, SNPs had to either be identified at both centers, or identified at only one center, with no contradictory data from the other. Phylogenetic analysis were carried using MEGA5 [26] using Maximum Likelihood (using the Tamura-Nei model), Maximum Parsimony (using Max-mini Branch-&-bound algorithm) and Neighbour-Joining (using the number of differences method) methods.

Availability of supporting data section

All raw read files have been submitted to the related trace archive of the National Center for Biotechnology Information Sequence Read Archive entry PRJNA200950. Alignments and phylogenetic trees have been submitted to TreeBASE with study ID 15537.