An ancient genome of Streptococcus pyogenes from a pre-Columbian Bolivian mummy

Domain: biology
Published: 2026-04-13
Editorial outcome: minor_revision

Abstract
Streptococcus pyogenes, or Group A Streptococcus (GAS), is a human pathogen responsible for a range of diseases, from mild infections to severe illnesses. Despite its significance in modern clinical settings, little is known about the pathogen’s evolutionary history or its presence in ancient human populations. Here, we present genomic evidence of S. pyogenes in the pre-Columbian Americas. We analysed a tooth from a naturally mummified individual dating to the Late Intermediate Period (1283–1383 cal AD), housed in the National Museum of Archaeology (MUNARQ) in La Paz, Bolivia. Mitochondrial DNA analysis confirmed the host's Native American ancestry. Shotgun metagenomic sequencing and de-novo assembly enabled the near-complete reconstruction of an ancient S. pyogenes genome displaying close similarity to contemporary strains linked to pharyngitis. The genome contains core virulence genes, but prophages lack streptococcal pyrogenic exotoxins. Phylogenetic analyses place the strain at the base of modern S. pyogenes diversity, and Bayesian analyses indicate that most extant lineages diversified globally within the past ~5,500 years. Our results push back the confirmed presence of S. pyogenes in the Americas by several centuries and suggest that the pathogen circulated among Indigenous populations prior to the European contact.

Introduction
Streptococcus pyogenes, also known as Group A Streptococcus (GAS), is a globally relevant human-adapted pathogen responsible for a wide spectrum of diseases – from mild infections such as pharyngitis (strep throat) to life-threatening conditions including necrotizing fasciitis, toxic shock syndrome, and serious post-infectious complications like rheumatic fever and glomerulonephritis 1, 2, 3. Despite its critical impact on public health, our understanding of its long-term evolutionary history, particularly its presence and virulence in ancient human populations, remains virtually non-existent.

This knowledge gap is especially striking given the pathogen´s current prominence. Recent decades have seen a resurgence in GAS-related illnesses, including scarlet fever outbreaks and rising rates of invasive infections 2, 4. Advances in genomic epidemiology have revealed the emergence of highly virulent lineages and underscored the importance of mobile genetic elements—particularly prophages—in driving GAS evolution and virulence 5, 6, 7, 8. Scarlet fever, once a leading cause of childhood mortality prior to the antibiotic era, is caused by toxigenic strains of S. pyogenes. The disease was mentioned in the London Bills of Mortality (1657–1758) 9, together with other fevers as the second most common cause of death at that time 10. Historically feared for its rapid transmission and high fatality rates, scarlet fever has reemerged in recent years. Since late 2022, several European countries have reported significant outbreaks 11, 12, 13, predominantly afflicting children and often resulting in higher-than-expected mortality rates. This resurgence raises urgent questions about the continued threat posed by this pathogen despite medical advances. The pathogenicity of scarlet fever is driven in part by potent exotoxins (superantigens) that bypass normal immune regulation and provoke excessive immune responses, manifesting in symptoms such as rash, fever, and sore throat 14. Not all S. pyogenes infections progress to scarlet fever, indicating that disease manifestation depends on a combination of host immune factors and the presence of specific bacterial virulence genes.

The closest phylogenetic relative of S. pyogenes is S. dysgalactiae subsp. equisimilis (SDSE), sharing extensive genomic homology and many homologous virulence determinants 15. In contrast to S. pyogenes, that belongs Lancefield serogroup A, SDSE commonly expresses Lancefield group C or G antigens, but still causes a similar spectrum of human diseases ranging from pharyngitis to invasive soft tissue infections 16. Comparative genomics indicates frequent horizontal gene transfer between these species, contributing to overlapping pathogenic potential 17. Mobile genetic elements such as bacteriophages may play a central role in the pathogen’s evolution, reshaping population structure and enhancing virulence by facilitating horizontal gene transfer 3, 5. Prophages in modern GAS strains carry modular genetic cassettes that recombine frequently, allowing for the dissemination of key virulence factors (e.g., pyrogenic exotoxins) across subpopulations 18.

These insights have been pivotal for modern surveillance and vaccine development. Yet, all our understanding of this pathogen is based exclusively on modern strains. Despite the pathogen’s modern worldwide prevalence and clinical importance, no genomic data from ancient S. pyogenes strains have yet been reported, leaving it unclear whether this pathogen has accompanied humans for millennia or emerged more recently. Was it already present in the Americas before European contact, possibly arriving with early human migrations across Beringia? Or did it spread globally only after the onset of European colonialism in the 15th century?

In this study, we provide genomic evidence of S. pyogenes in the pre-Columbian Americas. We analysed a tooth sample from a naturally mummified individual dating to the Late Intermediate Period (1283–1383 cal AD), housed in the National Museum of Archaeology – MUNARQ in La Paz, Bolivia. Through shotgun metagenomic sequencing, we reconstructed the host mitochondrial genome, confirming the Native American ancestry, and applied de-novo metagenomic assembly to recover microbial genomes without reference bias. This approach led to the near-complete reconstruction of an ancient S. pyogenes genome—marking the earliest confirmed occurrence of this pathogen in the Americas.

The ancient genome assembly allowed us to apply modern S. pyogenes typing methods, infer the presence of virulence genes and prophages, and calibrate the pathogen’s phylogeny in time to finally place this strain within a broader phylogenetic and functional framework.

Results
We analysed a tooth sample from a pre-Columbian Andean mummified individual from the anthropological collection of the National Museum of Archaeology – MUNARQ in La Paz, Bolivia. The tooth sample ID:2730 was extracted from a partially mummified human head from a male individual (approx. 18-25 years-old) (Figure 1a). Radiocarbon dating (C14) assigned the individual to the Late Intermediate Period, i.e., 1283–1383 cal AD (with 95.4% probability) (Figure1a, Figure S1, Supplementary Data 1). Stable isotope analysis has a measured C:N ratio of 3,2 showing well-preserved tooth collagen. Carbon δ13C showed -8,6‰ suggesting high maize consumption (e.g., C4-based diet, typical for agricultural populations), and Nitrogen δ15N showed 9,0 ‰ suggesting a low trophic level and a limited intake of animal protein (e.g., little meat or marine food) (Supplementary Data 1). We subjected the sample to DNA extraction and shotgun metagenomic sequencing, which resulted in ~ 6.8 million merged reads (Methods). Analysis of the retrieved DNA showed that 1.72% was human endogenous DNA (Figure 1a, Supplementary Data 2). Next, we analysed the uniparental markers, i.e., mitochondrial DNA (mtDNA) and Y-Chromosome haplogroups in order to confirm the geographic origin and DNA authenticity of the individual. We report a B2 mtDNA haplogroup, that falls within the genetic diversity expected in pre-Columbian Americas 19. Not enough sequencing reads impaired the Y-Chromosome assignment, although molecular sexing 20 successfully confirmed the anthropological sex estimation (male, XY). The overall de-amination of C to T was ~ 3.7% while of G to A was ~ 5.1% (Figure 1b). The low deamination level can be explained by the sample origin, coming from high-altitude environments with dry conditions which have been previously shown to be less prone to DNA degradation 21. Further microbial profiling using MetaPhlAn 4 displayed unusual high abundances of different opportunistic human pathogenic bacterial species, i.e., Clostridium tetani (37.5%), Morganella morganii (29.5%), Streptococcus pyogenes (11.6%), Clostridium botulinum (6.7%), Asaccharospora irregularis (6.4%), and Paeniclostridium sordellii (4.9%) (Figure 1c, Supplementary Data 3). This is in addition to other less abundant commensal oral species, e.g., Desulfobulbus oralis (0.08%) and Actinomyces dentalis (0.005%) (Figure 1c).

The detection of Streptococcus pyogenes was particularly significant; despite its prevalence in modern outbreaks, this pathogen has not yet been detected in ancient times. Hence, we opted to devote particular focus to this species to understand its evolutionary trajectories. Therefore, we decided to perform a metagenomic screening on the Ancient Metagenomics Directory AMDir 22 using MetaPhlAn 4, which revealed the presence of S. pyogenes in human samples from different time periods (ranging from 4000 to 200 years BP) in Europe and Africa, as well as in gorilla samples from museum collections (Figure 1d). However, in nearly all cases, the available genomic data from these samples were insufficient to reconstruct complete genomes and thus precluded further genomic analyses (Supplementary Data 4). Only a 200-year-old dental calculus sample from a gorilla (Gorilla beringei beringei) from Congo yielded sufficient data to reconstruct a high-quality Streptococcus dysgalactiae genome (completeness = 96.26%, contamination = 0.12%) 23 (highlighted in Supplementary Data 4).

We then performed de-novo metagenomic assembly on the Bolivian sample, to reconstruct a reference-free genome of Streptococcus pyogenes. The de-novo assembly resulted in the reconstruction of 7 high-quality bacterial metagenome-assembled genomes (MAGs), including Clostridium tetani, Morganella morganii, and Streptococcus pyogenes (Supplementary Data 5), in addition to other 3 medium-quality genomes including Asaccharospora irregularis (Supplementary Data 5). The completeness of the assembled S. pyogenes genome was 99.98%, with a 0.28% contamination as estimated by CheckM2 (Methods). The total size of the MAG was 1,703,496 bases distributed on 29 contigs ranging in size from 2,244 to 250,505 bases with GC content of 38.44% (Table 1). Considering the minimum information about metagenome-assembled genomes 24, the ancient, assembled genome is of high-quality, given that the completeness is > 90% and contamination is < 5%, and the presence of rRNA genes (23S, 16S, and 5S) and 24 tRNAs (Table 1). The overall coverage of the reconstructed genome is 8-16 X (Figure S2). Similar to human DNA, analysis of ancient DNA damage revealed low levels of deamination of C-to-T and G-to-A (Figure S2).

To systematically and functionally assign the ancient S. pyogenes strain we applied Multi-locus Sequence Typing (MLST), emm, and pilin typing schemes that are commonly used to classify modern S. pyogenes strains 25, 26. MLST of seven housekeeping genes assigned the ancient strain to a distinct sequence profile, that has no direct modern equivalents and possibly represents an ancestral lineage of contemporary strains (Supplementary Data 6). The most similar modern S. pyogenes isolates come from remote aboriginal island communities in tropical Australia sharing three out seven MLST alleles with the ancient strain 27. Additional S. pyogenes typing schemes target the emm gene encoding the surface-exposed streptococcal M protein, and pilin adhesin and backbone genes 28. Both M proteins and pili are key virulence factors and targets of host protective immunity 29. The emm gene is located in the mga regulon with genes encoding additional M-like proteins such as mrp and enn, and scpA encoding a C5a peptidase 30. In the assembled Bolivian genome, the mga regulon genes are distributed to two different contigs, contig 17 that carries the transcriptional activator gene mga and a partial emm gene, and contig 13 that encodes the emm-like gene enn and the C5a peptidase gene scpA (Figure 2b). Considering the size range of 220 to 513 amino acids in modern M protein sequences 31, the emm gene of the ancient Bolivian strain is nearly complete (1161 bp, 387 aa in length) lacking only the 3´ gene part that encodes the D-repeat region of the M protein (Figure S3). This gene part could not be assembled with the currently available data and assembly strategy. Nevertheless, emm typing targeting 5´ regions of the gene was possible and assigned the ancient S. pyogenes strain to the modern emm type 46 (79% nucleotide identity) (see Figure S4, Supplementary Data 7). However, based on the currently existing typing criteria, the ancient strain should be assigned to a unique emm type, having less than 92% nucleotide identity over the first 90 bases encoding the mature M protein to any other emm type 32. Next, we extended our strain typing to the emm pattern analysis that groups emm types into functional patterns (A–C, D, E) that correlate with tissue tropism, ecological niche, and pathogenic behaviour 25. Based on the chromosomal arrangements of emm and emm-like genes (mga-emm-enn-scpA, Figure 2b) and their subfamily (SF) forms (enn SF-1), the ancient Bolivian strain can be associated with the emm pattern group A-C, being most likely a S. pyogenes emm pattern B strain. The analysis of a global collection of Group A Streptococcus revealed that majority of modern strains belong to emm pattern groups D and E, whereas emm pattern B or C strains are nowadays only rarely observed 31 (Supplementary Data 8). Interestingly, the modern strains displaying highest sequence identity to the partial Bolivian emm gene and that share a similar configuration of the mga regulon mostly come from the Americas or Oceania (Supplementary Data 8 and S9, Figure S3).

Further analysis of the fibronectin-binding, collagen-binding, T antigen (FCT) region revealed that the ancient Bolivian S. pyogenes strain harbors a complete Rof regulon on contig 1, including the transcriptional regulator rofA, the fibronectin-binding protein prtF1, the backbone pilin pilB9 (allele 15), the linker pilin pilL9, two pilus biosynthesis enzymes (sortases), and a putative LPxTG-anchored surface protein (fctZ) (Figure 2b). Based on the chromosomal organization of the FCT genes, the absence of an adhesin pilin gene (pilA), and the presence of pilB9, this strain is classified as carrying an FCT-9 region typical of group A streptococci (GAS) 28.

In addition to the surface-associated M- and T- proteins, the ancient S. pyogenes strain holds a diverse array of virulence genes that may enable it to colonize, invade, and evade the host immune system (Figure 2, Supplementary Data 10, S11) 33. The hyaluronic acid capsule (hasABC operon) encoded on contig 28 mimics host tissue and contributes to immune evasion. The full set of genes (gacA-L) of the group A carbohydrate (GAC) gene cluster are present on contig 5 with high homology, i.e., > 99% proteins identity, to modern GAS strains (Supplementary Data 10). Among the virulence associated enzymes, the ancient strain carries the cell-wall anchored proteases C5a peptidase (ScpA, contig 13) and chemokine protease (SpyCEP, contig 20), and the secreted proteases hyaluronidase (HylA, contig 4) and pyrogenic exotoxin B (SpeB, contig 19). In addition, the ancient strain encodes the two important S. pyogenes exotoxins streptolysin O (SLO, contig 1) and streptolysin S (SLS, contig 5), that both promote host cell lysis. Interestingly, the ancient strain carries only two streptococcal superantigens, the streptococcal pyrogenic exotoxin (Spe) type G (SpeG, contig 1) and J (SpeJ, contig 20), that are known to be encoded on the core chromosome of modern S. pyogenes 34, 35. Other important streptococcal superantigens such as SpeA, SpeC or the streptococcal superantigen SSA, that are all encoded in contemporary strains on bacteriophage elements, are missing in the ancient S. pyogenes strain.

There were two high-quality prophages identified in the ancient S. pyogenes genome (Figure 2, Figure S4). Despite the absence of superantigens encoded on these two prophages, one prophage harboured a virulence determinant, Virulence-Associated Protein E (VapE), that was previously associated with virulence in a mouse model of Streptococcus pneumoniae infection 36. Comparative analyses revealed that both prophages are closely related to previously characterized prophages from a large-scale study of Streptococcus 28, 36; however, the prophages in the ancient strain appear to lack a specific region that was present in those near relatives (Figure S4). Follow-up searches recovered four additional prophage contigs corresponding to these missing segments, but it remained ambiguous which two belonged to which prophage, because each prophage shared highly conserved regions that caused assembly breaks. Inspection of the assembly graph confirmed multiple valid pathways for contig arrangement. Despite this uncertainty, both prophages likely represent complete genomes with all necessary machinery for excision and replication. Further screening of publicly available S. pyogenes genomes indicated that these two prophages, or close variants of them, are widely distributed across the modern S. pyogenes phylogeny. This pattern of presence or absence across distinct bacterial clades highlights the role that prophages play in shaping the accessory genome and potentially driving the evolutionary trajectory of S. pyogenes.

Comparative genomic analyses of prophage architecture (Figure S4) uncovered a conserved modular organization— lysogeny, DNA replication, head, tail, and lysis genes—punctuated by mosaicism and recombination events. Interestingly, in some related strains, the DNA replication region contained a larger number of shorter predicted CDS, suggesting partial degradation that could lead to non-infective or “cryptic” prophage states. The detection of structural rearrangements, evidence of circularization, and functional gene conservation in these prophages collectively underscores their historical and ongoing relevance to S. pyogenes evolution and virulence potential.

Additional comparative analysis identified two antibiotic resistance genes in the ancient S. pyogenes genome using ABRicate v1.0.0 with the CARD database (Supplementary Data 12). The lmrP gene, encoding a multidrug MFS transporter, exhibited 99.9% identity and 97.5% coverage, and is associated with resistance to macrolides, lincosamides, streptogramins, and tetracyclines. Furthermore, the mef(E) gene, another macrolide efflux pump, was recovered with 95.9% identity and 82.4% coverage. Both genes were complete, with no evidence of frameshifts or premature stop codons, suggesting they were intact in the ancient strain. Their functional properties, however, remain to be determined.

The expression of the above-mentioned virulence genes, the modulation of antibiotic resistance and stress-response genes, and the control of prophage and other mobile genetic element genes is mediated by specific two-component regulatory systems (TCS) 37. Group A S. pyogenes possesses up to 14 TCS to adapt the bacterial behaviour to host and environmental conditions. The ancient Bolivian strain encodes for 11 TCS including the classic complete covR/S TCS that acts as a global control for many virulence genes (Supplementary Data 13). The three TCS salKR, silAB, and srtRS that are missing in the ancient strain are described to be completely absent in other modern GAS genomes as well 37.

To further investigate the evolutionary placement of the ancient S. pyogenes strain, we performed a phylogenetic analysis with the HC254 clusters of modern S. pyogenes, S. dysgalactiae and S. canis assemblies from the EnteroBase database (Supplementary Data 14). This analysis was based on the core genome (Please refer to the Methods for details), of 605 conserved genes in >99% of the genomes at 95% similarity. The phylogenetic tree placed the pre-Columbian Bolivian strain basal to the modern S. pyogenes diversity (Figure 3a). The second Streptococcus genome from the 200-year-old gorilla was assigned to the closest phylogenetic relative of S. pyogenes, S. dysgalactiae subsp. equisimilis (SDSE) (Figure 3b). To further support the phylogenetic placement of both ancient strains, we aimed to identify species-specific KEGG modules using pangenome analysis (Methods). The final set of KEGG modules clearly indicate metabolic differences between the three species. Either complete KEGG modules or certain genes within modules are present in one species but not the other (Figure 3b and S5, Supplementary Data 15, S16). As previously reported by Xie and colleagues 17 in their pangenome comparison of modern S. pyogenes and SDSE strains, we could confirm in our dataset the absence of the modules encoding glycogen biosynthesis (M00854) and threonine biosynthesis (M00018) in both the ancient Bolivian and modern S. pyogenes strains. Furthermore, unique to the ancient and modern S. pyogenes strains were the presence of multiple genes encoding V/A-type ATPases (M00159). There are additional differences in the metabolic gene set between the three Streptococcus species, such as unique presence of possible cydAB genes that encode subunits of the Cytochrome bd ubiquinol oxidase in SDSE, S. canis and the gorilla strains. However, we do not see any indications for genes that occur only in the ancient Bolivian or gorilla strain.

In summary, the phylogenetic placement of the ancient Bolivian strain basal to the modern S. pyogenes and the gorilla strain into the SDSE diversity is highly supported on pangenome level by the presence of species-specific genes in both ancient genomes.

Finally, to obtain improved estimates of the evolutionary history of S. pyogenes and to infer divergence times among its closest related species, we applied a Bayesian phylogenetic framework using BEAST. Time estimates of the full dataset under most fitting analytical parameters (see Methods) returned an age for the split of S. pyogenes and S. dysgalactiae spp. at 137,264 years (95% HPD: 124,428 to 149,847 cal years BP) (Figure 3a). The timing of the Most Recent Common Ancestor (MRCA) of all S. pyogenes strains (ancient Bolivian strain and modern known S. pyogenes) is estimated to be 10,270 years (95% HPD: 9314 to 11,218 cal years BP). The diversification of most modern S. pyogenes lineages took place in the past ~5,500 years of human history (Figure 3C). Time estimates are robust to demographic model variation and outgroup exclusion but become more recent when subsampling the dataset (Figure S6c). In all cases the split between the ancient Bolivian strain and the remaining modern S. pyogenes post-date by more than 15,000 years the supposed migration of humans into the Americas from the Bering strait at circa 22,000 years ago.

Methods
Ethics & Inclusion statement
This research was conducted in accordance with ethical guidelines for the treatment of human remains 55. The mummified human remains analysed in this study are curated by the National Museum of Archaeology – MUNARQ in Bolivia. All necessary permits and authorizations for the study were obtained from the Bolivian Ministry of Cultures, Decolonization, and Depatriarchalization (AUTORIZACIÓN MDCyT – UDAM No. 017/2018). The mummified head is part of a historical collection composed of human remains and individuals recovered from diverse geographical regions. Based on current records and consultations, no specific Indigenous community has claimed affiliation with or recognized a direct connection to this individual.

Sample collection
The National Museum of Archaeology – MUNARQ in La Paz, houses the largest collection of mummified human remains and human skulls in Bolivia 56. Minimally invasive sample collection was performed. We collected a tooth sample from a partially mummified head belonging to the anthropological collection (Sample ID: 2730 / Museum ID: MUNARQ 2-87) (Figure 1a).

Burial context and palaeopathological overview
The mummified human remains housed in the MUNARQ were recovered from the Chullpas, or funerary towers, spread across the Bolivian Altiplano. During the Late Intermediate Period (1100 – 1450 AD), these individuals were interred within funerary bundles and placed inside these structures. While there are some accounts of anthropogenic mummification practices, most of the remains found in this region likely underwent a process of natural mummification over time.

The palaeopathological examination of individual ID:2730 suggests that this male individual (approx. 18-25 years-old) experienced intentional cranial modification. Dental assessment reveals dental chipping and localized-periapical infections involving both mandibular first molars (M1). A linear fracture is observed on the inferior portion of the occipital bone, extending posteriorly from the foramen magnum along the sagittal plane. Furthermore, post-mortem fractures of the cervical vertebrae, specifically involving the axis, show evidence of mechanical tension consistent with the deliberate removal of the cranium after death.

DNA extraction and sequencing
All laboratory work was performed at the ancient DNA facilities of the Eurac Research, Institute for Mummy Studies in Bolzano, Italy. DNA sampling preparation was conducted following the protocol developed specifically optimized for the recovery of ancient pathogens 57. The procedure involved the transverse sectioning of the dental crown, followed by the mechanical drilling of the pulp chamber, which serves as a reservoir for endogenous pathogen DNA. This area is targeted because the highly vascularized pulp tissue historically trapped circulating pathogens, preserving their genetic material within the tooth structure. DNA extraction was performed on 50 mg of drilled tooth powder using EDTA-based demineralization followed by DNA precipitation using the linear polyacrylamide (LPAA) precipitation method 58. The extracted DNA was then converted into genomic DNA libraries following a protocol specific for ancient DNA 59. The method focuses on "double-indexing" to maximize sample identification during high-throughput sequencing. It involves blunt-end repair, adapter ligation, and a specialized indexing PCR to create libraries ready for Illumina platforms. The prepared libraries were subjected to Illumina next generation sequencing using Illumina HiSeq 2500.

Sequencing quality controls
After demultiplexing of the reads, the sequencing adapters and low-quality reads were trimmed using fastp v0.23.1 60, allowing for minimum base quality of 20. The quality filtered reads were merged using fastp with minimum overlap between forward and reverse reads of 10 bases.

Human DNA analysis
For the analysis of the human DNA, we mapped the merged reads against the reference human genome (build hg19), using BWA-aln v0.7.17 61 (using the following parameters: -l 1000 -n 0.04 -o 1 -e -1), then filtered for minimum mapping quality of 25 using SAMtools v1.17 62. The mapping quality and coverage were then checked using Qualimap v2.2.2d 63. The ancient DNA damage pattern was checked using MapDamage 2 v2.2.1 64. The mitochondrial DNA contamination was checked using Schmutzi v1.5.7 65. The damaged sites of aDNA reads in the BAM file were then rescaled using MapDamage 2 and used to call a consensus sequence using Schmutzi script “log2fasta”, for mitochondrial haplogroup assignment using Haplogrep3 v2.4.0 66. Molecular sex identification was performed from shotgun sequence data 20.

Radiocarbon dates and stable Isotopic data
The remaining sample ID:2730 was sent for radiocarbon analysis to the Department of Anthropology, AMS Radiocarbon Lab, Institute of Energy and the Environment, Penn State University, USA. The date was then used to assign the individual to a cultural period (Supplementary Data 1, Figure S1). Radiocarbon dates were calibrated using the Oxcal software v4.4.4 67 applying the (SHCal20) for the Southern Hemisphere calibration curve 68. Additionally, stable isotopic analysis of carbon (δ13C) and nitrogen (δ15N) in the dental sample were analysed to reconstruct the individual's diet, environmental conditions, and potential reliance on specific crops or animal resources during their lifetime (Supplementary Data 1).

Microbial DNA analysis and metagenomic assembly
The general microbial taxonomic classification was carried out using MetaPhlAn 4 v4.0.6 69, with the database version of “June23”. Then, to screen the other publicly available ancient metagenomic datasets, we used the tool AMDirT v1.4.5 70 to retrieve the raw fastq files from AMDir 22, which were analysed as described above.

Then, we used MEGAHIT v1.2.9 71 and metaSPAdes v3.15.5 72 to perform de-novo metagenomic assembly on the Bolivian sample as well as other samples from the ancient metagenome directory that displayed presence of S. pyogenes. For each sample, the resulting contigs were filtered for minimum length of 1000 prior using multi-binning approach employing the metagenomic binners MetaBAT2 v2.12.1, MaxBin2 v2.2.7 73, and CONCOCT v1.1.0 74. The outputs of different binners were then combined using DAS Tool v1.1.6 to select the highest quality genomes 75. To assess the assembly quality in terms of completeness and contamination, we used CheckM2 v1.0.2 76. We further used Genome Database Taxonomy toolkit GTDBtk v2.1.1 77 for taxonomic assignment of the bins.

Improving the quality of ancient Bolivian S. pyogenes genome
The initial binning process yielded two bins identified as Streptococcus pyogenes: one assembled with MEGAHIT and the other with metaSPAdes. To refine these bins, we used Kraken2 v2.17.1 (with the standard database 78) to classify all contigs and removed any contigs not assigned to the genus Streptococcus. This decontamination step excluded five contigs in total (two were unclassified and three classified as Proteus spp.) Notably, CheckM2 quality metrics for the bins remained unchanged after this filtering. Following decontamination, the cleaned bins were merged and reassembled using the long-read assembler Flye v2.9.3 79 to maximize genome completeness and quality.

Functional annotation and typing of the Bolivian S. pyogenes genome
To perform functional annotation of the ancient S. pyogenes genome, we used prokka v1.14.6 80, using the genome Streptococcus pyogenes M1 GAS (Strain: SF370; assembly ID: ASM4323122v1) as a reference. For the specific virulence factors, the annotated proteins were further compared against the NCBI non-redundant protein database using BLASTp v2.12.0 and the highest matches based on bitscore, query coverage, and sequence similarity were reported (Supplementary Data 10).

Collection of modern comparative genomes
To perform a comprehensive phylogenetic and comparative analysis, we utilized EnteroBase (https://enterobase.warwick.ac.uk/), a platform for exploring genomic epidemiology 81. Specifically, we downloaded genomes of Streptococcus pyogenes and other two closely related species, i.e., Streptococcus dysgalactiae and Streptococcus canis. In EnteroBase, Hierarchical Clustering (HC) groups bacterial isolates based on core-genome Multi-locus Sequence Typing (cgMLST) allelic profiles, with cluster levels (e.g., HC5, HC10, HC20, HC245) defined by the maximum number of allowable allelic differences, lower thresholds representing tighter, more recent genetic relationships. These clusters facilitate outbreak detection, epidemiological tracking, and population structure analysis by organizing isolates into hierarchically nested groups of increasing genetic diversity 82. We opted to use the HC245 after exclusion of genomes with missing origin and isolation date in the metadata, which resulted in retrieving 149 genomes for S. pyogenes, 10 S. dysgalactiae, and 6 genomes for S. canis. A full list of genomes with metadata is available in Supplementary Data 14. These genomes were used further for phylogenetic placement, pangenome analysis, and functional comparison.

Typing, antibiotic resistance gene, and virulence factor analysis
Multi-locus sequence types (MLST) and emm-types were determined from the Bolivian strain as well as the modern genomes using MLST v2.23.0 (https://github.com/tseemann/mlst) with the PubMLST database 83, and the emmtyper command line tool (https://github.com/MDU-PHL/emmtyper), respectively. For emmtyper, we used the following parameters to allow detection of both emm and enn proteins: “--percent-identity 75 --culling-limit 1 -f verbose --mismatch 35 --gap 4”. The investigation of virulence factors and antibiotic resistance genes was carried out with ABRicate v1.0.0 (https://github.com/tseemann/abricate) tool by using the Virulence Factor Database (VFDB) 84, and the Comprehensive Antibiotic Resistance Database (CARD) 85 with 80% coverage (“--mincov”) and 80% identity (“—minid”) parameters.

Phylogenetic and molecular clocking analysis
Then, gene prediction and annotation were carried out using the prokka pipeline 80. To identify the core genome, we used Panaroo v1.5.2 86, a bioinformatics tool designed for bacterial pan-genome analysis. Panaroo identifies core genes (present in all strains), accessory genes (present in some but not all strains), and strain-specific genes. For this study, we set the similarity threshold at 95% for amino acid sequences to define the core genome. The core genes identified by Panaroo were aligned using MAFFT v7.525 87, a highly accurate multiple sequence alignment tool. The aligned sequences were then concatenated into a single alignment for phylogenetic analysis.

Divergence times have been estimated using BEAST.X. v1.10.5 88 and the BEAGLE library 89. We preliminary tested between HKY and GTR replacement model and between a Constant and a Bayes skyGrid coalescent model by comparing marginal likelihoods: GTR+G was favoured over HKY (D=+66292), and Constant was favoured over Bayesian SkyGrid (D=1865). According to model selection we therefore used the following priors for the main analyses: 1) GTR replacement model coupled with 4 categories of Gamma distribution; 2) A Constant Coalescent demographic/tree model; 3) The age of all samples for tip calibrations 90 by sampling uniformly from date uncertainty using as date the year of sample collection, an uncertainty of one years for all sample, and an uncertainty of 50 years for the Bolivian mummy (see suppl. Supplementary Data 14 for age of samples); 4) Strict clock with a normal prior rate distribution with mean 1.8x10-6 and sd 7x10-7: this allows covering between 4.3e-7 and 3.2e-6 at the 97.5 % a range that is compatible with the rates of 1.8x10-6 previously identified for S. pyogenes 91 and safely covers the posterior proposed for other Streptococcus, S. dysgalactiae (5.8x10-7) and S. canis (8.6x10-7) in respectively 92, 93. To test for outgroup effect we repeated the analysis using a dataset composed of only S. pyogenes and to test for taxon sampling we repeated the analysis using a reduced set of 50 samples including both outgroups and S. pyogenes; for both datasets we used the same combination of priors of the main analysis (GTR+G, Constant coalescent, all tip priors). All MCMC were run for a minimum of 50 million generations until convergence. Convergence was inspected using Tracer by checking that all parameters had an ESS>200 and setting a burn-in that maximised ESS. Analysis was run in double to assess further convergence, and trees were pulled from one of the two runs after burn-in using HIPSTER method in TreeannotaorXv10. The resulting phylogenetic tree was visualized using iTOL v7 (Interactive Tree of Life), a platform for annotating and exploring phylogenies 94.

Phage genome analysis
Further examination of all 29 assembled contigs revealed the gene composition (by contigs) including the presence of bacteriophages in contig 5 and in contig 26. Prophage predictions were made with PHASTEST 95, geNomad v1.7.4 96, VIBRANT v1.2.1 97, and PhageBoost v0.1.3 98. Manual inspection of the prophage predictions revealed little difference between the prophage miners and only the geNomad predictions were retained. Predictions were assessed with CheckV v1.0.3 99 and two near-complete prophages were retained. Annotation was first performed with Pharokka v1.7.5 100 with -g prodigal-gv 101, 102 that uses tRNAscan-SE 2.0 103, Aragorn 104, CRT 105, PHROGs 106, VFDB 84, CARD 107, MMseqs2 108, PyHMMER 109, INPHARED 110 and MASH 111. Subsequent annotation was performed with PHOLD v0.2.0 112 that implements Foldseek 113 and ProstT5 114, and Phynteny v0.1.12 115.

The two prophage predictions were compared to INPHARED (January, 2024) 110 using MASH v2.3 111. Comparison of the prophages to their nearest hits revealed a seemingly missing section on both of the prophages. We then used MASH to compare the nearest hits to other contigs in the MAG, revealing four other contigs that were likely part of the two prophage genomes, however it remains unclear which of these additional contigs corresponds to which of the two prophage regions.

Prophage sequences were extracted from the modern genomes from Enterobase (BLAST reference database) using geNomad v1.7.4 96. An all vs all comparison of prophages was made using MASH-dist to determine which reference S. pyogenes genomes contained similar prophages to those in the MAG. Similar prophages from bacterial genomes in the same clade as our MAG (Figure 3, Figure S4) were extracted and de-replicated following Minimum Information about an Uncultivated Virus Genome (MIUViG) standards (95 % ANI over 85 % length) using blast v2.14.0+ alongside the anicalc.py and aniclust.py scripts described in the CheckV documentation 99. The species representative reference prophages were compared to our prophage using Clinker v1.3.2 116.

Analysis of functional modules in the pangenome
The resulting protein sequences (faa files) from prokka from each genome were used as input for further functional annotation. Functional annotation was performed using KofamScan v1.3.0 117 to assign KEGG Orthology (KO) identifiers to predicted genes. KofamScan searches were conducted against the KEGG profile database using a threshold score of 1 and output in detailed TSV format. KO assignments were made based on profile Hidden Markov Models (HMMs) with pre-defined score thresholds to ensure high-confidence annotations.

Assigned KO identifiers were mapped to KEGG modules using the official KEGG ko_to_module mapping file (accessed January 2025). For each sample, unique KO identifiers were extracted from KofamScan results and cross-referenced with the ko_to_module database to determine module presence. Module abundance was quantified as the count of unique KO identifiers belonging to each module within each sample.

Statistical analysis and visualization
All statistical analyses and visualizations were performed in R version 4.5.1 using the tidyverse (v2.0.0) and pheatmap (v1.0.13) packages. Module count matrices were constructed with modules as rows and samples as columns, with values representing the number of unique KOs detected per module per sample. Modules present in fewer than 5 samples were excluded from downstream analysis to reduce noise.

To identify modules that differ significantly between the three Streptococcus species (S. pyogenes, S. dysgalactiae, and S. canis), we applied the Kruskal-Wallis rank sum test to each module independently. P-values were adjusted for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) correction, with modules showing FDR < 0.001 considered statistically significant.

Module abundance patterns were visualized using hierarchical clustering heatmaps with Euclidean distance and complete linkage clustering. Count data were log10-transformed (log10(count + 1)) and row-scaled (z-score normalization) for visualization. Sample metadata including species, geographic continent, and host type were displayed as color-coded annotation bars.

Data availability
Sequencing data and the assembled genome are available at the European Nucleotide Archive (ENA) under ENA: PRJEB91735. The sequencing reads are available at the Sequence Read Archive under accession ERR15308372 and the assembly under accession number GCA_982145515.1. Source data is provided with this article.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Noteworthy Results:\nThis work presents compelling findings for both the ancient DNA and modern microbiology/epidemiology communities. It\nsuccessfully identifies Streptococcus pyogenes in a pre-Columbian context while providing a detailed genomic\nreconstruction of the ancient strain, complete with virulence factors and prophage content. The study offers valuable insights\ninto disease dynamics among ancient Amerindian populations, including a plausible hypothesis about systemic infection or\nsepsis as a potential cause of death.\nSignificance and Originality:\nThe study represents a significant advancement as the first genomic evidence of S. pyogenes in pre-Columbian America,\nextending our knowledge of the pathogen's historical presence by centuries. It effectively bridges gaps in understanding the\nevolutionary history of this clinically important pathogen, while aligning with established literature on GAS virulence and\nphage-mediated evolution (e.g., Brouwer et al., 2023; McShan et al., 2019). The complete absence of prior ancient genomic\ndata for S. pyogenes underscores the originality of these findings.\nSupport for Conclusions:\nThe authors' conclusions are well-supported through comprehensive analyses including metagenomic assembly,\nphylogenetic placement, virulence profiling, and prophage characterization. The data establish convincing connections\nbetween the ancient strain and modern lineages while demonstrating remarkable conservation of pathogenic mechanisms\nover time.\nFlaws and Limitations:\nWhile no major flaws were identified in the data analysis and the methodology meets field standards, certain aspects merit\nconsideration. The phylogenetic analysis, while technically sound, requires cautious interpretation of geographic dispersal\npatterns due to the inherent weak phylogeographic signal in S. pyogenes. The shared prophages might reflect horizontal\ngene transfer rather than direct ancestral inheritance.\nThe mixed infection hypothesis presents a compelling case supported by metagenomic evidence of multiple\ncontemporaneous pathogens. This approach proves particularly valuable for detecting ancient systemic infections that\ntypically leave no skeletal markers. Future investigations incorporating biochemical markers, such as proteomic signatures\nof inflammation, could further validate such findings in cases of ancient sepsis.\nThe detection of mreA, lmrP, and mef(E) resistance genes in the ancient genome is intriguing, though the clinical relevance\nin this pre-antibiotic-era context remains uncertain without functional validation through phenotypic assays. Similarly, while\nthe prophage reconstructions suggest functional potential, assembly ambiguities and reliance on computational predictions\nsomewhat limit definitive conclusions about their activity in vivo.\nMethodology and Reproducibility:\nThe methods are sound and meticulously detailed, fully meeting the expected standards for ancient DNA and metagenomic\nresearch. The complete workflow—from DNA extraction through phylogenetic analysis—is presented transparently with\nappropriate controls for contamination and damage patterns, ensuring reproducibility.\nFinal Assessment:\nThis study delivers unprecedented insights into the ancient genomic landscape of S. pyogenes, supported by robust\nmethodology and high-quality data. While certain interpretations regarding phylogeography and prophage activity would\nbenefit from additional caveats, the work undoubtedly establishes a new benchmark in pathogen paleogenomics. It will\nundoubtedly inspire and inform future research into the evolutionary history of human pathogens.\nMinor comments:\nIn figure 1d, I would place a blue dot, to keep consistency, rather than colouring the whole country as blue.\nFigure S2 does not contain any black star as stated in the caption.", "reviewer_2": "The authors provide a high quality genome sequence of S. pyogenes from a Bolivian pre-Columbian mummified individual. I\ncannot comment on the accuracy of the radiocarbon dating as it is not within my area of expertise, yet the pathogen genome\nanalyses seems robust where a 1.7Mb S. pyogenes draft genome is reported. Thank you for providing access to the raw\ndata which appears of high quality and the assembly looks adequate (considering a <10x read coverage). The report is\nlargely descriptive in line with the nature of the study, which may be enhanced with some additional core genome analyses\nto improve temporal evolutionary estimates. This genome provides a framework for more informed understanding of the\narchitecture and evolutionary framework of the ‘modern’ S. pyogenes population. Comments:\n1. While reporting virulence genes, is the putative location of these genes consistent with a modern S. pyogenes genome?\nEg. Is the location of the emm gene within an archetypical Mga operon? And located in the same core genome location it is\nin modern genomes?\n2. One element lacking is a report of the classic regulatory elements. Does the ancient genome have a classic complete\ncovR/S regulator? The incorporation of other known regulatory genes would be of benefit.\n3. The annotated map of the pyogenes genome in Figure 2a lacks some information. Perhaps a mauve alignment with a\nmodern reference would help to contextualise the overall architecture. Why a question mark in the legend when reporting\n“biins’? The shaded heat map of phage mash distances is uninformative.\n4. Where are the two intact prophages located? Please use the McShan et al nomenclature for this. When reporting the\nprophage elements in the text, please add which virulence genes are carried by which phage (and also in figure 3). Can the\nauthors confirm whether the ssa gene (mentioned on page 8) was associated with a prophage (in line with modern\ngenomes)?\n5. S. pyogenes prophage typically are ordered by the integrase gene. In figure 3, it appears as though the assembly is\ndisordered and the genetic makeup of the two prophase are very similar across several phage modules. This figure appears\nto be a pairwise alignment, yet this is not clear with the shading applied for the coloured phage modules. Suggest that the\nauthors revise this figure to make the pairwise alignment more clear. The authors infer in their results that these assignments\nmay be ambiguous. Can the assembly graphs be provided as supplementary to support their inferences. This confusion is\nreiterating in the final results paragraph, where the authors suggest degradation in these prophage and as such, may not\nrepresent “represent complete genomes with all necessary machinery for excision and replication“ as inferred in the\npreceeding paragraphs.\n6. Could the authors provide some comparison on the overall gene functional categories in the ancient genome relative to\nmodern genomes? Ie are there any genes in the ancient genome that is NOT found in modern genomes (pan-genome\nanalyses)? Was the number and type of gene function categories the same? Eg metabolic genes. eggNOG or equivalent\nwould enable this functional categorisation of the strain relative to modern references. A reader may find this comparison\ninteresting considering that diet changes that have occurred in the human host over time (or as the authors suggest may in\nfact be multi-host).\n7. Generally in ancient genome studies, the presence of an ancient phylogenetic root facilitate temporal style phylogenetic\nanalysis. Perhaps the authors would consider Bayesian analysis to improve estimates regarding ancestral evolution of S.\npyogenes (core genome)? Eg. Within the “HC10 clusters” (unclear to this reviewer what this cluster designation refers to)\nsection where “evolutionary placement” phylogenetic inferences were undertaken.\n8. In the discussion, the authors extrapolate their findings to infer invasive potential “These findings also suggest that it\nwould belong to an invasive group A Streptococcus (GAS) lineage capable of causing severe infections in ancient human\npopulations.“ This is a loose inference, and should be tempered.\nMinor comments:\nLowercase emm when referring to the gene (even at the start of a sentence). Eg. Emm137 should be emm137\nPlease italicise “Streptococcus” throughout.\nPlease double check references. There appears to be some duplication (eg 39 and 5).", "reviewer_3": "This is an extremely interesting data set. My primary criticism (points #1-#5) is that the sequence data has not been analyzed to\nsufficient depth. There are some key questions that remain unanswered - however they could be readily addressed, and\nsuch an analysis is likely to elevate the impact and value of the study. A second concern lies in speculations on data\ninterpretation (points #6-#9).\nOverall, I find some of the summations of the S. pyogenes and related literature to be simplified and/or out of context (and\nprobably unnecessary). The manuscript could likely be strengthened by a more focused approach that includes a more\nthorough evolutionary analysis of this important genome sequence. Specific comments are listed below:\n1. One key piece of missing analysis is an in-depth comparison of the Bolivian isolate genome to Streptococcus equisimilis\nsubspecies dysgalactiae (SDSE) – SDSE are the closest known genetic relatives of S. pyogenes. SDSE and S. pyogenes\ndisplay an extensive history of horizontal gene exchange and share a very large number of virulence genes (emm, FCT-9\npili, streptolysins, streptokinase, many others), core genes (including MLST overlap) and mobile genetic elements\n(phage genes, R-genes). For example, please see Xie et al (Nature Comm, 2024), among several other publications on S.\npyogenes-SDSE gene flow. A deeper analysis of the Bolivian isolate versus closely related groups C and G strep genomes\nseems warranted. Strep. agalactiae was chosen for comparison as a related species, however, it is more genetically distant.\n2. A second key piece of missing analysis: Most S. pyogenes strains have an emm gene plus 1 or 2 emm-like (emmL)\ngenes; thus, one cannot necessarily assign emm-type based on a sequence hit alone. Evidence should be provided\nshowing that the so-called ‘emm137’ gene corresponds to the correct genome map position for the type-specific emm locus -\nand not the downstream (paralogous) emmL/enn locus. Most emm genes have a history of horizontal transfer and\nrecombination, and most emm types have been recovered from multiple divergent genetic backgrounds. The (rare) emm137\nFiji strain cited in reference 29 has a distinct FCT-region (not FCT-9) and therefore, does not appear share a recent/direct\ncommon ancestor with the Bolivian strain.\n3. Related to #2: An in-depth and granular Blast analysis that includes the entirety of the emm-region – spanning from the\nMga locus to the C5a peptidase locus (whether assembled or in fragments) – might address the emm/emmL ‘137’ issue, as\nwell as provide key evolutionary insights since the emm-region of the genome is subject to strong diversifying selection. It\nmay be insightful to know which sub-portions of the emm-region and emm/emmL gene product functional domains display\nhigh homology to contemporary S. pyogenes genes, and which parts are unique and non-extant.\n4. A map of the FCT region genes might also be a nice addition. It is stated that the Bolivian isolate displays the FCT-9\nregion form. Is the pilB gene highly homologous to any contemporary alleles? How about other FCT-9 region genes?\nSeveral FCT-9 loci of S. pyogenes are shared with group C and G streptococci.\n5. One of the benefits of a dated ancient genome is that the data can be used to refine the ‘molecular clock’ and mutation\nrates. Although often such analyses yield soft conclusions, a thoughtful attempt here might be worthwhile. A deeper analysis\nof select core genes, for their sequence divergence with modern-day strains, may also provide interesting new evolutionary\ninsights. Additional data supporting the novelty of the ancient isolate can help to minimize any lingering doubts that it is a\nmodern-day contaminant.\n6. Descriptive annotations derived from bioinformatics software do not by themselves equate to biological function. There\nare several instances in the manuscript where the authors delve into commentaries on annotated genes as virulence- or\nresistance-associated - yet those ascribed gene functions in this species (S. pyogenes) are not necessarily experimentally\nproven phenotypes. Unless biological functions are supported by (accurate) citations, the discussions seem way too\nspeculative; this concern extends to the Abstract.\n7. What is meant by “invasive” lineage? Most strains/emm-types of S. pyogenes can (more or less) cause invasive disease.\n8. The commentary on a systemic infection, with a possible mixed infection, seems weak and not well-connected to current\nclinical knowledge. The ‘co-isolates’ Clostridium tetani and botulinum (and Morganella) are primarily environmental (soil;\nand the clostridia form endospores, perhaps preserving their DNA); it is extremely rare for these species to be present in the\nbloodstream (unlike S. pyogenes). Perhaps the person had sepsis due to S. pyogenes or perhaps they had a\npolymicrobial tooth abscess (although a mixed infection would more likely include anaerobes from the oral microbiota), or maybe the\ncause of death was non-infectious – one can only speculate.\n9. The discussion on the population structure of S. pyogenes seems a bit out of place and unclear.\n10. There are a few reports of natural recovery of S. pyogenes from non-human primates, although there is no strong\nevidence that this bacterium is a pathogen of consequence (thereby ‘human-specific’). The authors may wish to consider\ntheir discussion within this context.\n11. The references lack journal titles. Some citations do not provide the evidentiary support that the authors claim.\n12. Line numbers would be helpful.", "reviewer_4": "In this manuscript, Valverde and colleagues present the ancient DNA analysis of a tooth sample from a pre-Columbian era\nmummy from Bolivia. The highlight of the study is the recovery and characterization of a near-complete genome of a\npathogen, Streptococcus pyogenes.\nThe paper is generally well-written. The analyses related to the phylogeny and genomic architecture of the S. pyogenes\nstrain are well done and interesting to read. However, the authors need to do a better job of convincing the readers regarding\nthe ancient origin of the pathogens presented in this paper. The following comments need to be addressed before I can\nrecommend the manuscript for publication:\nThe authors state “Further microbial profiling using MetaPhlAn4 displayed unusual high abundances of different human\npathogenic bacterial species - Clostridium tetani (37.5%), Morganella morganii (29.5%), Streptococcus pyogenes (11.6%),\nClostridium botulinum (6.7%), Asaccharospora irregularis (6.4%), and Paeniclostridium sordellii“. Typically when one is\nanalyzing ancient pathogen DNA from a tooth sample, the pathogen’s DNA is very low in abundance as compared to that\nfrom oral microbes. The profile given here strikes me as very unusual. Given the high abundance of these six taxa in the\nsample, one might hypothesize that these are contaminants that have leached into the tooth from the burial environment.\nThis possibility needs to be discussed and the authors need to make a better case for why they think these are ancient taxa,\nendogenous to the tooth sample, and why their presence in this sample indicates that the individual suffered from an\ninfection by these pathogens. For e.g. Clostridium tetani and C. botulinum can both be found in the soil. Asaccharospora\nirregularis was formerly considered to be a Clostridium, so it is likely found in the soil too, but I can find no reference of it\nbeing a human pathogen.\nFrom Table 1, it seems the S. pyogenes DNA shows 1.2% and 1.7% damage at the 5’ and 3’ ends, respectively, whereas\nthe human DNA from the same tooth shows 3.7% and 5.1% damage. The authors need to discuss this discrepancy between\nthe amount of damage for the host and pathogen DNA.\nThe argument for mixed infection seems weak to me, especially since I find it less likely that the Clostridium species are\nendogenous to the tooth sample. In any case, the authors will need to give damage patterns for the other pathogens in order\nto support an argument for mixed infection.\nIt is disappointing that there is no discussion of the burial context for this individual. Furthermore, was any paleopathological\nanalysis done that suggests the presence of any particular disease(s)? Also, please specify which tooth was sampled.\nThe DNA extraction and sequencing section need more details. Please briefly describe the extraction and library\npreparation protocol used."}, "author_response": "General response to the Editor and Reviewers\nWe sincerely thank the Editor and all four reviewers for their careful evaluation of our manuscript\nand for the insightful, constructive comments. We greatly appreciate the time and expertise invested\nin these reviews. The feedback has been extremely valuable and has helped us to substantially\nimprove the clarity, rigor, and scope of the study.\nIn response to the reviewers’ main concerns, we have extensively revised and restructured the\nmanuscript. In particular, we have refined our interpretations to avoid overstatement, expanded and\nstrengthened the evolutionary analyses, and focused the functional discussion on the most robust\nand well-supported genomic features of Streptococcus pyogenes.\nThe major changes implemented in the revised manuscript include the following:\n1) We have toned down the mixed-infection hypothesis, limiting our interpretation to reporting\nthe presence of additional microbial taxa detected in the metagenomic data, and clarifying\nthe uncertainties regarding their origin and biological relevance.\n2) We have redone the phylogenetic analyses, incorporating the closely related species\nStreptococcus dysgalactiae and Streptococcus canis, and applied molecular clock\napproaches to estimate divergence times and better contextualize the evolutionary\nplacement of the ancient S. pyogenes genome.\n3) We have restricted the virulence analysis to the most prominent and well-documented\nvirulence determinants of S. pyogenes, supported by established literature, and removed or\ndownplayed more speculative functional interpretations.\n4) We have toned down the prophage analyses, including simplifying the discussion of\nprophage similarity and reducing emphasis on MASH distance-based comparisons, while\nclarifying uncertainties related to assembly completeness and functional inference.\nIn addition, we have addressed numerous specific comments regarding genome architecture,\nregulatory elements, comparative genomics, figure presentation, terminology, formatting, and clarity\nof methods and discussion. All reviewer comments are addressed point by point, with references to\nspecific line numbers where changes have been made in the revised manuscript. Instead of\nhighlighting the changes with colours in the main manuscript we decided to copy the text changes\ninto the respective reviewer response sections. Our responses are in blue font.\nWe hope that these substantial revisions adequately address the reviewers’ concerns and\nstrengthen the manuscript. We are grateful for the opportunity to revise our work and look forward to\nyour further evaluation."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for their careful and thorough revision. My concerns have been adequately addressed, and the\nmanuscript has been substantially strengthened. I have no further comments and support publication in its current form.", "reviewer_2": "The authors have made substantial alternations to the manuscript including new analyses which address the major\nconcerns/queries. The revised version has substantially improved and the authors should be commended for their diligence.\nThe supplementary material is important for validation.\nWhile the raw data looks to have been deposited, it is strongly recommended that draft assemblies generated be made\npublic available.", "reviewer_3": "The extensive additional analyses and thoughtful responses to reviewer concerns are greatly appreciated. The revised\nmanuscript seems much better focused and importantly, seems well-poised to have a high impact on the S. pyogenes field.\nI have a few suggestions for further clarifications:\n1. Lines 79-80. “Like S. pyogenes, SDSE commonly expresses Lancefield group C or G antigens…..” Incorrect as written: S.\npyogenes express Lancefield group A carbohydrate.\nAlso, if the authors know that their isolate has high homology to the group A carbohydrate biosynthesis genes (gac genes), it\nmay be of added value to include that point in their analysis of other key genes (Results section)\n2. Line 83: Omit the “Thereby” – Exchange of many [virulence] genes is by mechanisms unknown (i.e., not yet proven for\nmobile genetic elements as the vehicle, except for [largely] pyrogenic exotoxins via transduction]). Also lines 85-87.\n3. Line 253 (and related Discussion): animal infection models for pneumococcus (not S. pyo)? Ref. 35\n4. Line 276-paragraph (and related Discussion): Is the lmrP gene product actually known to confer antibiotic resistance in\nthis species (S. pyo)? Citation? Also, a functional ‘mef[E]’ is far more typical for other streptococcal species. The ideas in this\nsection seem a bit overstated.\n5. Lines 396-399: It is not all that uncommon for modern strains to lack phage-associated virulence (superantigen) genes. So\nperhaps this idea is not quite correct as presented.\n6. Line 400: “S. pyogenes is predominantly a clonal pathogen” (ref. 43 cited). Authors should probably clarify or soften.\nGlobal strains of S. pyogenes taken as a group are nearly as ‘recombinogenic’ as pneumococci via numerous measures.\nThere exist some strongly ‘clonal’ lineages within S. pyogenes (e.g., M1T1) and many lineages/clusters (phylogroups) that\nare more diverse.", "reviewer_4": "The manuscript by Valverde and colleagues is much improved following the previous round of revision. I am satisfied with\ntheir response to the concerns I had previously raised. Below I have noted some additional (mostly minor) concerns that\nneed to be addressed.\nLines 319-324: The authors infer that the timing of the origin of S. pyogenes is the same as the timing of its split from S.\ndysgalactiae i.e. 137,264 years. I don’t quite agree with this inference; in my opinion, the timing of the MRCA of all S.\npyogenes strains (ancient and modern), which is estimated to be 10,270 years, is what one should consider as the origin of\nthis species. If the authors disagree with this, I encourage them to include a sentence stating why they think using the stem\nage (time at which S. pyogenes diverged from its sister taxon), as opposed to the crown age (time of MRCA for all S.\npyogenes), would be appropriate to consider as the time at which this species originated.\nLine 532: Please specify which bwa algorithm was used – mem or aln? Were default parameters used? If using aln, were\nrecommended aDNA parameters used (such as disabling the seed, changing mapping stringency to allow reads with\ndamaged nucleotides to map)? If yes, please specify the parameter values.\nLine 533: With regards to Qualimap, do you mean to say coverage instead of mapping quality?\nLine 536: It is unclear how the consensus sequence for the human DNA was called, since this cannot be done using either\nMapDamage2 or Haplogrep. Did you mean to say you rescaled the damaged sites using MapDamage2 and then called the\nconsensus using schmutzi?\nLine 631: Correct \"EES\" to \"ESS\". Also, specify what was considered as a “satisfactory” ESS value.\nApart from the above minor concerns, the manuscript needs to be thoroughly proof-read to correct instances of awkward\nphrasing and lingering typos. The instances I could find are given below.\nLine 30: Missing period at end of sentence. Add “the” before “European colonial expansion”.\nLine 37: Use “the” in place of “an”.\nLine 56: Add comma after “health”.\nLines 96, 124: Add “the” before “pre-Columbian”.\nLine 104: Missing period after the “S” of Streptococcus.\nLine 105: Rephrase to “… and calibrate the pathogen’s phylogeny…”\nLine 109: Rephrase to “…tooth sample from a pre-Columbian…”\nLine 129: Rephrase to “…environments with dry conditions which have been previously shown to be less prone to DNA\ndegradation.”\nLine 130: Remove space between MetaPhlAn and 4.\nLine 138: Rephrase to “…to this species to understand…”\nLine 140-142: Rephrase to “…which revealed the presence of S. pyogenes in human samples from different time periods\n(ranging from 4000 to 200 years BP) in Europe and Africa, as well as in gorilla samples from museum collections”\nLine 144: I recommend using the word “samples” instead of the word “findings”.\nLine 164: Rephrase to “The completeness of the assembled S. pyogenes genome was 99.98%..”\nTable 1: Needs to be properly formatted with respect to gray and white cell colors.\nLine 179: Redundant “typing” after MLST.\nLine 195: Remove the word “situation”. I recommend using “currently available” instead of “current available”.\nLine 197: Rephrase to “…pyogenes strain to the modern emm type 46 (79% nucleotide identity)”\nLine 234: Rephrase to “In addition to the surface-associated M- and T-proteins, the ancient S. pyogenes strain holds…”\nLine 252: Change to “…pathogenicity in animal infection…”\nLine 255: Specify here what you mean by “our prophages”. I recommend saying “however, the prophages in the ancient\nstrain…”.\nLine 290: Rephrase to “to be completely absent in other modern GAS genomes as well”.\nLine 297, 311, 312, 314, and elsewhere: Do not capitalize the word “gorilla”.\nLine 317: Remove the word “ancestral”. It is redundant in that sentence.\nLine 322 – 324: Rephrase to “The split of the Bolivian strain from all other S. pyogenes (divergence of known S. pyogenes)\nwas estimated to be 10,270 years ago (95% HPD: 9314 to 11,218 cal yr BP), and a diversification of most modern S.\npyogenes in the past 5,500 years of human history.”\nLine 370: Missing period after the word “infections”.\nLine 476: I recommend saying tuberculosis, instead of TB, as it is the first (and only) time the disease is mentioned in the\nmanuscript.\nLine 503: Remove comma after 2730.\nLine 534 and 536: Remove space between MapDamage and 2.\nLine 559: Give the full form of GTDBtk.\nLine 579: This line is missing a verb after “we” for e.g. downloaded or acquired.\nLines 622-625: Please ensure that numbers that need to be shown as superscript are correctly formatted. For eg. 1.8 10-6\nneeds to be corrected to 1.8x10^-7 or given as superscript."}, "author_response": "Point-by-point response to the reviewers' comments\nSecond round of review\nOur responses in blue\nREVIEWERS' COMMENTS\nReviewer #1 (Remarks to the Author):\nI thank the authors for their careful and thorough revision. My concerns have been\nadequately addressed, and the manuscript has been substantially strengthened. I have\nno further comments and support publication in its current form.\nWe thank the reviewer for this positive assessment\nReviewer #2 (Remarks to the Author):\nThe authors have made substantial alternations to the manuscript including new\nanalyses which address the major concerns/queries. The revised version has\nsubstantially improved and the authors should be commended for their diligence. The\nsupplementary material is important for validation.\nWe thank the reviewer for this positive assessment\nWhile the raw data looks to have been deposited, it is strongly recommended that draft\nassemblies generated be made public available.\nSequencing data and the assembled genome are available at the European Nucleotide\nArchive (ENA) under ENA: PRJEB91735.PRJEB91735. The sequencing reads are\navailable at the Sequence Read Archive under accession ERR15308372\nERR15308372 and the assembly under accession number\nGCA_982145515.1.GCA_982145515.1. Source data are provided with this paper.\nReviewer #3 (Remarks to the Author):\nThe extensive additional analyses and thoughtful responses to reviewer concerns are\ngreatly appreciated. The revised manuscript seems much better focused and\nimportantly, seems well-poised to have a high impact on the S. pyogenes field.\nI have a few suggestions for further clarifications:\n1. Lines 79-80. “Like S. pyogenes, SDSE commonly expresses Lancefield group C or\nG antigens…..” Incorrect as written: S. pyogenes express Lancefield group A\ncarbohydrate.\nThanks for pointing out this incorrect phrasing. We rewrote the sentence as follows:\n“In contrast to S. pyogenes, that belongs Lancefield serogroup A, SDSE commonly\nexpresses Lancefield group C or G antigens, but still causes a similar spectrum of\nhuman diseases ranging from pharyngitis to invasive soft tissue infections.”\nAlso, if the authors know that their isolate has high homology to the group A\ncarbohydrate biosynthesis genes (gac genes), it may be of added value to include that\npoint in their analysis of other key genes (Results section)\nWe followed the reviewer’s suggestion and subjected the twelve gac operon genes\nlocated on the contig 5 to BlastP homology analysis and included the results into the\nSupplementary Data 10 table. All twelve gac genes show high homology (>99%) to S.\npyogenes gac genes.\nWe added the following sentence to the results section of the revised manuscript:\n\"The full set of genes (gacA-L) of the group A carbohydrate (GAC) gene cluster are\npresent on contig 5 with high homology, i.e., > 99% proteins identity, to modern GAS\nstrains (Supplementary Data 10). \"\n2. Line 83: Omit the “Thereby” – Exchange of many [virulence] genes is by\nmechanisms unknown (i.e., not yet proven for mobile genetic elements as the vehicle,\nexcept for [largely] pyrogenic exotoxins via transduction]). Also lines 85-87.\nWe omitted the “Thereby” and used the pyrogenic exotoxins as an example.\n3. Line 253 (and related Discussion): animal infection models for pneumococcus (not\nS. pyo)? Ref. 35\nThank you for pointing this out. We changed the part in Line 253 accordingly and\nwrote: “Despite the absence of superantigens encoded on these two prophages, one\nprophage harboured a virulence determinant, Virulence-Associated Protein E (VapE),\nthat was previously associated with virulence in a mouse model of Streptococcus\npneumoniae infection.”\nThe Discussion part has been removed.\n4. Line 276-paragraph (and related Discussion): Is the lmrP gene product actually\nknown to confer antibiotic resistance in this species (S. pyo)? Citation? Also, a\nfunctional ‘mef[E]’ is far more typical for other streptococcal species. The ideas in this\nsection seem a bit overstated.\nThank you for pointing this out. We rephrased the last sentence in line 253 accordingly:\n“Both genes were complete, with no evidence of frameshifts or premature stop codons,\nsuggesting they were intact in the ancient strain. Their functional properties, however,\nremains to be determined.”\nIn addition, we removed the following sentence in the discussion:\n“The high sequence identity and broad coverage observed indicate that these genes\nwere likely functional at the time, supporting the hypothesis that resistance traits were\nalready circulating in microbial communities before the antibiotic era”\n5. Lines 396-399: It is not all that uncommon for modern strains to lack phage-\nassociated virulence (superantigen) genes. So perhaps this idea is not quite correct as\npresented.\nWe agree and removed the following sentence:\n“Their absence in the ancient genome supports the view that such superantigens\nrepresent later-acquired, mobile genetic elements rather than ancestral core virulence\nfactors, highlighting the role of horizontal gene transfer in increasing the pathogenic\npotential of modern S. pyogenes strains.”\n6. Line 400: “S. pyogenes is predominantly a clonal pathogen” (ref. 43 cited). Authors\nshould probably clarify or soften. Global strains of S. pyogenes taken as a group are\nnearly as ‘recombinogenic’ as pneumococci via numerous measures. There exist some\nstrongly ‘clonal’ lineages within S. pyogenes (e.g., M1T1) and many lineages/clusters\n(phylogroups) that are more diverse\nThank you for pointing this out! We rephrased the sentence accordingly:\n“S. pyogenes lineages can be grouped into related clusters (‘phylogroups’) defined by\nwhole genome clustering methods. These phylogroups are not genetically static, as\noccasional recombination occurs via mobile genetic elements such as….”\nReviewer #4 (Remarks to the Author):\nThe manuscript by Valverde and colleagues is much improved following the previous\nround of revision. I am satisfied with their response to the concerns I had previously\nraised. Below I have noted some additional (mostly minor) concerns that need to be\naddressed.\nLines 319-324: The authors infer that the timing of the origin of S. pyogenes is the\nsame as the timing of its split from S. dysgalactiae i.e. 137,264 years. I don’t quite\nagree with this inference; in my opinion, the timing of the MRCA of all S. pyogenes\nstrains (ancient and modern), which is estimated to be 10,270 years, is what one\nshould consider as the origin of this species. If the authors disagree with this, I\nencourage them to include a sentence stating why they think using the stem age (time\nat which S. pyogenes diverged from its sister taxon), as opposed to the crown age\n(time of MRCA for all S. pyogenes), would be appropriate to consider as the time at\nwhich this species originated.\nWe agree with the reviewer that our wording was unclear and that the crown is more\ncorrect for describing the age of S.pyogenes. Technically, the origin of S.pyogenes is\nanytime between the stem and the crown, with more confidence as we get closer to the\nCrown.\nWe rephrased the part accordingly:\n“Time estimates of the full dataset under most fitting analytical parameters (see\nMethods) returned an age for the split of S. pyogenes and S. dysgalactiae spp. at\n137,264 years (95% HPD: 124,428 to 149,847 cal years BP) (Figure 3a). The timing of\nthe Most Recent Common Ancestor (MRCA) of all S. pyogenes strains (ancient\nBolivian strain and modern known S. pyogenes) is estimated to be 10,270 years (95%\nHPD: 9314 to 11,218 cal years BP).The diversification of most modern S. pyogenes\ntook place in the past 5,500 years of human history (Figure 3C).”\nLine 532: Please specify which bwa algorithm was used – mem or aln? Were default\nparameters used? If using aln, were recommended aDNA parameters used (such as\ndisabling the seed, changing mapping stringency to allow reads with damaged\nnucleotides to map)? If yes, please specify the parameter values.\nWe specified accordingly:\n“For the analysis of the human DNA, we mapped the merged reads against the\nreference human genome (build hg19), using BWA-aln v0.7.17 59 (using the following\nparameters: -l 1000 -n 0.04 -o 1 -e -1), then filtered for minimum mapping quality of 25\nusing SAMtools v1.17 60.”\nLine 533: With regards to Qualimap, do you mean to say coverage instead of mapping\nquality?\nWe rephrased accordingly:\n“The mapping quality and coverage were then checked using Qualimap v2.2.2d 63.”\nLine 536: It is unclear how the consensus sequence for the human DNA was called,\nsince this cannot be done using either MapDamage2 or Haplogrep. Did you mean to\nsay you rescaled the damaged sites using MapDamage2 and then called the\nconsensus using schmutzi?\nWe rephrased accordingly:\n“The damaged sites of aDNA reads in the BAM file were then rescaled using\nMapDamage 2 and used to call a consensus sequence using Schmutzi script\n“log2fasta”, for mitochondrial haplogroup assignment using Haplogrep3 v2.4.0 64.”\nLine 631: Correct \"EES\" to \"ESS\". Also, specify what was considered as a\n“satisfactory” ESS value.\nWe rephrased accordingly:\n“All MCMC were run for a minimum of 50 million generations until convergence.\nConvergence was inspected using Tracer by checking that all parameters had an\nESS>200 and setting a burn-in that maximised ESS.”\nApart from the above minor concerns, the manuscript needs to be thoroughly proof-\nread to correct instances of awkward phrasing and lingering typos. The instances I\ncould find are given below.\nLine 30: Missing period at end of sentence. Add “the” before “European colonial\nexpansion”.\nChecked (new Abstract)\nLine 37: Use “the” in place of “an”.\nChecked (new Abstract)\nLine 56: Add comma after “health”.\nChecked and corrected\nLines 96, 124: Add “the” before “pre-Columbian”.\nChecked and corrected\nLine 104: Missing period after the “S” of Streptococcus.\nChecked and corrected\nLine 105: Rephrase to “… and calibrate the pathogen’s phylogeny…”\nChecked and corrected\nLine 109: Rephrase to “…tooth sample from a pre-Columbian…”\nChecked and corrected\nLine 129: Rephrase to “…environments with dry conditions which have been previously\nshown to be less prone to DNA degradation.”\nChecked and corrected\nLine 130: Remove space between MetaPhlAn and 4.\nChecked and corrected\nLine 138: Rephrase to “…to this species to understand…”\nChecked and corrected\nLine 140-142: Rephrase to “…which revealed the presence of S. pyogenes in human\nsamples from different time periods (ranging from 4000 to 200 years BP) in Europe and\nAfrica, as well as in gorilla samples from museum collections”\nChecked and corrected\nLine 144: I recommend using the word “samples” instead of the word “findings”.\nChecked and corrected\nLine 164: Rephrase to “The completeness of the assembled S. pyogenes genome was\n99.98%..”\nChecked and corrected\nTable 1: Needs to be properly formatted with respect to gray and white cell colors.\nChecked and corrected\nLine 179: Redundant “typing” after MLST.\nChecked and corrected\nLine 195: Remove the word “situation”. I recommend using “currently available” instead\nof “current available”.\nChecked and corrected\nLine 197: Rephrase to “…pyogenes strain to the modern emm type 46 (79% nucleotide\nidentity)”\nChecked and corrected\nLine 234: Rephrase to “In addition to the surface-associated M- and T-proteins, the\nancient S. pyogenes strain holds…”\nChecked and corrected\nLine 252: Change to “…pathogenicity in animal infection…”\nChecked and corrected\nLine 255: Specify here what you mean by “our prophages”. I recommend saying\n“however, the prophages in the ancient strain…”.\nChecked and corrected\nLine 290: Rephrase to “to be completely absent in other modern GAS genomes as\nwell”.\nChecked and corrected\nLine 297, 311, 312, 314, and elsewhere: Do not capitalize the word “gorilla”.\nChecked and corrected\nLine 317: Remove the word “ancestral”. It is redundant in that sentence.\nChecked and corrected\nLine 322 – 324: Rephrase to “The split of the Bolivian strain from all other S. pyogenes\n(divergence of known S. pyogenes) was estimated to be 10,270 years ago (95% HPD:\n9314 to 11,218 cal yr BP), and a diversification of most modern S. pyogenes in the past\n5,500 years of human history.”\nWe rephrased to:\n“The timing of the Most Recent Common Ancestor (MRCA) of all S. pyogenes strains\n(ancient Bolivian strain and modern known S. pyogenes) is estimated to be 10,270\nyears (95% HPD: 9314 to 11,218 cal years BP). The diversification of most modern S.\npyogenes took place in the past 5,500 years of human history (Figure 3C).”\nLine 370: Missing period after the word “infections”.\nChecked and corrected\nLine 476: I recommend saying tuberculosis, instead of TB, as it is the first (and only)\ntime the disease is mentioned in the manuscript.\nChecked and corrected\nLine 503: Remove comma after 2730.\nChecked and corrected\nLine 534 and 536: Remove space between MapDamage and 2.\nChecked and corrected\nLine 559: Give the full form of GTDBtk.\nChecked and corrected\nLine 579: This line is missing a verb after “we” for e.g. downloaded or acquired.\nCorrected to: “we downloaded...”\nLines 622-625: Please ensure that numbers that need to be shown as superscript are\ncorrectly formatted. For eg. 1.8 10-6 needs to be corrected to 1.8x10^-7 or given as\nsuperscript.\nChecked and corrected"}]

Decision letter