iScience ll OPEN ACCESS Article Salmonella enterica from a soldier from the 1652 siege of Barcelona (Spain) supports historical transatlantic epidemic contacts Toni de-Dios, Pablo Carrión, Iñigo Olalde, ..., François Balloux, Lucy van Dorp, Carles Lalueza- Fox lucy.dorp.12@ucl.ac.uk (L.v.D.) carles.lalueza.fox@gmail.com (C.L.-F.) Highlights Genome-wide data from two soldiers of the 1652 siege of Barcelona (Spain) Retrieval of an ancient Salmonella enterica serovar Paratyphi C strain Seventeenth century European S. enterica Paratyphi C clusters with colonial Mexican strains Further support for transatlantic dispersals of paratyphoid fever de-Dios et al., iScience 24, 103021 September 24, 2021 ª 2021 The Author(s). https://doi.org/10.1016/ j.isci.2021.103021 iScience ll OPEN ACCESS Article Salmonella enterica from a soldier from the 1652 siege of Barcelona (Spain) supports historical transatlantic epidemic contacts Toni de-Dios,1 Pablo Carrión,1 Iñigo Olalde,1 Laia Llovera Nadal,1 Esther Lizano,1 Dı́dac Pàmies,2 Tomas Marques-Bonet,1,3,4,5 François Balloux,6 Lucy van Dorp,6,* and Carles Lalueza-Fox1,7,* SUMMARY Ancient pathogen genomics is an emerging field allowing reconstruction of past epidemics. The demise of post-contact American populations may, at least in part, have been caused by paratyphoid fever brought by Europeans. We retrieved genome-wide data from two Spanish soldiers who were besieging the city of Bar- celona in 1652, during the Reapers’ War. Their ancestry derived from the Basque region and Sardinia, respectively, (at that time, this island belonged to the Span- ish kingdom). Despite the proposed plague epidemic, we could not find solid evidence for the presence of the causative plague agent in these individuals. However, we retrieved from one individual a substantial fraction of the Salmo- nella enterica serovar Paratyphi C lineage linked to paratyphoid fever in colonial period Mexico. Our results support a growing body of evidence that Paratyphi C enteric fever was more prevalent in Europe and the Americas in the past than it is today. INTRODUCTION Next generation sequencing techniques allow for the recovery of historic pathogens from a variety of sour- ces including medical collections (De-Dios, 2019; Van Dorp et al., 2020), ancient parchments (de-Dios et al., 1Institute 2020; Piñar et al., 2020), ancient ‘‘chewing gum’’ (Jensen et al., 2019), and human archaeological remains of Evolutionary Biology (CSIC-UPF), 08003 (Maixner et al., 2016; Mühlemann et al., 2018; Vågene et al., 2018; Zhou et al., 2018; Key et al., 2020). Barcelona, Spain One of the most widely studied historic pathogen is Yersinia pestis, the bacterial agent of plague, for which 2Antequem. Arqueologia- numerous ancient genomes—almost exclusively sequenced from teeth—have been analyzed spanning Patrimoni Cultural, 08301 different time periods and locations (Bos et al., 2011, 2016; Rasmussen et al., 2015; Andrades Valtueña Mataró, Spain 3Catalan et al., 2017; Rascovan et al., 2019). Institution of Research and Advanced Studies (ICREA), 08010 More recently, ancient genomics has recovered DNA attributable to the pathogen Salmonella enterica, a Barcelona, Spain species comprising over 2500 serovars. These include Salmonella enterica serotypes Typhi (or typhoid fe- 4CNAG-CRG, Centre for ver, which accounts for 75% of cases) and Paratyphi A, B, and C (or paratyphoid fevers) (Cash-Goldwasser Genomic Regulation, Barcelona Institute of Science and Barry, 2018). The former serovar is globally prevalent among paratyphoid fevers (Ochiai et al., 2005; and Technology (BIST), 08028 Stanaway et al., 2019), while Paratyphi B and C are relatively scarce (Achtman et al., 2012; Cash-Goldwasser Barcelona, Spain and Barry, 2018). Typhoid and paratyphoid fevers affect up to 14 million people and cause the death of 5Institut Català de 135,900 people annually (Stanaway et al., 2019). Although a non-life-threatening disease with appropriate Paleontologia Miquel Crusafont, Universitat antibiotic treatment, without it, mortality rates can reach 10–20% (Crump et al., 2008; World Health Orga- Autònoma de Barcelona, nisation, 2018). Typhoid and paratyphoid fevers are particularly prevalent in developing countries of sub- 08193 Cerdanyola del Vallès, Saharan Africa, South East Asia, and South Asia, where they represent one of the leading causes of death Spain 6UCL Genetics Institute, and disability (Naghavi et al., 2017; Vos et al., 2017). University College London, London WC1E 6BT, UK Genomic analyses of ancient strains of S. enterica support its prevalence as a pathogen of humans in the 7Lead contact past, with S. enterica directly observed in association with humans across Eurasia over at least the last 6500 *Correspondence: years (Key et al., 2020). Together with other Salmonella enterica serotypes, Paratyphi C belongs to a sub- lucy.dorp.12@ucl.ac.uk group of S. enterica serovars defined as the Ancient Eurasian Super Branch (AESB), a cluster that comprises (L.v.D.), carles.lalueza.fox@gmail.com phylogenetically close serotypes that infect different wild animals and livestock, as well as strictly human- (C.L.-F.) adapted lineages. The process of host adaptation of the AESB likely accompanied the Neolithic transition, https://doi.org/10.1016/j.isci. when changes in lifestyle and closer interactions with domesticated animals may have led to recurrent 2021.103021 iScience 24, 103021, September 24, 2021 ª 2021 The Author(s). 1 This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ll iScience OPEN ACCESS Article exposure and infections (Key et al., 2020). The characterization of the Paratyphi C genome ‘‘Ragna’’ iden- tified from the sequenced remains of both the bones and the teeth of an 800-year-old skeleton from Trond- heim, Norway, demonstrates the presence of Paratyphi C and invasive salmonellosis in past Medieval Europe (Zhou et al., 2018). More recent evidence for the past epidemic potential of this bacterium is the discovery of Paratyphi C in human burials from colonial Mexico associated with the Cocoliztli epidemics (Vågene et al., 2018). The Co- coliztli (disease or plague in Nahuatl) was a series of epidemics that devastated the native populations of New Spain after the arrival of Spanish colonizers during the early 16th century (Malvido and Viesca, 1985). These epidemics were believed to have been caused by the introduction of infectious diseases by Euro- peans with suggested agents including measles, smallpox, malaria, and unknown haemorrhagic fever (So- molinos d’Ardois, 1970; Malvido and Viesca, 1985; Marr and Kiracofe, 2000). Detection of Paratyphi C in association with victims of the Cocoliztli lends support to the hypothesis that S. enterica was a contributor to epidemics during recent historical times. S. enterica has also been suggested as the causative agent of the famous Plague of Athens (430–426 BCE) following the amplification of two DNA fragments from indi- viduals from the ancient Kerameikos mass grave, dating to that period (Papagrigorakis et al., 2006). How- ever, subsequent phylogenetic assessment of these sequences could not authenticate them as S. enterica (Shapiro et al., 2006). This discrepancy can be partly explained by the difficulties of retrieving significant portions of ancient microbial genomes before the advent of second-generation sequencing technologies. Therefore, although it has been suggested that Paratyphi C was a ‘‘globally’’ distributed pathogen of the past, more evidence based on genome-wide data from different locations and time periods is needed to assess its putative role as an epidemic pathogen. In this study, we generate and analyze DNA from the teeth of two individuals following archaeological ex- cavations at La Sagrera (Sant Martı́ de Provençals), the Spanish army encampment that besieged Barcelona (August 1651 to October 1652) during the Reapers’ War (Monguiló et al., 2012). As well as characterizing the ancestry of these two soldiers, we retrieve from one individual a large number of sequence reads mapping to Salmonella enterica Paratyphi C. We place the resulting Paratyphi C genome in the context of previously published modern and ancient strains including from medieval Norway and colonial Mexico. The accurate dating of this wartime site, coupled with the observations provided by this genome, extends our under- standing of the spatiotemporal breadth of Paratyphi C across Europe and the Americas. RESULTS Ancient DNA sequencing and characterization of La Sagrera soldiers Archaeological excavations at La Sagrera uncovered 576 skeletons of soldiers grouped in shallow graves. Most of these contained between one and ten individuals, but they also included two large pits with 69 and 79 individuals, respectively. Based on the archaeological context of the burials, the soldiers were hastily buried—sometimes still dressed and with boots on—without signs of war injuries and within a short period of time. Archaeologists tend to associate such burial procedures with disease outbreaks (Signoli et al., 2004; Raoult et al., 2006). In contemporaneous chronicles, it is mentioned that both the defenders—within the city—and the besiegers suffered a bout of ‘‘pestilence’’ that has been previously considered to be a plague outbreak (Parets, 1660; Luis and Moya, 1990). We sampled the teeth of two soldiers of the Spanish army, labeled F1691-1810 and F1364-1436, for ancient DNA sequencing (Figure 1A). From these samples, we obtained 246,473,297 and 99,018,404 DNA reads, respectively. The ends of the DNA reads obtained displayed typical ancient DNA (aDNA) damage patterns in a ratio consistent with the estimated age of the site (Figure S1). In both cases, the molecular sex was as- signed to male (Table S1). Based on the presence of different mitochondrial DNA (mtDNA) haplogroups and on heterozygosity levels of the X chromosome, we found that one of the libraries was partially contam- inated (around 9% of human reads). Therefore, we only considered the non-contaminated libraries for the subsequent human population genetics analysis. The mtDNA haplogroups were assigned as U5b1f1a and H2a5a for F1691-1810 and F1364-1436, respectively (Table S2). Y chromosome haplogroups were assigned with differing degrees of resolution to R1 and R1b1a1b1 for F1691-1810 and F1364-1436, respectively. Population genetic affinities of the analyzed individuals Little is known about the geographic origins of soldiers enlisted in the Spanish army at the time, and un- certainties exist over whether these predominately comprised locally or foreign raised troops. Projection 2 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS Figure 1. Population genetic affinities of the soldiers analyzed at the archaeological site at La Sagrera (A) Image of some of the soldiers’ remains found at the site. (B) Genetic data for both individuals (see legend at top right) are projected onto the two main principal components (PC) defined by 1431 present-day West-Eurasian individuals genotyped by human origins array (gray-points). Modern populations more closely related to the Siege of Barcelona soldiers are highlighted (present-day Basques (green); Spanish (yellow); Sardinians (cyan)). of the diversity detected both individuals into a principal component analysis (PCA) conducted on modern West Eurasians placed F1364-1436 neighboring the diversity observed in present day Basques and F1691- 1810 falling close to the genetic diversity observed in present day Sardinians (Figure 1B) (Álvarez-Iglesias et al., 2009; Cardoso et al., 2013). What is now the Basque Country and the island of Sardinia were under dominion of the Spanish crown during the 17th century, implying that all the regions of the Spanish realm contributed with manpower to the war effort. In fact, it is recorded that the fall of Barcelona in 1652 was iScience 24, 103021, September 24, 2021 3 ll iScience OPEN ACCESS Article celebrated across the Spanish possessions, including Cagliari in Sardinia (Caredda, 2015). The analysis of further individuals could provide additional information on the composition of professional Spanish armies during the 17th century; for example, contemporaneous chronicles explain there was a contingent of Irish mercenaries (Morales Ó, 2004). Ancient pathogens screening of La Sagrera soldiers Although contemporaneous assessments suggested that soldiers died of plague during the siege (Luis and Moya, 1990), metagenomic screening of the non-human DNA content from both individuals was unable to identify the presence of Yersinia pestis DNA reads; only 0.00003%–0.000003% of all DNA reads were as- signed to pathogenic Yersinia species (Table S3). The absence of Y. pestis in these individuals does not exclude outbreaks of plague during the siege, noting that identification of Y. pestis tends to require large-scale screening of individuals. Interestingly, 0.01% of sequencing reads were assigned to Salmonella in one of the individuals F1691-1810, with only 0.00001% in the other. After mapping all reads in F1691-1810 against a comprehensive set of Salmonella enterica reference genomes, Paratyphi C was determined as the best representative based on the mean coverage, the percentage of the genome covered by at least a sin- gle read, and the mean read edit distance (Table S4) (Liu et al., 2009). Mapping solely against the Paratyphi C reference chromosome (see STAR Methods), we obtained 24,225 uniquely mapped DNA reads, account- ing for 30% of the reference genome at an average depth of coverage of 0.38X (Table S5). In addition, we obtained 59.02% and an average depth of coverage of 0.91X over the 5 4kb Paratyphi C virulence plasmid (pSPCV), supporting its vertical inheritance and co-evolution with the Paratyphi C lineage (Zhou et al., 2018; Key et al., 2020). The expected coverage for both the chromosome and the plasmid matched the expected theoretical value (Figures 2A, 2B and 2D) (see STAR Methods). Postmortem damage at the end of the Sal- monella DNA reads (up to 15%) and read length distribution, as well as their absence in two extraction blanks (Table S5), supported these sequences as authentically old rather than modern contaminants (Fig- ures 2C, 2E, and 2F). Phylogenetic analysis of the Salmonella enterica Paratyphi C-retrieved sequences After filtering out sequences by depth, heterozygosity, and the presence of low coverage transitions (see STAR Methods), we generated 1,146,808 confidently covered positions, of which 10,250 are single-nucleotide poly- morphisms (SNPs). Finally, we explored the depth of coverage in virulence-related regions associated with the pSPCV (NC_012124.1) plasmid (Liu et al., 2009) and other virulence-associated genes. No apparent variation in coverage over these genes was detected, suggesting that virulence gene composition of pSPCV has been largely maintained from the time of our sample to the present day (Figure S2; Table S6). To place our strain into the wider context of Salmonella enterica, we built a maximum likelihood (ML) tree over the chromosomal alignment that included 415 isolates (both ancient and modern), encompassing the AESB clade and including S. enterica serovar Enteritidis as an outgroup. This placed our 17th-century Bar- celona strain within the diversity of sampled Paratyphi C. Repeating the ML tree using only Paratyphi C strains and Typhisuis as an outgroup, F1691-1810 was positioned in a clade (with 100% bootstrap support) falling basal to all modern Paratyphi C diversity and including the colonial Mexican Paratyphi C strains asso- ciated with the Cocoliztli epidemics (Figures S3A and S4) (Vågene et al., 2018). The only strain basal to this clade is the medieval sample from Trondheim (Norway) dated to 1,200 CE (Zhou et al., 2018). Under the expectation of constant accumulation of mutations with time, the terminal branch lengths of the ancient samples are longer than seen in modern strains, suggesting the possibility of recombination events or downstream incorporation of SNPs which may be erroneously called due to low coverage or postmortem damage (Figure S3B). Due to the high genetic affinity and close temporal range of the Paratyphi C genome obtained from the F1691-1810 La Sagrera individual and the ancient lineage associated with the Mexican Cocoliztli, we explored if the split time of these clades was more consistent with an introduction of S. enterica Paratyphi C into the Americas by Spanish colonizers or if the strain from Barcelona was imported into Spain from the American continent. To create a robust alignment suitable for the application of phylogenetic tip dating, we first set out to prune from the alignment all sites in conflict with the expectation under clonal evolution. In particular, we applied ClonalFrameML (Didelot et al., 2015) to identify all putative homologous recom- bination events between sets of donors and recipients in the Paratyphi C phylogeny. We identified putative recombinant tracts comprising a total of 104 different events totaling 30 kb in length (Figures S5 and S6; Table S7). 4 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS Figure 2. Characteristics of the S. enterica Paratyphi C recovered sequences (A and B) Coverage plot of the Paratyphi C chromosome (A) and plasmid (B); the outer blue circle provides the genomic position, the second outermost circle in gray provides the reference mappability, the third most outermost ring provides the presence of genes (blue), pseudogenes (orange), and RNAs (green) along the reference, and the innermost purple circle provides the mean depth (binned) of coverage. (C) Read length distribution of sequences mapped against Paratyphi C reference genome. (D) Comparison between the calculated and expected coverage and the actual observed coverage for both the pSPCV plasmid and chromosome. (E) aDNA damage patterns (site specific substitutions) at both ends of the DNA reads mapped against the Paratyphi C reference genome. (F) damage patterns at both ends of the DNA reads mapped to the pSPCV plasmid. For both (E) and (F), red provides the C to T substitution frequency; blue provides the G to A substitution frequency; and gray provides the frequency relative to all other substitutions. Temporal relationship of La Sagrera with other taxa in the Paratyphi C clade In addition, we applied a further strict filtering technique, identifying all homoplastic positions in the align- ment given the consensus relationships supported by the ML phylogeny. Such sites may derive from recombination events otherwise undetected by ClonalFrameML or correspond to spurious variant calls iScience 24, 103021, September 24, 2021 5 ll iScience OPEN ACCESS Article Figure 3. Maximum likelihood phylogeny of S. enterica falling within the Ancient Eurasian Super Branch (AESB). (A) ML tree of La Sagrera and 413 S. enterica strains representing the ancient and modern diversity of the AESB from previous studies (Vågene et al., 2018; Zhou et al., 2018; Key et al., 2020). The tree is rooted using S. enterica ser. Enteriditis from a previous study (Alikhan et al., 2018). Serovars are colored at the branch tip. (B) Maximum likelihood phylogeny of the Paratyphi C clade including 124 Paratyphi C modern and ancient strains. The tree includes La Sagrera strain, 119 Modern Paratyphi C strains (Key et al., 2020), 2 strains from Colonial Mexico (Vågene et al., 2018), a strain from Medieval Norway (Zhou et al., 2018), and an S. enterica ser. Typhisuis as an outgroup (Key et al., 2020). Tip points are colored according to their origin as given by the legends. In both cases, phylogenies were constructed from clonal alignments with recombination events and homoplasies removed. included as a result of low coverage or DNA damage. Applying HomoplasyFinder (Crispell et al., 2019), we identified homoplasies by calculating the consistency index, which, for each site, provides a measure of the observed number of changes divided by the minimum number of changes needed to achieve a certain state at tip of a tree. We detected 532 positions in the alignment with a consistency index <0.5 and subse- quently pruned these from the recombination-stripped Paratyphi C alignment (Crispell et al., 2019). The resulting filtered alignment comprised 3680 variant sites fully covered across all of the 124 Paratyphi C strains. We generated a further ML phylogeny using only those positions that passed the filtering criteria. The new tree was topologically congruent; however, the terminal branch length of the ancient tips in the tree was reduced, suggesting that the alignment filtering procedure removed recombinant and/or spurious sites (Figure 3). Previous studies have utilized ancient Paratyphi C genomes to calibrate the evolutionary rate for mem- bers of the AESB (Zhou et al., 2018; Key et al., 2020). Similarly, we estimate the age of the most recent common ancestor (MRCA) (the split times) of our sample at relevant nodes by calibrating the recombi- nation-filtered phylogenetic tree by time. To do so, we first confirmed the existence of a significant tem- poral signal over the alignment by computing the correlation between the root-to-tip genetic distance and the estimated time of sample collection. For our ancient Spanish sample, we set a prior on the tip date to 1652 to coincide with the timing of the siege of Barcelona and the archaeological dating of the site. We recovered a positive correlation (R2 = 0.3) between the root-to-tip distance and sampling time, which was significant following 1 3 104 randomizations of the tip sampling date (p < 1.00 3 104) (Figure 4A). We then applied BactDating to formally estimate the mutation rate over the alignment running until full convergence of the MCMC algorithm (Didelot et al., 2018) (Figure S7). We estimated a mutation rate for this clade over the clonal frame to 1.02 3 104 mutations per site per year. This mu- tation rate corresponds well with the rates reported for other S. enterica serovars (Okoro et al., 2012; Hawkey et al., 2013; Octavia et al., 2015) and is similar to recent estimates following equivalent approx- imations (Didelot et al., 2021). Based on the application of an additive clock model, the estimated rate of evolution leads to an inferred date of the split between the genome identified in medieval Trondheim (Ragna) and the rest of the diversity of Paratyphi C to 1058.43 CE (CI: 949.85 CE–1167.0 CE), a timing consistent with an equivalent analysis by 6 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS Figure 4. Temporal signal in the Paratyphi C clade (A) Regression of the root-to-tip distance and collection date of the samples present in the tree, with 1 3 107 permutations. R2 = 0.44, p < 1.00 3 104. Dotted lines provide the 95% CI. (B) Time-calibrated phylogenetic tree of the curated Paratyphi C dataset. X axis provides the time of collection date. Gray bars at each node denote the 95% CI of the estimated divergence times. Key divergence estimates are highlighted and discussed in the text. Zhou et al. (2018) (1162–1526 ybp, CI: 833–1975 ybp), thereby supporting the robustness of temporal esti- mates when including our low coverage sample. We estimate the divergence time between the 1652 La Sagrera Spanish genome and samples from colonial Mexico from sampled modern diversity to 1240.11 CE (CI: 1134.32 CE–1345.89 CE). Finally, we estimated the split between our 1652 sample and those linked solely to New World epidemics in Mexico to 1301.93 CE (CI: 1194.5 CE–1409.35 CE) (Figure 4B), suggesting a common ancestor of Mexican and the Spanish Paratyphi C around the turn of the 14th century. DISCUSSION Despite an expectation that the two soldiers whose remains we analyzed had died of the plague (Parets, 1660; Luis and Moya, 1990), no reliable traces of Yersinia pestis could be found. This does not preclude plague outbreaks during the Siege of Barcelona; a far larger number of skeletons will have to be screened before Y. pestis can be formally ruled out as the agent of the disease that affected Barcelona’s besiegers. However, we did recover the partial genome of the pathogen Salmonella enterica Paratyphi C in one soldier, which may have contributed to his demise, possibly as part of a wider epidemic affecting the Span- ish army. Indeed, outbreaks of enteric fever are plausible considering the sanitary conditions of a cramped military camp during a siege, as well as the surrounding marshes (that were drained in modern times (Mar- tı́n-Vide, 2015)). Both of these factors and the presence of contaminated water may favor the rapid trans- mission of disease (Gunn et al., 2014). Nonetheless, it is not possible to rule out the co-existence of Y. pestis given that so far we have only analyzed two individuals and the fact that Plague was a well-known disease at the time with clear symptoms (Parets, 1660; Luis and Moya, 1990). Recent studies have demonstrated the possible epidemic character of S. enterica Paratyphi C in historical times. For example, Salmonella Paratyphi C DNA was found in ancient remains from the New World (Vå- gene et al., 2018) associated with mass burials in Mexico dating back to the mid-16th century, plausibly attributed to the Cocoliztli epidemics. The discovery and characterization of the pathogen in a 12th century woman from Trondheim, Norway, further supports the notion that S. Paratyphi C may have noticeably contributed to past pandemics of enteric fever in recent human history (Zhou et al., 2018), despite being rare in Europe and the Americas today. iScience 24, 103021, September 24, 2021 7 ll iScience OPEN ACCESS Article Our work, recovering a 1652 S. enterica Paratyphi C genome from a Spanish soldier, extends the spatiotemporal distribution of S. Paratyphi C to Southern Europe to as recently as the mid-17th century and provides a further observation of the presence of this pathogen in epidemic burial context. The phylogenetic placement of our Paratyphi C genome from 1652 Barcelona basal to the diversity observed in post-contact Mexican strains supports one of two scenarios: (i) an introduction of S. Para- typhi C strains from Europe or (ii) an extensive geographic range of S. Paratyphi C prior to the colonial period. Given that other cases of introduction of pathogens in the Americas during the colonial period have been recently proposed, aided by the analysis of ancient microbial genomics, including parvovirus and hepatitis B virus, leprosy, syphilis, and malaria (Yalcindag et al., 2012; Schuenemann et al., 2013; Guzmán-Solı́s et al., 2020; Majander et al., 2020; Van Dorp et al., 2020), it is plausible that colonial mi- grations contributed to the distribution of past Paratyphi C observed. The continental origin of some infectious diseases putatively associated with colonial settlement, such as syphilis, has been debated for decades with little progress from the analysis of modern strains alone. Syphilis was once thought to have been introduced in Europe from the Americas; however, strains likely predating Columbian times have been recently discovered in Northern European skeletal remains (Schuenemann et al., 2018; Majander et al., 2020). The presence of S. Paratyphi C in mid-17th century Spain is particularly interesting in the context of these historical epidemics due to the fact that these strains are no longer present in the current world diversity of S. Paratyphi C. Our estimated split between the Spanish and Mexican strains to the early 14th century pre- dates Columbian times. Such an observation may be consistent with a possible European reservoir in Spain of the Salmonella lineages which are postulated to have caused widespread epidemics in the native pop- ulations of the Americas during the discovery epoch. Adding support to our estimates, the divergence dates we recover over the wider phylogeny, for instance, the split time of S. enterica Paratyphi C from S. enterica Typhisuis, are highly consistent with independent estimates obtained in studies providing the first characterizations of the timing of evolutionary events in this serovar using ancient DNA (Zhou et al., 2018; Key et al., 2020). This suggests the inclusion of even partial ancient genomes, as we obtained here, can be useful in temporal phylogenomic analyses. Despite uncertainties associated to phylogenetic tip-dating methods and to the low coverage data from ancient skeletal remains, our results reinforce the presence of S. enterica Paratyphi C in Europe for centuries, opening the possibility that this pathogen contributed to other historical pandemics of debated cause—including the one suffered by Spanish be- siegers of Barcelona in 1652. Nevertheless, the likely complex pattern of trans-Atlantic connectivity detected here may only be fully clar- ified with the study of additional ancient samples, both from the Americas, Africa, and Eurasia. One chal- lenge caused by the lack of osteological signals associated with typhoid fevers means that such a survey will need to be conducted blindly, analyzing, like we did in this study, historical mass graves potentially attrib- uted to other common past pandemics such as the plague. Limitations of the study It has not been possible to retrieve the full S. enterica Paratyphi C genome, and the low coverage may have resulted in the inclusion of erroneous sites despite our stringent filters (Figure S8). It is possible this could subtly impact some of our temporal estimates. In the future, we would like to generate DNA from remains of additional soldiers from this site in the hope to identify further, higher coverage, representatives of the clade. STAR+ METHODS Detailed methods are provided in the online version of this paper and include the following: d KEY RESOURCES TABLE d RESOURCE AVAILABILITY B Lead contact B Materials availability B Data and code availability d EXPERIMENTAL MODEL AND SUBJECTS DETAILS B The site d METHOD DETAILS 8 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS B DNA extraction and library preparation B Human DNA mapping B Sex determination and uniparental markers analysis B Contamination estimates B Human population genetics analysis B Pathogens’ screening and Salmonella sequences mapping B Coverage and mappability of Salmonella enterica Paratyphi C B Virulence gene presence B Variant calling and phylogenetics dataset creation B Maximum likelihood tree B Temporal signal exploration d QUANTIFICATION AND STATISTICAL ANALYSIS SUPPLEMENTAL INFORMATION Supplemental information can be found online at https://doi.org/10.1016/j.isci.2021.103021. ACKNOWLEDGMENTS We are grateful to the Museum of the History of Barcelona (MUHBA) for allowing access to these samples. C.L.-F. is supported by a PGC2018-0955931-B-100 grant (MCIU/AEI/FEDER, UE) of Spain; T.M.-B. is sup- ported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 864203), BFU2017-86471-P (MINECO/FEDER, UE), ‘‘Unidad de Excelencia Marı́a de Maeztu’’, funded by the AEI (CEX2018-000792-M), Howard Hughes International Early Career, Secretaria d’Universitats i Recerca, and CERCA Programme del Department d’Economia i Coneixement de la Generalitat de Catalunya (GRC 2017 SGR 880). L.v.D. is supported by a UCL Excellence Fellowship. AUTHOR CONTRIBUTIONS Conceptualization, C.L.-F.; investigation, D.P., L.L.N., E.L., and T.d.-D.; resources, T.M.-B. and C.L.-F.; su- pervision, F.B., L.v.D., T.M.-B., and C.L.-F.; writing-original draft, T.d.-D., L.v.D., and C.L.-F.; writing-review & editing, all authors. DECLARATION OF INTERESTS The authors declare no competing interests. Received: February 26, 2021 Revised: May 14, 2021 Accepted: August 19, 2021 Published: September 24, 2021 REFERENCES Achtman, M., Wain, J., Weill, F.-X., Nair, S., Zhou, Altschul, S.F., Gish, W., Miller, W., Myers, E.W., Andrews, R.M., Kubacka, I., Chinnery, P.F., Z., Sangal, V., Krauland, M.G., Hale, J.L., and Lipman, D.J. (1990). Basic local alignment Lightowlers, R.N., Turnbull, D.M., and Harbottle, H., Uesbeck, A., et al. (2012). search tool. J. Mol. Biol. 215, 403–410. http:// Howell, N. (1999). Reanalysis and revision of Multilocus sequence typing as a replacement for www.sciencedirect.com/science/article/pii/ the Cambridge reference sequence for serotyping in Salmonella enterica. PLoS Pathog. S0022283605803602. human mitochondrial DNA [5]. Nat. Genet. 8, e1002776. https://doi.org/10.1371/journal. 23, 147. ppat.1002776. Álvarez-Iglesias, V., Mosquera-Miguel, A., Cerezo, M., Quintáns, B., Zarrabeitia, M.T., Batzilla, J., Höper, D., Antonenka, U., Cuscó, I., Lareu, M.V., Garcı́a, Ó., Pérez-Jurado, Acman, M., van Dorp, L., Santini, J.M., and Heesemann, J., and Rakin, A. (2011). Complete L., Á, Carracedo, et al. (2009). New population Balloux, F. (2020). Large-scale network genome sequence of Yersinia enterocolitica and phylogenetic features of the internal analysis captures biological features subsp. palearctica serogroup O:3. J. Bacteriol. variation within mitochondrial DNA macro- of bacterial plasmids. Nat. Commun. 11, 193, 2065–2067. haplogroup R0. PLoS One 4, e5112. https:// 2452. https://doi.org/10.1038/s41467-020- pubmed.ncbi.nlm.nih.gov/19340307. 16282-w. Biagini, S.A., Solé-Morata, N., Matisoo-Smith, E., Andrades Valtueña, A., Mittnik, A., Key, F.M., Zalloua, P., Comas, D., and Calafell, F. (2019). Alikhan, N.F., Zhou, Z., Sergeant, M.J., and Haak, W., Allmäe, R., Belinskij, A., Daubaras, M., People from Ibiza: an unexpected isolate in the Achtman, M. (2018). A genomic overview of the Feldman, M., Jankauskas, R., Jankovic, I., et al. Western Mediterranean. Eur. J. Hum. Genet. 27, population structure of Salmonella. PLoS Genet. (2017). The stone age plague and its persistence 941–951. https://doi.org/10.1038/s41431-019- 14, e1007261. in Eurasia. Curr. Biol. 27, 3683–3691.e8. 0361-1. iScience 24, 103021, September 24, 2021 9 ll iScience OPEN ACCESS Article Bos, K.I., Schuenemann, V.J., Golding, G.B., Dabney, J., Knapp, M., Glocke, I., Gansauge, Hawkey, J., Edwards, D.J., Dimovski, K., Hiley, L., Burbano, H.A., Waglechner, N., Coombes, B.K., M.T., Weihmann, A., Nickel, B., Valdiosera, C., Billman-Jacobe, H., Hogg, G., and Holt, K.E. McPhee, J.B., Dewitte, S.N., Meyer, M., Garcı́a, N., Pääbo, S., Arsuaga, J.L., et al. (2013). (2013). Evidence of microevolution of Salmonella Schmedes, S., et al. (2011). A draft genome of Complete mitochondrial genome sequence of a Typhimurium during a series of egg-associated Yersinia pestis from victims of the Black Death. Middle Pleistocene cave bear reconstructed from outbreaks linked to a single chicken farm. BMC Nature 478, 506–510. https://doi.org/10.1038/ ultrashort DNA fragments. Proc. Natl. Acad. Sci. Genom. 14, 800. https://doi.org/10.1186/1471- nature10549. U S A 110, 15758–15763. https://www.pnas.org/ 2164-14-800. content/110/39/15758. Bos, K.I., Herbig, A., Sahl, J., Waglechner, N., Ihaka, R., and Gentleman, R. (1996). R: a language Fourment, M., Forrest, S.A., Klunk, J., de-Dios, T., van Dorp, L., Charlier, P., for data analysis and graphics. J. Comput. Gr. Schuenemann, V.J., Poinar, D., Kuch, M., et al. Morfopoulou, S., Lizano, E., Bon, C., Le Bitouzé, Stat. 5, 299–314. http://www.tandfonline.com/ (2016). Eighteenth century Yersinia pestis C., Alvarez-Estape, M., Marquès-Bonet, T., doi/abs/10.1080/10618600.1996.10474713. genomes reveal the long-term persistence of an Balloux, F., et al. (2020). Metagenomic analysis of historical plague focus. eLife 5, e12994. Edited by a blood stain from the French revolutionary Jean- Jensen, T.Z.T., Niemann, J., Iversen, K.H., Fotakis, R.A. Neher. https://doi.org/10.7554/eLife.12994. Paul Marat (1743–1793). Infect. Genet. Evol. 80, A.K., Gopalakrishnan, S., Vågene, Å.J., Pedersen, 104209. M.W., Sinding, M.H.S., Ellegaard, M.R., Allentoft, Broad Institute (2018). Picard, Available at: http:// M.E., et al. (2019). A 5700 year-old human broadinstitute.github.io/picard/. De-Dios, T., van Dorp, L., Gelabert, P., Carøe, C., genome and oral microbiome from chewed birch Sandoval-Velasco, M., Fregel, R., Escosa, R., pitch. Nat. Commun. 10, 5520. https://doi.org/10. Cardoso, S., Valverde, L., Alfonso-Sánchez, M.A., Aranda, C., Huijben, S., Balloux, F., et al. (2019). 1038/s41467-019-13549-9. Palencia-Madrid, L., Elcoroaristizabal, X., Algorta, Genetic affinities of an eradicated European J., Catarino, S., Arteta, D., Herrera, R.J., Plasmodium falciparum strain. Microb. Genom. 5, Jónsson, H., Ginolhac, A., Schubert, M., Johnson, Zarrabeitia, M.T., et al. (2013). The expanded e000289. https://www.microbiologyresearch. P.L.F., and Orlando, L. (2013). MapDamage2.0: mtDNA phylogeny of the Franco-Cantabrian org/content/journal/mgen/10.1099/mgen.0. Fast approximate Bayesian estimates of ancient region upholds the pre-neolithic genetic 000289. DNA damage parameters. Bioinformatics 29, substrate of Basques. PLoS ONE 8, e67835. 1682–1684. https://doi.org/10.1371/journal.pone.0067835. Depristo, M.A., Banks, E., Poplin, R., Garimella, K.V., Maguire, J.R., Hartl, C., Philippakis, A.A., Del Jun, G., Wing, M.K., Abecasis, G.R., and Kang, Caredda, S. (2015). Un agente de la Corona Angel, G., Rivas, M.A., Hanna, M., et al. (2011). A H.M. (2015). An efficient and scalable analysis hispánica en Cerdeña: Pedro Martı́nez Rubio framework for variation discovery and framework for variant extraction and refinement (1614-1667) y la relación de las fiestas calaritanas genotyping using next-generation DNA from population-scale DNA sequence data. por la rendición de Barcelona (1652). In Las sequencing data. Nat. Genet. 43, 491–501. Genome Res. 25, 918–925. relaciones de sucesos en los cambios polı́ticos y sociales de la Europa Moderna. Universitat Didelot, X., Wilson, D.J., Bryant, D., Quail, M., and Key, F.M., Posth, C., Esquivel-Gomez, L.R., Autònoma de Barcelona, Studia Aurea Cockfield, J. (2015). ClonalFrameML: efficient Hübler, R., Spyrou, M.A., Neumann, G.U., Monográfica, 6, J.G. López and S. Boadas, eds., inference of recombination in whole bacterial Furtwängler, A., Sabin, S., Burri, M., Wissgott, A., pp. 259–269. genomes. PLoS Comput. Biol. 11, e1004041. et al. (2020). Emergence of human-adapted Edited by A. Prlic. https://doi.org/10.1371/ Salmonella enterica is linked to the Neolithization journal.pcbi.1004041. process. Nat. Ecol. Evol. 4, 324–333. Carøe, C., Gopalakrishnan, S., Vinner, L., Mak, S.S.T., Sinding, M.H.S., Samaniego, J.A., Wales, Didelot, X., Croucher, N.J., Bentley, S.D., Harris, Korneliussen, T.S., Albrechtsen, A., and Nielsen, N., Sicheritz-Pontén, T., and Gilbert, M.T.P. S.R., and Wilson, D.J. (2018). Bayesian inference R. (2014). ANGSD: analysis of next generation (2018). Single-tube library preparation for of ancestral dates on bacterial phylogenetic sequencing data. BMC Bioinform. 15, 356. degraded DNA. Methods Ecol. Evol. 9, 410–419. trees. Nucleic Acids Res. 46, e134. https://doi. org/10.1093/nar/gky783. Lazaridis, I., Patterson, N., Mittnik, A., Renaud, G., Cash-Goldwasser, S., and Barry, M. (2018). CDC Mallick, S., Kirsanow, K., Sudmant, P.H., Yellow Book 2018: Health Information for Didelot, X., Siveroni, I., and Volz, E.M. (2021). Schraiber, J.G., Castellano, S., Lipson, M., et al. International Travel, Clinical Infectious Diseases, Additive uncorrelated relaxed clock models for (2014). Ancient human genomes suggest three Edited by G.W. Brunette and J.B. Nemhauser. the dating of genomic epidemiology ancestral populations for present-day Europeans. https://wwwnc.cdc.gov/travel/yellowbook/2020. phylogenies. Mol. Biol. Evol. 38, 307–317. Nature 513, 409–413. Church, D.M., Schneider, V.A., Graves, T., Auger, Fitch, W.M. (1971). Toward defining the course of Lazaridis, I., Nadel, D., Rollefson, G., Merrett, K., Cunningham, F., Bouk, N., Chen, H.C., evolution: minimum change for a specific tree D.C., Rohland, N., Mallick, S., Fernandes, D., Agarwala, R., McLaren, W.M., Ritchie, G.R.S., topology. Syst. Zool. 20, 406. http://www.jstor. Novak, M., Gamarra, B., Sirak, K., et al. (2016). et al. (2011). Modernizing reference genome org/stable/2412116. Genomic insights into the origin of farming in the assemblies. PLoS Biol. 9, e1001091. ancient Near East. Nature 536, 419–424. https:// Fortea, J., de la Rasilla, M., Garcı́a-Tabernero, A., doi.org/10.1038/nature19310. Connor, T.R., Owen, S.V., Langridge, G., Connell, Gigli, E., Rosas, A., and Lalueza-Fox, C. (2008). S., Nair, S., Reuter, S., Dallman, T.J., Corander, J., Excavation protocol of bone remains for Letunic, I., and Bork, P. (2019). Interactive Tree of Tabing, K.C., Le Hello, S., et al. (2016). ‘What’s in a Neandertal DNA analysis in El Sidrón cave Life (iTOL) v4: recent updates and new name? Species-wide whole-genome sequencing (Asturias, Spain). J. Hum. Evol. 55, 353–357. developments. Nucleic Acids Res. 47, 256–259. resolves invasive and noninvasive lineages of https://doi.org/10.1093/nar/gkz239. Salmonella enterica serotype Paratyphi B’. mBio Gunn, J.S., Marshall, J.M., Baker, S., Dongol, S., 7, 1–9. Charles, R.C., and Ryan, E.T. (2014). Salmonella Li, H., and Durbin, R. (2009). Fast and accurate chronic carriage: epidemiology, diagnosis, and short read alignment with Burrows–Wheeler Crispell, J., Balaz, D., and Gordon, S.V. (2019). gallbladder persistence. Trends Microbiol. 22, transform. Bioinformatics 25, 1754–1760. Homoplasyfinder: a simple tool to identify 648–655. https://pubmed.ncbi.nlm.nih.gov/ homoplasies on a phylogeny. Microb. Genomics 25065707. Li, H., Handsaker, B., Wysoker, A., Fennell, T., 5, e000245. https://pubmed.ncbi.nlm.nih.gov/ Ruan, J., Homer, N., Marth, G., Abecasis, G., and 30663960. Guzmán-Solı́s, A., Blanco-Melo, D., Villa-Islas, V., Durbin, R. (2009). The sequence alignment/map Bravo-López, M., Sandoval-Velasco, M., Wesp, J., format and SAMtools. Bioinformatics 25, 2078– Crump, J.A., Ram, P.K., Gupta, S.K., Miller, M.A., Gómez-Valdés, J., de la Luz Moreno-Cabrera, M., 2079. and Mintz, E.D. (2008). Part I. Analysis of data Meraz-Moreno, A., Solı́s-Pichardo, G., et al. gaps pertaining to Salmonella enterica serotype (2020). Ancient viral genomes reveal introduction Liu, W.Q., Feng, Y., Wang, Y., Zou, Q.H., Chen, F., Typhi infections in low and medium human of HBV and B19V to Mexico during the Guo, J.T., Peng, Y.H., Jin, Y., Li, Y.G., Hu, S.N., development index countries, 1984–2005. transatlantic slave trade. bioRxiv. http://biorxiv. et al. (2009). Salmonella paratyphi C: genetic Epidemiol. Infect. 136, 436–448. https://pubmed. org/content/early/2020/06/06/2020.06.05. divergence from Salmonella choleraesuis and ncbi.nlm.nih.gov/17686194. 137083.abstract. pathogenic convergence with salmonella typhi. 10 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS PLoS One 4, e4510. https://doi.org/10.1371/ Okonechnikov, K., Conesa, A., and Garcı́a- Rasmussen, S., Allentoft, M.E., Nielsen, K., journal.pone.0004510. Alcalde, F. (2016). Qualimap 2: advanced multi- Orlando, L., Sikora, M., Sjögren, K.G., Pedersen, sample quality control for high-throughput A.G., Schubert, M., Van Dam, A., Kapel, C.M.O., Luis, J., and Moya, B. (1990). Sociedad y peste en sequencing data. Bioinformatics 32, 292–294. et al. (2015). Early divergent strains of Yersinia la Barcelona de 1651. Manuscrits Revista Història https://pubmed.ncbi.nlm.nih.gov/26428292. pestis in Eurasia 5,000 years ago. Cell 163, Modern. 8, 255–282. 571–582. Okoro, C.K., Kingsley, R.A., Connor, T.R., Harris, Maixner, F., Krause-Kyora, B., Turaev, D., Herbig, S.R., Parry, C.M., Al-Mashhadani, M.N., Kariuki, S., Renaud, G., Slon, V., Duggan, A.T., and Kelso, J. A., Hoopmann, M.R., Hallows, J.L., Kusebauch, Msefula, C.L., Gordon, M.A., de Pinna, E., et al. (2015). Schmutzi: estimation of contamination U., Vigl, E.E., Malfertheiner, P., Megraud, F., et al. (2012). Intracontinental spread of human invasive and endogenous mitochondrial consensus (2016). The 5300-year-old Helicobacter pylori Salmonella Typhimurium pathovariants in sub- calling for ancient DNA. Genome Biol. 16, 1–18. genome of the Iceman. Science 351, 162–165. Saharan Africa. Nat. Genet. 44, 1215–1221. https://science.sciencemag.org/content/351/ https://doi.org/10.1038/ng.2423. Schmieder, R., and Edwards, R. (2011). Quality 6269/162. control and preprocessing of metagenomic Papagrigorakis, M.J., Yapijakis, C., Synodinos, datasets. Bioinformatics (Oxford, England) 27, Majander, K., Pfrengle, S., Kocher, A., Neukamm, P.N., and Baziotopoulou-Valavani, E. (2006). DNA 863–864. J., du Plessis, L., Pla-Dı́az, M., Arora, N., Akgül, G., examination of ancient dental pulp incriminates Salo, K., Schats, R., et al. (2020). Ancient bacterial typhoid fever as a probable cause of the Plague Schubert, M., Ginolhac, A., Lindgreen, S., genomes reveal a high diversity of Treponema of Athens. Int. J. Infect. Dis. 10, 206–214. Thompson, J.F., AL-Rasheid, K.A.S., Willerslev, E., pallidum strains in early modern Europe. Curr. Krogh, A., and Orlando, L. (2012). Improving Biol. 30, 3788–3803.e10. Paradis, E., and Schliep, K. (2019). Ape 5.0: an ancient DNA read mapping against modern environment for modern phylogenetics and reference genomes. BMC Genom. 13, 178. http:// Malvido, E., and Viesca, C. (1985). Epidemia de evolutionary analyses in R. Bioinformatics 35, bmcgenomics.biomedcentral.com/articles/10. cocoliztli de 1576. In Historias, pp. 26–33. https:// 526–528. https://doi.org/10.1093/bioinformatics/ 1186/1471-2164-13-178. www.estudioshistoricos.inah.gob.mx/ bty633. revistaHistorias/wp-content/uploads/ Schubert, M., Lindgreen, S., and Orlando, L. historias_11_24-33.pdf. Parets, M. (1660). De molts successos que han (2016). AdapterRemoval v2: rapid adapter succeı̈t dins Barcelona i molts altres llocs de trimming, identification, and read merging. BMC Marr, J.S., and Kiracofe, J.B. (2000). Was the Huey Catalunya dignes de memòria (1626–1660), Res. Notes 9, 88. https://www.ncbi.nlm.nih.gov/ cocoliztli a haemorrhagic fever? Med. Hist. 44, Barcelona. pubmed/26868221. 363–388. https://pubmed.ncbi.nlm.nih.gov/ 10954969. Patterson, N., Price, A.L., and Reich, D. (2006). Schuenemann, V.J., Singh, P., Mendum, T.A., Population structure and eigenanalysis. PLoS Krause-Kyora, B., Jäger, G., Bos, K.I., Herbig, A., Martı́n-Vide, J.P. (2015). Restauración del rı́o Genet. 2, 2074–2093. Economou, C., Benjak, A., Busso, P., et al. (2013). Besòs en Barcelona. Historia y lecciones. Ribagua Genome-wide comparison of medieval and 2, 51–60. https://doi.org/10.1016/j.riba.2015.07. modern Mycobacterium leprae. Science 341, Patterson, N., Moorjani, P., Luo, Y., Mallick, S., 001. 179–183. Rohland, N., Zhan, Y., Genschoreck, T., Webster, T., and Reich, D. (2012). Ancient admixture in Monguiló, E., Hernandez, J., and Molinas, R. Schuenemann, V.J., Kumar Lankapalli, A., human history. Genetics 192, 1065–1093. http:// (2012). Intervenció arqueològica a l’espai Barquera, R., Nelson, E.A., Iraı́z Hernández, D., www.genetics.org/content/192/3/1065. delimitat pels carrers d’Espronceda, ronda de Acuña Alonzo, V., Bos, K.I., Márquez Morfı́n, L., Sant Martı́, carrer de Josep Soldevila i passeig de Herbig, A., and Krause, J. (2018). Historic la Verneda (Ave Sagrera), Mataró. Piñar, G., Tafer, H., Schreiner, M., Miklas, H., and Treponema pallidum genomes from Colonial Sterflinger, K. (2020). Decoding the biological Mexico retrieved from archaeological remains. Morales Ó, R. (2004). ‘‘Una nación inclinada al information contained in two ancient Slavonic PLoS Negl. Trop. Dis. 12, e0006447. https:// parchment codices: an added historical value. pubmed.ncbi.nlm.nih.gov/29927932. ruido de las armas’’ La presencia irlandesa en los ejércitos españoles, 1580-1818: ¿La historia de un Environ. Microbiol. 22, 3218–3233. https:// éxito? Tiempos Modern. 10, 1–15. pubmed.ncbi.nlm.nih.gov/32400083. Shapiro, B., Rambaut, A., and Gilbert, M.T.P. (2006). No proof that typhoid caused the Plague Mühlemann, B., Jones, T.C., De Barros Plummer, M., et al. (2006). CODA: convergence of Athens (a reply to Papagrigorakis et al.). Int. J. Damgaard, P., Allentoft, M.E., Shevnina, I., diagnosis and output analysis for MCMC. R. News Infect. Dis. 10, 334–335. Logvin, A., Usmanova, E., Panyushkina, I.P., 6, 7–11. http://cran.r-project.org/doc/Rnews/. Boldgiv, B., Bazartseren, T., et al. (2018). Ancient Signoli, M., Ardagna, Y., Adalian, P., Devriendt, hepatitis B viruses from the bronze age to the Price, A.L., Patterson, N.J., Plenge, R.M., W., Lalys, L., Rigeade, C., Vette, T., Kuncevicius, medieval period. Nature 557, 418–423. https:// Weinblatt, M.E., Shadick, N.A., and Reich, D. A., Poskiene, J., Barkus, A., et al. (2004). doi.org/10.1038/s41586-018-0097-z. (2006). Principal components analysis corrects for ‘Découverte d’un charnier de la Grande Armée stratification in genome-wide association studies. en Lituanie (Vilnius, 1812)’. C. R. Palevol 3, Naghavi, M., Abajobir, A.A., Abbafati, C., Abbas, Nat. Genet. 38, 904–909. 219–227. http://www.sciencedirect.com/science/ K.M., Abd-Allah, F., Abera, S.F., Aboyans, V., article/pii/S163106830400065X. Adetokunboh, O., Ärnlöv, J., Afshin, A., et al. Quinlan, A.R., and Hall, I.M. (2010). BEDTools: a (2017). Global, regional, and national age-sex flexible suite of utilities for comparing genomic Skoglund, P., Storå, J., Götherström, A., and specifc mortality for 264 causes of death, 1980- features. Bioinformatics 26, 841–842. https://doi. Jakobsson, M. (2013). Accurate sex identification 2016: a systematic analysis for the Global Burden org/10.1093/bioinformatics/btq033. of ancient human remains using DNA shotgun of Disease Study 2016. Lancet 390, 1151–1210. sequencing. J. Archaeol. Sci. 40, 4477–4482. Raoult, D., Dutour, O., Houhamdi, L., Jankauskas, Ochiai, R.L., Wang, X., Von Seidlein, L., Yang, J., R., Fournier, P.E., Ardagna, Y., Drancourt, M., Skoglund, P., Northoff, B.H., Shunkov, M.V., Bhutta, Z.A., Bhattacharya, S.K., Agtini, M., Deen, Signoli, M., La, V.D., Macia, Y., et al. (2006). Derevianko, A.P., Pääbo, S., Krause, J., and J.L., Wain, J., Kim, D.R., et al. (2005). Salmonella ‘Evidence for louse-transmitted diseases in Jakobsson, M. (2014). Separating endogenous paratyphi A rates, Asia. Emerg. Infect. Dis. 11, soldiers of Napoleon’s Grand army in Vilnius’. ancient DNA from modern day contamination in 1764–1766. https://pubmed.ncbi.nlm.nih.gov/ J. Infect. Dis. 193, 112–120. https://doi.org/10. a Siberian Neandertal. Proc. Natl. Acad. Sci. U S A 16318734. 1086/498534. 111, 2229–2234. http://www.pnas.org/content/ 111/6/2229.abstract. Octavia, S., Wang, Q., Tanaka, M.M., Sintchenko, Rascovan, N., Sjögren, K.-G., Kristiansen, K., V., and Lan, R. (2015). Genomic variability of serial Nielsen, R., Willerslev, E., Desnues, C., and Somolinos d’Ardois, G. (1970). La epidemia de human isolates of Salmonella enterica serovar Rasmussen, S. (2019). Emergence and spread of Cocoliztli de 1545 señalada en un codice. Tribuna Typhimurium associated with prolonged basal lineages of Yersinia pestis during the Med. 15, 85. carriage. J. Clin. Microbiol. 53, 3507–3514, Edited Neolithic decline. Cell 176, 295–305.e10. http:// by N.A. Ledeboer. http://jcm.asm.org/content/ www.sciencedirect.com/science/article/pii/ Sprague, L.D., and Neubauer, H. (2014). Genome 53/11/3507.abstract. S0092867418314648. sequence of Yersinia similis Y228T, a member of iScience 24, 103021, September 24, 2021 11 ll iScience OPEN ACCESS Article the Yersinia pseudotuberculosis complex. (2020). Plasmodium vivax malaria viewed through Epidemiol. Rec. 13, 153–172. www.who.int/ Genome Announc. 2, e00216-14. the lens of an eradicated European strain. Mol. immunization/position_papers/. Biol. Evol. 37, 773–785. https://academic.oup. Stamatakis, A. (2014). RAxML version 8: a tool for com/mbe/article/37/3/773/5614438. Yalcindag, E., Elguero, E., Arnathau, C., Durand, phylogenetic analysis and post-analysis of large P., Akiana, J., Anderson, T.J., Aubouy, A., Balloux, phylogenies. Bioinformatics 30, 1312–1313. Vos, T., Abajobir, A.A., Abbafati, C., Abbas, K.M., F., Besnard, P., Bogreau, H., et al. (2012). Multiple Abate, K.H., Abd-Allah, F., Abdulle, A.M., Abebo, independent introductions of Plasmodium Stanaway, J.D., Reiner, R.C., Blacker, B.F., T.A., Abera, S.F., Aboyans, V., et al. (2017). falciparum in South America. Proc. Natl. Acad. Goldberg, E.M., Khalil, I.A., Troeger, C.E., Global, regional, and national incidence, Sci. U S A 109, 511–516. https://pubmed.ncbi. Andrews, J.R., Bhutta, Z.A., Crump, J.A., Im, J., prevalence, and years lived with disability for 328 nlm.nih.gov/22203975. et al. (2019). The global burden of typhoid and diseases and injuries for 195 countries, 1990– paratyphoid fevers: a systematic analysis for the 2016: a systematic analysis for the Global Burden Global Burden of Disease Study 2017. Lancet of Disease Study 2016. Lancet 390, 1211–1259. Yu, G., Smith, D.K., Zhu, H., Guan, Y., and Lam, Infect. Dis. 19, 369–381. T.T.Y. (2017). Ggtree: an R package for Weissensteiner, H., Pacher, D., Kloss- visualization and annotation of phylogenetic Vågene, Å.J., Herbig, A., Campana, M.G., Robles Brandstätter, A., Forer, L., Specht, G., Bandelt, trees with their covariates and other associated Garcı́a, N.M., Warinner, C., Sabin, S., Spyrou, H.J., Kronenberg, F., Salas, A., and Schönherr, S. data. Methods Ecol. Evol. 8, 28–36. M.A., Andrades Valtueña, A., Huson, D., Tuross, (2016). HaploGrep 2: mitochondrial haplogroup N., et al. (2018). Salmonella enterica genomes classification in the era of high-throughput Zhou, Z., Lundstrøm, I., Tran-Dien, A., Duchêne, from victims of a major sixteenth-century sequencing. Nucleic Acids Res. 44, W58–W63. S., Alikhan, N.F., Sergeant, M.J., Langridge, G., epidemic in Mexico. Nat. Ecol. Evol. 2, 520–528. Fotakis, A.K., Nair, S., Stenøien, H.K., et al. (2018). Wickham, H. (2016). ggplot2: Elegant Graphics Pan-genome analysis of ancient and modern Van der Auwera, G.A., Carneiro, M.O., Hartl, C., for Data Analysis (Springer-Verlag). https:// Salmonella enterica demonstrates genomic Poplin, R., del Angel, G., Levy-Moonshine, A., ggplot2.tidyverse.org. stability of the invasive Para C lineage for Jordan, T., Shakir, K., Roazen, D., Thibault, J., Millennia. Curr. Biol. 28, 2420–2428.e10. et al. (2013). From fastQ data to high-confidence variant calls: the genome analysis toolkit best Wood, D.E., Lu, J., and Langmead, B. (2019). practices pipeline. Curr. Protoc. Bioinform. 43, Improved metagenomic analysis with Kraken 2. Zhou, Z., Alikhan, N.F., Mohamed, K., Fan, Y., and 11.10.1–33. Genome Biol. 20, 257. https://doi.org/10.1186/ Achtman, M. (2020). ‘The EnteroBase user’s s13059-019-1891-0. guide, with case studies on Salmonella Van Dorp, L., Gelabert, P., Rieux, A., De Manuel, transmissions, Yersinia pestis phylogeny, and M., De-Dios, T., Gopalakrishnan, S., Carøe, C., World Health Organisation (2018). Typhoid Escherichia core genomic diversity’. Genome Sandoval-Velasco, M., Fregel, R., Olalde, I., et al. vaccine: WHO position paper - March 2018. Wkly. Res. 30, 138–152. 12 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS STAR+ METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Biological Samples Tooth from a 17th century soldier La Sagrera site (Barcelona, Spain) F1364-1436 Tooth from a 17th century soldier La Sagrera site (Barcelona, Spain) F1691-1810 Deposited Data Sequencing data This paper ENA Project: PRJEB43276 Human reference genome NCBI build 37, GRCh37 Genome Reference Consortium http://www.ncbi.nlm.nih.gov/projects/genome/ assembly/grc/human/ Human Origins Dataset (Lazaridis et al., 2016) https://reich.hms.harvard.edu/datasets Y chromosome haplogroup database version 15.73 International Society of Genetic http://www.isogg.org Genealogy Yersinia enterocolitica assembly ASM25317v1 (Batzilla et al., 2011) https://www.ncbi.nlm.nih.gov/assembly/ GCF_000253175.1/ Yersinia pestis assembly ASM379822v1 Food and Drug Administration https://www.ncbi.nlm.nih.gov/assembly/ GCF_003798225.1/ Yersinia pseudotuberculosis assembly 51108_B01 Wellcome Sanger Institute https://www.ncbi.nlm.nih.gov/assembly/ GCF_900637475.1/ Yersinia similis assembly ASM58251v1 (Sprague and Neubauer, 2014) https://www.ncbi.nlm.nih.gov/assembly/ GCF_000582515.1/ Salmonella enterica serovar Paratyphi C assembly (Liu et al., 2009) https://www.ncbi.nlm.nih.gov/assembly/ RKS4594 GCF_000018385.1/ Plasmid Database (Acman et al., 2020) https://github.com/macman123/plasmid_ network_analysis Salmonella enterica AESB strains whole genome data (Key et al., 2020) http://enterobase.warwick.ac.uk/species/senterica/ search_strains?query=workspace:12971 Salmonella enterica Enteriditis whole genome data (Alikhan et al., 2018) NCBI Sample: ERS217420 Medieval Norway Salmonella enterica Paratyphi C (Zhou et al., 2018) ENA Project: PRJEB19916 whole genome data Colonial Mexico Salmonella enterica Paratypic C whole (Vågene et al., 2018) ENA Project: PRJEB23438 genome data Phylogenetic tree of ‘‘la Sagrera’’ Salmonella This Paper https://microreact.org/project/ sequences within the Paratyphi C diversity G7XdKYogeqhqkWBQLu34t Recombination tracks of ‘‘la Sagrera ‘‘– Paratyphi C This Paper N/A phylogeny Dation of ‘‘la Sagrera ‘‘– Paratyphi C phylogeny This Paper N/A Software and algorithms AdapterRemoval v.2.2.2 (Schubert et al., 2016) https://github.com/MikkelSchubert/adapterremoval BWA v.0.1.17 (Li and Durbin, 2009) http://bio-bwa.sourceforge.net/ SAMtools v.1.12 (Li et al., 2009) http://www.htslib.org/ picard v.2.25.5 (Broad Institute, 2018) https://broadinstitute.github.io/picard/ Qualimap v.2.2 (Okonechnikov et al., 2016) http://qualimap.conesalab.org/ pmdtools v.0.60 (Skoglund et al., 2014) https://github.com/pontussk/PMDtools mapDamage2 v.2.0.8 (Jónsson et al., 2013) https://ginolhac.github.io/mapDamage/ Ry_compute v.0.4 (Skoglund et al., 2013) https://github.com/pontussk/ry_compute/blob/ master/ry_compute.py (Continued on next page) iScience 24, 103021, September 24, 2021 13 ll iScience OPEN ACCESS Article Continued REAGENT or RESOURCE SOURCE IDENTIFIER Haplogrep2 v.2.1.21 (Weissensteiner et al., 2016) https://haplogrep.i-med.ac.at/app/index.html schmutzi v.1.5.4 (Renaud et al., 2015) https://github.com/grenaud/schmutzi/tree/master angsd v.0.925 (Korneliussen et al., 2014) http://www.popgen.dk/angsd EIGENSOFT v.7.2.1 (Patterson et al., 2006; https://www.hsph.harvard.edu/alkes-price/software/ Price et al., 2006) prinseq v.0.20.4 (Schmieder and Edwards, 2011) https://sourceforge.net/projects/prinseq/files/ Kraken2 v.2.1.1 (Wood et al., 2019) http://www.ccb.jhu.edu/software/kraken2/ bedtools v.2.24.0 (Quinlan and Hall, 2010) https://bedtools.readthedocs.io/en/latest/ GATK v.3.7 (Van der Auwera et al., 2013) https://gatk.broadinstitute.org/hc/en-us BamUtil v.1.0.13 (Jun et al., 2015) https://github.com/statgen/bamUtil RAxML v.8.2.4 (Stamatakis, 2014) https://cme.h-its.org/exelixis/web/software/raxml/ R v.4.1.0 (Ihaka and Gentleman, 1996) https://cran.r-project.org/ Ggtree v.1.14.6 (Yu et al., 2017) https://guangchuangyu.github.io/software/ggtree/ iTOL v.6 (Letunic and Bork, 2019) https://itol.embl.de/ ClonalFrameML v.1.2 (Didelot et al., 2015) https://github.com/xavierdidelot/ClonalFrameML HomoplasyFinder v. 0.0.0.9 (Crispell et al., 2019) https://github.com/JosephCrispell/homoplasyFinder BactDating v.1.0.11 (Didelot et al., 2018) https://github.com/xavierdidelot/BactDating ggplot2 v.3.0.0 (Wickham, 2016) https://cran.r-project.org/web/packages/ggplot2/ index.html Ape v.5.4 (Paradis and Schliep, 2019) https://cran.r-project.org/web/packages/ape/index. html CODA v.0.19-4 (Plummer et al., 2006) https://cran.r-project.org/web/packages/coda/index. html RESOURCE AVAILABILITY Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Carles Lalueza-Fox carles.lalueza.fox@gmail.com. Materials availability This study did not generate new unique reagents. Data and code availability d All bacterial and human shotgun sequencing data have been deposited at ENA (project PRJEB43276) and are publicly available as of the date of publication. Accession numbers are listed in the key resources table. d This paper does not report original code. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. EXPERIMENTAL MODEL AND SUBJECTS DETAILS The site A recently excavated archaeological site located in La Sagrera -currently a north-eastern quarter of the city of Barcelona- dates from 1651-1652 CE, during the conflict known as the Reapers’ War (Guerra dels Sega- dors, 1640-1659), the last stage of the Thirty Years’ War (Monguiló et al., 2012). During 1651, the city was under siege by Spanish forces commanded by Juan José de Austria. Despite the efforts of the local garrison and the arrival of French reinforcements, the city capitulated in the spring of 1652. The site consists of multiple mass graves containing mostly male individuals (118 out of 140 studied by physical 14 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS anthropologists so far), with ages ranging from 16 to 40 years and skeletal hallmarks of having undertaken intense physical activities. The individuals analyzed in this study come from the burial F1691 (individual 1810) and F1364 (individual 1436). An upper canine was extracted from each of the skulls and used for ge- netic analyses. The chronology was confirmed by a direct radiocarbon dating of one individual (F1334) that yielded a date of 360+/-30 years BP (Beta-573178) and by the analysis of 152 coins found across the site. Of these, 104 coins correspond to coinage minted in the first half of the 17th century from Philip IV of Spain and also Louis XIII and XIV of France. The remaining 48 coins are a diverse assembly of coins minted during the second half of the 16th century and the first half of the 17th century. METHOD DETAILS DNA extraction and library preparation All DNA extraction and initial library preparation steps (prior to amplification) were performed in a dedi- cated ancient DNA laboratory, physically isolated from the laboratory used for post-PCR analyses. Strict protocols were followed to minimize the amount of human DNA in the ancient DNA laboratory, including the existence of positive air pressure in the clean rooms, the wearing of a full body suit, sleeves, shoe covers, clean shoes, facemask, hair net and double gloving. All lab surfaces, consumables, disposables, tools and instruments were wiped with bleach and ethanol, and UV irradiated before and after use. First, the teeth samples were UV irradiated (245 nm) for 10 minutes and the outermost surface of the teeth was scrapped off with a drill engraving cutter, followed by another UV irradiation for 10 more minutes, in order to exclude surface DNA contamination. Second, approximately 30 mg of tooth cementum were drilled into a fine powder by a Dremel drilling machine at low speed (5000 rpm). The fact that despite these precautions one initial library was contaminated suggests this took place at the excavation itself, probably during handling. Such an observation would suggest it is advisable to excavate, even recent sites such as this one, with anti-contamination protocols such as those implemented in far more ancient samples (Fortea et al., 2008). DNA extraction from tooth powder was performed following the method proposed by Dabney et al. (2013). The teeth-powder samples, including an extraction blank, were added to 1ml of extraction buffer (final con- centrations: 0.45 M EDTA, 0.25 mg/mL Proteinase K, pH 8.0), resuspended by vortexing and incubated at 37 C overnight (24h) on rotation (750-900 rpm). Remaining tooth powder was then pelleted by centrifuga- tion in a bench-top centrifuge for 2 min at maximum speed (16,100 3 g). The supernatant was added to 10mL of binding buffer (final concentrations: 5 M guanidine hydrochloride, 40% (vol/vol) isopropanol, 0.05% Tween-20, and 90mM sodium acetate (pH 5.2)) and purified on a High Pure Extender column (Roche). DNA samples were eluted with 45ml of low EDTA TE buffer (pH 8.0). A total of 35ml of each DNA extract were converted into Illumina sequencing libraries following the BEST protocol (Carøe et al., 2018). Each library was amplified by PCR using two uniquely barcoded primers. After index PCR, libraries were purified with a 1.5x AMPure clean (Beckman Coulter) and eluted in 25ml of low EDTA TE buffer (pH 8.0). Libraries were quantified using BioAnalyzer and sequenced by HiSeq 4000 (Illumina). Human DNA mapping Reads were trimmed of sequencing adapters, filtered for reads of less than 30bp and merged using Adapt- erRemoval with default parameters (Schubert et al., 2016). Clipped reads were then mapped against the human reference genome hg37/19 and the Revised Cambridge mitochondrial DNA reference sequence using BWA aln/samse with the seeding option disabled (Andrews et al., 1999; Li and Durbin, 2009; Li et al., 2009; Church et al., 2011; Schubert et al., 2012). Next, duplicated reads were removed using Picard and only reads with a mapping quality equal or above 30 were considered for downstream analysis (Broad Institute, n.d.). Mapping statistics were calculated using SAMtools and Qualimap2 (Li et al., 2009; Okonech- nikov et al., 2016). Finally, to assess the ancient authenticity of the DNA reads, post-mortem associated DNA damage was estimated using mapDamage2 and PMDtools (Jónsson et al., 2013; Skoglund et al., 2014). Sex determination and uniparental markers analysis Molecular sex was assigned using Ry_compute, a script designed to determine the sex of individuals sequenced at low coverage based on the ratio of reads mapping to each sex chromosome (Skoglund iScience 24, 103021, September 24, 2021 15 ll iScience OPEN ACCESS Article et al., 2013). Mitochondrial haplogroups were determined using haplogrep2 (Weissensteiner et al., 2016). Y chromosome haplogroup determination was performed by manually annotating variants from the Interna- tional Society of Genetic Genealogy (http://www.isogg.org). Contamination estimates Modern human contamination was estimated using two approaches. For mitochondrial contamination, we applied Schmutzi, which calculates the degree of modern contamination by considering the profile of aDNA associated deamination in the sample (Renaud et al., 2015). For assessing nuclear DNA contamina- tion, and considering that both individuals are compatible with being males, we estimated the exogenous DNA contamination based on the heterozygosity of the X chromosome sites using angsd (Korneliussen et al., 2014) To discard contamination in the reagents, two extraction blanks generated along the project were also sequenced. Human population genetics analysis To analyze these 17th-century individuals in the context of present-day human genetic diversity, their genomic data was merged with 1134 West-Eurasian individuals genotyped in the Human origins (HO) array (Patterson et al., 2012; Lazaridis et al., 2014, 2016; Biagini et al., 2019). Principal Component Analysis (PCA) was computed using the modern HO individuals, and the ancient samples were projected onto the first two components (PC1 and PC2) using options ‘lsqproject: YES’ and ‘shrinkmode: YES’ of smartpca built-in module of EIGENSOFT (Patterson et al., 2006; Price et al., 2006). Pathogens’ screening and Salmonella sequences mapping To explore the presence of relevant microbial organisms in the samples we collapsed unique reads from the human-free sequences and removed from the dataset low complexity sequences using Prinseq (Schmieder and Edwards, 2011). Afterward, we applied kraken2 to assign reads against a standard data- base (bacteria, archaea, fungi, protozoa and viral) (Wood et al., 2019). We found almost no evidence of Yer- sinia pestis DNA reads, but one of the samples indicated a significant presence of Salmonella enterica. To validate this signal, and to identify the closest modern representative to our sample amongst a diverse set of S. enterica, we downloaded 839 Salmonella published assemblies from NCBI (as for 03/04/2020) and created a custom database. Human-free reads were then mapped against this database using the local alignment algorithm mem of BWA (Li and Durbin, 2009). We mapped the sequences resultant of the local alignment of the free DNA reads against each Salmonella assembly downloaded independently using BWA’s global alignment algorithm backtrack. The settings used were then optimized for mapping ancient genomes by disabling the seeding option, setting an edit distance value of 0.01 and a gap open penalty value of 2. After that, duplicated sequenced were removed using Picard tools and DNA reads with mapping quality >25 were retained (Schubert et al., 2012). Likewise, the absence of Y. pestis sequences in the sam- ples was checked by mapping against the Y. pestis, Y. pseudotuberculosis, Y. similis and Y. enterocolitica reference genomes, and afterward, possible candidate sequences were explored for taxonomic origin us- ing blast to further asses their authenticity (Altschul et al., 1990). The mapping statistics and presence of post-mortem aDNA damage were determined as described for the human sequences. We determined that the most suitable assembly to map against was the Salmonella enterica subsp. enter- ica ser. Paratyphi C reference genome (Liu et al., 2009) (NCBI accession GCF_000018385.1). The number of mapped reads, the mean coverage, the fraction of the genome recovered, and the mean edit distance are provided in Table S4. The presence of Salmonella related virulence plasmids was searched by mapping the raw sequences against a comprehensive plasmid database (Acman et al., 2020). Coverage and mappability of Salmonella enterica Paratyphi C We determined the mappability of the reads to the Paratyphi C reference genome by mapping k-mers of length 40 to 100 base pairs. The coverage of the sample was determined using bedtools specifying win- dows of 1kb for the chromosome and 100 bases for the plasmid (Quinlan and Hall, 2010). To ascertain if the fraction of the reference recovered with our sample matches the expected coverage given the se- quences’ characteristics we assumed a random distribution matching the actual coverage, assuming sam- ple authenticity, following a previously described approximation (Rasmussen et al., 2015) which calculates the probability of a position being covered given the presence of r reads of length l using the following 16 iScience 24, 103021, September 24, 2021 iScience ll Article OPEN ACCESS Q formula: c = 1  N ri i = 1 ð1  li =gÞ , in which N are the different l1 to lN, read lengths with counts r1 to rN. This value must be corrected by multiplying c for the mappability of the reference genome, otherwise the true coverage may not match the expected. An additional scenario, leading to the expected and actual coverage not matching, is when the reference has a region with no coverage in the ancient sample (Rasmus- sen et al., 2015). Virulence gene presence The sequences of a set of previously described virulence-related genes in Paratyphi C were concatenated to create a reference panel (Connor et al., 2016; Vågene et al., 2018). We mapped our sequences against this reference using the parameters described in the previous section. The percentage of recovered posi- tions for each gene was calculated using bedtools and visualized using ggplot2 (Quinlan and Hall, 2010; Wickham, 2016). Variant calling and phylogenetics dataset creation Additionally, a dataset of 407 modern Salmonella enterica representative of the AESB (Alikhan et al., 2018; Key et al., 2020; Zhou et al., 2020), as well as five historical samples (Vågene et al., 2018; Zhou et al., 2018; Key et al., 2020) were curated. For all modern samples raw read data processing was as described for F1691-1810 except for employing the BWA mem algorithm with default settings instead of the algorithm backtrack, the latter being advised for ancient DNA mapping pipelines. Historical samples were processed as described for the ancient samples. Variants were called from processed published sequences using GATK algorithm UnifiedGenotyper (Depristo et al., 2011; Van der Auwera et al., 2013). For each sample, positions were filtered for minimum coverage of 5X and indels were excluded. This resulted in a high-qual- ity dataset containing 230,597 SNPs. The Colonial strains from Mexico, Tepos-20 and Tepos-37 were not considered due to their low coverage; and Tepos-10 was also discarded due to its high heterozygosity levels. Although it did not reach the coverage thresholds, the la Sagrera genome was used in the subsequent analysis after undergoing strict filtering. We compared different filtering parameters to find a balance between covered positions and possible erroneous calling (Figure S8). Accordingly, we decided to trim 10 bases from the ends of each of the DNA reads using TrimBam (Jun et al., 2015). Following this step all positions in the la Sagrera sample were called with the reference allele considered as valid independent of their coverage. For variant posi- tions, only SNPs found in the high-quality dataset were considered independently again of their coverage; all other positions were discarded. Maximum likelihood tree All positions satisfying the filtering criteria were used to create a consensus genome fasta using the refer- ence genome as template. Those positions that were filtered out in the previous steps were masked in the resultant fasta file. This resulted in a phylogenetic dataset comprising 414 Salmonella enterica AESB and an outgroup (ERS217420). An initial Maximum Likelihood (ML) tree was built using RAxML with the nucleotide substitution model GTRCAT and a strain of Salmonella enterica ser. Enteriditis as an outgroup. The resultant tree was visual- ized using the R package ggtree and iTOL (Stamatakis, 2014; Yu et al., 2017; Letunic and Bork, 2019). Following this, we applied ClonalFrameML, with default parameters, to correct the subsequent tree for the presence of homologous recombination events (Didelot et al., 2015). ClonalFrameML can infer the location of a recombination event in each branch of a ML tree, which can impact the resultant terminal branch lengths. Additionally, ClonalFrameML imputes missing data from the input sequences. We partic- ularly analyzed the presence of recombination events in the Paratyphi C subtree. The transition / transver- sion ratio used (k) was 2.097. After removing the recombination events present in the Paratyphi C lineage, we proceeded to remove all remaining homoplasies from the dataset using HomoplasyFinder (Crispell et al., 2019). A homoplasy is defined as a substitution event which has arisen independently in different separate lineages, and which can alter the tree topology, hence these may arise from undetected recombination or spurious variant calls which may appear in low coverage samples. Homoplasies are determined by their consistency index value which ranges between 0 and 1, with values closer to zero denoting more homoplasic variants (Fitch, 1971). A final ML tree, based on the filtered alignment of 3680 variant sites, was generated using RAxML with the iScience 24, 103021, September 24, 2021 17 ll iScience OPEN ACCESS Article parameters described above. As before, the Paratyphi C phylogeny was visualized using ggtree (Yu et al., 2017). Temporal signal exploration In order to perform phylogenetic tip-calibration it is necessary to first confirm the presence of significant temporal signal in the alignment. Specifically, we tested the Paratyphi C branch of the phylogeny. A tem- poral regression of sampling date (years) against root-to-tip distance was performed using the R package BactDating, which additionally conducts a date randomization test of significance (Didelot et al., 2018). Samples without an assigned collection date were treated as missing. We established the most suitable model by running all models and comparing their DIC and ESS value (Plummer et al., 2006; Didelot et al., 2018) (Table S8). We performed the analysis using the model arc inference model and 1$107 itera- tions. The resulting time calibrated phylogeny was visualized using ggtree (Table S9). QUANTIFICATION AND STATISTICAL ANALYSIS Quantification of model statistics were performed in R v4.1 using the packages CODA and BactDating. Lowest DIC values and highest ESS values were prioritized. A p-value of % 0.05 was considered statistically significant. 18 iScience 24, 103021, September 24, 2021