As infecções treponêmicas, causadas pela bactéria T. pallidum , estão aumentando a taxas alarmantes em todo o mundo 7,8,9,10,11. Increasing evidence suggests that many treponemal strains have developed antibiotic resistance, which is expected to facilitate their spread12. This re-emerging threat has led to many modern genetic and medical studies8,13,14,15. intimamente relacionadas de T. pallidum As subespécies T. pallidum pallidum (TPA), T. pallidum pertenue (TPE) e T. pallidum endemicum (TEN) - responsáveis ​​​​pela sífilis, bouba e bejel, respectivamente - têm sequências genômicas altamente semelhantes que diferem em aproximadamente 0,03% 16 , 17 . 

Hoje, o bejel está geograficamente concentrado em ambientes áridos e quentes, especialmente no Mediterrâneo oriental e no oeste da Ásia, enquanto a bouba é encontrada principalmente nos trópicos úmidos e quentes, como a África ou a América do Sul. 18 . Entre as treponematoses, a sífilis é a mais distribuída globalmente; é comum mesmo em populações ocidentais ricas com fácil acesso a cuidados de saúde 7 , 13 . Por outro lado, a bouba e o bejel afetam principalmente os países em desenvolvimento e permanecem menos extensivamente estudados 18 .

Historicamente, a sífilis venérea é conhecida por ter causado um surto devastador na Europa no final do século XV. Os sintomas que podem desenvolver-se na ausência de tratamento eficaz incluem desfiguração física grave, cegueira e deficiência mental. 19 . Como manifestações semelhantes podem aparecer em todas as treponematoses 20 , 21 , 22 , a sua distinção ao nível da subespécie é muitas vezes pouco fiável e baseia-se principalmente na localização de úlceras cutâneas características (nos órgãos genitais ou noutros locais), especialmente em países em desenvolvimento com recursos médicos limitados 23 , 24 , 25 . Os diagnósticos de casos históricos são igualmente difíceis: embora as treponematoses possam deixar alterações patológicas nos ossos, estas aparecem em apenas aproximadamente 5–30% dos casos avançados. 26 , 27 , resultando em provável subestimação da prevalência passada de treponematoses.

The early presence and potential origin of syphilis in Europe was proposed in the pre-Columbian hypothesis, based on osteological analyses of treponemal lesions, whereas the Columbian hypothesis associates its emergence with Columbus’ first American expedition and considers the contradicting palaeopathological evidence to be unreliable2. Before the distinctions among the subspecies could be genetically defined28,29, the unitarian hypothesis claimed that all treponematoses were the same disease, which only manifested differently under different environmental and cultural factors. Attempts to identify subspecies using palaeopathology have had ambiguous results and require DNA evidence as confirmation: previous ancient DNA studies have, for example, revealed that some cases of presumed syphilis instead correspond to yaws4,5, and recovered at least one previously unknown T. pallidum strain4. Since treponemes possess an impressive ability to adjust to various environments and are known to have previously occupied geographical regions outside their present distributions4,5, apenas evidências inequívocas de DNA treponêmico pré-colombiano podem iluminar as origens da sífilis e também desvendar aspectos importantes da história evolutiva de todos os treponemas.

Aqui, apresentamos evidências de doença treponêmica pré-colombiana no Novo Mundo a partir de um cemitério de sambaqui brasileiro com quase 2.000 anos de idade, Jabuticabeira II, através de quatro genomas reconstruídos de T. pallidum com cobertura de até 33,6x, filogeneticamente basal ao diversidade moderna da subespécie causadora de bejel, T. pallidum endemicum .

Origens geográficas e paleopatologia

Noventa e nove exemplares de Jabuticabeira II da região de Laguna de Santa Catarina, na costa brasileira, com e sem patologias, foram incorporados neste estudo. Análises osteológicas anteriores revelaram patologias relacionadas com infecções, sugerindo potenciais infecções treponémicas 30 , como periostite (24 casos), remodelação óssea (4 casos) e marcas de traça no crânio (4 casos). Das 37 amostras consideradas preliminarmente positivas para DNA treponêmico após a triagem inicial, 12 eram de indivíduos com patologias e o restante era proveniente de espécimes não patológicos (Tabela Suplementar 1 ). Quatro amostras ósseas, de quatro indivíduos diferentes, produziram dados genômicos suficientes para análises abrangentes. A amostra ZH1390 (Tabela 1 e Fig. 1a ) representa um fragmento de tíbia apresentando periostite. A amostra ZH1540 veio de um conjunto de ossos misturados de um esqueleto incompleto, nomeadamente de uma fíbula com lesões patológicas (Tabela 1 e Fig. 1a ). As amostras ZH1541 e ZH1557 originaram-se de ossos longos sem nenhuma patologia identificada (Tabela 1 e Figura 1a ). Todas as amostras foram datadas por radiocarbono e testadas quanto ao efeito de reservatório marinho. Os dados brutos, calibrados e corrigidos de 14 A datação C é apresentada em Dados Estendidos Fig. 1a, b , Tabela Suplementar 2 e Métodos , 'Informações Arqueológicas'. Os indivíduos positivos para DNA treponêmico não foram enterrados separadamente dos demais indivíduos em Jabuticabeira II, sugerindo que foram tratados igualmente.

Table 1 Summary of the sample information and central statistics
Full size table
Fig. 1: The archaeological site and the T. pallidum-positive samples that yielded the reconstructed genomes.
figura 1

a , Mapa mostrando a localização do sítio arqueológico Jabuticabeira II, no litoral sul do estado de Santa Catarina, Brasil, e as amostras ZH1390, ZH1540, ZH1541 e ZH1557, cujos genomas foram reconstruídos. b , Uma árvore filogenética de máxima verossimilhança das cepas modernas e antigas de T. pallidum usando GTR + G + I (ver Métodos ) como modelo evolutivo e 1.000 repetições de bootstrap. Todos os genomas antigos utilizados neste estudo (genomas antigos recentemente reconstruídos e publicados anteriormente; ver Tabela Suplementar 3 ) estão marcados em negrito. Os pontos rosa representam nós com valores de bootstrap superiores a 70%. A barra de escala está em unidades por substituições por site.

Imagem em tamanho real

Preliminary pathogen screening

Na triagem inicial a partir de dados de sequenciamento shotgun, 37 de 99 amostras mostraram entre 7 e 133 ocorrências para táxons da família Treponema no banco de dados Kraken e foram incluídas no processo de enriquecimento de alvo (Tabela Suplementar 1 e Métodos , 'Processamento de amostra') . Destas amostras, 9 tiveram mais de 5.000 leituras mapeadas para 3 genomas de referência de T. pallidum (BosniaA, CDC2 e Nichols) pós-captura e foram, portanto, consideradas positivas para infecção treponêmica (Tabela Suplementar 1 ). Para estas amostras positivas, três bibliotecas adicionais de cadeia dupla foram produzidas para uma segunda rodada de enriquecimento de todo o genoma 4 , 31 (para uma metodologia detalhada, consulte Métodos , 'Processamento de amostras'). Após o enriquecimento adicional, os dados de todas as bibliotecas feitas a partir dos mesmos extratos originais e todos os dados de sequenciamento produzidos nas duas rodadas de captura foram combinados para cada amostra. Quatro amostras, ZH1390, ZH1540, ZH1541 e ZH1557, tiveram leituras cobrindo 9,2–99,4% do genoma de referência da BósniaA em 1×, com uma cobertura média entre 2× e 33,6× (Tabela 1 ). Essas quatro amostras foram consideradas como tendo maior potencial para reconstrução do genoma completo e análises posteriores.

Estimativa de autenticidade de DNA antigo

A autenticidade do DNA antigo foi confirmada examinando a desaminação de bases nas extremidades das leituras: 21%, 10%, 12% e 7% nas extremidades 5' e 17%, 12%, 14% e 6% nas extremidades 3. ′ termina para as amostras ZH1390, ZH1540, ZH1541 e ZH1557, respectivamente (Extended Data Fig. 1c – f ). As amostras tinham comprimentos médios de fragmentos 32 , 33 , 34 variando de 64 pb a 74 pb (Tabela 1 ). Além disso, a atribuição dos cromossomos sexuais dos dados do shotgun foi consistente com XX para as amostras ZH1540 e ZH1541. Embora os indivíduos tenham sido anteriormente considerados prováveis ​​​​do sexo masculino nas análises osteológicas, as amostras ZH1390 e ZH1557 produziram dados insuficientes para a determinação molecular do sexo (Tabela 1 e Métodos , 'Informações arqueológicas').

Reconstrução do genoma

After high-throughput Illumina sequencing of the enriched DNA from the 4 selected samples, the resulting 20–100 million raw reads were merged sample-wise and duplicate reads were removed (Table 1). Genomes were reconstructed by mapping each sample to three representative high-quality reference genomes of T. pallidum subspecies: CDC2 for TPE, BosniaA for TEN, and Nichols for TPA (Methods, ‘Sample processing’ and ‘Dataset selection’). We filtered positions on the basis of read coverage, variant allele frequency, P value and base quality, and obtained three different consensus sequences for each sample, each with a different number of covered bases, as well as SNPs. The number of SNPs in each sequence, along with the phylogenetic analyses consistently supported a placement of all four samples within the TEN clade (Figs. 1b and 2b, Table 1 and Supplementary Table 3). Although the consensus sequences from three samples, ZH1390, ZH1541 and ZH1557, were assigned to T. pallidum endemicum (Fig. 1b), read coverage was below the threshold required for downstream analyses (for details, see Methods ‘Sample processing’ and ‘Read processing and multiple reference-based genome alignment generation’).

Figura 2: Análise do genoma de alta cobertura ZH1540.
Figura 2

a , Gráfico circular do genoma ZH1540. Os círculos indicam (de dentro para fora): posição genômica, conteúdo de GC (preto) e cobertura (azul). A borda externa (cinza) mostra um conjunto de 60 genes candidatos associados à virulência e descritos em estudos anteriores 3,4. b , Árvore filogenética de máxima verossimilhança recolhida obtida usando o alinhamento do genoma completo e usando Treponema paraluiscuniculi para enraizar a árvore usando GTR + G + I (ver Métodos ) como modelo evolutivo e 1.000 repetições de bootstrap. O novo genoma antigo ZH1540 está destacado em negrito e com uma estrela. As sub-linhagens dentro de cada subespécie ou linhagem são recolhidas, com exceção dos genomas antigos e do clado TEN. Os pontos rosa representam nós com valores de bootstrap superiores a 70%.

Full size image

A sequência final obtida para a amostra ZH1540 resultou em 99,38% de cobertura em relação ao genoma de referência TEN (BósniaA), uma profundidade de cobertura mínima de 5× e uma profundidade mediana de 33,6× (Tabela 1 e Fig. 2a ). A chamada de variante resultou na identificação de 123 SNPs, cada um dos quais foi verificado individualmente (detalhes fornecidos na Tabela Suplementar 3 e em Métodos , 'Processamento de amostra' e 'Processamento de leitura e geração de alinhamento de genoma baseado em múltiplas referências'). Das referências modernas disponíveis, o novo genoma antigo da TEN exibe uma diferença de 123 SNPs em comparação com as amostras da BósniaA e do IraqueB. No entanto, o número de posições diferentes é muito maior em comparação com os 4 genomas RTE japoneses (205 SNPs) e os genomas RTE cubanos (504 SNPs).

Alinhamento do genoma baseado em múltiplas referências

The new ancient genome ZH1540 was analysed together with an additional 98 publicly available genomes, including 8 modern TEN strains, 30 TPE strains (including 9 genomes from primates and 4 ancient genomes), 30 Nichols-lineage and 30 SS14-lineage TPA strains (including 4 ancient genomes) (Supplementary Table 3). Assembly files for 33 of these 98 genomes were available and downloaded directly from the public databases European Nucleotide Archive (ENA) and National Center for Biotechnology Information (NCBI). For the remaining 65 genomes, we mapped the raw sequencing data to the closest of four representative reference genomes (CDC2, BosniaA, Nichols and SS14), to obtain new assembly files. The genome reference selected for each sample was based on the subspecies and/or lineage classification of each sample from the original publications (Supplementary Table 3). A multiple reference-based genome alignment of 98 sequences from several sources was generated according to the previously published methodology35. O alinhamento resultante abrangeu um total de 1.141.812 nucleotídeos com 6.149 SNPs detectados (consulte Disponibilidade de dados e métodos , 'Processamento de amostras' e 'Processamento de leitura e geração de alinhamento de genoma baseado em referências múltiplas').

Phylogenetic and recombination analyses

Uma reconstrução filogenética confiável exigiu a remoção de regiões genômicas herdadas não verticalmente, como regiões recombinantes ou loci com conversão intra ou intergênica. Em uma análise de recombinação com o método de incongruência filogenética 36 (PIM), we detected 34 recombinant regions across 27 genes, encompassing a total of 957 SNPs (15.56% of the total SNPs) (Supplementary Table 3 and Supplementary Table 4). Owing to the exclusion of the highly passaged Seattle-81 strain, 3 of the previously detected recombinant genes were not detected here, and 11 detected genes were novel in relation to the previously published results. The average length of the recombinant regions was 368 bp, with a minimum length of 4 bp and a maximum of 2,209 bp. Notably, all the recombination events detected here correspond to inter-subspecies transfers with the exception of an intra-subspecies recombination event found in the tp0117 gene and three additional genes for which the putative donors are unidentified external sources (Supplementary Table 4 and Methods, ‘Recombination analysis using PIM’ and ‘PIM procedure for likelihood mapping and topology tests’).

Para construir um alinhamento de herança estritamente vertical, removemos os 27 genes recombinantes detectados aqui junto com três genes, tp0316 , tp0317 e tp0897 , que são conhecidos por serem hipervariáveis ​​e/ou sujeitos à conversão gênica 37 , 38 , desde o alinhamento inicial (ver Disponibilidade de dados e métodos , 'Análise filogenética'). O alinhamento final livre de recombinação abrangeu 1.103.436 pb com 3.718 SNPs. Árvores de máxima verossimilhança foram construídas usando alinhamentos múltiplos de genoma (Fig. 2b e Dados Estendidos Figs. 2 e 3 ). Na Fig. 2 de dados estendidos , as topologias das duas árvores de máxima verossimilhança com e sem os loci recombinantes ou hipervariáveis ​​são comparadas.

A eliminação de genes herdados não verticalmente teve um efeito menor na reconstrução da filogenia de T. pallidum (Extended Data Figs. 2b e 3 e Methods , 'Phylogenetic Analysis'). Os resultados após a remoção dos sítios recombinantes detectados no PIM foram confirmados com dois outros programas de detecção de recombinação, Gubbins e ClonalFrameML ( Methods , 'Recombination analysis using Gubbins and ClonalFrameML'). Além disso, foram avaliadas mutações genéticas (A2058G e A2059G) relacionadas à resistência a antibióticos macrólidos 12 , e foram encontrados ausentes no genoma antigo ZH1540 (Methods, 'Caracterização exploratória dos genes 16S-23S').

Molecular clock dating

A datação do relógio molecular foi realizada no mesmo conjunto de dados acima, com 27 genes recombinantes, tp0316 , tp0317 e tp0897 removidos. Na filogenia calibrada no tempo estimada, todas as três subespécies (TEN, TPE e TPA), bem como as linhagens SS14 e Nichols de TPA receberam alto suporte para formar clados distintos (probabilidade posterior> 0,97; Fig. 3a e Tabela de Dados Estendidos 1 ). Tal como na filogenia de máxima verossimilhança, o novo genoma antigo, ZH1540, ocupa uma posição basal dentro do clado TEN, com todas as estirpes modernas de TEN formando um subgrupo monofilético (probabilidade posterior 0,96; Fig. 3a ). A maioria das cepas SS14 se enquadra no subclado SS14-Ω previamente definido 4 , 16 , que também recebe alto suporte posterior. De acordo com os resultados das análises de regressão raiz-a-ponta (Extended Data Fig. 4 ), um modelo de relógio relaxado não correlacionado log-normalmente distribuído (UCLD) e um modelo de relógio relaxado não correlacionado distribuído exponencialmente (UCED) foram escolhidos para a análise de datação do relógio molecular, ambos com um prior lognormal estreito com uma média (no espaço real) de 1 × 10 −7 substituições por site por ano e dp de 0,25 na taxa média de clock. Consistente com os resultados relatados anteriormente de datação por relógio molecular 4 , 35 , descobrimos que todas as cepas históricas de TPA são basais a todas as cepas modernas de SS14 e, juntas, formam um clado bem suportado (probabilidade posterior de 0,97). Portanto, consideramos que as cepas históricas se enquadram no clado SS14.

Figura 3: Resultados da datação por relógio molecular.
Figura 3

a , árvore de máxima credibilidade do clado (MCC) de genomas antigos e modernos publicados anteriormente, e o genoma antigo ZH1540 deste estudo ( n = 99). As barras azuis indicam os intervalos HPD de 95% das idades dos nós e o texto em vermelho a probabilidade posterior de um grupo ser monofilético (mostrado apenas para nós com probabilidade posterior> 0,8). b , Árvore MCC recolhida simplificada da análise com ZH1540 incluído (parte superior) e excluído (parte inferior). c , Densidades posteriores dos tempos dos tMRCAs da subespécie T. pallidum e linhagens principais estimadas pela datação do relógio molecular com ZH1540 incluída (densidades coloridas, correspondentes a b , superior) e excluídas (densidades cinza, correspondentes a b , inferior ). As linhas verticais dentro das curvas de densidade indicam os limites superior e inferior dos intervalos de 95% HPD.

Imagem em tamanho real

A idade de ZH1540, que é parametrizada pelos resultados da datação por radiocarbono, atua como uma restrição no tempo do ancestral comum mais recente (tMRCA) do clado TEN. A inclusão de uma amostra muito mais antiga resulta em tempos de divergência anteriores, com intervalos de credibilidade mais amplos para todos os clados principais da árvore (Fig. 3b e Dados Estendidos Figs. 5 – 7 ). Este efeito é mais pronunciado para TEN, onde o intervalo de densidade posterior (HPD) 95% mais alto do tMRCA se estende de 780 aC a 449 dC (236-1845 dC para o subclado compreendendo apenas cepas modernas de TEN), mas é limitado a 1077– de 1855 Anúncio ao excluir ZH1540 (Dados Estendidos Figs. 5 – 7 e Tabela de Dados Estendidos 1 ). Para todas as outras linhagens principais, o efeito é mais moderado e, embora o limite inferior do intervalo HPD de 95% possa ser várias centenas de anos mais antigo ao incluir ZH1540 (cerca de 400 anos no caso de TPE), o limite superior nunca é muito maior. mais de 50 anos (Fig. 3c , Tabela de Dados Estendidos 1 e Dados Estendidos Figs. 6 e 7 ).

Embora as estimativas medianas dos tempos de divergência de linhagem sejam mais antigas do que aquelas relatadas anteriormente (TEN, 47 ad ; TPE, 835 ad ; TPA, 844 ad ; Nichols, 1238 ad ; SS14, 1127 ad ; e SS14-Ω, 1738 ad ), o Os intervalos de HPD de 95% se sobrepõem em grande parte às estimativas relatadas em outros lugares 5 , 35 (Tabela de dados estendida 1 ). As duas exceções são TEN e SS14-Ω, que estimamos ter uma origem possivelmente muito mais antiga do que se pensava anteriormente, independentemente de ZH1540 estar incluído. Isto provavelmente se deve ao conjunto de dados mais diversificado usado aqui, que representa com mais precisão toda a diversidade genética da linhagem SS14-Ω. Da mesma forma, ao incluir ZH1540, o tMRCA geral de T. pallidum é estimado como muito mais antigo do que o estimado anteriormente 4 , 5 , 35 .

Advertimos que, embora tenhamos realizado uma análise de relógio relaxado, não modelamos explicitamente as taxas de substituição específicas da linhagem ou dependentes do tempo. Ambos os fenômenos poderiam explicar a idade mais avançada das linhagens de TPA estimadas aqui em comparação com estudos anteriores, e as taxas dependentes do tempo também poderiam empurrar as subespécies e os tMRCAs gerais de T. pallidum ainda mais para o passado. Como tal, os resultados aqui apresentados devem ser interpretados como limites inferiores nos tempos de divergência dos clados de T. pallidum , deixando aberta a possibilidade de estimar tempos de divergência mais antigos com a recuperação de genomas antigos de maior qualidade e o desenvolvimento de modelos de relógio molecular melhorados. .

Discussão

Muitas hipóteses anteriores – baseadas apenas em evidências paleopatológicas – sugeriram infecções treponêmicas precoces entre as populações pré-históricas nas Américas 2 , 30 , 39 . Aqui apresentamos evidências antigas de DNA de uma treponematose pré-colombiana do Novo Mundo, reconstruindo um genoma de T. pallidum de alta cobertura recuperado de restos humanos indígenas brasileiros de quase 2.000 anos de idade, juntamente com três genomas de baixa cobertura do mesmo contexto espaço-temporal. Inesperadamente, estes genomas são notavelmente semelhantes aos do agente causador do bejel moderno, T. pallidum endemicum . Como a sífilis tem sido o foco central da pesquisa sobre o Treponema , as treponematoses endêmicas têm recebido menos atenção 2 , 40 . Ao contrário da bouba e da sífilis, ambas encontradas anteriormente em contextos do Velho e do Novo Mundo desde o início do período moderno 3 , 4 , 5 , 6 , este genoma recentemente reconstruído representa o primeiro patógeno semelhante ao TEN isolado de vestígios arqueológicos e a única comparação antiga com o conjunto atual de oito genomas de bejel publicados 17 , 38 , 41 , 42 . Casos recentes mostraram que as treponematoses são ocasionalmente transmitidas de forma atípica para as suas subespécies geneticamente confirmadas e desafiaram as suas categorizações geográficas e baseadas em sintomas. 43 , 44 , 45 . Nossas descobertas neste estudo apenas reforçam esta visão: um antigo agente semelhante ao NET, identificado longe do nicho geográfico moderno da doença, em uma região costeira brasileira úmida, atesta a capacidade dos treponemas de se adaptarem a vários climas e localizações geográficas. Excluindo os danos ósseos observados em alguns dos restos estudados, os sintomas clínicos, a gravidade e a história evolutiva do antigo patógeno semelhante ao NET recém-descoberto permanecem desconhecidos. Na verdade, as descobertas de outros patógenos antigos, como a presença da peste na Eurásia desde o final do Neolítico, 46 , e Salmonella paratyphi C, com a possível conexão com a grande epidemia de cocoliztli no México em meados do século XVI, mostraram que doenças historicamente devastadoras podem ter representado sorovares inesperados e exibir hoje uma distribuição muito alterada 47 , 48 , 49 , 50 . O DNA treponêmico de alta qualidade recuperado de uma fonte pré-histórica valida o uso de técnicas antigas de DNA no estabelecimento de uma hipótese inteiramente nova e mais informada sobre os eventos que levaram à disseminação do Treponema pallidum pelo mundo.

Bejel em foco

Atualmente, a doença tropical negligenciada bejel é encontrada principalmente nas regiões áridas da África, da Ásia Ocidental e do Mediterrâneo, tornando-a uma candidata improvável a uma potencial treponematose sul-americana em um contexto costeiro. 23. Although genetically unconfirmed, palaeopathological cases of potential treponematoses found worldwide2,51,52,53 may indicate that bejel was previously more widespread and possibly associated with different environmental habitats. Our genomic investigation, together with the radiocarbon dating of both human remains and stratigraphy, places the newly found treponematosis in South America long before the European contact in the 15th century, even predating the Viking expeditions to the North American coast—firmly attesting to the presence of bejel-like treponemal infections in the pre-contact New World. Phylogenetically, this prehistoric form belongs indisputably to the TEN clade, basal to all of its modern strains. Overall, the TEN genomes are highly similar to each other, which may indicate a slow evolution of the lineage as a whole, at least until recently. Regardless of the improved genomic representation of the modern TEN genomes and the newly reconstructed pre-Columbian genome in this study, a larger representation of this lineage would be needed to draw robust conclusions about the evolution and diversification of the subspecies.

Consequências para do T. pallidum a evolução

Os dados apresentados aqui incluem um genoma treponêmico antigo excepcionalmente de alta qualidade e alta cobertura 4 , 5 , 6 , 35 , e adiar as datas das mais antigas cepas de T. pallidum reconstruídas em mais de 1.000 anos. Embora os antigos fragmentos de DNA recuperados neste estudo não fossem adequados para montagem de novo ou análises pan-genômicas, a alta cobertura obtida permitiu-nos realizar uma análise detalhada ao nível do gene e do SNP. Nossas descobertas revelam numerosos eventos de recombinação inter-subespécies, que são conhecidos por serem um mecanismo chave na evolução bacteriana que resulta, por exemplo, na aquisição de fatores de virulência ou outras características adaptativas. Uma vez que os eventos de recombinação aqui identificados envolvem estirpes antigas e modernas, pelo menos uma das formas endémicas permaneceu em proximidade geográfica e num conjunto de hospedeiros comum com as estirpes de TPA após a sua divergência inicial. Quando e onde exatamente esses eventos de recombinação ocorreram é desconhecido. No entanto, a divergência de clados pode ser estimada através da datação por relógio molecular para os diferentes ramos da árvore filogenética de T. pallidum . A calibração deste método é baseada nas idades conhecidas dos genomas utilizados, tornando genomas antigos datados por radiocarbono com segurança indispensáveis ​​para a análise. Nosso novo genoma antigo de alta cobertura fornece um ponto de calibração pré-histórico sem precedentes para a datação do relógio molecular e nos permite concluir que todas as três subespécies já haviam divergido umas das outras antes das viagens de Colombo. As novas estimativas para o tMRCA de todos T. pallidum (12,006–545 bc) and the emergences of the modern clades (TEN: 780 bc–449 ad, TPE: 28–1299 ad and TPA: 42 bc–1376 ad) are much earlier than previous estimates that relied on modern and previously published historical genomes. Yet, these are only the lower bounds of the divergence times, and the subspecies could have originated even earlier. Only genetically ancestral forms of treponemes could illuminate whether the early American strains spread with the early human dispersals—some 15,000–23,000 years ago29,54—ou resultou de um evento local, talvez zoonótico. Finalmente, como a descoberta revolucionária de uma treponematose pré-colombiana aqui é o resultado de uma combinação de genómica de patógenos antigos e da seleção cuidadosa de amostras arqueológicas, podemos esperar que descobertas futuras iluminem os eventos que levaram ao surgimento e propagação da sífilis venérea, e ajudar a resolver os fatores evolutivos responsáveis ​​pelo sucesso global da Treponema . família

Métodos

Inclusão e ética

Os estudos genéticos de doenças humanas antigas esclarecem como as populações do passado prosperaram e lidaram com problemas de saúde, o que pode desencadear preocupações como a estigmatização devido a doenças ou direitos e questões legais entre as pessoas que vivem hoje. As injustiças históricas, a colonização e a expropriação complicaram frequentemente a capacidade das comunidades indígenas de afirmarem e manterem os seus direitos territoriais num quadro jurídico ou administrativo. É, portanto, crucial considerar, além dos aspectos científicos, também as perspectivas das comunidades e pessoas vivas (indígenas) ao realizar este trabalho. 55 .

Aqui estudamos restos humanos de indivíduos totalmente anônimos que morreram há mais de 1.000 anos e foram sepultados no sítio arqueológico Jabuticabeira II, no município de Tubarão, no estado de Santa Catarina, Brasil. Este local foi escavado por P. de Blasis e equipe 56 , financiado pela Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP). e foi obtida autorização de pesquisa junto ao Instituto do Patrimônio Histórico e Artístico Nacional (IPHAN), conforme correspondência 1793/2019 GAB PRESI-IPHAN do processo 01506.000720/2019-65 de K. Santos Bogia. A utilização das amostras dos restos mortais para este estudo também foi aprovada por P. de Blasis, guardião das coleções Jabuticabeira II do Museu de Arqueologia e Etnologia da Universidade de São Paulo. Os restos mortais foram curados, estudados e amostrados por SE e equipe da Universidade de São Paulo até 2016 e, posteriormente, no Museu de História Natural de Viena.

Os territórios e sítios que abrangem o Rio Grande do Sul, Santa Catarina e Paraná são inerentes à herança ancestral das comunidades Kaingang, Guarani e Xokleng (também chamados de 'povo Sol' ou 'grupo costeiro'), que ainda vivem hoje na região. Estas sociedades não só utilizaram a região para mobilização e migração em busca de abastecimento alimentar, mas também percorreram tradicionalmente distâncias significativas, deixando um rasto de marcas culturais, particularmente no domínio das práticas funerárias. Em um estudo anterior 57, samples of five of the individuals exhumed at Jabuticabeira II were studied, revealing some genetic affinity with the Kaingang (a Ge-speaking group of Southern Brazil). However, to the best of our knowledge, the Kaingang are not seen as direct descendants of sambaqui societies, nor do they identify with the people who once dwelled at Jabuticabeira II or request their remains. Finally, research in the Instituto Socioambiental (https://www.socioambiental.org/; for the defence of Brazilian socio-environmental diversity, including Indigenous Rights) states that the region around Jabuticabeira II is not part of any Indigenous reserve, nor are there claims of groups for territorial rights of this region or for the archaeological remains of this site (P. de Blasis, personal communication).

Os processos degenerativos, muitas vezes resultantes de contextos de marginalização, conflito e deslocamento, testemunham o impacto das relações históricas dos grupos indígenas com os colonizadores e invasores. As aflições e doenças vividas por estes grupos carregam ramificações históricas e ambientais de notável significado, garantindo reconhecimento e exame explícitos. Em relação à possível estigmatização das comunidades locais (indígenas) e das pessoas afetadas pelo bejel, deve-se ressaltar que esta doença contagiosa é uma doença endêmica, principalmente não sexualmente transmissível, comum em regiões quentes onde as pessoas vivem em contato próximo umas com as outras, não têm necessidade de cobrir especialmente roupas e compartilhar utensílios. Hoje, o bejel, que pode levar à estigmatização devido a feridas desfigurantes, ocorre especialmente em comunidades do Mediterrâneo Oriental e da África Ocidental com acesso limitado a cuidados médicos modernos. Embora a Organização Mundial da Saúde perceba a importância das ações tomadas para erradicar o bejel em todo o mundo desde 1949 (WHA2.36 Bejel e outras treponematoses ( https://www.who.int/publications/i/item/wha2.36) ), a doença não é vista como um problema atual de saúde pública no Brasil, como é para alguns outros países 58 , 59 . Isto contrasta com as altas prevalências de doenças sexualmente transmissíveis, como o HIV e a sífilis venérea, que afetam as comunidades indígenas no Brasil (Em São Paulo, ação em aldeias promove debate e testagem rápida de HIV e Sífilis — Fundação Nacional dos Povos Indígenas ( http://www.gov.br/funai/pt-br/assuntos/noticias/2019/em-sao-paulo-acao-em-aldeias-promove-debate-e-testagem-rapida-de-hiv-e- sífilis )). É notável, entretanto, que, no contexto arqueológico, nada implicava que aqueles povos pré-históricos de Jabuticabeira II portadores da treponematose local teriam sido discriminados em sua época e cultura.

Além disso, as descrições culturalmente insensíveis em artigos de investigação paleogenómica são uma questão ética de preocupação 60 . Para garantir a discrição, selecionamos expressões potencialmente insensíveis ou discriminatórias no manuscrito. É importante ressaltar que tivemos a ajuda inestimável de E. Krenak, líder da Cultural Survival no Brasil, ativista indígena e doutorando na Universidade de Viena, para analisar criticamente nossos textos e fornecer conselhos sobre o uso eticamente correto e justo da terminologia.

Informação arqueológica

Os sambaquis da região de Laguna

Um sambaqui é o tipo predominante de sítio arqueológico na costa brasileira: um monturo ou montículo de conchas construído pelo homem, de dimensões variadas, localizado em áreas ricas em recursos, como lagoas, manguezais ou estuários. Sambaquis consistem em sedimentos inorgânicos, conchas de moluscos, restos de alimentos e matéria orgânica misturados em intrincadas estratigrafias associadas a funções domésticas e/ou funerárias. 61 . Mais de 1.000 sambaquis estão mapeados ao longo dos 7.500 km de extensão da costa brasileira e são datados entre 7.500 e 1.000 anos bp 61 , 62 . Pesquisas arqueológicas recentes sugerem que essas populações construtoras de sambaquis eram sedentárias, com uma abundante e estável subsistência marinha, horticultura e alto crescimento populacional. 61 , 63 , rituais funerários elaborados 64 e apropriação da paisagem 56 .

Jabuticabeira II excavation site

Jabuticabeira II (UTM 22 J − 0699479E; 6835488 S) é um sambaqui de médio porte (400 × 250 × 10 m de altura), assentado sobre uma paleoduna e localizado na região de Laguna, área de maior densidade de sambaquis do sul Costa brasileira, a 3 km da Laguna do Camacho, uma das diversas fontes de água associadas a um sistema geológico barreira-lagoa formado durante o Holoceno (Fig. 1 ). Jabuticabeira II, construída ao longo de quase mil anos, é um dos 65 sambaquis mapeados ao redor do sistema lagunar. Este grande número de assentamentos e sua história de ocupação cronologicamente sobreposta atestam uma ocupação bastante densa e intensas interações dos construtores de sambaqui entre 7.500 e 900 anos calibrados (ano cal) bp 56. According to stratigraphic studies, Jabuticabeira II is the result of incremental funerary rituals accumulated over centuries. Although Jabuticabeira II was not completely excavated, 204 burials containing the remains of 282 individuals were exhumed from a 373 m2 area64. Radiocarbon dates of Jabuticabeira II stratigraphy56,64 suggest a long occupation period between 1214–830 cal bc and 118–413 cal ad or 3137–2794 to 1860–1524 cal yr bp (2σ), roughly in line with the radiocarbon datings from bone material of the four individuals in this study, ranging from 350 cal bc to 573 cal ad.

Os restos mortais do sambaqui Jabuticabeira II foram encontrados em sepultamentos simples, duplos e múltiplos, dispersos em aglomerados. Os esqueletos recuperados estavam em sua maioria incompletos, evitando estimativas categóricas de idade e sexo ou outros achados osteológicos. O padrão de sepultamento era bem flexionado e sugeria tratamento intencional do corpo antes do internamento. O pequeno tamanho das sepulturas sugeria que os corpos sofreram dessecação ou decomposição prévia de tecidos moles, mas não o suficiente para produzir desarticulação completa (os ossos das mãos e dos pés foram encontrados articulados). Muitos enterros vêm de perfis e estão incompletos. Os ossos de vários indivíduos estão manchados de ocre vermelho 65 , uma prática comum em sítios arqueológicos do estado de Santa Catarina 66 , 67 . As oferendas são comuns em contextos funerários de sambaqui e incluem adornos confeccionados com material faunístico e ferramentas líticas em diversas formas, desde entulhos até ferramentas polidas e zoólitos, com diferenças na frequência de ocorrência entre diferentes locais e estratos 68 . A oferta mais comum na Jabuticabeira II era o peixe.

Ao todo, 99 indivíduos de Jabuticabeira II, com e sem alterações ósseas sugestivas de infecção, foram triados quanto ao conteúdo de DNA do patógeno. 37 amostras foram consideradas positivas para DNA treponêmico na triagem inicial e quatro amostras produziram dados suficientes para a reconstrução do genoma de T. pallidum (Tabela Suplementar 1 ).

Palaeopathological analysis of treponematoses

Bioarchaeological analyses showed results compatible with increasing population growth and high population density in Jabuticabeira II, including high frequencies of nonspecific stress markers69 and occasional infant stress70, but no evidence of trauma associated with interpersonal conflicts over resources or territory69.

Há, no entanto, evidências de doenças sistêmicas transmissíveis em Jabuticabeira II e outros sambaquis locais brasileiros 30. Eleven 14C accelerator mass spectrometry dates obtained directly from the presumably treponematosis-affected individuals suggest that these diseases are very old on the east coast of South America, with a time-range between 6,300 and 500 yr BP. Among the possible treponemal cases based on osteological analysis, three came from Jabuticabeira II. However, these did not overlap with the individuals yielding the detected genetic evidence in this study.

Informações sobre indivíduos

Indivíduo 41A-L2.05-E4, amostra ZH1390

O indivíduo é um homem adulto de constituição robusta, com estatura estimada em 150,49 ± 2,6 cm (ref. 70 ). Embora fragmentados, os ossos deste indivíduo constituíam um esqueleto quase completo (80%), articulado e enterrado em uma matriz oval rica em conchas em posição hiperflexionada. Os ossos do indivíduo apresentavam sinais de doença infecciosa sistêmica em membros inferiores. Fêmures, tíbias e fíbulas apresentam periostite generalizada discreta e osteoartrose. Também foi observado alargamento na porção lateral das clavículas. Segundo Filippini et al. 30 , aplicando o método ESPIRAL 71 , a doença desse indivíduo poderia ser classificada de forma não conclusiva como sífilis, bouba ou bejel. A amostragem foi realizada em lesão ativa no fragmento da tíbia.

FS9–L3–T2 individual, amostra ZH1540

A amostra vem de um conjunto de ossos misturados, provavelmente de mais de um indivíduo. Os ossos atribuídos a este indivíduo consistem em vários elementos esqueléticos, alguns com alterações patológicas, como osteomielite grave no terço distal do úmero direito, periostite grave na ulna esquerda, periostite na diáfise da fíbula e dois corpos vertebrais com osteofitose. A amostra foi retirada do fragmento da fíbula, na área com periostite.

FS3B-L3-T4 individual, amostra ZH1541

A amostra vem de um de três indivíduos separados, encontrados misturados. Os elementos esqueléticos pertencentes a este adulto robusto de idade e sexo desconhecidos incluem um rádio esquerdo com artrite, um fragmento da ulna esquerda (muito robusto), um fragmento do úmero esquerdo, fragmentos de um fêmur, uma tíbia e uma fíbula e um primeiro metatarso. A amostra foi retirada de um fragmento de fêmur, sob a superfície imediata do osso, para melhor evitar a possível introdução de contaminantes externos.

Individual 2B-L6-E3, ZH1557

A amostra provém de um indivíduo provavelmente adulto do sexo masculino. O indivíduo estava articulado e em posição flexionada com outro indivíduo adulto do sexo feminino enterrado em cima. Os achados osteopatológicos nos ossos do indivíduo amostrado incluíram sinais de doença articular degenerativa, osteoartrite intervertebral lombar grave, escoliose e possíveis lesões na patela. No entanto, não foram observadas lesões típicas sugestivas de infecção treponêmica. A amostra foi retirada de um pequeno pedaço de osso longo, sob a superfície imediata do osso, para melhor evitar a possível introdução de contaminantes externos.

Correção do efeito do reservatório marinho para 14 C namoro

A datação por radiocarbono foi realizada pelo Laboratório de Física de Feixes de Íons da ETH Zurich (laboratório número: ETH-127328) utilizando colágeno ósseo purificado por um método de ultrafiltração modificado 72 . A calibração dos dados foi feita com OxCal v4.4.4. A dieta dos habitantes de Jabuticabeira II, constituída substancialmente por fontes alimentares marinhas, produz efeito reservatório nas datas de radiocarbono calculadas como idade média de 247,8 ( σ = 103,7) anos 73 . Considerando a alta contribuição de carbono marinho para o colágeno ósseo dos indivíduos de Jabuticabeira II, as datas de radiocarbono foram modeladas com Calib Rev 8.20 74 ( http://calib.org/calib/calib.htm ) usando a curva de calibração Mixed Marine SHCal20 75 , 76 e aplicando o valor médio estimado de correção de reservatórios de radiocarbono marinhos locais (Δ R ) de -126 ± 29 para a costa sul do Brasil (banco de dados de Correção de Reservatórios Marinhos) 73,77. We considered the average relative contribution of marine carbon to collagen derived from Bayesian Mixing Models for Jabuticabeira II individuals, calculated at a mean value of 42.5%78,79. Para as estimativas individuais para as amostras, consulte a Tabela Suplementar 2 .

Sample processing

Samples were documented and carried through sampling, DNA extraction, library preparation and library indexing in facilities dedicated to ancient DNA work at the University of Zurich, including decontamination of samples, laboratory equipment and reagents with UV irradiation and using protective clothing and minimum contamination-risk working methods.

Todas as etapas pós-amplificação foram realizadas nas instalações laboratoriais regulares disponíveis para o Grupo de Paleogenética do Instituto de Medicina Evolutiva (IEM) da Universidade de Zurique (UZH). O sequenciamento de DNA foi realizado nas instalações de Sequenciamento de Próxima Geração das Instalações Principais do BioCenter de Viena (VBCF) ou no Centro de Genômica Funcional da Universidade de Zurique (FGCZ).

Extração de DNA antigo

Todas as superfícies das amostras foram irradiadas com luz ultravioleta para minimizar a contaminação potencial do DNA moderno. O pó ósseo foi obtido utilizando-se broca odontológica e brocas com cabeça diamantada. A extração de DNA foi realizada em cerca de 50-100 mg de pó de osso, de acordo com um protocolo de extração bem estabelecido para DNA antigo. 80 . Os controles negativos para processos de extração e biblioteca foram processados ​​paralelamente em todos os experimentos, um controle para cada dez amostras, sequenciados e comparados bioinformaticamente com seus correspondentes lotes de amostras, como precaução contra possível contaminação.

Preparação da biblioteca

Bibliotecas de DNA de cadeia dupla foram produzidas para triagem inicial com sequenciamento shotgun, sem tratamento UDG (ou seja, tratamento químico com o objetivo de limitar os danos relacionados à idade no DNA). Duas bibliotecas adicionais para cada uma das amostras potencialmente positivas da primeira rodada de captura foram produzidas para maximizar a complexidade do DNA. Para a preparação de bibliotecas de DNA, 20 µl de extrato de DNA foram convertidos em bibliotecas de DNA de fita dupla 31 . Códigos de barras (índices) específicos da amostra foram adicionados a ambas as extremidades dos fragmentos de DNA nas bibliotecas 81 . As bibliotecas indexadas foram então amplificadas para atingir uma concentração mínima de DNA de aproximadamente 90 ng ml −1 . A amplificação foi realizada utilizando 1 × tampão Herculase II, 0,4 mM IS5 e 0,4 mM IS6 primer 81 , DNA polimerase de fusão Herculase II (Agilent Technologies), dNTPs 0,25 mM (100 mM; 25 mM cada dNTP) e 5 ml de biblioteca indexada como modelo de DNA. Foram preparadas quatro reações por biblioteca e o volume total da reação de amplificação foi de 100 ml. O perfil térmico incluiu uma desnaturação inicial por 2 min a 95 °C e 3-18 ciclos, dependendo da concentração de DNA após a indexação das bibliotecas, desnaturação por 30 s a 95 °C, 30 s de recozimento a 60 °C e 30 s de recozimento a 60 °C. alongamento a 72°C, seguido por uma etapa de alongamento final por 5 min a 72°C. Todas as divisões de uma biblioteca indexada foram reunidas e purificadas utilizando o kit de purificação QIAGEN MinElute PCR. As bibliotecas de DNA foram então quantificadas com D1000 ScreenTape em um Agilent 2200 TapeStation (Agilent Technologies) e combinadas em pools equimolares para sequenciamento.

Pathogen screening

Os dados do Shotgun foram usados ​​para uma triagem inicial das 99 amostras candidatas, com o software Kraken2 82 , e 41 amostras que tiveram mais de 7 ocorrências para T. pallidum foram selecionadas para enriquecimento alvo. As amostras selecionadas foram submetidas a um processo de enriquecimento alvo e posteriormente processadas pelo FastQ Screen v0.15.1 83 para verificar o número de leituras mapeadas em relação a três genomas de referência representativos de alta qualidade da subespécie de T. pallidum (CDC2, BosniaA e Nichols). As nove amostras mais promissoras (> 5.000 acessos de Kraken a T. pallidum após a primeira rodada de captura em solução) foram transformadas em duas bibliotecas extras e recapturadas conforme explicado em detalhes nas seções a seguir.

Enriquecimento alvo para T. pallidum DNA de

O enriquecimento de bibliotecas de cadeia dupla em todo o genoma foi realizado através de kits de enriquecimento de alvo personalizados (Arbor Bioscience). Iscas de RNA com comprimento de 60 nucleotídeos e densidade de ladrilhos de 4 pb foram projetadas com base em três genomas de referência: Nichols ( CP004010.2 ), SS14 ( CP000805.1 ), Fribourg-Blanc ( CP003902 ). Conjuntos de bibliotecas de 500 ng foram enriquecidos de acordo com as instruções do fabricante. As bibliotecas capturadas foram amplificadas em reações de 100 µl contendo 1 unidade de DNA polimerase de fusão Herculase II (Agilent), 1 × tampão de reação Herculase II, dNTPs 0,25 mM, iniciadores IS5 e IS6 0,4 mM. 81 e modelo de biblioteca de 15 µl, com o seguinte perfil térmico: desnaturação inicial a 95 °C por 2 min, 14 ciclos de desnaturação a 95 °C por 30 s, recozimento a 60 °C por 30 s e alongamento a 72 °C por 30 s, seguido de um alongamento final a 72°C por 5 min. As bibliotecas capturadas foram purificadas com colunas de rotação MinElute (QIAGEN) e quantificadas com uma ScreenTape de alta sensibilidade D1000 em um Agilent 2200 TapeStation.

Sequencing

For both shotgun data retrieval and after the capture processing, the samples were pooled in unimolar quantity (for SG sequencing up to 50 samples per pool, and for the capture process 2–8 samples per pool), and sequenced on an Illumina NextSeq500 with 2 × 75  +  8  +  8 cycles using the manufacturer’s protocols for multiplex sequencing at the Functional Genomics Center in Zurich or at the Vienna BioCenter Core Facilities.

Statistical analyses

Dataset selection

Montamos um conjunto de dados genômicos compreendendo 98 genomas de T. pallidum disponíveis publicamente (8 TEN, 30 TPE e 60 TPA) de estudos publicados anteriormente (incluindo 8 genomas antigos) e o recém-gerado genoma ZH1540. Os genomas representam a variação genética das três subespécies conhecidas de T. pallidum (TPA, TPE e TEN) disponíveis até dezembro de 2022, e foram selecionados com foco em TEN e TPE, devido à sua proximidade com o novo genoma antigo classificado como TEN .

Os dados publicados para o conjunto de dados do genoma moderno neste estudo estão disponíveis no banco de dados European Nucleotide Archive (ENA): PRJNA313497 (números de acesso: SRR3268682 , SRR3268724 , SRR3268715 , SRR3268694 , SRR3268696 , SRR3268709 , SRR3268710 ), 11481 (números de acesso: ERR1470343 , ERR3596780 , ERR3596747 , ERR3596783 ), PRJEB28546 (números de acesso: ERR4045394 , ERR3684452 , ERR3684456 , ERR3684465 , SRR13721290 , ERR4853530 , ERR4993349 , ERR485358 7 , ERR4899206 , ERR5207017 , ERR5207018 , ERR5207019 , ERR4899215 , ERR4853623 , ERR4853625 ), PRJNA508872 (números de acesso: SRR8501165 , SRR8501164 , SRR8501167 , SRR8501166 , SRR8501168 , SRR8501171 ), PRJNA723099 (números de acesso: SRR14277267 , SRR14277266 , SRR14277458 , SRR14277444 ), PRJEB11481 (número de acesso: ERR 1470331 ), PRJDB9408 (números de acesso: DRR213712 , DRR213718 ), PRJNA588802 (números de acesso: SRR10430858 , SRS5636328 ), PRJNA322283 (número de acesso: SRR3584843 ), PRJNA754263 (números de acesso: SRR15440297 , SRR15440150 , SRR15440451 , SRR15440240 ), PRJEB40752 (números de acesso: ERR4690809 , ERR4690806 , ERR4690810 , ERR4690812 , ERR4690811 ). Arquivos de montagem foram usados ​​para 9 genomas do banco de dados do National Center for Biotechnology Information (NCBI): CP002375.1 , CP002376.1 , NC_016842.1 , NC_017268.1 , NC_018722.1 , NC_021490.2 , NC_021508.1 , GCA_000813285.1 , CP035193.1 e para 24 genomas modernos do European Nucleotide Archive (ENA): CP021113.1 , CP073572.1 , CP073557.1 , CP073553.1 , CP073536.1 , CP073526.1 , CP073490.1 , CP073487.1 , CP07347 0 .1 , CP073447.1 , CP073446.1 , CP073399.1 , CP040555.1 , LT986433.1 , LT986434.1 , CP032303.1 , CP020366.1 , CP024088.1 , CP024089.1 , CP07812 , 1.1 , CP078090.1 , CP081507.1 , CP051889.1 e CP003902 .1 . Os dados de sequência bruta (arquivos fastq) usados ​​para 6 genomas modernos estão disponíveis no banco de dados NCBI: PRJEB20795 (números de acesso: ERS1724928 , ERS1724930 , ERS1884567 ) e PRJNA343706 (números de acesso: SRR4308604 , SRR4308606 , SRR4308597 ). Genomas treponêmicos antigos publicados anteriormente aqui usados ​​estão disponíveis na ENA: PRJEB37490 (número de acesso: ERR4065503 ), PRJEB37633 (número de acesso: ERR4000645 ), PRJEB35855 , PRJEB21276 (números de acesso: ERS2470995 , ERS2470994 ) e PRJEB62102 . Informações detalhadas sobre a fonte do conjunto de dados de referência estão documentadas na Tabela Suplementar 3 .

Selecionamos todos os oito genomas TEN disponíveis publicamente, todos com mais de 99,4% de cobertura genômica, com exceção de C77 17 (81,4%). Selecionamos 30 genomas de TPE (Tabela Suplementar 3 ). Para representar cada linhagem ou sub-linhagem, selecionamos pelo menos um genoma, preferindo aqueles com maior profundidade de sequenciamento e cobertura genômica. Todos os genomas TPE incluídos têm mais de 95,3% de cobertura genômica, exceto os quatro antigos genomas TPE: SJN003, AGU007, 133 e CHS119, exibindo 97,4%, 92,7%, 57% e 62% de cobertura genômica, respectivamente. Além disso, foram incluídos 60 genomas de TPA das principais linhagens e sub-linhagens descritas em estudos anteriores (Tabela Suplementar 3 ). Todos esses genomas tinham mais de 90% de cobertura, exceto os quatro genomas antigos, PD28, W86, SJ219 e 94B, todos com cobertura genômica de 30% ou mais. Todos os genomas no conjunto de dados são separados uns dos outros por pelo menos 5 SNPs. A cepa TPA Seattle-81 foi excluída do conjunto de dados final devido a mutações provavelmente acumuladas durante passagens extensas em coelhos que podem causar posicionamento ambíguo em filogenias 4 , 16 , 36 .

Os dados brutos e/ou arquivos de montagem para cada genoma em nosso conjunto de dados foram baixados dos bancos de dados públicos: European Nucleotide Archive (ENA) 84 e Centro Nacional de Informações sobre Biotecnologia (NCBI) 85 . Os números de acesso são fornecidos na Tabela Suplementar 3 .

Processamento de leitura e geração de alinhamento de genoma baseado em múltiplas referências

To reconstruct the individual genomes from the raw data, we carried out raw read quality control and preprocessing, removing duplicates, variant calling and filtering using the default parameters when not otherwise specified. After processing the de-multiplexed sequencing reads, sample sequencing quality was analysed with FastQC version 0.11.983, filtering reads with a QC value < 25. Following processing by cutadapt version 4.186 to remove the sequencing adapters, in order to reduce the reference bias, and improve the posterior phylogenetic inference and assignment87, the genome reference selection for mapping each sample was determined according to the results from the original manuscript where the genomes were published (see Supplementary Table 3). The mapping was carried out by BWA mem88(usando parâmetros: -k 19, -r 2,5). Foram utilizados quatro genomas de referência; os bem estudados genomas TEN e TPE BosniaA (NZ_CP007548.1) e CDC2 ( NC_016848.1 ), bem como os genomas Nichols ( NC_021490.2 ) e SS14 ( NC_010741.1 ), representando as duas linhagens principais dentro do TPA. No entanto, para as novas amostras antigas obtidas aqui, os genomas de cada amostra foram reconstruídos mapeando três genomas de referência de alta qualidade, representando as três subespécies de T. pallidum (CDC2, BosniaA e Nichols).

CleanSam, do Picard Toolkit versão 2.18.29 ( http://broadinstitute.github.io/picard ), foi usado para limpar os arquivos SAM ou BAM fornecidos. As leituras duplicadas foram removidas usando MarkDuplicates, do kit de ferramentas Picard versão 2.18.29. AddOrReplaceReadGroups, do Picard Toolkit versão 2.18.29, foi usado para atribuir todas as leituras em um arquivo a um único novo grupo de leitura antes de usar mapDamage versão 2.2.0-86-g81d0aca 89 para estimar os parâmetros de dano ao DNA e redimensionar os índices de qualidade de posições provavelmente danificadas nas leituras (usando o parâmetro: --rescale).

Depois de gerar uma saída de pilha de texto para os arquivos BAM com a ferramenta mpileup do Samtools versão 1.7 90 , os SNPs foram chamados usando o VarScan versão 2.4.3 91 (usando parâmetros: -p-valor 0,01, -min-reads2 1, -min-coverage 1, -min-freq-for-hom, 0,4 -min-var-freq 0,05, -output-vcf 1). Em seguida, uma filtragem SNP também foi realizada com VarScan (usando para as amostras modernas os parâmetros: -p-value 0,01, -min-reads2, 5 -min-coverage 10, -min-avg-qual 30 -min-freq-for -hom 0,4, -min-var-freq 0,9, -output-vcf 1 e modificação de alguns parâmetros para as amostras antigas devido à sua menor cobertura e qualidade de leitura: -p-value 0,01 -min-reads2 3, -min-coverage; 5, -min-avg-qual 30, -min-freq-for-hom 0,4, -min-var-freq 0,9 -output-vcf 1). Além disso, todas as posições com menos de 3 leituras mapeadas foram mascaradas com Genomecov do Bedtools versão 2.26.0 92 para amostras modernas e antigas. Todas as etapas de geração do genoma foram visualizadas e confirmadas manualmente com Tablet versão 1.21.02.08 93 , verificando cada SNP um por um e descartando os possíveis SNPs espúrios do novo genoma antigo ZH1540. As sequências finais resultantes foram obtidas por maskfasta do Bedtools v2.26.0.

Além disso, utilizamos metodologias testadas de sequenciamento e análise posterior 17 , 42 para obter maior cobertura e genomas modernos de T. pallidum mais confiáveis . Sempre que possível, foram obtidos arquivos de montagem em vez de dados brutos (Tabela Suplementar 3 ). Um alinhamento de genoma baseado em referência múltipla para todas as sequências foi gerado no MAFFT v7.467 94 (usando parâmetros: --adjustdirection --auto --fastaout --reorder). Porém, devido ao uso de diferentes referências genômicas, regiões com baixa cobertura para alguns genomas, correspondendo principalmente aos genes tpr e arp , foram revisadas e alinhadas manualmente com o Aliview versão 1.25. 95 .

As amostras ZH1390, ZH1541 e ZH1557 tinham dados suficientes para tentar uma reconstrução do genoma e foram determinadas como tendo o maior número de SNPs em comum com a referência TEN, mas foram excluídas das análises downstream devido à cobertura limitada adquirida para cada uma delas, o que fez os SNPs obtidos são menos confiáveis. A amostra ZH1540, no entanto, rendeu uma cobertura genômica notável de 33,6× e foi selecionada para análises aprofundadas subsequentes.

Proteinorto versão 6.0b 96 (usando parâmetros: -p=blastn -singles -keep) foi usado para conduzir um estudo de ortologia a fim de encontrar genes ortólogos nos quatro genomas de referência utilizados 96 . Cada gene presente em pelo menos um dos quatro genomas de referência teve suas coordenadas genômicas determinadas com base em sua localização no alinhamento final mesclado (ver Tabela Suplementar 3 ).

Para verificar a precisão do alinhamento final do genoma múltiplo, e que nenhum gene codificador de proteína foi inadvertidamente truncado, as traduções de proteínas para cada gene presente em pelo menos um genoma de referência foram comparadas com os arquivos gff3 originais de cada uma das quatro referências (Suplemento Tabela 3 ). O genoma ZH1540 reconstruído e suas principais características foram representados graficamente usando BRIG versão 0.95-dev.0003 97 .

Análise de recombinação usando PIM

Como observado anteriormente 36 , a presença de recombinação nos genomas de T. pallidum pode interferir nas topologias das árvores filogenéticas inferidas. Para analisar a potencial recombinação genética, usamos o pipeline PIM 36 para detectar recombinação gene por gene. Em resumo, o processo envolveu as seguintes etapas:

  1. (1)

    Using IQ-TREE version 1.6.10, a maximum-likelihood tree was created for the multiple genome alignment98. All maximum-likelihood trees for the remaining steps were obtained using GTR99 + G100 + I101 as an evolutionary model and 1,000 bootstraps replications.

  2. (2)

    The 1,161 genes found in at least one of the reference genomes were extracted, and the number of SNPs for each gene was calculated. Genes with less than three SNPs were excluded.

  3. (3)

    The phylogenetic signal in each gene alignment for each of the remaining genes was evaluated by likelihood mapping102 in IQ-TREE (using parameters: -lmap 10000 -n 0), retaining only those genes that showed a phylogenetic signal.

  4. (4)

    A maximum-likelihood tree was generated for each of the remaining genes using IQ-TREE.

  5. (5)

    Para cada gene incluído, testamos a congruência filogenética entre árvores usando IQ-TREE (usando parâmetros: -m GTR + G8 -zb 10000 -zw), comparando a árvore de máxima verossimilhança obtida a partir do alinhamento do gene e a árvore de máxima verossimilhança obtida do alinhamento do genoma completo usando dois métodos diferentes: Shimodaira – Hasegawa 103 e pesos de verossimilhança esperados (ELW) 104 . Genes para os quais pelo menos um teste rejeitou a topologia da árvore de referência com o alinhamento do gene adotando uma abordagem conservadora ( P <0,2, valor de peso próximo a 0, para os testes Shimodaira-Hasegawa e ELW, respectivamente) e o alinhamento completo do genoma rejeitou a topologia de as árvores construídas utilizando o alinhamento gênico (incongruência recíproca, P < 0,2 e valor de peso próximo a 0) em pelo menos um deles foram selecionadas e examinadas mais de perto na etapa seguinte.

  6. (6)

    Usando MEGAX 105 , os genes selecionados que apresentavam incongruência recíproca foram posteriormente examinados para avaliar e descrever potenciais eventos de recombinação. Um gene deve ter pelo menos três SNPs homoplásicos próximos – SNPs que são compartilhados por vários grupos (TPE, TEN, TPA-Nichols ou TPA-SS14) e produzir uma distribuição polifilética – para ser rotulado como recombinante. Os SNPs homoplásicos encontrados no alinhamento dos genes serviram como limites das áreas recombinantes.

  7. (7)

    Utilizando um critério de parcimônia na distribuição de estados alternativos dos SNPs homoplásticos, foram inferidos os potenciais clados ou cepas doadoras e receptoras de cada evento de recombinação.

DNA sections, a number of genes have a high percentage of sites with missing data. The majority of these genes are members of the tpr and arp families, which include collections of paralogous genes. In order to continue analysing these intriguing genes with the PIM pipeline, strains that had a high percentage of missing data in each of these genes were eliminated. Following previous findings35,36, the hypervariable gene tprK (tp0897), with seven hypervariable regions that undergo intrastrain gene conversion17,37,106,107,108,109, and the tp0316 and tp0317 genes, also under gene conversion, were completely excluded from the recombination analysis.

PIM procedure for likelihood mapping and topology tests

Um teste de mapeamento de verossimilhança foi realizado usando IQ-TREE para determinar quais genes (Tabela Suplementar 4 ) mostraram um sinal filogenético (dos 382 genes para os quais> 3 SNPs foram encontrados em comparação pareada com pelo menos um genoma de referência). Para cada quarteto (subconjunto de quatro sequências) nos dados, o teste cria árvores filogenéticas não enraizadas. As probabilidades do quarteto são então traçadas dentro de um triângulo, onde a posição denota a “semelhança em árvore” do quarteto em questão. Os quartetos de canto estão completamente resolvidos, os quartetos nas laterais estão parcialmente resolvidos e os quartetos no centro não estão resolvidos. Dos 382 genes, 29 tinham muitos valores faltantes para serem testados usando o método de mapeamento de verossimilhança. Para incluir esses genes nas próximas etapas do pipeline PIM e nas comparações de topologia, as sequências problemáticas com mais de 50% de posições com dados faltantes foram removidas.

Após o teste de mapeamento de verossimilhança, 9 genes pertencentes à zona central do triângulo foram descartados (Tabela Suplementar 4 ). Em seguida, usando os testes de topologia Shimodaira-Hasegawa e ELW, comparamos as árvores gênicas dos genes restantes com a árvore de referência preliminar do alinhamento do genoma completo para avaliar sua congruência filogenética (Tabela Suplementar 4 ). Dos 373 genes que testaram positivo para incongruência filogenética, 27 continham pelo menos três SNPs consecutivos, apoiando um evento de recombinação. A estes adicionamos tp0859 , que foi detectado como recombinante em um estudo anterior 35 , resultando em um total de 27 genes recombinantes.

Análise de recombinação usando Gubbins e ClonalFrameML

Gubbins version 2.3.1110 and ClonalFrameML version 1.11-1111 are frequently used tools for the genome-wide identification of recombinant positions in bacterial genomes. To test the robustness of our recombination analysis using PIM, we also ran these two programs, with default parameters and the same whole-genome alignment used with PIM. Gubbins identified 301 distinct recombination events associated with 103 genes, ranging in size from 5 bp to 13,866 bp. Similarly, ClonalFrameML detected 656 events, with 32 of them being 1 or 2 bp long, and the longest event spanning 782 bp. Notably, all the genes identified by PIM as having a recombinant region were also detected by both ClonalFrameML and Gubbins, except for gene tp0558, which was missed by ClonalFrameML but detected by Gubbins. Additionally, genes tp0164 and tp0179 were detected by ClonalFrameML but missed by Gubbins.

Phylogenetic analysis

Uma árvore de máxima verossimilhança baseada no alinhamento incluindo todos os genes foi construída com IQ-TREE, usando GTR + G + I como modelo evolutivo e 1.000 replicações de bootstrap (Extended Data Fig. 2a ). Em seguida, os genes identificados como recombinantes pelo PIM foram removidos do alinhamento do genoma múltiplo. Três genes adicionais ( tp0897 , tp0316 e tp0317 ), que contêm regiões repetitivas e foram identificados como hipervariáveis ​​e/ou sob conversão gênica no passado, também foram removidos para evitar a introdução de um viés potencial. Como o gene tp0317 está aninhado dentro do gene tp0316 e as coordenadas do genoma de referência da BósniaA para tp0316 cobriram uma área maior do que as dos outros genomas de referência, tp0316 e tp0317 foram removidos de acordo com as coordenadas tp0316 do genoma de referência da BósniaA. Uma árvore filogenética de referência foi então construída empregando o novo alinhamento do genoma de herança vertical, também com IQ-TREE usando GTR + G + I como modelo evolutivo e 1.000 replicações de bootstrap (Extended Data Fig. 2b ). Ambas as árvores obtidas foram comparadas e são mostradas na Fig. 2 de dados estendidos .

A linhagem SS14 foi descrita anteriormente como um grupo amplamente epidêmico e resistente a macrolídeos que surgiu após, e possivelmente foi motivado pelo, uso clínico de antibióticos após sua descoberta. 12 , 16 . Com base nos resultados da nossa análise filogenética e expandindo as classificações filogenéticas anteriores e a nomenclatura da linhagem SS14 12 , 16 , definimos o clado que contém quase todos os genomas SS14 de amostras clínicas e contemporâneas como a sub-linhagem SS14-Ω. No entanto, duas amostras clínicas contemporâneas (MD18Be e MD06B) não foram classificadas como sub-linhagem SS14-Ω, porque essas amostras se agrupam com o genoma do MéxicoA, em linha com resultados publicados anteriormente 42.

To compare the PIM-based analysis with other widely used recombination detection methods, Gubbins and ClonalFrameML, we followed a similar procedure of removing the recombinant positions detected by these tools and inferred maximum-likelihood trees with the retained positions in the corresponding multiple genome alignments. All the phylogenetic trees with recombination events removed exhibit general congruence with each other, whether the events were identified by PIM, Gubbins or ClonalFrameML. Furthermore, the placement of the ZH1540 genome remained consistent in the phylogenetic trees, regardless of the recombination detection method employed, and despite the elimination of recombinant genes to generate the vertically inherited alignment.

Exploratory characterization of the 16S-23S genes

T. pallidum contains two rRNA (rrn) operons, each of which encodes the 16S-23S-5S rRNA genes and intergenic spacer regions (ISRs). There is evidence that the random distribution of rrn spacer patterns in T. pallidum may be generated by reciprocal translocation of rrn operons mediated by a recBCD-like system found in the intergenic spacer regions (ISRs)112 . De acordo com estudos anteriores 112 , 113 , 114 , 115 , descobrimos que os ISRs 16S-23S das cepas de TPA contêm os genes tRNA-Ile (tRNA-Ile-1; tp0012 ) e tRNA-Ala (tRNA-Ala-3; tp00t15 ) dentro dos rrn1 e rrn2 operons , respectivamente . Por outro lado, os genomas TPE mostram um padrão espaçador Ala/Ile, onde os tp0012 e tp00t15 ortólogos estão localizados dentro dos operons rrn2 e rrn1 , respectivamente.

Identificamos 68 SNPs nos genes r0001 , r0002 , r0004 e r0005 , codificando os genes 16S-23S rRNA do novo e antigo genoma ZH1540, colocando-os entre os genes mais variáveis ​​em nosso alinhamento e aumentando o potencial de que incluí-los no alinhamento poderia resultar em uma reconstrução filogenética tendenciosa. Embora os SNPs encontrados pareçam ser bem suportados pelas leituras obtidas a partir do mapeamento de sequências (Tabela Suplementar 3 ), a sua origem de possível contaminação não pode ser completamente descartada e seriam necessárias análises adicionais para os confirmar.

A exclusão destes genes do alinhamento, além dos genes recombinantes e tp0316, tp0317 e tp0897 , não resultou em quaisquer alterações na topologia (Extended Data Figs. 2 b e 3 ), embora os comprimentos dos ramos tenham sido alterados. Como se sabe que esses genes possuem regiões conservadas, além de regiões variáveis, usadas para explorar as relações evolutivas entre bactérias patogênicas 116 , 117 , 118 , decidimos mantê-los no alinhamento para todas as análises subsequentes. Finalmente, notamos que o genoma ZH1540 não possuía nenhuma das duas mutações do gene de RNA ribossômico 23 S de T. pallidum conhecidas por conferir resistência a macrólidos (A2058G e A2069G). Em contraste, quatro estirpes modernas de TEN do Japão possuem a mutação A2048G, sugerindo uma recente pressão de selecção para mutações de resistência a antibióticos.

Molecular clock dating

We used the Bayesian phylogenetics package BEAST2 v2.6.7119 to estimate a time-calibrated phylogeny of the context dataset of 98 T. pallidum genomes along with our new ancient genome, ZH1540. We removed hypervariable and recombining genes from the alignment, as described above, reduced it to variable sites and used an ascertainment bias correction to account for constant sites.

Análises de regressão da raiz às pontas (Dados Estendidos Fig. 4 ) mostram que, embora haja uma correlação positiva entre o ano de amostragem e a divergência da raiz às pontas entre todas as cepas clínicas modernas, indicando uma evolução semelhante a um relógio, a correlação é muito fraca quando incluindo também cepas passadas e negativo quando incluindo cepas antigas. Dentro dos clados TPE, TEN e SS14 existe uma correlação positiva entre todas as cepas clínicas modernas e passadas. Por outro lado, a correlação é negativa para cepas de Nichols, mesmo quando se olha apenas para cepas clínicas. A fim de levar em conta a variação da taxa e os longos ramos terminais em algumas cepas (provavelmente devido a uma infinidade de efeitos, incluindo erros de sequenciamento, contaminação e mutações introduzidas durante a passagem do coelho), usamos um modelo de relógio UCLD e UCED para a datação do relógio molecular. análise 120 . Para ambos os modelos, colocamos um prior lognormal estreito com uma média (no espaço real) de 1 × 10 −7 substituições por local por ano e desvio padrão de 0,25 na frequência média do clock. Este forte anterior foi usado para compensar o fraco sinal temporal entre os genomas de T. pallidum e foi calibrado em estimativas anteriores da taxa de substituição 4 , 35 . Usamos ainda um modelo de substituição GTR + G + I 118 e um gráfico do horizonte bayesiano 121 modelo demográfico (árvore anterior) com 10 grupos. Para todos os genomas onde as datas de amostragem não são conhecidas exatamente, usamos anteriores uniformes ao longo dos intervalos de datas relatados nos estudos originais para levar em conta a incerteza 4 , 5 , 6 , 16 , 122 . Para ZH1540, definimos o intervalo de datas para 364–573 dC , de acordo com os resultados da datação por radiocarbono corrigidos pelo efeito do reservatório marinho acima. Priores padrão foram usados ​​para todos os outros parâmetros do modelo. A mesma análise foi repetida sem ZH1540, a fim de avaliar o efeito do nosso novo genoma antigo nas datas de divergência. Repetimos ainda mais a análise usando um lognormal anterior amplo com uma média (no espaço real) de 1 × 10 −7 substitutions per site per year and standard deviation 1 on the mean clock rate and using both constant-size and exponential growth coalescent models to assess the impacts of the mean clock rate prior and demographic models on divergence time estimates.

Para cada análise, executamos quatro cadeias de Markov Monte Carlo (MCMC) de 5 × 10 8 steps each, sampling parameters and trees every 10,000 steps. After assessing convergence in Tracer v1.7123 and confirming that all four chains converged to the same posterior distribution, we combined the chains after discarding the first 10% of samples as burn-in. In the resulting combined chains all parameters have effective sample size (ESS) values > 150. TreeAnnotator v2.6.7 was used to compute MCC trees and the results were visualized using ggplot2124, ggtree125 and custom scripts. The 95% HPD of the coefficient of variation estimated under the UCLD model excluded 0 (median = 1.46, 95% HPD 1.08–1.9), indicating that a strict clock model is not appropriate for our dataset. Robustness analyses show that under a narrow mean clock rate prior both the UCED and UCLD clock models result in similar divergence time estimates (Extended Data Fig. 5a–f), with the UCED model estimates tending to be more recent and the UCLD model estimates usually having longer tails. Under a wide mean clock rate prior, estimates with the UCED are broadly similar, albeit wider, while the UCLD model estimates very wide posterior distributions for divergence times, indicating little information under this model. Divergence time estimates were not sensitive to the demographic model used. The MCC trees under the UCED model with a narrow prior, both with and without ZH1540 included in the analysis are shown in Extended Data Figs. 6 and 7, respectively.

Finalmente, realizamos um teste bayesiano de randomização de datas 126 , 127 , 128 (DRT) para avaliar melhor a força do sinal temporal em nosso conjunto de dados, permutando datas de amostragem entre genomas e realizando 50 análises replicadas. Para as análises, utilizou-se o conjunto de dados completo, um modelo de relógio UCED com um anterior estreito e o modelo demográfico Bayesiano do horizonte, fixando as datas de amostragem de cepas antigas às médias dos intervalos de datas de radiocarbono para simplificar. Cadeias MCMC foram executadas por 1 × 10 8 passos, amostragem de parâmetros a cada 10.000 passos. A convergência foi avaliada usando a coda 129 pacote para garantir que todos os parâmetros em todas as cadeias tenham valores de ESS > 150. Os resultados do DRT mostram que os intervalos HPD de 95% da taxa de clock média em réplicas com datas de amostragem permutadas são muito menores do que o esperado se todas as informações vierem da taxa de clock média anterior (Dados Estendidos Fig. 5g ). Em geral, os intervalos HPD não se sobrepõem ao intervalo HPD de 95% da taxa média de clock estimada com as datas de amostragem verdadeiras.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.