Genetic affinities of an eradicated European Plasmodium falciparum strain.

Malaria was present in most of Europe until the second half of the 20th century, when it was eradicated through a combination of increased surveillance and mosquito control strategies, together with cross-border and political collaboration. Despite the severe burden of malaria on human populations, it remains contentious how the disease arrived and spread in Europe. Here, we report a partial Plasmodium falciparum nuclear genome derived from a set of antique medical slides stained with the blood of malaria-infected patients from Spain's Ebro Delta, dating to the 1940s. Our analyses of the genome of this now eradicated European P. falciparum strain confirms stronger phylogeographical affinity to present-day strains in circulation in central south Asia, rather than to those in Africa. This points to a longitudinal, rather than a latitudinal, spread of malaria into Europe. In addition, this genome displays two derived alleles in the pfmrp1 gene that have been associated with drug resistance. Whilst this could represent standing variation in the ancestral P. falciparum population, these mutations may also have arisen due to the selective pressure of quinine treatment, which was an anti-malarial drug already in use by the time the sample we sequenced was mounted on a slide.

Data Summary

Plasmodium falciparum reads from Ebro-1944 have been deposited in the European Nucleotide Archive under accession number ERP114811. All modern P. falciparum samples used for the population genomics analyses were reported by the MalariaGEN Plasmodium falciparum Community Project in 2016 [1] and are provided in Table S1 (available in the online version of this article). All software used in the bioinformatic analyses are publicly available. Positions screened for anti-malarial drug resistance are available in Table S2.

Malaria is a serious infectious disease affecting over 200 million people annually. The disease is caused by species of parasitic protozoans from the genus Plasmodium, which are transmitted by several species of mosquitoes from the genus Anopheles. Today, Plasmodium is restricted to tropical and subtropical latitudes. However, malaria was historically present in most of Europe, spanning from southern Britain and the Mediterranean to as far north as Finland. Spain represented one of its last footholds, where it persisted until the 1960s. Here, we report a substantial fraction of the genome of a 20th century European Plasmodium falciparum strain isolated from slides stained with the blood of malaria-infected patients in the 1940s. We analyse this genome in the context of worldwide modern strains to trace the historical dispersal of P. falciparum in Europe. We find evidence supporting a longitudinal spread from Asia into Europe, over a latitudinal spread from Africa, as well as variants in anti-malarial resistance genes predating the use of most common anti-malarial drugs. Our work highlights the potential of collections of antique medical slides to open new possibilities in the study of ancient microbial genomics, including malaria.

Introduction

Classical Greek accounts in the fourth and fifth centuries BCE describe people with intermittent fevers and infectious symptoms characteristic of malaria [2]. The Roman author Celsius was able to accurately differentiate the clinical symptoms of Plasmodium vivax versus Plasmodium malariae infections [3]. In contrast, it is unclear whether Plasmodium falciparum, the deadliest form of the pathogen, was already present in classical times. While some authors argue that P. falciparum only spread to southern Europe with the dawn of the Roman Empire, historical accounts suggest it may have affected western and central Italy as early as 400–100 BCE [4], before reaching the Po Delta region of northern Italy around 1000 years later [3]. The spread of this pathogen through Italy in historical times is supported by the discovery of a large infant cemetery at Lugnano (Teverina, Umbria, Italy), dating to approximately 450 CE, and the report of P. falciparum ribosomal DNA sequences obtained from one of these skeletons by means of traditional PCR [5]. The recent retrieval of larger amounts of P. falciparum genetic data, by means of capture baits and second-generation sequencing technologies, from teeth sampled in Velia and Vagnari cemeteries (Italy), directly places this parasite in southern Italy by the beginning of the Roman Imperial period [6].

Genetic analyses of DNA recovered from a unique set of antique microscopy slides (1942–1944), stained with the blood of malaria-infected patients in the Ebro Delta of Spain for immediate diagnostic purposes, allowed us to report the complete mitochondrial DNA (mtDNA) genomes of historical P. falciparum and P. vivax [7]. In this study, we analyse data generated through merging the sequence information derived from four different slides, allowing us to reconstruct the partial nuclear genome of this eradicated European P. falciparum strain.

Methods

Sample collection

Four microscope slides, dated between 1942 and 1944, were selected for this study. All were obtained from Dr Ildefonso Canicio’s family collection. Dr Canicio was in charge of the anti-malarial hospital at Sant Jaume d’Enveja (Ebro Delta, Spain), inaugurated in 1925, until his death in 1961. Patients were predominantly local people who worked in the Ebro rice fields and had no history of international travel. The samples consist of four labelled microscopy slides stained with Giemsa – probably not previously fixed – in which parasites were still visible under the microscope. DNA extraction was performed in dedicated ancient DNA (aDNA) laboratories at the Institute of Evolutionary Biology in Barcelona (Spain) and the Centre for GeoGenetics in Copenhagen (Denmark) in 2015 and 2017, respectively, as described in the Supplementary Material.

Ancient sample mapping and assembly

We first analysed the sequenced reads obtained from our eradicated European P. falciparum samples using FastQC (v0.11.7) [8] in order to determine their quality before and after trimming of the adapter sequences. We removed sequencing adapters and reads shorter than 30 bp using cutadapt (v1.3) [9]. We then mapped our reads against the P. falciparum 3D7 and P. vivax Sal1 reference genomes using BWA (v0.7.17) [10] aln, specifying no seeding, a gap open penalty of 2, an edit distance of 0.01 and no trimming; parameters shown to optimize mapping of ancient microbial samples [11]. We removed duplicated reads using Picard (v2.18.6) ‘MarkDuplicates’ and retained all mapped reads with a map quality of at least 30 in SAMtools (v1.6) [12]. As the blood samples had known co-infection with P. vivax, we extracted the sequencing reads mapping more confidently to P. falciparum by selecting those reads that had a lesser edit distance between P. falciparum 3D7 than with P. vivax Sal1 (Fig. S1). The obtained G+C content of the reads matched the expected and extremely low G+C content characteristic of the P. falciparum genome (Fig. S2).

Post-mortem aDNA damage patterns were determined using MapDamage (v2.0.8) [13] (Fig. S3). Most of the P. falciparum reads (0.6294× mean coverage, representing 88 % of the total genome) come from a single slide (labelled CA) sequenced in 2015. The remaining reads come from three other slides with mean genome-wide coverage of 0.0213×, 0.0285× and 0.0314×. The reduced coverage obtained from these slides relates to the overall endogenous sequence content and quality, which may also vary depending on the stage of the parasite’s life cycle when the slides were prepared. To increase the overall coverage, we therefore merged data from all four slides, always calling the dominant allele. We call our resultant composite genome Ebro-1944.

We also generated a reference panel of modern P. falciparum samples from publicly available sequence data [1]. In all cases, we mapped reads using BWA (v0.7.17) [10] mem, before removing duplicated sequence reads and filtering by mapping quality using thresholds as described above. The same procedure was followed for mapping the genome of Plasmodium praefalciparum, which was used as an outgroup. Qualimap (v2.2) [14] was used to generate the mapping metrics for all samples.

Variant calling and dataset creation

To compare Ebro-1944 to current strains in global circulation, we selected modern strains with a mean depth of coverage equal or above 50× and with at least 90 % of the reference genome covered. This filtering strategy resulted in a dataset comprising 206 global samples. We used GATK (v3.7) [15], algorithm UnifiedGenotyper, to call variants, specifying EMIT_ALL_CONFIDENT_SITES and a standard confidence threshold over 50. We subsequently used VCFtools (0.1.14) [16] to filter out positions with less than 40× coverage, heterozygous positions, multi-allelic positions, indels and recombinant sub-telomeric regions (see the Supplementary Material). These samples were later merged using bcftools (v1.3.1) [12] and all variants that were not present in at least three samples were removed leaving 681 486 single nucleotide positions (SNPs). We used these positions to call additional modern samples, which mapped with lower coverage (<50×), applying the same filters. The resultant calls were merged using bcftools. SNPs were called for P. praefalciparum using the same parameters.

For Ebro-1944, we used a different approach to overcome the post-mortem damage associated with aDNA samples. We generated pseudo-haploid calls at the positions identified in the population genetics dataset using SAMtools to call a random base drawn from all possible reads at each site [17]. We then merged all the samples (including Ebro-1944) using bcftools. The complete dataset was filtered for a minimal minor allele frequency of 0.01 and multi-allelic SNPs were removed using VCFtools. All SNPs that were not present in at least 75 % of the samples considered were also removed. This procedure resulted in a final dataset comprising 435 samples and 14 346 bi-allelic SNPs, with 50.85 % of these covered in Ebro-1944.

Population genetics analyses

This dataset was compiled by filtering all modern samples for sites called in Ebro-1944 and retaining all samples that had at least 50 % of these sites covered. This resulted in a dataset of 306 samples and 6 755 SNPs. A principal component analysis (PCA) was performed on this dataset using SmartPCA within the EIG (v.6.0.1) suite of tools [18].

To formally test the relative affinity of Ebro-1944 to South Asian, Oceanian and African strains, we used the full 14 346 biallelic SNP dataset to generate f4 statistics of the form f4(P. praefalciparum, Ebro-1944; X, Y), where we use P. praefalciparum as an outgroup and iterate through all combinations of geographical samples (X and Y) included in our global dataset. This statistic is designed to quantify the covariance in allele frequency differences between P. praefalciparum and Ebro-1944, and X and Y. If P. praefalciparum and Ebro-1944 form a clade with respect to X and Y, then their allele frequency differences should be uncorrelated and f4 will have a value of 0. Deviations from 0, thus, provide the relative affinity of Ebro-1944 to X and Y; positive values indicating a closer relationship to Y relative to X and negative values indicating a closer relationship to X relative to Y. f4 statistics were generated in qpDstat of AdmixTools (v.5.0) [19] and statistical significance was assessed through Z-scores following block jack-knife resampling (Fig. S4).

We additionally sub-sampled the global population genetics dataset more stringently, retaining only the strains with no missingness across the entire alignment, which led to a dataset of 30 strains over 8195 sites (Table S1). We generated patterns of pairwise haplotype sharing across this dataset using Chromopainter v2 [20], specifying a uniform recombination rate of 9.6 kb cM−1 [21] (Fig. S5) and using default estimates of mutation and switch parameters.

We also considered the mitochondrial relationships by placing our Ebro-1944 mtDNA genome (16.1×) in a minimum spanning network together with 435 global strains (Fig. S6). The mtDNA genome was also screened for specific variants restricted to South Asian strains.

Resistance variants screening

A set of variants previously reported in the literature as associated with resistance to different anti-malarial drugs (Table S2) was screened in Ebro-1944 [22-46]. The position of the derived allele was determined using Jvarkit (v.2018.04.05) [47].

Results

Our sequencing of the four antique microscopy slide samples yielded 218 952 P. falciparum reads producing a nuclear genome of 0.67× mean depth and covering 40.99 % of the P. falciparum reference genome 3D7 (Fig. S1). In addition, 1398 reads mapped to the mtDNA genome, generating a mean coverage of 16.1× over 99.43 % of the reference genome. Sequence reads showed damage patterns typical of aDNA (C to T and G to A substitutions at the 5′ and 3′ ends) in ratios of 3.8 and 2.2 %, respectively (Fig. S3). The characteristic aDNA damage of miscoding lesions indicates that our P. falciparum reads are authentically old, deriving from our historic specimen rather than from modern contamination where we would expect no such damage profile.

PCA (Fig. 1a), f4 allele sharing statistics (Figs 1b, c, S4) and haplotype sharing analyses (Fig. S5) indicate that the closest genetic affinity of our European strain is to contemporary samples from central south Asia, including those currently in circulation in Laos, Myanmar and Vietnam. We also detect a shared genetic component with samples from Papua New Guinea, in particular PN0008-C (Fig. 1a). Whilst Ebro-1944 is more closely related to current strains from central south Asia than those from Africa, we found a higher proportion of African haplotypes in Ebro-1944 than in any central south Asian samples (two sample t-test P<2.2×10−16). Such a pattern would be consistent with a common origin of the European and central south Asian P. falciparum populations, followed by secondary introgression of African strains into the European population after the split between the European and central south Asian lineages (Fig. S5).

Fig. 1.

(a) PCA with worldwide sample locations (inset map). (b) f4-statistics under the test relationship f4(P. praefalciparum, Ebro-1944, X, Democratic Republic of Congo), where X iterates through the geographical sampling locations of our included modern P. falciparum strains. A more negative f4 value indicates a closer relationship of Ebro-1944 to X relative to strains sampled from the Democratic Republic of Congo. (c) Z-scores under the test relationship f4(P. praefalciparum, Ebro-1944, X, Colombia) assessed through block jack-knife resampling. All possible f4 topologies were tested (Fig. S4), with Ebro-1944 showing a consistent closer affinity to strains sampled from central south Asia.

The close genetic affinity of the nuclear genome to strains in circulation in central south Asia is further confirmed at the mtDNA level (Fig. S6). Ebro-1944 carries two Indian-specific mtDNA substitutions at positions 276 and 2763 [48]. One additional mutation (position 725) is also present in our sample and is shared with nine contemporary Indian samples and two African samples included in our modern reference dataset (Table S1). We observe mutations at these three positions across all four of our sampled slides, suggesting the strains combined in Ebro-1944 are phylogenetically similar (Table S3). Unfortunately, these three mutations are not covered in a previously published partial strain retrieved from Roman cemetery sites [6]. Additionally, none of our reads overlapped with the ribosomal DNA fragment of European P. falciparum retrieved in a previous work [5].

We screened the Ebro-1944 nuclear genome for the presence of 117 variants that have been previously associated with anti-malarial drug resistance (Table S2). We achieved sequence coverage of at least one read at 62 positions. Of these, only two mutations (H191Y and I876V) displayed the resistance-associated allele (Table S2). Both alleles are located in the multidrug resistance protein 1 encoding gene (pfmrp1), which encodes an ABC family transporter, and has been associated to alterations in quinine, chloroquine, artemisinin, piperaquine and primaquine sensitivity [49]. Although one of these mutations, H191Y, is a C to T substitution that is characteristic of post-mortem DNA deamination, and might, thus, be expected in ancient samples, we believe this to be authentic as the transition is present in two overlapping reads, and is not located at the end of the reads where damage tends to accumulate [50]. While the second variant, I876V, is only covered by a single read, it is an A to G substitution and is, thus, not a common DNA damage motif. To provide additional support for the presence of these derived alleles, we imputed these positions in our partial genome by using a reference panel of modern strains carrying the H191Y and I876V variants (see the Supplementary Material). Interestingly, while the latter variant is distributed worldwide, the former is currently restricted to Asia and Oceania, consistent with an Asian dispersal of P. falciparum into Europe.

Discussion

We showed that an eradicated European P. falciparum strain from the 20th century is most closely related to extant strains from central south Asia, such as those currently in circulation in Laos, Myanmar and Vietnam. Although we detect some evidence for secondary introgression of African P. falciparum into the extinct European population, the significantly stronger affinity to Asian strains argues against a direct African origin [4] of European P. falciparum, and points instead to a migration event between Europe and Asia.

We cannot infer the directionality of the migration between Europe and Asia from the genetic evidence alone. A recent expansion of P. falciparum from Europe to Asia and Oceania, coinciding with European colonial expansion, might be conceivable given the well-characterized role of colonialism in the widespread dissemination of other major infectious diseases, such as lineage 4, the globally distributed agent of tuberculosis [51].

However, a migration of P. falciparum from Europe to Asia does not sit well with historical evidence of the arrival of the parasite in Europe during antiquity [2-6]. Given that our extinct European strain shares significantly more alleles with extant Asian rather than African strains, the arrival of P. falciparum malaria parasites into Europe likely took place from the Asian sub-continent rather than spreading from Africa via the Mediterranean during the Roman Empire.

Plausible historical migrations responsible for a spread into Europe from the East include the extensive commercial exchanges and movement between people of various ethnicities in the Achaemenid Empire (550-330 BCE). Alternatively, P. falciparum might have reached Europe during the subsequent Hellenistic period, connecting India with the Mediterranean, following the conquest of the Achaemenid Empire by Alexander the Great.

The availability of a strain pre-dating most of the currently used anti-malarial drugs allowed us to look for the presence of resistance variants that may inform on the spread of such resistances in the future. The two resistance variants we observed in the pfmrp1 gene could be explained by standing variation for drug-resistance mutations in P. falciparum or may have arisen following the use of quinine for over three centuries; chloroquines were not introduced in Spain until 1948 and initially only in the African colonies.

Our results provide novel insights into the evolution and past demography of one of the world’s deadliest pathogens, which could not have been reached by studying the genomes of extant strains alone. Additional genomic evidence from both medical collections and ancient remains will be needed to reconstruct more precise timings and routes for the spread of malaria into Europe, and could also help in determining the emergence and drivers of resistance to anti-malarial drugs.

Data bibliography

1. The Plasmodium falciparum reference genome used was P. falciparum 3D7, assembly ASM276v2.

2. The Plasmodium vivax reference genome used was P. vivax Sal-1, assembly ASM241v2.

3. Details of the Plasmodium falciparum samples used in the population genetics analyses can be found under the study accession number PRJEB2136. Accession numbers for each sample are described in Table S1.

4. The sequence of Plasmodium praefalciparum can be found under accession number SAMEA2464702.

Click here for additional data file.

Click here for additional data file.