Diverse variola virus (smallpox) strains were widespread in northern Europe in the Viking Age
Viking smallpox diversity
Humans have a notable capacity to withstand the ravages of infectious diseases. Smallpox killed millions of people but drove Jenner's invention of vaccination, which eventually led to the annihilation of this virus, declared in 1980. To investigate the history of smallpox, Mühlemann et al. obtained high-throughput shotgun sequencing data from 1867 human remains ranging from >31,000 to 150 years ago (see the Perspective by Alcamí). Thirteen positive samples emerged, 11 of which were northern European Viking Age people (6th to 7th century CE). Although the sequences were patchy and incomplete, four could be used to infer a phylogenetic tree. This showed distinct Viking Age lineages with multiple gene inactivations. The analysis pushes back the date of the earliest variola infection in humans by ∼1000 years and reveals the existence of a previously unknown virus clade.
Structured Abstract
INTRODUCTION
Variola virus (VARV), the causative agent of smallpox, is estimated to have killed between 300 million and 500 million people in the 20th century and was responsible for widespread mortality and suffering for at least several preceding centuries. Humans are the only known host of VARV, and smallpox was declared eradicated in 1980. The timeline of the emergence of smallpox in humans is unclear. Based on sequence data up to 360 years old, the most recent common ancestor of VARV has been dated to the 16th or 17th century. This contrasts with written records of possible smallpox infections dating back at least 3000 years and mummified remains suggestive of smallpox dating to 3570 years ago.
RATIONALE
Ancient virus sequences recovered from archaeological remains provide direct molecular evidence of past infections, give detail of genetic changes that have occurred during the evolution of the virus, and can reveal viable virus sequence diversity not currently present in modern viruses. In the case of VARV, ancient sequences may also reduce the gap between the written historical record of possible early smallpox infections and the dating of the oldest available VARV sequences. We therefore screened high-throughput shotgun sequencing data from skeletal and dental remains of 1867 humans living in Eurasia and the Americas between ~31,630 and ~150 years ago for the presence of sequences matching VARV.
RESULTS
VARV sequences were recovered from 13 northern European individuals, including 11 dated to ~600–1050 CE, overlapping the Viking Age, and we reconstructed near-complete VARV genomes for four of them. The samples predate the earliest confirmed smallpox cases by ~1000 years. Eleven of the recovered sequences fall into a now-extinct sister clade of the modern VARVs in circulation prior to the eradication of smallpox, while two sequences from the 19th century group with modern VARV. The inferred date of the most recent common ancestor of VARV is ~1700 years ago.
The number of functional genes is generally reduced in orthopoxviruses with narrow host ranges. A comparison of the gene content of the Viking Age sequences shows great contrast with that of modern VARV. Three genes that are active in all modern VARV sequences were inactive over 1000 years ago in some or all ancient VARV. Among 10 genes inactive in modern and Viking Age VARV, the mutations causing the inactivations are different and the genes are predicted to be active in the ancestor of both clades, suggesting parallel evolution. Fourteen genes inactivated in modern VARV are active in some or all of the ancient sequences, eight of which encode known virulence factors or immunomodulators. The active gene counts of the four higher-coverage Viking Age viral genomes provide snapshots from an ~350-year period, showing the reduction of gene content during the evolution of VARV. These genomes support suggestions that orthopoxvirus species derive from a common ancestor containing all genes present in orthopoxviruses today, with the reduction in active gene count conjectured to be the result of long-term adaptation within host species.
CONCLUSION
The Viking Age sequences reported here push the definitive date of the earliest VARV infection in humans back by ~1000 years. These sequences, combined with early written records of VARV epidemics in southern and western Europe, suggest a pan-European presence of smallpox from the late 6th century. The ancient viruses are part of a previously unknown, now-extinct virus clade and were following a genotypic evolutionary path that differs from modern VARV. The reduction in gene content shows that multiple combinations of active genes have led to variola viruses capable of circulating widely within the human population.
Sample locations and genome evolution of Viking Age variola virus.
(Top left) Sample VK533 with dagger, in situ. (Bottom left) Positive sample locations. Higher- and lower-coverage samples have solid and open circles, respectively. (Right) Diverse gene inactivation patterns in the four higher-coverage 600–1000 CE Viking Age sequences, within clade and as compared to sequences from other human VARVs and the phylogenetically nearest animal poxviruses, camelpox and taterapox. Open circles show intact genes, circle colors (independent for each gene) indicate different inactivations. Gene names at the bottom.
Abstract
Smallpox, one of the most devastating human diseases, killed between 300 million and 500 million people in the 20th century alone. We recovered viral sequences from 13 northern European individuals, including 11 dated to ~600–1050 CE, overlapping the Viking Age, and reconstructed near-complete variola virus genomes for four of them. The samples predate the earliest confirmed smallpox cases by ~1000 years, and the sequences reveal a now-extinct sister clade of the modern variola viruses that were in circulation before the eradication of smallpox. We date the most recent common ancestor of variola virus to ~1700 years ago. Distinct patterns of gene inactivation in the four near-complete sequences show that different evolutionary paths of genotypic host adaptation resulted in variola viruses that circulated widely among humans.
Variola virus, the causative agent of smallpox, is estimated to have killed between 300 million and 500 million people in the 20th century and was responsible for widespread mortality and suffering for at least several preceding centuries (1) (see the supplementary text section of the supplementary materials). Smallpox drove the 1796 development of vaccination by Edward Jenner, who later declared that “it now becomes too manifest to admit of controversy, that the annihilation of the Small-pox, the most dreadful scourge of the human species, must be the final result of this practice” (2). His prophecy was fulfilled in 1980, when the disease was declared eradicated by the World Health Organization, following a worldwide vaccination campaign. The virus is now exclusively stored in laboratories in the United States and Russia. Despite the eradication of smallpox, there are ongoing concerns regarding the reemergence of a smallpox-like disease via accidental or deliberate reintroduction of variola virus, adaptation of monkeypox virus to humans, or zoonosis or genetic engineering of another orthopoxvirus (3–5). Specific information regarding the evolutionary history of variola virus is therefore of substantial interest.
Orthopoxviruses are a genus of the Poxviridae and have large, linear, double-stranded DNA genomes (6). They differ in the range of mammalian host species they infect and the severity of the disease they cause: variola virus (VARV) and camelpox virus (CMLV) have a narrow host range and can be highly virulent, whereas taterapox virus (TATV), which infects gerbils, is only known to cause unapparent infections (7, 8). Cowpox virus (CPXV), monkeypox virus (MPXV), and vaccinia virus (VACV) infect several host species, with varying severity of disease. VACV, whose origin and reservoir host are unknown (9), generally has low virulence in humans, elicits an immune response that is protective against VARV, and was used as the vaccine during the smallpox eradication campaign (8). VACV infection of mice is a useful model for studying orthopoxvirus gene function in vivo. Orthopoxvirus genomes are typically between ~186,000 and ~228,000 nucleotides (nt) long (6). Conserved genes, necessary for virus transcription and replication, are located in the central ~100,000-nt region of the genome and are commonly used in phylogenetic inference (10). In vitro and in vivo studies show that orthopoxvirus genomes also contain many genes that are important for modulating host innate immunity and determining host range but whose inactivation does not prevent virus replication (11, 12). Those genes are predominantly located near the genome termini and vary between orthopoxvirus species. The factors governing host range and virulence of orthopoxviruses are complex and not fully understood. For example, some genes that promote VACV virulence are inactivated in the more virulent VARV (13–15), whereas in other cases the loss or inactivation of host immune system–modulating genes in VACV can result in increased virulence (12). In general, orthopoxvirus species with a narrow host range have fewer genes than those with broad host ranges, for example, compare VARV (~162 functional genes) with viruses of the CPXV species (~209 functional genes) (6). No orthopoxvirus species has a gene exclusive to that species, and all genes present across all orthopoxvirus species are also present in CPXV. It has therefore been proposed that, by a process of gene inactivation and loss, the different orthopoxvirus species arose from a common ancestor containing a full set of CPXV genes (6).
The timeline of the emergence of smallpox in humans is unclear. Analysis of sequences from a 17th-century Lithuanian mummy (VARV-VD21) (16), two Czech museum specimens from the 19th and 20th centuries (17), and VARV from the 20th century sampled during the eradication campaign dates their most recent common ancestor (MRCA) to the 16th or 17th century (10, 16–18). However, written records of possible smallpox infections date back to at least 3000 years ago (ya), and mummified remains suggestive of smallpox date to 3570 ya (1, 19). Ancient virus sequences recovered from archaeological remains provide direct molecular evidence for past infections and can reconcile the discrepancy between the written historical record of possible early infections and the time to the oldest available genetic sequences.
Recovering ancient variola virus DNA
To investigate the evolutionary history of VARV, we screened high-throughput shotgun sequencing data from skeletal and dental remains of 1867 humans living in Eurasia and the Americas between ~31,630 and ~150 ya. Using Kraken (20) on the microbial, protozoan, and viral genomes from the National Center for Biotechnology Information Reference Sequence (NCBI RefSeq) databases (21–23), we identified a total of 26 samples with between 1 and 730 reads assigned to VARV (fig. S1 and table S1). On the basis of the availability of additional sampling material, 13 of those 26 samples were then enriched for viral sequences using a targeted capture approach, as described in (24). Eleven of the 13 individuals are dated to 603–1050 CE, overlapping the Viking Age (793–1066 CE), and are from northern Europe, western Russia, or the United Kingdom, and two are from the 19th century from western Russia (Fig. 1A, Table 1, and supplementary text). All 13 samples showed evidence of VARV infection after capture, with the following characteristics supporting authenticity of the recovered reads: elevated rates of mismatching nucleotides (due to postmortem DNA damage) toward the read termini when aligned against a reference sequence (25), consistent with their respective ages and sequencing depth (fig. S2); mapped reads distributed evenly across the reference genome (fig. S3); lowest edit distances (the number of nucleotide differences in a read alignment) to VARV (19th-century samples) or to either CMLV or TATV (Viking Age samples) when comparing read mapping against eight diverse orthopoxvirus reference genomes (fig. S4); and 100% orthopoxvirus specificity when mapping reads against the NCBI nucleotide database (table S2 and supplementary text).
Fig. 1 Geographic location and phylogenetic placement of samples with reads matching VARV.
(A) Map showing geographic locations of samples. The locations of the higher-coverage and lower-coverage samples are indicated with solid and open red circles, respectively. The location of the published VD21 sample is shown with a blue circle. Gray dots indicate locations of screened samples in which VARV was not detected. (B) Condensed maximum likelihood tree showing the higher-coverage samples and the placement of the lower-coverage samples. Lower-coverage samples were placed onto the tree containing the higher-coverage samples using EPA-ng (27). Clades that do not have any lower-coverage samples placed on them are collapsed and indicated with black circles. Color is used to indicate the branch with the highest likelihood weight ratio assigned by EPA-ng for each lower-coverage sample. Support values are based on 1000 bootstrap replicates. The full tree is shown in fig. S6. Complete phylogenetic trees, showing the placement of each sample, individually and combined, are shown in figs. S7 and S8, respectively. s/s, substitutions per site.
| Sample | Date (CE) |
Site | Sex | Age (years) |
Sample type |
TATV reads |
TATV genome coverage (%) |
TATV mean coverage depth |
VARV reads |
VARV genome coverage (%) |
VARV mean coverage depth |
|---|---|---|---|---|---|---|---|---|---|---|---|
| VK382 | 640–770 | Öland, SWE | M | >15 | Tooth | 85,036 | 91.079 | 40.723 | 88,324 | 99.564 | 45.188 |
| VK388 | 603–653 | Nordland, NOR | M | 12–17 | Tooth | 20,119 | 90.396 | 7.082 | 20,629 | 98.999 | 7.755 |
| VK470 | 950–1000 | Gnezdovo, RUS | F | 25–35 | Tooth | 18,894 | 89.522 | 6.944 | 19,117 | 97.114 | 7.500 |
| VK281 | 885–990 | Zealand, DEN | M | 30–35 | Tooth | 10,801 | 88.358 | 4.594 | 11,006 | 95.885 | 5.012 |
| VK168 | 880–1000 | Oxford, UK | M | 16–25 | Petrous | 8,876 | 77.065 | 2.397 | 8,798 | 84.008 | 2.553 |
| VK533 | 800–1050 | Öland, SWE | * | 50–60 | Tooth | 590 | 18.694 | 0.211 | 610 | 20.712 | 0.233 |
| VK515 | 664–768 | Nordland, NOR | M | 18–22 | Tooth | 648 | 12.143 | 0.149 | 639 | 12.975 | 0.157 |
| VK255 | 950–1000 | Gnezdovo, RUS | F | 20–25 | Tooth | 109 | 2.572 | 0.027 | 101 | 2.583 | 0.027 |
| FIN1 | 1800–1900 | Varzino, RUS | F | 40–45 | Tooth | 83 | 1.869 | 0.021 | 95 | 2.275 | 0.025 |
| KHA1 | 1800–1900 | Yamalo-Nenetskiy, RUS | N/D | 3–4 | Petrous | 58 | 1.763 | 0.018 | 58 | 1.932 | 0.019 |
| VK443 | 800–1050 | Öland, SWE | M | 20–23 | Tooth | 88 | 1.719 | 0.018 | 85 | 1.767 | 0.018 |
| VK108 | 800–1000 | Malmö, SWE | F | 55–75 | Tooth | 69 | 1.413 | 0.015 | 70 | 1.537 | 0.017 |
| VK138 | ~1000 | Fyn, DEN | M | 25–35 | Tooth | 49 | 1.327 | 0.013 | 49 | 1.389 | 0.014 |
Table 1 Overview of samples with reads matching orthopoxviruses.
From left to right, columns give the sample name, approximate sample date based on either carbon-14 (14C) dating (VK382, VK388, and VK515) or archaeological context (all other samples), the site where the sample was found (DEN, Denmark; UK, United Kingdom; NOR, Norway; RUS, Russia; and SWE, Sweden), the sex of the individual, the age of the individual at the time of death, the type of tissue that was sequenced, the number of reads that aligned to TATV (GenBank accession no. NC_008291.1) using Bowtie2 local alignment and with mapping quality of at least 30, the percentage of the TATV genome covered (fig. S3), the mean TATV genome coverage depth, the number of reads that aligned to VARV (GenBank accession no. NC_001611.1) using Bowtie2 local alignment with mapping quality of at least 30, the percentage of the VARV genome covered, and the mean VARV genome coverage depth (fig. S3). Rows are ordered by decreasing mean coverage depth (table S1). The aVARV sequences do not correspond exactly to VARV (185,578 nt) or TATV (198,050 nt), these are just the most logical similar modern reference sequences to align against (table S4). In fact, the aVARV-VK382 consensus contains 192,255 nt and is 99.98% complete, with only 41 unresolved nucleotides. Eleven of the samples date from the Viking Age (sample names starting with “VK”) or just before it and are from Scandinavian countries, western Russia, or the United Kingdom. The 14C (2σ) age ranges for VK382, VK388, and VK515 were calibrated with IntCal13 (February 2020) (82) with no correction made for possible marine reservoir effect. Site descriptions, with additional sample information, are given in the supplementary text. N/D, not determined.
*VK533 was identified as female by osteology but is genetically male.
In four samples (VK281, VK382, VK388, and VK470; hereafter referred to as higher-coverage samples), alignment against VARV (accession no. NC_001611.1) resulted in a mean coverage depth of 5.01× to 45.19× (Table 1 and fig. S3) and consensus sequences covering 95.89 to 99.56% of the VARV genome, with the sequence for VK382 being 99.98% complete (Table 1 and table S3). On the basis of nucleotide sequence identity, the consensus sequences for the four higher-coverage samples were collectively most similar first to each other and then to TATV (table S4). For the remaining seven “lower-coverage” Viking Age samples, read edit distances were lowest against the assembled VK382 consensus sequence, further supporting their authenticity (fig. S5).
Phylogenetic analysis reveals an unknown variola virus clade
Before inferring phylogenetic trees on the basis of the conserved central region of the genomes from the four higher-coverage sequences, we tested whether modern VARVs show evidence of genomic regions that descend from a virus population containing sequences resembling these four higher-coverage ancient VARV (aVARV) sequences and whether those sequences are themselves recombinants. We examined the higher-coverage sequences and seven representative sequences from related modern orthopoxvirus species for evidence of recombination using seven algorithms from the RDP4 toolset (26). A 125-nt noncoding region of aVARV-VK470 was flagged as possibly attributable to recombination with an unknown sequence, but the evidence of recombination was uncertain. No other recombination event was suggested in the higher-coverage aVARV sequences, and no modern VARV sequence in the selection was suggested to be a recombinant derived from any aVARV sequence (Materials and methods and table S5).
To establish the phylogenetic placement of the four sequences from the higher-coverage samples, we inferred a maximum likelihood (ML) tree including 84 sequences chosen to represent the full orthopoxvirus diversity. In the resulting tree (Fig. 1B and fig. S6), the sequences form a now-extinct monophyletic clade in sister relationship with the clade consisting of all modern VARV and VARV-VD21 (collectively referred to as mVARV). We then placed partial consensus sequences of the lower-coverage samples onto the inferred tree using the evolutionary placement algorithm EPA-ng (27). Four of the seven Viking Age samples (VK108, VK168, VK515, and VK533) are confidently placed within the aVARV clade (Fig. 1B and figs. S7 to S12). Their placement position suggests a common ancestor with the more recent VK281 and VK470 lineages, supporting the existence of distinct Viking Age lineages (fig. S8). The placements of the other samples show some uncertainty but are consistent with their respective sample ages: the Viking Age samples (VK138, VK255, and VK443) in the vicinity of the aVARV clade, and the 19th-century samples (FIN1 and KHA1) within the mVARV clade. The placement of the lower-coverage sequences provides further evidence of their authenticity and independently confirms the presence of a separate aVARV clade, as already shown by the higher-coverage samples alone.
The existence of a temporal signal in the sequence data is a precondition for inferring dated phylogenetic trees (28). A linear regression of sample ages and root-to-tip distances from an ML tree that includes sequences from the four higher-coverage samples and mVARV (fig. S13) shows clear evidence of clocklike evolution (fig. S14), in agreement with earlier results (16, 18). Consequently, we proceeded to infer dated coalescent trees using BEAST2 (29) under six different evolutionary models (table S6). With the best-fitting model, the MRCA of the combined mVARV and aVARV sequences is dated to 1.7 thousand years ago (kya) [95% highest priority density interval (HPD95): 2.2 to 1.4 kya], and the MRCA of the aVARV clade is dated to 1.5 kya (HPD95: 1.7 to 1.4 kya). The MRCA of the mVARV clade is dated to 0.44 kya (HPD95: 0.6 to 0.35 kya), and to 0.31 kya (HPD95: 0.43 to 0.21 kya) for the mVARV clade without VARV-VD21, both of which are in line with previous molecular dating analyses (10, 16–18) (table S7 and fig. S15).
The ancient variola viruses have a diverse pattern of gene inactivation
Many genes that are active in other orthopoxviruses are inactivated (truncated, fragmented, or completely absent) in mVARV, possibly owing to host adaptation (6, 13). We therefore compared and contrasted the gene status of sequences from the four higher-coverage aVARV samples, mVARV, plus CMLV and TATV (the orthopoxviruses most closely related to aVARV) (Fig. 2, fig. S16, table S8, and supplementary text). Gene status analysis allows a detailed comparison of interclade and intraclade variation at the level of possible phenotypic differences. Gene differences are of interest individually and in combination, and the overall pattern of differences in gene status complements the whole-genome diversity summarized by phylogenetic analysis.
Fig. 2 Gene-inactivating mutations in mVARV, aVARV, CMLV, and TATV sequences, with inferred gene status for internal nodes.
The cladogram on the left shows the topology of the phylogenetic relationship between mVARV, the aVARV sequences we report in this paper, TATV, and CMLV (see trees in Fig. 1B and figs. S13 and S15). Each row to the right indicates gene status (present, inactivated) in the virus at that level in the cladogram (bottom eight rows), or in an inferred internal node or the root (top seven rows). The 40 columns represent genes that are either absent or that have an inactivating mutation in at least one of the mVARV or aVARV sequences, excluding 19 genes that are absent in both mVARV and aVARV. Genes are sorted left-to-right by category and then by position in the genome. Empty circles indicate genes assumed to be present and functional. Colored circles indicate genes with a gene-inactivating mutation, genes that are absent, or genes that do not have coverage in the ancient sequences (also indicated by an X). Within a column, genes with identical gene-inactivating mutations are shown in the same color. Color correspondence between columns carries no meaning. Genes with partial coverage (in the aVARV sequences) are marked with a black dot. A horizontal bar across a filled circle indicates a gene-inactivating mutation considered uncertain, because of less than three reads covering the position of the mutation. A horizontal bar across an empty circle indicates a gene that may be functional (with a length intermediate between its length in viruses where it is functional and viruses where it is not). Gene labels indicate either the name of the VACV-Cop homolog or the final offset of the gene in CPXV-Gri or CPXV-Ger (6). We inferred gene status at internal nodes using the DELTRAN maximum parsimony algorithm implemented in PastML (47). Genes are organized into four groups: group A, inactivated in all mVARV and all aVARV sequences; group B, inactivated in mVARV, but present in some or all aVARV sequences; group C, present in mVARV, but inactivated in some aVARV sequences; and group D, present in mVARV and aVARV (not shown in this figure). See supplementary text for further information and the description of subgroups.
Interclade differences in gene status can be logically organized into four groups (fig. S17): group A, inactive in mVARV and aVARV clades; group B, inactive in mVARV but active in at least one aVARV sequence; group C, active in mVARV but inactive in at least one aVARV sequence; and group D, active in both clades. Figure 2 shows details of groups A to C, highlighting gene differences between clades.
In group A (genes inactive in both clades), comparing the exact position of gene-inactivating mutations allows us to infer the likely status of the genes in the common ancestors of mVARV and aVARV. Five genes in group A [B23R;B24R/C17L;C18L, 26765, A25L, A37R, and B16R; nonitalic numeric values correspond to unnamed genes with offsets as given in (6)] have identical inactivating mutations in both clades (Fig. 2 and supplementary text). Identical inactivating mutations among orthopoxviruses not in sister relationship are rare: The 17 representative orthopoxviruse sequences included in our gene-loss analysis (Materials and methods) have 53 genes with gene-inactivating mutations in viruses not in sister relationship, but only two of these have the mutation in the same location. This suggests that the five genes listed above were already inactivated in the common ancestor of mVARV and aVARV (Fig. 1B). An additional 10 genes (20088, A39R, A44L, A53R, B2R;B3R, 215231, C9L, A9L, A52R, and 202711) are inactivated in both clades, although in different ways, and are thus projected to be active in the ancestor of mVARV and aVARV (fig. 2). This suggests parallel virus evolution after the divergence of the clades ~1.7 kya (fig. S15).
A gene inactivated in both mVARV and aVARV may have induced a similar outcome after infection of humans. The VACV gene encoding a soluble interleukin-1β (IL-1β) receptor is inactive in VACV strain Copenhagen (VACV-Cop) (gene B16R) but functional in VACV strain Western Reserve (VACV-WR; where it is denoted B15R) (30). Infection of mice with either VACV-Cop or VACV-WR in which B16R was inactivated resulted in a febrile response (15). Conversely, infection with either virus expressing a functional IL-1β receptor prevented the induction of fever (15). Human infections with VARV, which all have a disrupted B16R gene, are accompanied by high fever (8). The identical inactivation of B16R in the aVARV clade may suggest that the symptoms of the disease caused by those viruses also included high fever and that the inactivation occurred before the divergence of the clades.
Group B (genes inactivated in mVARV but present or with undetermined status in some or all of the aVARV sequences) contains 14 genes. Eight of these (B7R, CrmE, C2L, K1L, F3L, A35R, A40R, and A55R) encode known virulence factors or immunomodulators in VACV, three of which (C2L, F3L, and A55R) are members of the kelch-like family (supplementary text) (12). Deletion of the kelch-like proteins C2L or A55R in VACV-WR leads to increased lesion size in intradermally infected mice (31, 32), while the virulence of VACV lacking F3L is reduced in the same model (33). Deletion of kelch-like proteins has also been proposed to reduce the host range and virulence of CPXV-Gri in vitro (34). Also in contrast with mVARV, aVARV-VK382 and aVARV-VK388 encode a predicted functional guanylate kinase (A57R). This gene is otherwise only active in CPXV and appears to have been lost independently in CMLV, TATV, mVARV, aVARV-VK281, and aVARV-VK470 (supplementary text), suggesting that the gene was functional in their common ancestors.
Group C (genes inactivated in at least one aVARV but not in any mVARV sequence) contains seven genes (CrmB, C1L, E5R, F10L, I8R, B20R, and 210863, although only the first three have unequivocal gene-inactivating mutations). The effect of the inactivation of C1L and E5R in VACV is unknown. CrmB is a tumor necrosis factor receptor homolog and chemokine binding protein (35). Interestingly, mVARVs have a functional CrmB, whereas it is identically inactivated in aVARV-VK281 and -VK470, and the start codon in aVARV-VK388 is replaced by a valine. However, unlike in mVARV, CrmE, a tumor necrosis factor receptor homolog, and the chemokine binding protein B7R (35) are functional in aVARV-VK281, -VK470, and -VK388 and together may provide the same functionality as CrmB. The older aVARV-VK382 has functional versions of CrmB, CrmE, and B7R. The seven group C genes, already inactivated in the aVARV clade, show that these viruses had taken a different genetic (although perhaps not functional) evolutionary path than that of the mVARV sequences, in which none of these genes were inactivated during an additional ~1000 years of evolution.
The gene status analysis and molecular clock dating of the internal nodes allow us to investigate the temporal trend of gene inactivation within the combined mVARV and aVARV clades. We observe an inactivation of genes over time (Fig. 3), supporting suggestions that orthopoxvirus species may derive from a common ancestor containing all genes found in orthopoxviruses today (6). A simple linear model of gene loss suggests that an ancestor of aVARV and mVARV with a gene content similar to CPXV would have existed ~4 kya (Fig. 3). But any such model will necessarily involve unsupported and currently untestable assumptions, some of which will almost certainly be invalid, so caution is warranted (supplementary text). Loss of host-range genes has been correlated with an increase in host specificity (11), suggesting that the aVARVs and the common ancestor of mVARV and aVARV may have had a broader host range than mVARV. However, host range in poxviruses is likely to be a multigenic trait that depends on the cooperation of several genes, as is the case in myxoma virus (36) and VACV (37–41). The higher number of active genes thus does not necessarily imply that the aVARVs had a host tropism for animals other than humans, although this is a possibility.
Fig. 3 Gene inactivation over time.
The number of absent genes and genes with inactivating mutations is shown on the y axis, plotted against time (in years into the past) on the x axis. The evolutionary relationships between internal nodes and the tips of the phylogenetic tree (blue dots), as given by the ML trees in Fig. 1B and fig. S13, are indicated by black. Point estimates for dates for the aVARV sequences and the internal nodes are taken from the dated coalescent tree in fig. S15, as estimated using BEAST2 (29) with a relaxed log-normal clock and a Bayesian skyline population prior. The x-axis placement of the VARV dot corresponds to the inferred ancestor of all modern VARV. Gene counts for internal nodes are inferred on the basis of gene status in their descendant viruses (see Fig. 2). Orange dots plotted vertically below VK281 and VK470 indicate the number of gene-inactivating mutations that can be identified with certainty, as well as genes with no coverage in the ancient samples, and genes where the presence of a gene-inactivating mutation is uncertain (Materials and methods). The dashed dark-gray line and surrounding gray area show a linear least squares regression and 95% confidence region, with gene counts from the inferred internal nodes included in the calculation. The x-axis intercept (not shown) is 4017 years into the past. The other parameters of the regression are: slope = 0.014, R (correlation coefficient) = 0.91, R2 (coefficient of determination) = 0.84, y intercept = 58.8, Stderr (standard error of estimated gradient) = 0.002, and P = 0.0001. For conciseness, “aVARV-” has been omitted from the Viking Age sample names (aVARV-VK281, aVARV-VK382, aVARV-VK388, and aVARV-VK470), and “VARV-” has been omitted from VARV-VD21. See supplementary text for comments regarding modeling gene inactivation over time.
Within the aVARV clade, we observe a diverse pattern of gene inactivation (Figs. 2 and 3). Eleven genes are unequivocally inactivated in some, but not all, of the sequences from the four higher-coverage samples (in group B: 22682, C2L, K1L, F3L, A35R, A40R, 168145, and A57R; in group C: CrmB, C1L, and E5R). An additional four genes are inactivated in all four viruses, but with differing gene-inactivating mutations (in group A: A39R, 193435*, C13L;C14L, and 202711), suggesting that these genes were lost independently (supplementary text). By comparison, among 50 available mVARV sequences, the median number of pairwise differences in gene activation status is zero (mean: 0.6). On the basis of gene status, aVARV-VK281 and aVARV-VK470 can be grouped together, as can aVARV-VK382 and aVARV-VK388, differing unequivocally in status in one (A35R) and five (A35R, A40R, CrmB, C1L, and E5R) genes, respectively (fig. S18), in agreement with the grouping based on their overall gene inactivation counts (Fig. 3), sample ages (Table 1 and Fig. 3), and phylogenetic positions (Fig. 1B and figs. S13 to S15). Our data suggest that VARV followed at least two independent trajectories of gene inactivation during its evolution: one leading to aVARV-VK281 and aVARV-VK470 and one to mVARV (Fig. 3). The geographical and temporal spread of the viruses in the aVARV clade (Fig. 1A and figs. S8 and S15) indicates that variola viruses with a gene composition pattern different from that of mVARV existed over a period of at least ~450 years and circulated widely among humans during the Viking Age.
Historical perspective
We found evidence of VARV in 11 of 525 individuals (~2%) predating and overlapping the Viking Age. If aVARV, like mVARV, did not lead to persistent infections (8) and if no viral particles remained in teeth or bone after an infection, these individuals all died during an acute infection. The 11 samples with reads matching VARV are likely not the only infections among those 525 individuals because: (i) The presence of viral DNA in ancient teeth and bones is certainly an imperfect diagnostic test of viral infection, with an elevated false negative rate due to differences in DNA preservation, source tissue, or sequencing depth in the screening; (ii) VARV is not thought to be viremic for the entire duration of the illness (8); (iii) On the basis of mVARV case fatality rates (up to ~30%), many individuals survived infection; (iv) Five Viking Age samples with reads assigned to VARV by Kraken were not available for capture, because of a lack of sample material. Moreover, the 525 individuals are drawn from an ~450-year period, during which the overall population of Scandinavia rose from ~750 thousand to ~1 million (42). If the overall population increased linearly during this period, and life expectancy at birth was ~30 years throughout, a total of ~13.1 million people would have lived in Scandinavia during this period, so a sample of 525 individuals is only ~0.004% of the estimated total population. Given this considerable uncertainty, it would not be prudent to use the 2% detection rate to estimate smallpox prevalence or case fatality rates during the Viking Age. The likely inaccurate (low) detection rate and the very small sample size argue against regarding the figure as a prevalence. Even if all other sources of uncertainty were resolved, it would still not be clear how the figure relates to case fatality rates during that period, because we cannot be sure that the individuals died as a result of their infections.
Hypotheses regarding the earliest history of VARV have been derived exclusively from often ambiguous historical accounts and from the visual examination of mummies dating from as early as 3570 ya (1, 19). The dates of the earliest likely presence of smallpox in specific regions are difficult to determine from written records and have in some cases been a matter of dispute (1, 8, 19, 43). Ancient DNA sequences can resolve such questions by providing clear evidence of VARV presence in time and space. The dating of the aVARV samples, from as early as 603 CE, matches that of multiple written accounts of likely smallpox infections in southern and western Europe from the late 6th century onward (1, 8, 19, 44, 45). Our finding of the virus in northern Europe at these times disproves various suggestions of first introductions involving later dates. For example, the introduction (or later establishment) of smallpox into Europe via returning Crusaders (8, 9) or by way of Spain during the Moorish invasion of 710 CE (9), or into England by Normans in 1241–1242 CE (46), as well as the claim that the virus had not reached northern Europe by 1000 CE (8). We did not find evidence for infections in samples from earlier time periods, despite screening 619 individuals predating the Viking Age (fig. S1). Together, the written record and the confirmed aVARV infections suggest a pan-European presence of the virus by the end of the Viking Age at the latest, as suggested for parts of southern and western Europe based purely on written records (1, 19, 44).
The Viking Age sequences reported here push the definitive date of the earliest VARV infection in humans back by ~1000 years and reveal the existence of a previously unknown, now-extinct virus clade. The ancient viruses were following a genotypic evolutionary path that differs from modern VARV. This is highlighted by at least three genes inactivated in some aVARV sequences but still active in the mVARV clade, and by 10 genes inactivated in both clades but with gene-inactivating mutations that differ between clades. The ancient viral genomes show reduction of gene content during the evolution of VARV and that multiple combinations of gene inactivations have led to viruses capable of circulating widely within the human population.
Materials and methods
Software package versions
The following software package versions were used, with default arguments unless otherwise noted. BEAST2 2.4.7 (29), PastML 1.9.20 (47), BLASTn 2.2.29+ (48), bModelTest 1.0.4 (49), Bowtie2 2.3.5.1 (50), Picard 2.20.2 (51), BWA 0.7.15-r1140 (52), DIAMOND 0.9.23 (53), EPA-ng 0.3.6 (27), Gappa 0.5.0 (54), gargammel (released without a version number, October 2017) (55), Geneious 9 (56), ggtree 1.16.6 (57), Kraken 1.0 (20), MAFFT 7.419 (58), mapDamage 2.0.8-2 (59), Mosdepth 0.2.6 (60), R 3.6.0 (61), RAxML-NG 0.7.0 BETA (62), RDP4 4.9.5 (26), SAMtools 1.9 (63), SciPy 1.4.1 (64), TempEst 1.5 (65), Tracer 1.6.0 (66), and TreeAnnotator 2.2.4 prerelease (29).
Screening and sample preparation
Using Kraken on a database of all microbial, protozoan, and viral genomes in the NCBI RefSeq databases (21–23), we analyzed whole-genome shotgun sequencing data from 1867 individuals for the presence of reads that best matched VARV. We identified a total of 26 samples with between 1 and 730 reads assigned to VARV (fig. S1) as candidates for enrichment of viral target sequences in library by capture. Of these, sufficient sample materials were available for capture processing for 13 samples. Uniquely dual-indexed sequencing libraries (2 to 6 libraries per sample) were prepared from 1 to 3 ancient DNA (aDNA) extracts per sample. Libraries were pooled into capture reactions containing 1183 to 1880 ng per reaction. Hybridization and sequencing on HiSeq4000 were performed as described previously (24). An additional round of sequencing was performed on pools of enriched libraries, using the Illumina NovaSeq6000 (S4 flowcell, XP workflow, 2 × 100 cycles). Read length and insert size distributions are shown in fig. S19.
Authenticity
The following points argue for the authenticity of the recovered DNA sequences.
1) Standard precautions (67) and recommended methodologies (68) for working with ancient DNA were applied.
2) Clear damage patterns characteristic of DNA sequences recovered from archaeological remains are present in all samples with sufficient viral reads available. Note that the absence of a damage pattern in a sample with few reads does not (by itself) indicate that the sample reads are not ancient. The human reads all show clear damage patterns in all samples, showing that the nucleic acids present in the samples are indeed ancient.
3) The ancient samples included in the study have reads matching across the entire TATV reference genome (fig. S3). This is in stark contrast to false positives that regularly occur with human DNA sequences in which there are matches of only one or two poxvirus proteins out of ~160 to 220, in which cases the reads in fact are better matches for similar-but-unrelated sequences, e.g., human mRNA.
4) Reads matching TATV were matched against the entire NCBI nucleotide database using BLASTn. Matches were found in the nucleotide database for 118,049 of 118,118 reads (99.94%), all of which had a poxvirus as the best match, ignoring matches against artificial constructs (table S2). A small number of reads with no BLASTn match in the nucleotide database is expected (see Read poxvirus specificity section).
5) The edit distances of the reads are fully consistent with the placement of (full or partial) consensus sequences in phylogenetic trees. Reads from all nine aVARV clade samples very clearly have lowest edit distances against TATV and CMLV, the intermediate samples (VK138 and VK255) have read edit distances that are closest to VARV, TATV, or CMLV, and the two modern samples (FIN1 and KHA1) have read edit distances very clearly closest to VARV. As would also be expected, edit distances against our aVARV-VK382 sequence are the lowest of all for the Viking Age samples.
6) The fact that the almost-complete genome consensus sequences from the four higher-coverage ancient samples form a self-consistent separate basal clade indicates that these sequences are very unlikely to be due to four independent modern contaminations. If, for example, 20th-century smallpox sampling entirely missed a still-circulating clade that these samples actually came from, those lineages would need to have evolved with a far lower substitution rate than other variola viruses because of the relatively very short branch lengths leading to the ancient sequences in the ML tree (fig. S6).
7) The DNA sequenced in these samples was drilled out of the interior of well-preserved teeth (11 samples) and petrous bone (two samples). Pox virus DNA was therefore very likely circulating in some form in the bloodstream of the individual at the moment of death. Subsequent contamination of any form of virus or viral DNA before aDNA extraction in our dedicated clean laboratory is unlikely.
8) The EPA-ng placements from the lower-coverage samples provide important information via the phylogenetic placement of the (incomplete) consensus sequences from those samples. The fact that all ancient samples are placed consistently (according to sample ages) within or basal to the aVARV or mVARV clades is further evidence of the overall authenticity and consistency of the reads. We also show that reads simulated in silico from a different orthopoxvirus genome that was not included in the tree are consistently placed on their expected branch, not with the aVARV or mVARV clades (fig. S12).
9) Present-day analysis of purported VARV reads has the advantage that contamination by a modern variola virus can be definitively ruled out, given the nature of the associated disease. No member of the analysis team has any poxvirus infection and the facilities where the sample extraction, library preparation, and sequencing was done is not a virology laboratory (rather it is focused on ancient nonmicrobial host genome recovery) and has never been involved in processing modern viral material. The paper author (L.V.) who carried out all the targeted laboratory work received a routine VARV vaccination before 1980.
Read poxvirus specificity
To check that sample reads found to match poxvirus sequences were not in fact better matches for other sequences, reads matching TATV (GenBank accession no. NC_008291.1) using Bowtie2 end-to-end matching and with mapping quality (MAPQ) ≥ 30 were checked for specificity by matching them against the NCBI nucleotide database (downloaded during the week of 13 January 2020) using BLASTn (-task blastn -evalue 0.01). The best-matching database sequence for each read was checked to see if it corresponded to a poxvirus. Matches against bacterial artificial chromosomes and synthetic constructs [matching the case-insensitive regular expression (BAC cloning vector|synthetic construct)] were ignored. Table S2 shows that, apart from 69 reads (from a total of 118,118) that did not match anything in the nucleotide database, all reads had a poxvirus as their best match, as determined by the nucleotide database subject title matching the case-insensitive regular expression (pox|ectromelia|variola|akhmeta|abatino|vaccinia).
The end-to-end option to Bowtie2 and the MAPQ ≥ 30 filtering were applied to examine only reads that mapped in full and with a high confidence regarding their genome match position. This is necessary because orthopox viruses contain repetitive low-complexity regions near their termini, and such sequences often match many hundreds of host chromosome sequences from a diverse range of unrelated organisms. These matches sometimes exceed the default 500 matches returned by BLASTn, meaning no orthopox match may be reported at all, despite a perfect match with an orthopox genome. The stricter matching (in particular, only considering reads matching with MAPQ ≥ 30) eliminates low-complexity viral reads with these ambiguous placements from consideration.
It is expected that some orthopox reads matched by Bowtie2 will not match anything with BLASTn. This is due to search algorithm variations in seeding, match scoring values and thresholds, and reporting thresholds. For example, despite nucleotide identity exceeding 90%, BLASTn (-task blastn) will not match a sequence against a copy of itself with every 11th nucleotide mutated, because a perfectly matching seed region of 11 nucleotides is required before a region is further considered as a potential match. A different search algorithm will identify such matches. These algorithmic idiosyncrasies are the reason we also examined read matches found by DIAMOND and BWA as part of the consensus and gene-level analysis, with both analyses subject to manual visual inspection of aligned reads.
Investigation of gene-inactivating mutations
To assemble the reads used for the generation of consensus sequences and the investigation of gene-inactivating mutations, reads from the capture were mapped against orthopoxvirus reference genomes using DIAMOND, matching against version 13.0 of the RVDB protein database (69), and BWA [using both the mem and aln (with and without the −l option set to 1000) algorithms, with other arguments left at their defaults], matching against a collection of orthopox genomes from NCBI (GenBank accession nos. AY009089.1, AF438165.1, AF482758.2, X94355.2, AF012825.2, AF380138.1, DQ437594.1, AM501482.1, M35027.1, AY243312.1, L22579.1, X69198.1, Y16780.1, and KY358055.1). Reads from all mapping algorithms were combined and de-duplicated using Picard after mapping to a TATV genome (GenBank accession no. NC_008291.1) using Bowtie2.
To assess gene content of the higher-coverage aVARV samples in the absence of a reference genome, we investigated the presence of gene-inactivating mutations in the ancient samples by comparison of alignments made to 17 orthopoxviruses. Information about the offset locations of 219 homologous genes in 17 orthopoxvirus reference genomes (dataset 3; see supplementary text) is presented in supplementary table 1 of (6). These offsets were used to extract the sequences of the homologous genes in the 17 reference sequences. Excluding cases where a gene is absent in a reference sequence, we extracted 219 × 17 − 301 = 3422 gene sequences. Supplementary table 1 in (6) was also used to acquire information on status (presence, absence, fragmentation, and truncation) of genes in the 17 reference sequences. We aligned reads independently at the gene level with sequences from this broad set of references to make it possible to correctly deal with cases where, for example, two genes from an ancient sequence matched one gene best in one modern reference (e.g., CMLV) and another gene best in a different modern reference (e.g., TATV), or where a gene is present in the ancient virus but completely absent in the modern reference. Aligning against just one of the modern full-genome references would necessarily risk potential error. We considered this the only defensible approach, given the fact that the gene composition of the ancient sequences was unknown.
For each gene, all reads from an ancient sample were aligned to one or more of the 17 reference gene sequences using Geneious, medium sensitivity, fast parameters, and then visually inspected for the presence of gene-inactivating mutations. The reference(s) to align against were chosen in the following way: Reads for each ancient sample were matched against each reference gene sequence using BLASTn, and a consensus sequence for each gene was generated using BWA (aln algorithm) and SAMtools. The reference gene sequence with >50% coverage of the gene sequence and the highest median bit score was chosen as the reference. If this was inconclusive, additional references were used, usually one of CPXV-Gri, -Ger, or -BR. For each ancient sample, gene status, coverage, and reference sequence are given in table S8. Gene-inactivating mutations in the ancient sequences were classified as either “certain,” “uncertain,” or “uncertain if functional” according to the following criteria.
A “certain” gene-inactivating mutation: There is a mutation that leads to a stop codon (insertion, deletion, or point mutation) that is covered by at least three reads, and the reference against which the ancient reads are aligned is also truncated or fragmented in an identical fashion. An “uncertain” gene-inactivating mutation: There is a mutation that leads to a stop codon that is covered by less than three reads.
An “uncertain if functional” gene: A gene of intermediate length—shorter than a gene that is functional in some references but substantially longer than inactivated versions of the same gene in others (see gene A9L).
The ancient gene status was further confirmed by mapping the reads, consensus sequences resulting from the Geneious alignments against the gene reference(s), and the contigs from a de novo assembly of the reads (also performed in Geneious, default parameters) against the CPXV-Ger reference genome (GenBank accession no. DQ437593.1) and inspecting the resulting consensus for the expected genes. The status of each gene in the resulting consensus sequence agreed with the analysis performed above.
Ten genes [1393 and 222274 (C23L/B29R), 2286 and 221381 (CrmB), 3420 and 220247 (C19;C20L/B25/B27R), 5392 and 218275 (C17L;C18L/B23R;B24R), 7652 and 216015 (C16L/B22R)] in the orthopoxvirus genome are occasionally diploid. While both versions of these genes are present in CPXV-Ger, -Gri, and -BR, only one version of the gene is present in VARV. Given low genome coverage and short DNA template fragments, it is not clear how to determine whether these genes are present at the 3′ or 5′ (or both) ends of the aVARV genomes. The (potentially) two copies of the gene were thus treated as one in the analyses described above.
To infer the gene status at internal nodes of the phylogeny of mVARV, aVARV, CMLV, and TATV, we used the DELTRAN algorithm implemented in PastML, to perform a maximum parsimony ancestral state reconstruction. Gene-inactivating mutations were determined to be identical in two tips if the same mutation leading to an inactivation of the gene was present, whether or not a different upstream mutation may have previously ablated the gene in one of them.
We compared differences in gene activation status using a larger set of modern VARV sequences. The Virus Orthologous Cluster database contains (among other orthopoxviruses) the gene annotations of 50 modern VARV (70). Each gene is assigned to a cluster on the basis of BLASTp results and manual curation, genes within each cluster are orthologous. We determined correspondence between the assignments in (6) and the orthologous clusters on the basis of sequence identity. Only clusters that had corresponding genes assigned by Hendrickson et al. (6) were included. Genes in the Virus Orthologous Cluster database were deemed functional if they were at least 80% of the length of the corresponding functional gene in (6). To compare the gene content of the 50 modern VARV sequences, we counted the number of differences in gene status between two sequences. Gene status was considered different if the gene was present in one sequence and had a gene-inactivating mutation or was absent in the other.
Generation of consensus sequences
Reads from aVARV-VK382, the sample with the highest coverage, were mapped against VARV strain SLN68_258 (GenBank accession no. DQ441437.1) in Geneious using Medium sensitivity/fast and iterate up to five times. The alignment was visually inspected, and sections not present in the reference sequence were determined by iteratively mapping the reads against the new consensus sequence from either side of the insertion, until the consensus sequences from the mapped reads overlapped. The resulting consensus had 41 undetermined positions out of 192,255 nucleotides. Gene status and position in the genome order were both fully consistent with the analysis performed in the previous section. Consensus sequences for the remaining samples were generated by mapping their reads against the aVARV-VK382 consensus using Geneious using Medium sensitivity/fast and iterate up to five times. Manual inspection showed that all expected genes were present in the ancient consensus sequences, consistent with the gene loss analysis (table S3). The consensus sequences were also consistent with contiguous sequences (contigs) resulting from a de novo assembly of the reads (performed in Geneious, using default parameters). We called majority rule consensus sequences at a depth of one read, taking into account residue quality.
Recombination analyses
Seven algorithms, as implemented in the RDP4 package, were used to search for evidence of recombination in a selection of 11 orthopoxvirus sequences. The algorithms were: 3Seq (71), BootScan (72), Chimaera (73), GENECONV (74), MaxChi (75), RDP (76), and SiScan (77). The sequences (with GenBank accession numbers in parentheses) were: Akhmeta (MH607143), CMLV (AY009089), CPXV (KC813492), CPXV (KC813508), TATV (NC_008291), VARV-VD21 (BK010317), VARV (NC_001611), aVARV-VK281, aVARV-VK382, aVARV-VK388, and aVARV-VK470. A region of 125 nucleotides of aVARV-VK470 was identified as being of possible recombinant origin by three of the seven RDP4 algorithms (3Seq, BootScan, and GENECONV), with aVARV-VK281 as the major parent and the minor contribution from an unknown parent. The RDP4 program warned that “the proposed recombination signal may be attributable to a process other than recombination” and of a “possible misidentification of recombinant.” No recombination signals were detected with aVARV-VK281, aVARV-VK382, or aVARV-VK388 as recombinant sequence. No modern sequence was identified as being the possible descendant of a recombination involving any ancient sequence. Details on the algorithms, P values, and the possible breakpoints for aVARV-470 are shown in table S5.
Phylogenetic analyses
1) Maximum likelihood trees: Separate trees were made from datasets 1 and 2 (see supplementary text). The sequences were aligned using MAFFT. Maximum likelihood (ML) trees were generated using RAxML-NG with a K81uf+G+I (selected using bModelTest) substitution model and an ML estimate of tree topology, branch lengths, substitution rates, and nucleotide frequencies, with support estimated using 1000 bootstrap replicates.
2) Lower-coverage samples trees: The EPA-ng algorithm was used to infer the positions of consensus sequences from the nine lower-coverage samples (VK108, VK138, VK168, VK255, VK443, VK515, VK533, FIN1, and KHA1) on the ML tree made from dataset 2, including the sequences from the four higher-coverage samples. For all lower-coverage samples, consensus sequences were obtained by extracting the majority allele from BAM files made from reads mapped against the TATV reference genome (GenBank accession no. NC_008291.1) using the local alignment algorithm of Bowtie2, to reduce the impact of genotyping errors due to postmortem DNA damage (fig. S2). The resulting consensus sequence was aligned to the multisequence alignment of dataset 2 from which the original ML tree was inferred (see Maximum likelihood trees paragraph, above) using MAFFT (version 7.407) with the “–add” and “–keeplength” options. EPA-ng was run using the same model parameters that were used for the original ML tree (K81uf+G+I). Placements of lower-coverage samples were further analyzed using Gappa (54), and the uncertainty in placements for each sample was visualized by plotting the inferred likelihood weight ratios (LWRs) for each branch on the reference tree, using the ggtree package in R (figs. S7 and S8).
To assess the EPA-ng placements of the lower-coverage samples, BAM files with mapped reads from the higher-coverage sample VK388 (aVARV) and the published sample VD21 (mVARV) were randomly subsampled to sets containing 30, 50, 100, 500, 1000, and 2000 reads, with 100 replicates each. Consensus sequences for each sample and replicate were made and placed back on the reference tree, using the same approach as for the lower-coverage samples (described above). Placement uncertainty was visualized using the normalized placement edge mass (78) across all 100 replicates, for each branch on the reference tree (fig. S9).
To further assess the reliability of the placements, two additional metrics were used (79). First, the distributions of LWRs for the first, second, and third most-likely placement for each subsample set were inferred. If many replicates of the placements for a particular subsample set have high confidence, a large fraction of the most-likely placements will show high LWRs (close to the maximum LWR of 1), and if second or third most-likely placements are observed they will have very low LWRs (fig. S11A). Second, the expected distance between placement locations (EDPL) was calculated. EDPL is defined as the average distance on the reference tree between different placements for a query sequence, weighted by their respective LWR. A particular query sequence can have low confidence in the actual placement (i.e., many placements with low LWR) but high confidence for a local neighborhood of the reference tree, for example, in the case of sets of highly similar sequences within a subclade of the reference tree. The EDPL thus gives a measure of the locality of the query placements across the reference tree, with low values indicating high locality (fig. S11B).
We also tested the ability of EPA-ng to correctly place a sequence not currently in the tree. For that, we chose the Abatino orthopoxvirus (GenBank accession no. MH816996.1) sequence, which diverges basal to Ectromelia virus (80). We simulated a total of 1,000,000 sequencing reads from the Abatino reference genome using gargammel. Sequencing reads were simulated from HiSeq 2500 Illumina single-end runs of length 81 base pairs (bp), using an empirical read length distribution derived from sample VK388. After removal of sequencing adapters, read mapping, consensus sequence building, and phylogenetic placement was carried out, as previously described, for random subsamples of 30, 50, 100, 500, 1000, and 2000 reads (fig. S12). As expected, the placement edge masses were concentrated on the branch leading to the Ectromelia virus clade, even for as few as 30 simulated reads.
To complement and confirm the EPA-ng placements, branch-defining single-nucleotide polymorphisms (SNP) were inferred from the reference alignment. A SNP is considered branch-defining if all sequences descending from the branch share the same nucleotide at that location and all other sequences (on the alternate branches) do not. For each lower-coverage sample and each branch of the reference tree, the number and fraction of SNPs where the allele observed in the sample matched the branch-defining allele is shown in fig. S20.
3) Dated coalescent trees using BEAST2: We performed a linear regression of root-to-tip dates against sampling dates. Root-to-tip distances were extracted using TempEst from the ML tree inferred from dataset 1 and the sequences from the four higher-coverage samples in fig. S6 (dataset 1). The regression analysis was performed using SciPy. Dated coalescent trees were inferred using BEAST2. Forty-eight modern VARV sequences were used, as well as two published ancient sequences from the Czech Republic, and VARV-VD21 (dataset 1). Using bModelTest, we selected a TPM1 substitution model with unequal base frequencies, invariant sites, and gamma distributed rate heterogeneity among sites. The clock rate was constrained using a uniform [1 × 10−9 to 1 × 10−3 substitutions per site per year (s/s/y)] prior. Proper priors were used throughout. Trees were inferred using strict and relaxed log-normal molecular clocks, as well as constant and exponential coalescent and Bayesian skyline population priors. The Markov chain Monte Carlo analysis was run for 100M (for the strict clock model and coalescent constant and exponential population priors) or 200M (for the strict clock model and Bayesian skyline population prior, and the log-normal relaxed clock model and coalescent constant, coalescent exponential, and Bayesian skyline population priors) generations, sampling every 2000 generations. Convergence and mixing were assessed using Tracer, and an effective sample size >200 was reached for all relevant parameters. Final tree files were subsampled to contain 10,000 trees, and maximum clade credibility trees were summarized using TreeAnnotator. Path sampling, as implemented in BEAST2, was used to select the best-fitting molecular clock model and population prior. Per path sampling run, 50 steps with a chain length of 1,000,000 generations were run. Likelihood values were compared using a Bayes factor test (81).
Acknowledgments
E.W. thanks St. John’s College, Cambridge, for providing an excellent environment for scientific discussions. We thank S. Rankin and the staff of the University of Cambridge High Performance Computing service, the National High-throughput DNA Sequencing Centre (Copenhagen; seqcenter.ku.dk), L. Drenzel (curator at the Historical Museum in Stockholm, Sweden), and S. Thelin and C. Ödman (archaeologists and curators, Malmö Museum, Sweden). We thank the three reviewers for their suggestions and constructive criticism. Funding: Funding was provided by the Danish National Research Foundation, the Danish National Advanced Technology Foundation (Genome Denmark platform, grant 019-2011-2), the Lundbeck Foundation, the Novo Nordisk Foundation, the Carlsberg Foundation, the Villum Kann Rasmussen Foundation (grant KU2016), and Deutsches Zentrum für Infektionsforschung. L.M.S. is funded by the Norwegian University of Science and Technology, the Royal Norwegian Society of Sciences and Letters, and Sparebank 1 SMN. G.L.S. is a Wellcome Trust Principal Research Fellow. E.W. was supported by the Lundbeck Foundation, the Danish National Research Foundation, the Wellcome Trust (214300/B/18/Z), and the Villum Kann Rasmussen Foundation (grant KU2016). B.M., C.D., R.A.M.F., D.J.S., and T.C.J. were supported by the European Union Horizon 2020 Research and Innovation Programme, COMPARE (grant agreement 643476). T.C.J. was supported by the Wellcome Trust (214300/B/18/Z). Author contributions: B.M.: Initiation of the study, sample identification and data processing, phylogenetic analysis, gene inactivation analysis, visualizations, consensus sequences, manuscript preparation (initial draft, review, and editing). L.V.: Design of capture array, all virus sequence–targeted data production, and contribution to data analysis and writing of the manuscript. A.M., C.d.l.F.C., P.d.B.D., and H.S.: Sampling, sequence data generation, and manuscript review and editing. H.W., A.B., T.P., C.F., V.K., V.M., M.L.S.J., and P.Ø.S.: Osteological analysis, manuscript review and editing, and writing of site descriptions. M.E.A.: Supervision of data production and manuscript review and editing. A.J.H.: Design of capture array. S.H.N.: Initial data collection, processing, and identification of positive samples. L.M.S.: Osteological analysis and writing of site descriptions. J.B.: Selection of material, manuscript review and editing, and writing of site descriptions. Y.M. and I.G.: Contribution to site descriptions. G.S., G.L.S., C.D. R.A.M.F., and D.J.S.: Analysis of results and manuscript review and editing. E.W.: Initiation of the study; leading of efforts related to sampling and data production; and participation in discussion, interpretation of results, and writing of the manuscript. T.C.J.: Initiation of the study, sample identification and data processing, phylogenetic analysis, recombination analysis, read specificity analysis, and manuscript preparation (initial draft, review, and editing). M.S.: Initiation of the study, sample identification and data processing, ancient DNA data authentication, lower-coverage sample phylogenetic analyses, edit distance analysis, coverage, damage figures, and manuscript review and editing. Competing interests: E.W. is involved in scientific collaboration with Illumina, Inc. All other authors declare no competing interests. Data and materials availability: Sequencing data and consensus genomes are available in the European Nucleotide Archive under project number PRJEB38129. Data for sequence alignments, BEAST2 XML files, RAxML-NG output for the ML tree of the higher-coverage sequences, consensus sequences of the lower-coverage samples used in EPA analysis, Python code and data for producing figures and gene loss analysis, and SNP tables for each lower-coverage sample, showing status (derived, ancestral, or missing) at all branch informative SNPs are available via Figshare (83).
Supplementary Material
Summary
Supplementary Text
Figs. S1 to S27
Tables S1 to S17
MDAR Reproducibility Checklist
Resources
File (aaw8977_mdar_reproducibility_cheklist.pdf)
- Download
- 100.16 KB
File (aaw8977_muhlemann_sm.pdf)
- Download
- 43.96 MB
File (aaw8977_tables1.xlsx)
- Download
- 193.96 KB
File (aaw8977_tables14.xlsx)
- Download
- 11.14 KB
File (aaw8977_tables3.xlsx)
- Download
- 81.41 KB
File (aaw8977_tables8.xlsx)
- Download
- 81.45 KB
References and Notes
1
D. R. Hopkins, Princes and Peasants: Smallpox in History (Univ. of Chicago Press, 1983).
2
E. Jenner, On the origin of the vaccine inoculation. Med Phys J. 5, 505–508 (1801).
3
G. L. Smith, G. McFadden, Smallpox: Anything to declare? Nat. Rev. Immunol. 2, 521–527 (2002).
4
C. B. Coyne, Horsepox: Framing a dual use research of concern debate. PLOS Pathog. 14, e1007344 (2018).
5
N. Sklenovská, M. Van Ranst, Emergence of monkeypox as the most important orthopoxvirus infection in humans. Front. Public Health 6, 241 (2018).
6
R. C. Hendrickson, C. Wang, E. L. Hatcher, E. J. Lefkowitz, Orthopoxvirus genome evolution: The role of gene loss. Viruses 2, 1933–1967 (2010).
7
S. Parker, R. Crump, H. Hartzler, R. M. Buller, Evaluation of taterapox virus in small animals. Viruses 9, 203 (2017).
8
F. Fenner, D. A. Henderson, I. Arita, Z. Jezek, I. D. Ladnyi, Smallpox and Its Eradication (World Health Organization, 1988).
9
D. Baxby, Jenner’s Smallpox Vaccine: The Riddle of Vaccinia Virus and Its Origin (Heinemann Educational Publishers, 1981).
10
C. Smithson, J. Imbery, C. Upton, Re-assembly and analysis of an ancient variola virus genome. Viruses 9, 253 (2017).
11
K. A. Bratke, A. McLysaght, S. Rothenburg, A survey of host range genes in poxvirus genomes. Infect. Genet. Evol. 14, 406–425 (2013).
12
G. L. Smith, C. T. O. Benfield, C. Maluquer de Motes, M. Mazzon, S. W. J. Ember, B. J. Ferguson, R. P. Sumner, Vaccinia virus immune evasion: Mechanisms, virulence and immunogenicity. J. Gen. Virol. 94, 2367–2392 (2013).
13
B. Aguado, I. P. Selmes, G. L. Smith, Nucleotide sequence of 21.8 kbp of variola major virus strain Harvey and comparison with vaccinia virus. J. Gen. Virol. 73, 2887–2902 (1992).
14
J. B. Moore, G. L. Smith, Steroid hormone synthesis by a vaccinia enzyme: A new type of virus virulence factor. EMBO J. 11, 3490 (1992).
15
A. Alcamí, G. L. Smith, A mechanism for the inhibition of fever by a virus. Proc. Natl. Acad. Sci. U.S.A. 93, 11029–11034 (1996).
16
A. T. Duggan, M. F. Perdomo, D. Piombino-Mascali, S. Marciniak, D. Poinar, M. V. Emery, J. P. Buchmann, S. Duchêne, R. Jankauskas, M. Humphreys, G. B. Golding, J. Southon, A. Devault, J.-M. Rouillard, J. W. Sahl, O. Dutour, K. Hedman, A. Sajantila, G. L. Smith, E. C. Holmes, H. N. Poinar, 17th century variola virus reveals the recent history of smallpox. Curr. Biol. 26, 3407–3412 (2016).
17
P. Pajer, J. Dresler, H. Kabíckova, L. Písa, P. Aganov, K. Fucik, D. Elleder, T. Hron, V. Kuzelka, P. Velemínsky, J. Klimentova, A. Fucikova, J. Pejchal, R. Hrabakova, V. Benes, T. Rausch, P. Dundr, A. Pilin, R. Cabala, M. Hubalek, J. Stríbrny, M. H. Antwerpen, H. Meyer, Characterization of two historic smallpox specimens from a Czech museum. Viruses 9, 200 (2017).
18
A. F. Porter, A. T. Duggan, H. N. Poinar, E. C. Holmes, Comment: Characterization of Two Historic Smallpox Specimens from a Czech Museum. Viruses 9, 276 (2017).
19
C. W. Dixon, Smallpox (J & A Churchill Ltd., 1962).
20
D. E. Wood, S. L. Salzberg, Kraken: Ultrafast metagenomic sequence classification using exact alignments. Genome Biol. 15, R46 (2014).
21
N. A. O’Leary, M. W. Wright, J. R. Brister, S. Ciufo, D. Haddad, R. McVeigh, B. Rajput, B. Robbertse, B. Smith-White, D. Ako-Adjei, A. Astashyn, A. Badretdin, Y. Bao, O. Blinkova, V. Brover, V. Chetvernin, J. Choi, E. Cox, O. Ermolaeva, C. M. Farrell, T. Goldfarb, T. Gupta, D. Haft, E. Hatcher, W. Hlavina, V. S. Joardar, V. K. Kodali, W. Li, D. Maglott, P. Masterson, K. M. McGarvey, M. R. Murphy, K. O’Neill, S. Pujar, S. H. Rangwala, D. Rausch, L. D. Riddick, C. Schoch, A. Shkeda, S. S. Storz, H. Sun, F. Thibaud-Nissen, I. Tolstoy, R. E. Tully, A. R. Vatsan, C. Wallin, D. Webb, W. Wu, M. J. Landrum, A. Kimchi, T. Tatusova, M. DiCuccio, P. Kitts, T. D. Murphy, K. D. Pruitt, Reference sequence (RefSeq) database at NCBI: Current status, taxonomic expansion, and functional annotation. Nucleic Acids Res. 44, D733–D745 (2016).
22
T. Tatusova, M. DiCuccio, A. Badretdin, V. Chetvernin, E. P. Nawrocki, L. Zaslavsky, A. Lomsadze, K. D. Pruitt, M. Borodovsky, J. Ostell, NCBI prokaryotic genome annotation pipeline. Nucleic Acids Res. 44, 6614–6624 (2016).
23
J. R. Brister, D. Ako-Adjei, Y. Bao, O. Blinkova, NCBI viral genomes resource. Nucleic Acids Res. 43, D571–D577 (2015).
24
B. Mühlemann, T. C. Jones, P. B. Damgaard, M. E. Allentoft, I. Shevnina, A. Logvin, E. Usmanova, I. P. Panyushkina, B. Boldgiv, T. Bazartseren, K. Tashbaeva, V. Merz, N. Lau, V. Smrčka, D. Voyakin, E. Kitov, A. Epimakhov, D. Pokutta, M. Vicze, T. D. Price, V. Moiseyev, A. J. Hansen, L. Orlando, S. Rasmussen, M. Sikora, L. Vinner, A. D. M. E. Osterhaus, D. J. Smith, D. Glebe, R. A. M. Fouchier, C. Drosten, K. G. Sjögren, K. Kristiansen, E. Willerslev, Ancient hepatitis B viruses from the Bronze Age to the Medieval period. Nature 557, 418–423 (2018).
25
E. Cappellini, A. Prohaska, F. Racimo, F. Welker, M. W. Pedersen, M. E. Allentoft, P. de Barros Damgaard, P. Gutenbrunner, J. Dunne, S. Hammann, M. Roffet-Salque, M. Ilardo, J. V. Moreno-Mayar, Y. Wang, M. Sikora, L. Vinner, J. Cox, R. P. Evershed, E. Willerslev, Ancient biomolecules and evolutionary inference. Annu. Rev. Biochem. 87, 1029–1060 (2018).
26
D. P. Martin, B. Murrell, M. Golden, A. Khoosal, B. Muhire, RDP4: Detection and analysis of recombination patterns in virus genomes. Virus Evol. 1, vev003 (2015).
27
P. Barbera, A. M. Kozlov, L. Czech, B. Morel, D. Darriba, T. Flouri, A. Stamatakis, EPA-ng: Massively Parallel Evolutionary Placement of Genetic Sequences. Syst. Biol. 68, 365–369 (2019).
28
A. J. Drummond, O. G. Pybus, A. Rambaut, R. Forsberg, A. G. Rodrigo, Measurably evolving populations. Trends Ecol. Evol. 18, 481–488 (2003).
29
R. Bouckaert, J. Heled, D. Kühnert, T. Vaughan, C.-H. Wu, D. Xie, M. A. Suchard, A. Rambaut, A. J. Drummond, BEAST 2: A software platform for Bayesian evolutionary analysis. PLOS Comput. Biol. 10, e1003537 (2014).
30
A. Alcamí, G. L. Smith, A soluble receptor for interleukin-1β encoded by vaccinia virus: A novel mechanism of virus modulation of the host response to infection. Cell 71, 153–167 (1992).
31
M. Pires de Miranda, P. C. Reading, D. C. Tscharke, B. J. Murphy, G. L. Smith, The vaccinia virus kelch-like protein C2L affects calcium-independent adhesion to the extracellular matrix and inflammation in a murine intradermal model. J. Gen. Virol. 84, 2459–2471 (2003).
32
P. M. Beard, G. C. Froggatt, G. L. Smith, Vaccinia virus kelch protein A55 is a 64 kDa intracellular factor that affects virus-induced cytopathic effect and the outcome of infection in a murine intradermal model. J. Gen. Virol. 87, 1521–1529 (2006).
33
G. C. Froggatt, G. L. Smith, P. M. Beard, Vaccinia virus gene F3L encodes an intracellular protein that affects the innate immune response. J. Gen. Virol. 88, 1917–1921 (2007).
34
G. Kochneva, I. Kolosova, T. Maksyutova, E. Ryabchikova, S. Shchelkunov, Effects of deletions of kelch-like genes on cowpox virus biological properties. Arch. Virol. 150, 1857–1870 (2005).
35
A. Alejo, M. B. Ruiz-Argüello, Y. Ho, V. P. Smith, M. Saraiva, A. Alcami, A chemokine-binding domain in the tumor necrosis factor receptor from variola (smallpox) virus. Proc. Natl. Acad. Sci. U.S.A. 103, 5995–6000 (2006).
36
J. Liu, S. Wennier, L. Zhang, G. McFadden, M062 is a host range factor essential for myxoma virus pathogenesis and functions as an antagonist of host SAMD9 in human cells. J. Virol. 85, 3270–3282 (2011).
37
G. Sivan, P. Ormanoglu, E. C. Buehler, S. E. Martin, B. Moss, Identification of restriction factors by human genome-wide RNA interference screening of viral host range mutants exemplified by discovery of SAMD9 and WDR6 as inhibitors of the vaccinia virus K1L−C7L− mutant. mBio 6, e01122-15 (2015).
38
S. Gillard, D. Spehner, R. Drillien, A. Kirn, Localization and sequence of a vaccinia virus gene required for multiplication in human cells. Proc. Natl. Acad. Sci. U.S.A. 83, 5573–5577 (1986).
39
M. E. Perkus, S. J. Goebel, S. W. Davis, G. P. Johnson, K. Limbach, E. K. Norton, E. Paoletti, Vaccinia virus host range genes. Virology 179, 276–286 (1990).
40
J. O. Langland, B. L. Jacobs, The role of the PKR-inhibitory genes, E3L and K3L, in determining vaccinia virus host range. Virology 299, 133–141 (2002).
41
B. Liu, D. Panda, J. D. Mendez-Rios, S. Ganesan, L. S. Wyatt, B. Moss, Identification of poxvirus genome uncoating and DNA replication factors with mutually redundant roles. J. Virol. 92, e02152-17 (2018).
42
C. McEvedy, R. Jones, Atlas of World Population History (Puffin, 1978).
43
C. Creighton, A History of Epidemics in Britain: From A.D. 664 to the Extinction of Plague (Cambridge Univ. Press, 1891).
44
R. Willan, Miscellaneous Works of the Late Robert Willan, M.D. (T. Cadell, 1821).
45
J. Moore, The History of the Small Pox (Longman, Hurst, Rees, Orme, and Brown, 1815).
46
E. Paschen, in G. Jochmann’s Lehrbuch der Infektionskrankheiten (Springer Verlag, ed. 2, 1924).
47
S. A. Ishikawa, A. Zhukova, W. Iwasaki, O. Gascuel, A fast likelihood method to reconstruct and visualize ancestral scenarios. Mol. Biol. Evol. 36, 2069–2085 (2019).
48
C. Camacho, G. Coulouris, V. Avagyan, N. Ma, J. Papadopoulos, K. Bealer, T. L. Madden, BLAST+: Architecture and applications. BMC Bioinformatics 10, 421 (2009).
49
R. R. Bouckaert, A. J. Drummond, bModelTest: Bayesian phylogenetic site model averaging and model comparison. BMC Evol. Biol. 17, 42 (2017).
50
B. Langmead, S. L. Salzberg, Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012).
51
Broad Institute, Picard Toolkit; http://broadinstitute.github.io/picard/.
52
H. Li, R. Durbin, Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 26, 589–595 (2010).
53
B. Buchfink, C. Xie, D. H. Huson, Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60 (2015).
54
L. Czech, P. Barbera, A. Stamatakis, Genesis and Gappa: Processing, analyzing and visualizing phylogenetic (placement) data. Bioinformatics 36, 3263–3265 (2020).
55
G. Renaud, K. Hanghøj, E. Willerslev, L. Orlando, Gargammel: A sequence simulator for ancient DNA. Bioinformatics 33, 577–579 (2017).
56
M. Kearse, R. Moir, A. Wilson, S. Stones-Havas, M. Cheung, S. Sturrock, S. Buxton, A. Cooper, S. Markowitz, C. Duran, T. Thierer, B. Ashton, P. Meintjes, A. Drummond, Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649 (2012).
57
G. Yu, D. K. Smith, H. Zhu, Y. Guan, T. T.-Y. Lam, ggtree: An R package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol. Evol. 8, 28–36 (2016).
58
K. Katoh, D. M. Standley, MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).
59
H. Jónsson, A. Ginolhac, M. Schubert, P. L. F. Johnson, L. Orlando, mapDamage2.0: Fast approximate Bayesian estimates of ancient DNA damage parameters. Bioinformatics 29, 1682–1684 (2013).
60
B. S. Pedersen, A. R. Quinlan, Mosdepth: Quick coverage calculation for genomes and exomes. Bioinformatics 34, 867–868 (2018).
61
R Core Team, R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2008); www.R-project.org.
62
A. M. Kozlov, D. Darriba, T. Flouri, B. Morel, A. Stamatakis, RAxML-NG: A fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics 35, 4453–4455 (2019).
63
H. Li, B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin; 1000 Genome Project Data Processing Subgroup, The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).
64
P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt; SciPy 1.0 Contributors, SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272 (2020).
65
A. Rambaut, T. T. Lam, L. Max Carvalho, O. G. Pybus, Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). Virus Evol. 2, vew007 (2016).
66
A. Rambaut, M. Suchard, A. Drummond, Tracer; http://tree.bio.ed.ac.uk/software/tracer/.
67
E. Willerslev, A. Cooper, Ancient DNA. Proc. R. Soc. London Ser. B 272, 3–16 (2005).
68
C. Warinner, A. Herbig, A. Mann, J. A. Fellows Yates, C. L. Weiß, H. A. Burbano, L. Orlando, J. Krause, A robust framework for microbial archaeology. Annu. Rev. Genomics Hum. Genet. 18, 321–356 (2017).
69
T. Bigot, S. Temmam, P. Pérot, M. Eloit, RVDB-prot, a reference viral protein database and its HMM profiles. F1000 Res. 8, 530 (2019).
70
C. Upton, S. Slack, A. L. Hunter, A. Ehlers, R. L. Roper, Poxvirus orthologous clusters: Toward defining the minimum essential poxvirus genome. J. Virol. 77, 7590–7600 (2003).
71
M. F. Boni, D. Posada, M. W. Feldman, An exact nonparametric method for inferring mosaic structure in sequence triplets. Genetics 176, 1035–1047 (2007).
72
D. P. Martin, D. Posada, K. A. Crandall, C. Williamson, A modified bootscan algorithm for automated identification of recombinant sequences and recombination breakpoints. AIDS Res. Hum. Retroviruses 21, 98–102 (2005).
73
D. Posada, K. A. Crandall, Evaluation of methods for detecting recombination from DNA sequences: Computer simulations. Proc. Natl. Acad. Sci. U.S.A. 98, 13757–13762 (2001).
74
M. Padidam, S. Sawyer, C. M. Fauquet, Possible emergence of new geminiviruses by frequent recombination. Virology 265, 218–225 (1999).
75
J. M. Smith, Analyzing the mosaic structure of genes. J. Mol. Evol. 34, 126–129 (1992).
76
D. Martin, E. Rybicki, RDP: Detection of recombination amongst aligned sequences. Bioinformatics 16, 562–563 (2000).
77
M. J. Gibbs, J. S. Armstrong, A. J. Gibbs, Sister-scanning: A Monte Carlo procedure for assessing signals in recombinant sequences. Bioinformatics 16, 573–582 (2000).
78
L. Czech, A. Stamatakis, Scalable methods for analyzing and visualizing phylogenetic placement of metagenomic samples. PLOS ONE 14, e0217050 (2019).
79
F. Mahé, C. de Vargas, D. Bass, L. Czech, A. Stamatakis, E. Lara, D. Singer, J. Mayor, J. Bunge, S. Sernaker, T. Siemensmeyer, I. Trautmann, S. Romac, C. Berney, A. Kozlov, E. A. D. Mitchell, C. V. W. Seppey, E. Egge, G. Lentendu, R. Wirth, G. Trueba, M. Dunthorn, Parasites dominate hyperdiverse soil protist communities in Neotropical rainforests. Nat. Ecol. Evol. 1, 0091 (2017).
80
C. E. M. Gruber, E. Giombini, M. Selleri, S. H. Tausch, A. Andrusch, A. Tyshaieva, G. Cardeti, R. Lorenzetti, L. De Marco, F. Carletti, A. Nitsche, M. R. Capobianchi, G. Ippolito, G. L. Autorino, C. Castilletti, Whole genome characterization of Orthopoxvirus (OPV) Abatino, a zoonotic virus representing a putative novel clade of Old World orthopoxviruses. Viruses 10, 546 (2018).
81
R. E. Kass, A. E. Raftery, Bayes Factors. J. Am. Stat. Assoc. 90, 773–795 (1995).
82
P. J. Reimer, E. Bard, A. Bayliss, J. W. Beck, P. G. Blackwell, C. B. Ramsey, C. E. Buck, H. Cheng, R. L. Edwards, M. Friedrich, P. M. Grootes, T. P. Guilderson, H. Haflidason, I. Hajdas, C. Hatté, T. J. Heaton, D. L. Hoffmann, A. G. Hogg, K. A. Hughen, K. F. Kaiser, B. Kromer, S. W. Manning, M. Niu, R. W. Reimer, D. A. Richards, E. M. Scott, J. R. Southon, R. A. Staff, C. S. M. Turney, J. van der Plicht, IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 Years cal BP. Radiocarbon 55, 1869–1887 (2013).
83
B. Mühlemann, T. C. Jones, M. Sikora, Processed data for phylogenetic analysis of ancient variola virus samples, Figshare (2020); https://doi.org/10.6084/m9.figshare.12185466.
84
W. M. Bass, Human Osteology: A Laboratory and Field Manual (Missouri Archaeological Society, 1987).
85
M. B. Brickley, Cribra orbitalia and porotic hyperostosis: A biological approach to diagnosis. Am. J. Phys. Anthropol. 167, 896–902 (2018).
86
P. L. Walker, R. R. Bathurst, R. Richman, T. Gjerdrum, V. A. Andrushko, The causes of porotic hyperostosis and cribra orbitalia: A reappraisal of the iron-deficiency-anemia hypothesis. Am. J. Phys. Anthropol. 139, 109–125 (2009).
87
E. Naumann, T. D. Price, M. P. Richards, Changes in dietary practices and social organization during the pivotal late iron age period in Norway (AD 550-1030): Isotope analyses of Merovingian and Viking Age human remains. Am. J. Phys. Anthropol. 155, 322–331 (2014).
88
D. Avdussin, Gnezdovo–der Nachbar von Smolensk. Zeitschrift für Archdologie, 263–290 (1977).
89
T. Puškina, V. Muraševa, N. Eniosova, in Die Rus’ im 9. – 10. Jahrhundert. Ein archaologisches Panorama, N. A. Makarov, Ed. (Wachholtz Murmann, 2017) pp. 250–281.
90
M. Müller-Wille, Boat-graves in northern Europe. Int. J. Naut. Archaeol. 3, 187–204 (1974).
91
M. H. Beskow-Sjöberg, K.-H. Arnell, Ölands järnåldersgravfält. Vol. I, Alböke, Köpings, Räpplinge, Löts, Egby, Bredsätra och Gärdslösa socknar (Riksantikvarieämbetet och Statens historiska museer, 1987).
92
M. Beskow-Sjöberg, U. E. Hagberg, Ölands järnåldersgravfält, Vol. II (Högskolan i Kalmar, 1991).
93
U. E. Hagberg, M. Beskow-Sjöberg, Ölands järnåldersgravfält. Vol. III (Riksantikvarieämbetet och Statens historiska museer, 1996).
94
J.-H. Fallgren, M. Rasch, Ölands järnåldersgravfält. Vol. IV (Riksantikvarieämbetet och Statens historiska museer, 2001).
95
J. E. Anderbjörk, Brev daterat 17/2 1939. Till Riksantikvarieämbetet. Dnr 0975 (1939).
96
T. D. Price, J. Peets, R. Allmäe, L. Maldre, E. Oras, Isotopic provenancing of the Salme ship burials in Pre-Viking Age Estonia. Antiquity 90, 1022–1037 (2016).
97
J. E. Anderbjörk, “Nabberör i Böda-undersökningen av en märklig öländsk båtgrav” in Meddelanden XXVII (Kalmar Läns Fornminnesförening, 1939) pp. 63–72.
98
H. Wilhelmson, Perspectives from a Human-Centred Archaeology: Iron Age People and Society on Öland, vol. 2 of Studies in Osteology (Lund University, 2017).
99
H. Wilhelmson, T. D. Price, Migration and integration on the Baltic island of Öland in the Iron Age. J. Archaeol. Sci. Rep. 12, 183–196 (2017).
100
H. Wilhelmson, Shifting diet, shifting culture? A bioarchaeological approach to island dietary development on Iron-Age Öland, Baltic Sea. Am. J. Phys. Anthropol. 163, 264–284 (2017).
101
A. Margaryan, thesis, University of Copenhagen (2017).
102
P. Ø. Sørensen, Markildegård – En tidligneolitisk samlingsplads (Kulturhistoriske Studier I, Sydsjællands Museum, 1995), pp. 13–45.
103
P. Ø. Sørensen, “Markildegård – mødested gennem årtusinder” in 5000 år under motorvejen J. Hertz, Ed.(Vejdirektoratet og Rigsantikvarens Arkæologiske Sekretariat, 1995), pp. 32–33.
104
S. Wallis, The Oxford Henge and Late Saxon Massacre: With Medieval and Later Occupation at St John’s College, Oxford (Thames Valley Archaeological Services Ltd, 2014).
105
T. D. Price, M. Naum, P. Bennike, N. Lynnerup, K. M. Frei, H. Wagnkilde, T. Pind, F. O. Nielsen, Isotopic investigation of human provenience at the eleventh century cemetery of Ndr. Grødbygård, Bornholm, Denmark. Dan. J. Archaeol. 1, 93–112 (2012).
106
M. Naum, Homelands Lost and Gained: Slavic Migration and Settlement on Bornholm in the Early Middle Ages, vol. 9 of Lund Studies in Historical Archaeology (Lund University, 2008).
107
Y. Bäckström, T. D. Price, Social identity and mobility at a pre-industrial mining complex, Sweden. J. Archaeol. Sci. 66, 154–168 (2016).
108
I. Gustin, “Contacts, identity, and hybridity:Objects from South-western Finland in the Birka Graves” in Identity Formationand Diversity in the Early Medieval Baltic and Beyond, vol. 75 of TheNorthern World, J. Callmer, I. Gustin, M. Roslund, Eds. (Brill, 2017), pp.205–258.
109
S. Brooks, J. M. Suchey, Skeletal age determination based on the os pubis: A comparison of the Acsádi-Nemeskéri and Suchey-Brooks methods. Hum. Evol. 5, 227–238 (1990).
110
C. Cunningham, L. Scheuer, S. Black, Developmental Juvenile Osteology (Academic Press, ed. 2, 2017).
111
U. E. Hagberg, Archive report of excavations during 1965 in Sörby-Störlinge gravefield, fornl 74,Gärdslösa sn, Öland (1966).
112
U. E. H. Näsman, “The Setting” in Eketorp. Fortification and Settlement on Öland/Sweden, U. E. H. Näsman, E. Wegraeus, Eds. (Royal Academy of Letters, History and Antiquities, 1979), pp. 120–121.
113
B. Arrhenius, Knivar från Helgö och Birka. Fornvannen 1970, 40–51 (1970).
114
T. Waldron, Paleopathology (Cambridge Univ. Press, 2008).
115
N. N. Gutsol, S. N. Vinogradova, A. G. Samorukova, Kola Saami Relocated Groups (Apatity: Kola Science Center RAS, 2007).
116
V. I. Khartanovich, in Sbornik Muzeya antropologii i ehtnografii, I. Gokhman, Ed. (Leningrad, 1980), vol. 36, pp. 35–47.
117
M. Krzewińska, G. Bjørnstad, P. Skoglund, P. I. Olason, J. Bill, A. Götherström, E. Hagelberg, Mitochondrial DNA variation in the Viking age population of Norway. Philos. Trans. R. Soc. London Ser. B 370, 20130384 (2015).
118
S. Niinimäki, What do muscle marker ruggedness scores actually tell us? Int. J. Osteoarchaeol. 21, 292–299 (2009).
119
W. K. Seow, Developmental defects of enamel and dentine: Challenges for basic science research and clinical management. Aust. Dent. J. 59 (suppl. 1), 143–154 (2014).
120
O. Murashko, N. Krenke, Aboriginal Culture of the Obdorsk North in the XIX Century (According to the Archaeological and Ethnographic Collections of the Museum of Anthropology of Moscow State University) (Nauka, 2001).
121
D. J. Ortner, Identification of Pathological Conditions in Human Skeletal Remains (Academic Press, 2003).
122
B. Å. Samuelsson, “Ljungbacka, Lockarp 8 Archive Report” (Malmö Museum, 1976).
123
B. Å. Samuelsson, Ljungbacka, a Late Iron Age Cemetery in Southwest Scania. Lund Archaeological Review, 89–108 (2001).
124
J. E. Buikstra, D. H. Ubelaker, Eds., Standards for Data Collection from Human Skeletal Remains, vol. 44 of Arkansas Archaeological Survey Research Series (Arkansas Archaeological Survey, 1994).
125
M. K. Moore, “Sex estimation and assessment” in Research Methods in Human Skeletal Biology, E. DiGangi, M. Moore, Eds. (Academic Press, 2013), pp. 91–116.
126
M. A. Kelley, Phenice’s visual sexing technique for the os pubis: A critique. Am. J. Phys. Anthropol. 48, 121–122 (1978).
127
D. H. Ubelaker, C. G. Volk, A test of the phenice method for the estimation of sex. J. Forensic Sci. 47, 19–24 (2002).
128
N. C. Lovell, Test of Phenice’s technique for determining sex from the os pubis. Am. J. Phys. Anthropol. 79, 117–120 (1989).
129
J. Bruzek, A method for visual determination of sex, using the human hip bone. Am. J. Phys. Anthropol. 117, 157–168 (2002).
130
P. L. Walker, Greater sciatic notch morphology: Sex, age, and population differences. Am. J. Phys. Anthropol. 127, 385–391 (2005).
131
J. L. Buckberry, A. T. Chamberlain, Age estimation from the auricular surface of the ilium: A revised method. Am. J. Phys. Anthropol. 119, 231–239 (2002).
132
S. M. Hens, M. G. Belcastro, Auricular surface aging: A blind test of the revised method on historic Italians from Sardinia. Forensic Sci. Int. 214, 209.e1–209.e5 (2012).
133
C. E. Merritt, The influence of body size on adult skeletal age estimation methods. Am. J. Phys. Anthropol. 156, 35–57 (2015).
134
K. Moraitis, E. Zorba, C. Eliopoulos, S. C. Fox, A test of the revised auricular surface aging method on a modern European population. J. Forensic Sci. 59, 188–194 (2014).
135
C. Rissech, J. Wilson, A. P. Winburn, D. Turbón, D. Steadman, A comparison of three established age estimation methods on an adult Spanish sample. Int. J. Legal Med. 126, 145–155 (2012).
136
M. San Millán, C. Rissech, D. Turbón, A test of Suchey–Brooks (pubic symphysis) and Buckberry–Chamberlain (auricular surface) methods on an identified Spanish sample: Paleodemographic implications. J. Archaeol. Sci. 40, 1743–1751 (2013).
137
C. A. Kunos, S. W. Simpson, K. F. Russell, I. Hershkovitz, First rib metamorphosis: Its possible utility for human age-at-death estimation. Am. J. Phys. Anthropol. 110, 303–323 (1999).
138
L. Nilsson, “Osteologisk Analys, Ljungbacka, Lockarp 8” (Malmö Museum, 1976).
139
J. M. Suchey, Instructions for use of the Suchey-Brooks system for age determination of the female os pubis. Instructional materials accompanying female pubic symphysial models of the Suchey-Brooks system (France Casting, 1988).
140
C. Arcini, Prone burials. Current Archaeology 2009, 30–35 (2009).
141
M. S. Toplak, Prone burials and modified teeth at the Viking Age cemetery of Kopparsvik: The changing of social identities at the threshold of the Christian Middle Ages. Analecta Archaeologica Ressoviensia 10, 77–97 (2015).
142
E. Arwill-Nordbladh, Nyckelsymbolik i järnålderns kvinnogravar. Fornvannen 1990, 85 (1990).
143
T. Zachrisson, in I Att befolka det förflutna. Fem artiklar om hur vi kan synliggöra människan och hennes handlingar i arkeologiskt material, A. Carlie, Ed. (Riksantikvarieämbetet, 2014), pp. 73–91.
144
T. D. Price, K. Prangsgaard, M. Kanstrup, P. Bennike, K. M. Frei, Galgedil: Isotopic studies of a Viking cemetery on the Danish island of Funen, AD 800–1050. Dan. J. Archaeol. 3, 129–144 (2014).
145
M. Kanstrup, in Forhistorisk Arkeologi (Aarhus Universitet, 2006).
146
L. Melchior, T. Kivisild, N. Lynnerup, J. Dissing, Evidence of authentic DNA from Danish Viking Age skeletons untouched by humans for 1,000 years. PLOS ONE 3, e2214 (2008).
147
J. D. Gardner, D. C. Tscharke, P. C. Reading, G. L. Smith, Vaccinia virus semaphorin A39R is a 50-55 kDa secreted glycoprotein that affects the outcome of infection in a murine intradermal model. J. Gen. Virol. 82, 2083–2093 (2001).
148
M. P. Holgado, J. Falivene, C. Maeto, M. Amigo, M. F. Pascutti, M. B. Vecchione, A. Bruttomesso, G. Calamante, M. P. Del Médico-Zajac, M. M. Gherardi, Deletion of A44L, A46R and C12L vaccinia virus genes from the MVA genome improved the vector immunogenicity by modifying the innate immune response generating enhanced and optimized specific T-cell responses. Viruses 8, 139 (2016).
149
P. C. Reading, J. B. Moore, G. L. Smith, Steroid hormone synthesis by vaccinia virus suppresses the inflammatory response to infection. J. Exp. Med. 197, 1269–1278 (2003).
150
P. C. Reading, A. Khanna, G. L. Smith, Vaccinia virus CrmE encodes a soluble and cell surface tumor necrosis factor receptor that contributes to virus virulence. Virology 292, 285–298 (2002).
151
K. Dai, Y. Liu, M. Liu, J. Xu, W. Huang, X. Huang, L. Liu, Y. Wan, Y. Hao, Y. Shao, Pathogenicity and immunogenicity of recombinant Tiantan Vaccinia Virus with deleted C12L and A53R genes. Vaccine 26, 5062–5071 (2008).
152
J. B. Eaglesham, Y. Pan, T. S. Kupper, P. J. Kranzusch, Viral and metazoan poxins are cGAMP-specific nucleases that restrict cGAS-STING signalling. Nature 566, 259–263 (2019).
153
A. H. Banham, G. L. Smith, Characterization of vaccinia virus gene B12R. J. Gen. Virol. 74, 2807–2812 (1993).
154
R. Liu, B. Moss, Vaccinia virus C9 ankyrin repeat/F-box protein is a newly identified antagonist of the type I interferon-induced antiviral state. J. Virol. 92, e00053-18 (2018).
155
M. T. Harte, I. R. Haga, G. Maloney, P. Gray, P. C. Reading, N. W. Bartlett, G. L. Smith, A. Bowie, L. A. J. O’Neill, The poxvirus protein A52R targets Toll-like receptor signaling complexes to suppress host defense. J. Exp. Med. 197, 343–351 (2003).
156
A. Bowie, E. Kiss-Toth, J. A. Symons, G. L. Smith, S. K. Dower, L. A. O’Neill, A46R and A52R from vaccinia virus are antagonists of host IL-1 and toll-like receptor signaling. Proc. Natl. Acad. Sci. U.S.A. 97, 10162–10167 (2000).
157
S. C. Graham, M. W. Bahar, S. Cooray, R. A.-J. Chen, D. M. Whalen, N. G. A. Abrescia, D. Alderton, R. J. Owens, D. I. Stuart, G. L. Smith, J. M. Grimes, Vaccinia virus proteins A52 and B14 Share a Bcl-2-like fold but have evolved to inhibit NF-kappaB rather than apoptosis. PLOS Pathog. 4, e1000128 (2008).
158
K. H. Martin, D. W. Grosenbach, C. A. Franke, D. E. Hruby, Identification and analysis of three myristylated vaccinia virus late proteins. J. Virol. 71, 5218–5226 (1997).
159
V. I. Chernos, T. S. Vovk, O. N. Ivanova, T. P. Antonova, V. N. Loparev, [Insertion mutants of the vaccinia virus. The effect of inactivating E7R and D8L genes on the biological properties of the virus]. Mol. Gen. Mikrobiol. Virusol. 1993, 30–34 (1993).
160
N. Price, D. C. Tscharke, M. Hollinshead, G. L. Smith, Vaccinia virus gene B7R encodes an 18-kDa protein that is resident in the endoplasmic reticulum and affects virus virulence. Virology 267, 65–79 (2000).
161
S. N. Shchelkunov, Orthopoxvirus genes that mediate disease virulence and host tropism. Adv. Virol. 2012, 524743 (2012).
162
J. L. Shisler, X.-L. Jin, The vaccinia virus K1L gene product inhibits host NF-κB activation by preventing IκBα degradation. J. Virol. 78, 3553–3560 (2004).
163
K. E. Rehm, R. F. Connor, G. J. B. Jones, K. Yimbu, R. L. Roper, Vaccinia virus A35R inhibits MHC class II antigen presentation. Virology 397, 176–186 (2010).
164
D. Wilcock, S. A. Duncan, P. Traktman, W.-H. Zhang, G. L. Smith, The vaccinia virus A40R gene product is a nonstructural, type II membrane glycoprotein that is expressed at the cell surface. J. Gen. Virol. 80, 2137–2148 (1999).
165
D. C. Tscharke, P. C. Reading, G. L. Smith, Dermal infection with vaccinia virus reveals roles for virus proteins not seen using other inoculation routes. J. Gen. Virol. 83, 1977–1986 (2002).
166
M. A. Pallett, H. Ren, R.-Y. Zhang, S. R. Scutts, L. Gonzalez, Z. Zhu, C. Maluquer de Motes, G. L. Smith, Vaccinia virus BBK E3 ligase adaptor A55 targets importin-dependent NF-κB activation and inhibits CD8+ T-cell memory. J. Virol. 93, e00051-19 (2019).
167
G. L. Smith, Y. S. Chan, S. T. Howard, Nucleotide sequence of 42 kbp of vaccinia virus strain WR from near the right inverted terminal repeat. J. Gen. Virol. 72, 1349–1376 (1991).
168
A. Murcia-Nicolas, G. Bolbach, J. C. Blais, G. Beaud, Identification by mass spectroscopy of three major early proteins associated with virosomes in vaccinia virus-infected cells. Virus Res. 59, 1–12 (1999).
169
World Health Organization Smallpox Fact Sheet (2012); https://web.archive.org/web/20120901000000*/http://www.who.int/mediacentre/factsheets/smallpox/en/index.html.
170
World Health Organization, WHO Hepatitis B Fact Sheet (2019); www.who.int/en/news-room/fact-sheets/detail/hepatitis-b.
171
A. G. Carmichael, A. M. Silverstein, Smallpox in Europe before the seventeenth century: Virulent killer or benign disease? J. Hist. Med. Allied Sci. 42, 147–168 (1987).
172
Y. Li, D. S. Carroll, S. N. Gardner, M. C. Walsh, E. A. Vitalis, I. K. Damon, On the origin of smallpox: Correlating variola phylogenics with historical smallpox records. Proc. Natl. Acad. Sci. U.S.A. 104, 15787–15792 (2007).
(0)eLetters
eLetters is a forum for ongoing peer review. eLetters are not edited, proofread, or indexed, but they are screened. eLetters should provide substantive and scholarly commentary on the article. Embedded figures cannot be submitted, and we discourage the use of figures within eLetters in general. If a figure is essential, please include a link to the figure within the text of the eLetter. Please read our Terms of Service before submitting an eLetter.
Log In to Submit a ResponseNo eLetters have been published for this article yet.
Information & Authors
Information
Published In
Science
Volume 369 | Issue 6502
24 July 2020
24 July 2020
Copyright
Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works.
This is an article distributed under the terms of the Science Journals Default License.
Submission history
Received: 13 February 2019
Accepted: 29 May 2020
Published in print: 24 July 2020
Acknowledgments
E.W. thanks St. John’s College, Cambridge, for providing an excellent environment for scientific discussions. We thank S. Rankin and the staff of the University of Cambridge High Performance Computing service, the National High-throughput DNA Sequencing Centre (Copenhagen; seqcenter.ku.dk), L. Drenzel (curator at the Historical Museum in Stockholm, Sweden), and S. Thelin and C. Ödman (archaeologists and curators, Malmö Museum, Sweden). We thank the three reviewers for their suggestions and constructive criticism. Funding: Funding was provided by the Danish National Research Foundation, the Danish National Advanced Technology Foundation (Genome Denmark platform, grant 019-2011-2), the Lundbeck Foundation, the Novo Nordisk Foundation, the Carlsberg Foundation, the Villum Kann Rasmussen Foundation (grant KU2016), and Deutsches Zentrum für Infektionsforschung. L.M.S. is funded by the Norwegian University of Science and Technology, the Royal Norwegian Society of Sciences and Letters, and Sparebank 1 SMN. G.L.S. is a Wellcome Trust Principal Research Fellow. E.W. was supported by the Lundbeck Foundation, the Danish National Research Foundation, the Wellcome Trust (214300/B/18/Z), and the Villum Kann Rasmussen Foundation (grant KU2016). B.M., C.D., R.A.M.F., D.J.S., and T.C.J. were supported by the European Union Horizon 2020 Research and Innovation Programme, COMPARE (grant agreement 643476). T.C.J. was supported by the Wellcome Trust (214300/B/18/Z). Author contributions: B.M.: Initiation of the study, sample identification and data processing, phylogenetic analysis, gene inactivation analysis, visualizations, consensus sequences, manuscript preparation (initial draft, review, and editing). L.V.: Design of capture array, all virus sequence–targeted data production, and contribution to data analysis and writing of the manuscript. A.M., C.d.l.F.C., P.d.B.D., and H.S.: Sampling, sequence data generation, and manuscript review and editing. H.W., A.B., T.P., C.F., V.K., V.M., M.L.S.J., and P.Ø.S.: Osteological analysis, manuscript review and editing, and writing of site descriptions. M.E.A.: Supervision of data production and manuscript review and editing. A.J.H.: Design of capture array. S.H.N.: Initial data collection, processing, and identification of positive samples. L.M.S.: Osteological analysis and writing of site descriptions. J.B.: Selection of material, manuscript review and editing, and writing of site descriptions. Y.M. and I.G.: Contribution to site descriptions. G.S., G.L.S., C.D. R.A.M.F., and D.J.S.: Analysis of results and manuscript review and editing. E.W.: Initiation of the study; leading of efforts related to sampling and data production; and participation in discussion, interpretation of results, and writing of the manuscript. T.C.J.: Initiation of the study, sample identification and data processing, phylogenetic analysis, recombination analysis, read specificity analysis, and manuscript preparation (initial draft, review, and editing). M.S.: Initiation of the study, sample identification and data processing, ancient DNA data authentication, lower-coverage sample phylogenetic analyses, edit distance analysis, coverage, damage figures, and manuscript review and editing. Competing interests: E.W. is involved in scientific collaboration with Illumina, Inc. All other authors declare no competing interests. Data and materials availability: Sequencing data and consensus genomes are available in the European Nucleotide Archive under project number PRJEB38129. Data for sequence alignments, BEAST2 XML files, RAxML-NG output for the ML tree of the higher-coverage sequences, consensus sequences of the lower-coverage samples used in EPA analysis, Python code and data for producing figures and gene loss analysis, and SNP tables for each lower-coverage sample, showing status (derived, ancestral, or missing) at all branch informative SNPs are available via Figshare (83).
Authors
Funding Information
Villum Fonden: KU2016
Wellcome: 214300/Z/18/Z
Wellcome: 214300/Z/18/Z
Horizon 2020 Framework Programme: 643476
Horizon 2020 Framework Programme: 643476
Horizon 2020 Framework Programme: 643476
Horizon 2020 Framework Programme: 643476
Danish National Advanced Technology Foundation: 019-2011-2
Horizon 2020 Framework Programme: 643476
Royal Norwegian Society of Sciences and Letters and Sparebank
Metrics & Citations
Metrics
Article Usage
Altmetrics
Citations
Cite as
- Barbara Mühlemann et al.
Export citation
Select the format you want to export the citation of this publication.
Cited by
- Molecular Mechanisms of Poxvirus Evolution, mBio, 14, 1, (2023).https://doi.org/10.1128/mbio.01526-22
- APOBEC3 deaminase editing in mpox virus as evidence for sustained human transmission since at least 2016, Science, 382, 6670, (595-600), (2023)./doi/10.1126/science.adg8116
- Ancient chicken remains reveal the origins of virulence in Marek’s disease virus, Science, 382, 6676, (1276-1281), (2023)./doi/10.1126/science.adg2238
- Analysis of variola virus molecular evolution suggests an old origin of the virus consistent with historical records, Microbial Genomics, 9, 1, (2023).https://doi.org/10.1099/mgen.0.000932
- Benchmarking metagenomics classifiers on ancient viral DNA: a simulation study, PeerJ, 10, (e12784), (2022).https://doi.org/10.7717/peerj.12784
- Detection of Ancient Viruses and Long-Term Viral Evolution, Viruses, 14, 6, (1336), (2022).https://doi.org/10.3390/v14061336
- An Update of Orthopoxvirus Molecular Evolution, Viruses, 14, 2, (388), (2022).https://doi.org/10.3390/v14020388
- In the Search of Marine Pestiviruses: First Case of Phocoena Pestivirus in a Belt Sea Harbour Porpoise, Viruses, 14, 1, (161), (2022).https://doi.org/10.3390/v14010161
- Preventive Measures against Pandemics from the Beginning of Civilization to Nowadays—How Everything Has Remained the Same over the Millennia, Journal of Clinical Medicine, 11, 7, (1960), (2022).https://doi.org/10.3390/jcm11071960
- Historic and Prehistoric Epidemics: An Overview of Sources Available for the Study of Ancient Pathogens, Epidemiologia, 3, 4, (443-464), (2022).https://doi.org/10.3390/epidemiologia3040034
- See more
Loading...
View Options
View options
PDF format
Download this article as a PDF file
Download PDFCheck Access
Log in to view the full text
AAAS login provides access to Science for AAAS Members, and access to other journals in the Science family to users who have purchased individual subscriptions.
- Become a AAAS Member
- Activate your AAAS ID
- Purchase Access to Other Journals in the Science Family
- Account Help
Log in via OpenAthens.
Log in via Shibboleth.
More options
Register for free to read this article
As a service to the community, this article is available for free. Login or register for free to read this article.
Buy a single issue of Science for just $15 USD.
RE: Diverse variola virus (smallpox) strains were widespread in northern Europe in the Viking Age
Mühlemann et al. ("Diverse variola virus (smallpox) strains were widespread in northern Europe in the Viking Age," Science 369 (2020), eaaw8977) report the capture, sequence and reconstruction of 4 Variola virus (VARV) genomes, the oldest yet drafted (600-1050 CE). While advancing markedly our knowledge of the evolutionary history and genomic content of VARV, the explanation and contextualization offered for their findings requires further clarification. The origins of smallpox have been debated for centuries. Before the disease was understood to be pathogenic, physicians and polymaths employed vague descriptions of pustular epidemics and pox-like endemic illness to comment on the location and age of smallpox's cradle.(1,2,3) Like purported soft tissue and skeletal evidence for the disease,(4,5) the interpretation of early 'poxes' is complicated and has precluded definitive diagnosis. At the same time, recent attempts to use sequence data and a 'molecular clock' to date VARV's emergence,(6,7,8) are subject to ascertainment bias and very distant outgroups (camelpox virus and taterapox virus), complicating VARV's origins.(12) The use of equivocal evidence for historical smallpox as calibration points has likewise confounded emergence estimates.(9) Lack of DNA-based unambiguous VARV detections pre-dating the 17th century has further problematized attempts to better characterize VARV evolution. For this reason, the new study adds considerable time depth and context to the evolutionary history of VARV. The question remains, however, did the VARV strains Mühlemann et al. report cause smallpox?
Although coverage depth of 3 of 4 genomes is low enough (<8x) that identified gene (in)activations remain tentative, the reconstructed VARV strains very clearly sit, as Mühlemann et al. observe, in a previously unknown and extinct sister clade. These VARV strains are not the direct ancestor of the VARV known to have caused smallpox, and importantly their gene content and active/inactive state, as the highest-coverage sample (VK382, 45x) confirms, contrasts sharply with both VARV eradicated in the twentieth century and VARV drafted using 19th- and 17th-century samples.(9,10) Consequently, the disease caused by the newly reconstructed VARV would have differed, perhaps significantly, both clinically and epidemiologically, from the disease we know and recognize as smallpox. The authors' repeated assertion, therefore, that they identified 'smallpox' in Northern Europe, 600-1050, is an overstatement of the evidence they present and regrettably misleading. They do not address the origins of smallpox sensu stricto. The data they present neither put smallpox in Europe more than a millennium ago nor 'disprove' earlier claims regarding when smallpox first appeared in Europe. In fact, the phylogenetic placement of the distinct clade Mühlemann et al. identify, suggests that the clade currently comprising 17th-20th century VARV strains, the only VARV strains unquestionably associated with smallpox, emerged sometime after these two clades diverged. Mühlemann et al. estimate the lower bound of that divergence to date to ca. 1,700-to-1,400 years ago (95% highest priority density). Importantly, VARV that lead to a disease clinically and epidemiologically recognizable as smallpox may have emerged at any point along the branch leading to the 17th-20th century VARV clade. Previous studies have proposed that that emergence took place centuries later.(9,11) Establishing when and where it occurred are critical questions.
Interdisciplinarity will be vital in reconstructing the evolutionary history of VARV as well as the history of smallpox. That the 'Viking' VARV study would have benefitted from more collaboration with historians and archaeologists is apparent in the supposition that the extinct VARV was widespread. Mühlemann et al. impressively recovered VARV sequences from 11 out of 525 individuals spanning 600-1050 and much of northern Europe. False negatives are the norm in ancient DNA and the sample size Mühlemnann et al. describe represents a vast undersampling of the lives lived across those centuries and thousands of kilometres. Yet, with thousands of premodern European skeletons now analyzed for pathogen DNA (12) widespread and persistent VARV would likely have been detected already.
Further, and contrary to what the authors suggest, the written evidence for European pustular epidemics in the 570s and 580s (13,14), which they do not consult, cannot confirm the reconstituted VARV was present in the 6th century or that it was at an any point 'pan-European'. The reported VARV did not originate from known epidemic contexts (Mühlemann et al. "Diverse variola virus (smallpox)", SI: Site Descriptions), its clinical manifestation is completely uncertain, confounding attempts to draw linkages with documented disease, and there were then no settlements in northern Europe with anywhere near the >200,000 inhabitants required to ensure long-term smallpox maintenance.(15) If this VARV sister-clade was indeed widespread and long persistent, and if 6th-century authors did observe the disease that these 4 Northern Europeans would later suffer, there would be, in fact, more reason to argue that the VARV Mühlemann et al. report was epidemiologically distinct from smallpox as we know it. How else to explain the persistence of VARV for centuries in a sparsely populated rural region? That some of the samples the authors present are described as 'predating' the Viking Age, but all are ultimately characterized as 'Viking' is, in this regard, also misleading. There was no singular 'Viking' Age. Bookends of the era vary region-to-region and disciplines delineate the period differently. VK388 died upwards of 190 years, and VK382 possibly 150, before the Viking Age as historians define it. VK281 and VK470 fit the periodization, but identifying these or the other nine VARV sequences as 'Viking' speciously reinforces the conjecture that the drafted VARV was widespread. 'Viking' to many conjures up images of frequent raids in foreign lands, long-distance trade, and exploration, but the great majority of northern Europeans then led far less eventful lives and were not superspreaders of exceptional potential.
Crucially, the identification of a pathogen must not be conflated with the identification of a disease and not all VARV would have caused smallpox as we understand it. As a pathogen evolves over centuries so too can the disease it causes, both clinically and epidemiologically. This is particularly true for viral pathogens, including DNA viruses like VARV. With this in mind, it is worth considering whether the novel 'Viking' virus should be labelled VARV at all. How we dub evolutionary relatives of contemporary pathogens from the past requires more reflection. There is no doubting, though, that as remarkable as the VARV Mühlemann et al. identify is, it was not smallpox.
1) J. Hahn, Variolarum antiquitates, nunc primum e Graecis erutae (Brieg, 1733).
2) J.-J. Paulet, Histoire de la petite vérole I (Paris, 1768).
3) A. Hirsch, "Blattern" in A. Hirsch, Handbuch der historisch-geographischen pathologie I (Stuttgart, 1860).
4) A. McCollum et al., "Poxvirus viability and signatures in historical relics" EID 20 (2014), 177-184.
5) Z. Patterson Ross et al., "The paradox of HBV evolution as revealed from a 16th century mummy" PLOS Pathogens 14 (2018), e1006750.
6) Y. Li et al., "On the origins of smallpox: Correlating variola phylogenetics with historical smallpox records" PNAS 104 (2007), 15787-15792.
7) I. Babkin, I. Babkina, "The origin of the variola virus" Viruses 7 (2015), 1100-1112.
8) G. Zehender et al., "Bayesian reconstruction of the evolutionary history and cross-species transition of variola virus and othropoxviruses" J Med Virology 90 (2018), 1134-1141.
9) A. Duggan et al., "17th century variola virus reveals the recent history of smallpox" Current Biology 26 (2016), 3407-3412.
10) A. Porter et al., "Comment: Characterization of Two Historic Smallpox Specimens from a Czech Museum" Viruses 9 (2017), 276.
11) A. Carmichael, A. Silverstein, "Smallpox in Europe before the seventeenth century: Virulent killer or benign disease?" J Hist Medicine Allied Sciences 42 (2987), 147-168.
12) S. Duchêne et al., "The recovery, interpretation and use of ancient pathogen genomes" Current Biology 30 (2020), R1-R17.
13) Gregory of Tours, Libri historiarum decem, eds. B. Krusch and W. Levison Monumenta Germaniae Historica Scriptores rerum merovingicarum I.1 (Hanover, 1851).
14) Marius of Avenches, Chronica, ed. T. Mommsen Monumenta Germaniae Historica Auctores antiquissimi XI (Berlin, 1894).
15) F. Fenner et al., Smallpox and its eradication (Geneva, 1988).