Amostragem, extração de DNA, preparação e sequenciamento de biblioteca de glicosilase de DNA de uracila parcial
Obtivemos permissão do Kunstkamera, Peter the Great Museum of Anthropology and Ethnography em São Petersburgo para a amostragem e análise de DNA antigo de 7 espécimes de dentes, escavados entre 1885 e 1892 nos cemitérios medievais de Kara-Djigach e Burana (Informações Complementares 2 ) . Não foram utilizados métodos estatísticos para predeterminar o número de amostras utilizadas neste estudo. Todos os procedimentos laboratoriais foram realizados nas instalações dedicadas de aDNA do Instituto Max Planck para a Ciência da História Humana e do Instituto Max Planck para Antropologia Evolutiva. Os procedimentos detalhados usados para amostragem de dentes podem ser encontrados na ref. 52. Resumidamente, os dentes foram seccionados na junção dentina-esmalte usando uma serra elétrica com lâmina diamantada. Após o seccionamento do dente, aproximadamente 50 mg de pó foram removidos da superfície da câmara pulpar de cada dente usando brocas odontológicas arredondadas.
O pó do dente recuperado foi utilizado para extrações de DNA utilizando um protocolo previamente estabelecido otimizado para a recuperação de fragmentos curtos de DNA 53 . As etapas exatas e as modificações do procedimento utilizado foram disponibilizadas na ref. 54 . Em resumo, o pó dental foi incubado durante a noite (12–16 h) a 37°C em 1 ml de tampão de lise de DNA contendo EDTA (0,45 M, pH 8,0) e proteinase K (0,25 mg ml- 1). Após a incubação, a ligação e o isolamento do DNA foram realizados usando um tampão de ligação à base de GuHCl personalizado e a purificação usando o Kit de Grande Volume de Ácido Nucleico Viral de Alta Pureza (Roche). Finalmente, o DNA foi eluído em 100 μl de Tris-EDTA-Tween contendo Tris-HCl (10 mM), EDTA (1 mM, pH 8,0) e Tween-20 (0,05%). Para o monitoramento do procedimento, os brancos de extração e os controles de extração positivos foram incluídos em todas as etapas de processamento do laboratório.
Todos os extratos de DNA foram convertidos em uma a duas bibliotecas de DNA de fita dupla para sequenciamento Illumina, usando 25 μl de extrato de entrada por biblioteca com uma glicosilase de DNA de uracila parcial inicial (UDG) e tratamento com endonuclease VIII (enzima USER; New England Biolabs) conforme protocolos estabelecidos 55 , 56 . O procedimento detalhado de preparação da biblioteca, incluindo as etapas de reparo da extremidade cega, ligação do adaptador e reação de preenchimento do adaptador, podem ser encontrados na ref. 57 . Após a preparação da biblioteca, cada biblioteca foi quantificada usando um sistema de PCR quantitativo (LightCycler 96 Instrument) usando os primers IS7 e IS8 55 . Para sequenciamento multiplex, realizamos dupla indexação de todas as bibliotecas usando procedimentos publicados anteriormente 58, descrito em detalhes na ref. 59 . Uma combinação de primers de índice únicos contendo identificadores de 8 pares de bases (bp) foram atribuídos a cada biblioteca. Para ajudar na eficiência da amplificação, as bibliotecas foram então divididas em várias reações de PCR para a etapa de indexação com base em sua quantificação inicial de IS7/IS8. O número de reações de PCR de indexação realizadas para cada biblioteca foi determinado de modo que cada reação recebeu uma entrada de não mais que 1,5 × 10 10cópias de DNA. Cada reação foi configurada usando a polimerase de DNA Pfu Turbo Cx Hotstart (Agilent Technologies) e foi executada por 10 ciclos usando as seguintes condições: desnaturação inicial a 95 °C por 2 min seguido por um ciclo de 95 °C por 30 s, 58 °C por 30 s e 72 °C por 1 min, bem como uma etapa final de alongamento a 72 °C por 10 min. Todos os produtos de PCR foram purificados com o MinElute DNA Purification Kit (QIAGEN), com algumas modificações no protocolo do fabricante 59 . Finalmente, todos os produtos de PCR de indexação foram quantificados por qPCR (LightCycler 96 Instrument) usando a combinação de primers IS5 e IS6 58 , 59 . Para evitar a formação de heteroduplex, as bibliotecas indexadas foram amplificadas para 10 13Cópias de DNA por reação com a Herculase II Fusion DNA Polymerase (Agilent Technologies) e quantificadas usando um instrumento 4200 Agilent TapeStation usando um sistema D1000 ScreenTape (Agilent Technologies). As bibliotecas foram diluídas para 10 nM e reunidas equimolarmente para sequenciamento. Realizamos o sequenciamento de DNA shotgun em uma plataforma Illumina HiSeq 4000 usando um kit de 76 ciclos (1 × 76 + 8 + 8 ciclos).
Processamento de leitura de sequenciamento de próxima geração e triagem metagenômica Shotgun
Após a demultiplexação, as leituras sequenciadas de shotgun brutas foram pré-processadas no pipeline EAGER v.1.92.58 usando AdapterRemoval v.2.2.0 (ref. 60 ), que foi usado para remover adaptadores Illumina (sobreposição mínima de 1 bp), bem como para filtragem de leitura de acordo com a qualidade de sequenciamento (qualidade de base mínima de 20) e comprimento (retenção de leituras ≥30 bp). Posteriormente, todos os conjuntos de dados foram selecionados para a presença de traços de DNA de patógenos usando o pipeline metagenômico HOPS 29 . Primeiro, as leituras pré-processadas foram alinhadas com um banco de dados RefSeq personalizado 61 (novembro de 2017) contendo todos os conjuntos completos de genomas bacterianos e virais, um subconjunto de conjuntos de patógenos eucarióticos e o GRCh38genoma de referência humano. As montagens do genoma que continham a palavra 'desconhecido' foram removidas do banco de dados, mantendo um total de 15.361 entradas. O banco de dados reteve várias entradas de espécies de Yersinia : Yersinia aldovae ( n = 1), Yersinia aleksiciae ( n = 1), Yersinia enterocolitica ( n = 16), Yersinia entomophaga ( n = 1), Yersinia frederiksenii ( n = 3), Yersinia intermedia ( n = 1), Yersinia kristensenii ( n = 2), Y. pestis (n = 39), fago Yersinia ( n = 17), Yersinia pseudotuberculosis ( n = 13), Yersinia rohdei ( n = 1), Yersinia ruckeri ( n = 4), Yersinia similis ( n = 1) e Yersinia sp. FDA-ARGOS ( n = 1). MALT v0.4 62foi executado usando os seguintes parâmetros: -id 90 -lcaID 90 -m BlastN -at SemiGlobal -topMalt 1 -sup 1 -mq 100 -verboseMalt 1 -memoryMode load -additionalMaltParameters. Os arquivos de alinhamento resultantes foram pós-processados com MALTExtract para uma avaliação qualitativa em relação a uma lista predefinida de 356 entradas taxonômicas de destino ( https://github.com/rhuebler/HOPS/blob/external/Resources/default_list.txt ). Especificamente, as leituras foram avaliadas de acordo com sua distância de edição em relação a uma sequência de patógeno específica no banco de dados e a possível ocorrência de incompatibilidades que poderiam significar a presença de danos no DNA 29 . Nos casos em que ambos os parâmetros foram atendidos, o alinhamento do patógeno correspondente foi considerado um forte candidato. As leituras pré-processadas foram mapeadas contra a Y. pestisCO92 (NC_003143.1) e genomas de referência humanos ( hg19 ) com o alinhador Burrows-Wheeler (BWA). Os parâmetros de mapeamento foram definidos para 0,01 para a distância de edição (-n) e o comprimento da semente foi desabilitado (-l 9999). Posteriormente, usamos SAMtools v.1.3 (ref. 63 ) para remover reads com qualidade de mapeamento inferior a 37 (para CO92) ou 30 (para hg19 ); As duplicatas de PCR foram removidas com MarkDuplicates v1.140 ( http://broadinstitute.github.io/picard/ ). Finalmente, os padrões de dano de aDNA foram avaliados com mapDamage v.2.0 (ref. 64 ).
Preparação de biblioteca de DNA de fita simples e captura de hibridização
Para os espécimes BSK001 e BSK003, bibliotecas extras de DNA de fita simples foram construídas a partir de um extrato de DNA de entrada de 30 μl. Realizamos a preparação da biblioteca no Instituto Max Planck de Antropologia Evolutiva usando um protocolo automatizado que está disponível publicamente 65 . Bibliotecas de fita simples e fita dupla de indivíduos BSK001, BSK003 e BSK007 foram enriquecidas usando sondas de DNA cobrindo todo o genoma de Y. pestis , bem como 1,24 milhão de sítios SNP de todo o genoma do genoma humano 66 , 67. Para a preparação da captura, todas as bibliotecas foram amplificadas para o número necessário de ciclos de PCR para atingir 1–2 μg de DNA de entrada. As reações de PCR foram realizadas usando a Herculase II Fusion DNA Polymerase. Eles foram então purificados usando o MinElute DNA Purification Kit e eluídos em tampão de eluição EB contendo 0,05% de Tween 20. Finalmente, as concentrações da biblioteca (ng μl −1 ) foram quantificadas usando um espectrofotômetro NanoDrop (Thermo Fisher Scientific). Para as capturas de Y. pestis em solução , o projeto do conjunto de sondas foi baseado em um conjunto de genomas modernos disponíveis publicamente, especificamente o cromossomo Y. pestis CO92 (NC_003143.1), plasmídeo CO92 pMT1 (NC_003134.1), plasmídeo CO92 pCD1 (NC_003131.1), cromossomo KIM10 (NC_004088.1), cromossomo F Pestoides (NC_009381.1) e oCromossomo Y. pseudotuberculosis IP32953 (NC_006155.1). Para as capturas de DNA humano em solução, o design do conjunto de sondas foi criado para atingir 1.237.207 variantes em todo o genoma que são informativas para estudar a história genética de populações humanas em todo o mundo 28 , 67 . As capturas de hibridização de DNA humano e Y. pestis foram realizadas por duas rodadas, conforme descrito anteriormente 28 , 69 , 68 , 67 , 66 , em que bibliotecas parcialmente tratadas com UDG do mesmo indivíduo foram agrupadas em razões equimolares para captura e fita simples bibliotecas foram capturadas separadamente.
Processamento de dados pós-captura de Y. pestis
Após a captura do genoma completo de Y. pestis , as bibliotecas foram sequenciadas em uma plataforma HiSeq 4000 (1 × 76 + 8 + 8 ciclos ou 2 × 76 + 8 + 8 ciclos) em uma profundidade de aproximadamente 11–27 milhões de leituras brutas. O pré-processamento de leituras desmultiplexadas brutas foi realizado conforme descrito na seção 'Processamento de leitura de sequenciamento de próxima geração do Shotgun e triagem metagenômica'. Nesta fase, os conjuntos de dados produzidos a partir de bibliotecas parcialmente tratadas com UDG do mesmo indivíduo foram agrupados e as bases terminais foram cortadas usando fastx_trimmer (FASTX Toolkit 0.0.14, http://hannonlab.cshl.edu/fastx_toolkit/ ) para evitar danos ao site interferência com a chamada SNP durante o processamento posterior. As etapas a seguir para mapeamento de leitura, remoção de duplicatas de PCR e cálculo de danos de aDNA foram realizadas no pipeline EAGER70 . Realizamos mapeamento de leitura com BWA v.0.7.12 contra o genoma de referência de Y. pestis CO92 (NC_003143.1). Para as bibliotecas tratadas com UDG parcial agrupadas e aparadas, os parâmetros BWA foram definidos para 0,1 para a distância de edição (-n) e o comprimento da semente foi desativado (-l 9999). Dado que as bibliotecas de fita simples construídas para este estudo retiveram danos associados ao DNA, os parâmetros BWA foram definidos para 0,01 para a distância de edição (-n) para permitir um número maior de incompatibilidades que poderiam derivar da desaminação; o comprimento da semente foi desabilitado (-l 9999). Realizamos mapeamento de leitura contra os plasmídeos usando os mesmos parâmetros contra uma referência concatenada de todos os três Y. pestisplasmídeos (pMT1: NC_003134.1; pPCP1: NC_003132.1; e pCD1: NC_003131.1), mascarando a região pPCP1 problemática entre os nucleotídeos 3000 e 4200 que mostrou ter alta similaridade com vetores de expressão usados em reagentes de laboratório 71 . SAMtools v.1.3 (ref. 63 ) foi usado para remover todas as leituras com qualidade de mapeamento inferior a 37 (-q), enquanto MarkDuplicates foi usado para remover duplicatas de PCR. Padrões de desaminação associados com dano de aDNA foram recuperados com mapDamage v.2.0 (ref. 64 ). Usamos MALT 62 para uma classificação taxonômica de leituras mapeadas, para tentar uma retenção de leituras que são mais prováveis de serem endógenas de Y. pestis. O MALT foi executado no mesmo banco de dados descrito na seção 'Processamento de leitura de sequenciamento de próxima geração do Shotgun e triagem metagenômica', usando os seguintes parâmetros: -m BlastN -at SemiGlobal -top 1 -sup 1 -mq 100 -memoryMode load -ssc -sp. O parâmetro de identidade de porcentagem mínima foi definido como padrão (-id 0,0), em oposição a um filtro de identidade de 90% usado para executar o HOPS 29 , para evitar qualquer viés de referência que possa surgir da remoção de leituras endógenas com um número maior de incompatibilidades. Após a conclusão da execução, para manter o número máximo de leituras contabilizando o algoritmo ingênuo do ancestral comum mais baixo, extraímos as leituras que foram atribuídas ao nó do gênero Yersinia ou resumidas sob o Y. pseudotuberculosisnó complexo. As leituras foram extraídas no formato FASTA de MEGAN v.6.4.12 (ref. 72 ). Posteriormente, os arquivos FASTA foram convertidos para o formato FASTQ com o script reformat.sh no BBMap da suíte BBtools (versão 38.86, https://sourceforge.net/projects/bbmap/ ). Os arquivos FASTQ foram então remapeados contra o genoma de referência CO92 usando os mesmos parâmetros descritos anteriormente nesta seção. Para bibliotecas de fita simples, mapDamage v.2.0 (ref. 64) foi usado para redimensionar as pontuações de qualidade em posições de leitura nas quais foram identificadas possíveis incompatibilidades associadas à desaminação para a referência. Posteriormente, os arquivos BAM correspondentes ao mesmo indivíduo foram concatenados após filtragem de qualidade do mapeamento e remoção de duplicatas de PCR. Realizamos a concatenação usando o comando 'merge' do SAMtools e com a ferramenta AddOrReplaceReadGroups no Picard ( http://broadinstitute.github.io/picard/ ) para atribuir um único grupo de leitura a todas as leituras em cada novo arquivo.
Chamada de SNP, estimativas de heterozigosidade e filtragem de SNP
A chamada de variantes foi realizada para BSK001 e BSK003, antes e depois da filtragem MALT 62 usando o UnifiedGenotyper no Genome Analysis Toolkit (GATK) v.3.5 (ref. 73 ). O GATK foi executado usando a opção EMIT_ALL_SITES, que produziu uma chamada para cada posição no genoma de referência cromossômico CO92. Os perfis genômicos resultantes de BSK001 e BSK003 foram comparados com um conjunto de 233 genomas modernos e 46 históricos de Y. pestis , bem como com o genoma de referência de Y. pseudotuberculosis IP32953 (NC_006155.1), usando a ferramenta Java MultiVCFAnalyzer v.0.85 ( https://github.com/alexherbig/MultiVCFAnalyzer). O MultiVCFAnalyzer v.0.85 foi executado com os seguintes parâmetros. Os SNPs foram chamados com uma cobertura mínima de três vezes e em casos de posições heterozigóticas, as chamadas foram feitas com um limite mínimo de suporte de 90%. Além disso, SNPs foram chamados com uma qualidade mínima de genotipagem de 30. Além disso, regiões não-core e repetitivas previamente definidas, bem como regiões contendo homoplasias, RNAs ribossômicos, RNAs de transferência-mensageiro e RNAs de transferência foram excluídos da chamada SNP comparativa 16 , 32 . Um conjunto de 6.567 locais variantes totais foi identificado no presente conjunto de dados.
Para investigar a extensão da possível contaminação exógena nos conjuntos de dados BSK001 e BSK003, estimamos o número de variantes heterozigóticas ambíguas além do limite de chamada do SNP. Para isso, foi utilizado o MultiVCFAnalyzer v.0.85 (ref. 74 ) para gerar uma tabela SNP de frequências alélicas alternativas variando entre 10 e 90%. Os resultados foram então usados para criar gráficos de histograma de 'heterozigosidade' das frequências estimadas em R v.3.6.1 (ref. 75 ). Gráficos de heterozigosidade foram criados antes e depois da filtragem MALT (consulte ' Processamento de dados de Y. pestis pós-captura ') para investigar se a filtragem informada por taxonomia poderia ajudar na eliminação de sequências contaminantes nos conjuntos de dados investigados (Fig. 7 suplementar ).
Uma tabela SNP criada com MultiVCFAnalyzer v.0.85, contendo todas as posições variantes no conjunto de dados presente, foi filtrada para identificar diferenças SNP entre os genomas BSK001 e BSK003. As diferenças identificadas ( n = 20) foram então avaliadas com a ferramenta Java SNP_Evaluation 30 (data de compilação 13 de agosto de 2018; https://github.com/andreasKroepelin/SNP_Evaluation ). A tabela de variantes e os arquivos VCF de cada genoma foram usados como entrada para SNP_Evaluation. Além disso, cada variante privada identificada foi avaliada dentro de uma janela de 50 pb e foi considerada 'verdadeira' ao preencher os seguintes critérios estabelecidos em estudos publicados anteriormente 17 , 21 , 30 , 76: (1) nenhum sítio multi-alélico foi permitido dentro da janela avaliada, a menos que fossem consistentes com desaminação de aDNA (significado como substituições espúrias de C-para-T ou G-para-A); (2) a própria posição do SNP avaliada não era consistente com o dano de aDNA (nenhuma base sobreposta ao SNP foi reduzida por mapDamage v.2.0 (ref. 64 )); (3) não foram identificadas lacunas na cobertura genômica na janela avaliada; (4) leituras sobrepostas aos sítios SNP mostraram especificidade para o complexo Y. pseudotuberculosis quando rastreadas com BLASTn ( https://blast.ncbi.nlm.nih.gov/Blast.cgi ).
Finalmente, para obter resolução filogenética, os conjuntos de dados BSK001 e BSK003 Y. pestis foram concatenados. Realizamos concatenação de arquivos BAM, filtragem MALT 62 e reescalonamento de danos de aDNA (com mapDamage v.2.0 (ref. 64 )) conforme descrito na seção ' Processamento de dados pós-captura de Y. pestis '. Além disso, o conjunto de dados foi incluído na análise comparativa de SNP usando MultiVCFAnalyzer v.0.85 (ref. 74 ) conforme descrito acima. Finalmente, SNPs únicos foram avaliados com SNP_Evaluation 30 de acordo com os quatro critérios listados acima.
Reconstrução filogenética e estimativas de diversidade
A análise filogenética foi usada para explorar 233 genomas de Y. pestis como parte do moderno conjunto de dados comparativos. Um alinhamento de SNP produzido por MultiVCFAnalyzer v.0.85 (ref. 74 ) foi usado para construir uma árvore filogenética em MEGA7, usando a abordagem de máxima parcimônia com 95% de deleção parcial (6.032 SNPs). Dos 233 genomas modernos de Y. pestis no conjunto de dados atual, 30 exibiram extensos comprimentos de ramos privados (Fig. 12 suplementar). Tal efeito em filogenias bacterianas pode resultar tanto da verdadeira diversidade biológica quanto de artefatos técnicos associados à falsa incorporação de SNP durante a reconstrução computacional do genoma. Embora não possamos excluir a presença de várias cepas com taxas de mutação extremamente mais altas no conjunto de dados atual, estudos anteriores mostraram que cepas modernas de Y. pestis com perfis 'mutadores' são incomuns 16 , 36. Neste estudo, 27 dos 30 genomas que mostraram disparidades em suas contagens de SNP privados em comparação com o restante do conjunto de dados foram derivados de montagens para as quais a qualidade das chamadas de SNP não pôde ser avaliada (dados brutos indisponíveis). Como possíveis erros de montagem ou chamadas de SNP falso-positivos podem afetar inferências evolutivas e estimativas de diversidade, esses genomas foram excluídos de análises adicionais. Portanto, realizamos análise filogenética usando um subconjunto de 203 genomas modernos de Y. pestis (Tabela Complementar 13 ). A lista de genomas excluídos é a seguinte: 2.MED1_139 (ref. 19 ), 2.MED1_A-1809 (ref. 18 ), 2.MED1_A-1825 (ref. 19 ), 2.MED1_A-1920 (ref. 19 ) , 2.MED0_C-627 (ref.19 ), 2.MED1_M-1484 (ref. 19 ), 2.MED1_M-519 (ref. 19 ), 0.ANT5_A-1691 (ref. 18 ), 0.ANT5_A-1836 (ref. 18 ), 0.PE2_C -678 (ref. 77 ), 0.PE2_C-370 (ref. 77 ), 0.PE2_C-700 (ref. 77 ), 0.PE2_C-746 (ref. 77 ), 0.PE2_C-535 (ref. 77 ) ), 0.PE2_C-824 (ref. 77 ), 0.PE2_C-712 (ref. 77 ), 0.PE2b_G8786 (ref. 16 ), 0.PE4_I-3446 (ref. 78 ), 0.PE4_I-3517 ( ref. 78 ), 0.PE4t_A-1815 (ref. 18 ), 0.PE4_I-3447 (ref. 78 ), 0.PE4_I-3518 (ref.78 ), 0.PE4_I-3443 (ref. 78 ), 0.PE4_I-3442 (ref. 78 ), 0.PE4_I-3519 (ref. 78 ), 0.PE4_I-3516 (ref. 78 ), 0.PE4_I -3515 (ref. 78 ), 0.PE4_Microtus91001 (ref. 79 ), 0.PE5_I-2238 (ref. 80 ) e 0.PE7b_620024 (ref. 16 ).
Um alinhamento de SNP em todo o genoma consistindo em 203 genomas modernos e 48 históricos de Y. pestis (Tabela Complementar 14 ), bem como o genoma de Y. pseudotuberculosis IP32953, foi usado como entrada para construir uma filogenia de máxima probabilidade, incluindo 2.960 SNPs e até para 4% de dados perdidos. Realizamos análise filogenética com RAxML 81 v.8.2.9 usando o modelo de substituição reversível no tempo generalizado (GTR) com 4 categorias de taxa gama. Finalmente, 1.000 réplicas de bootstrap foram usadas para estimar o suporte do nó para a topologia de árvore resultante. Após a conclusão da execução, as filogenias de máxima verossimilhança foram visualizadas com FigTree v.1.4.4 ( http://tree.bio.ed.ac.uk/software/figtree/ ) e GrapeTree (v1.5.0) 50.
Para estimar a proporção da diversidade moderna de Y. pestis descendente de BSK001/003, usamos o pacote R picante v1.8.2 82 para calcular o MPD e FPD 83 da árvore de substituição de máxima verossimilhança reconstruída. As medidas feitas em um subconjunto da árvore correspondente ao subclado descendente de BSK001/003 (ramos 1–4) foram comparadas com a filogenia completa de Y. pestis . Em ambos os casos, apenas cepas modernas foram incluídas no cálculo. Usamos uma abordagem bootstrapping para avaliar a sensibilidade de nossos resultados em relação à amostragem e incerteza filogenética 84. Para cada uma das 1.000 árvores bootstrap RAxML, reamostramos aleatoriamente as cepas modernas com reposição e apenas mantivemos os galhos da árvore correspondentes às cepas amostradas. Medidas de diversidade foram realizadas para cada uma das árvores bootstrap reamostradas obtidas, a partir das quais foram derivadas estimativas medianas e intervalos de percentis de 95%.
Para avaliar o impacto potencial da amostragem desigual entre os ramos (os ramos 1-4 continham 130 cepas modernas, enquanto o ramo 0 continha apenas 73), repetimos a mesma análise, mas acrescentando uma etapa inicial destinada a equalizar o número de genomas em ambas as partes da árvore . Subamostramos os ramos 1–4 para o mesmo número de cepas que no ramo 0 usando agrupamento de sequências nos ramos 1–4 para obter subamostras representativas. Realizamos agrupamento hierárquico com base em distâncias filogenéticas pareadas (derivadas da árvore filogenética de máxima verossimilhança) e a árvore resultante foi cortada para definir 73 agrupamentos (funções hclust 85e cutree em R v.4.0.3). Para cada árvore bootstrap, os agrupamentos foram aleatoriamente reduzidos para uma linhagem, resultando em um número igual de linhagens entre o ramo 1-4 e o ramo 0. A reamostragem com substituição foi então aplicada como anteriormente para cada uma das árvores reduzidas antes de calcular as medidas de diversidade.
Análise de plasmídeo SNP
Para investigar uma possível variação genética entre os plasmídeos de genomas históricos, realizamos o mapeamento de leitura de BSK001, BSK003 e BSK001/003 com BWA, bem como a chamada de SNP com GATK v.3.5, conforme descrito na seção acima 'chamada de SNP, estimativas de heterozigosidade e SNP filtrar' contra cada um dos três plasmídeos de Y. pestis (pMT1: NC_003134.1; pPCP1: NC_003132.1; e pCD1: NC_003131.1). Em seguida, realizamos a chamada SNP comparativa usando o MultiVCFAnalyzer v0.85 (ref. 74 ) contra um conjunto de 46 Y. pestis históricosgenomas, bem como as cepas de referência modernas CO92, KIM5 e 0.PE4-Microtus91001. As variantes foram filtradas em genomas individuais usando SNP_Evaluation de acordo com critérios previamente definidos (consulte a seção 'Chamada de SNP, estimativas de heterozigosidade e filtragem de SNP'). No presente conjunto de dados, identificamos dez variantes em pCD1, oito em pMT1 e duas em pPCP1 (Tabela Complementar 15 ).
Análise filogenética calibrada no tempo
Para estimar o tempo para a divergência dos ramos 1-4 de Y. pestis usando os genomas BSK001/003 como um novo ponto de calibração, usamos um conjunto de dados compreendendo todos os genomas modernos dos ramos 1-4 usados para análise filogenética ( n = 130), genomas da linhagem ramificada ancestral 0.ANT3 ( n = 8) e todos os 29 genomas históricos (séculos XIV-XVIII) em nosso conjunto de dados representando genótipos únicos. Em casos de genomas idênticos, o genoma de maior cobertura foi escolhido para esta análise. Aplicamos um teste de relógio molecular usando um método de máxima verossimilhança no MEGA7 (ref. 86), usando um modelo de substituição GTR no qual as diferenças nas taxas evolutivas entre os sítios foram estimadas usando uma distribuição gama discreta com quatro categorias de taxa. Com base nesse teste do relógio molecular, a hipótese nula de taxas evolutivas iguais entre os ramos filogenéticos testados foi rejeitada, o que é consistente com estudos anteriores mostrando variação da taxa de substituição entre linhagens de Y. pestis 16 , 17 . Portanto, um modelo de relógio relaxado log-normal foi usado para todas as análises de datação molecular subsequentes.
Para a análise de datação molecular, usamos o framework estatístico Bayesiano BEAST2 v.6.6 (ref. 87 ). As idades de todos os isolados antigos foram usadas como pontos de calibração para construir uma filogenia calibrada no tempo com suas faixas etárias de radiocarbono ou contexto arqueológico definidas como anteriores uniformes (consulte a Tabela Suplementar 19 para todas as faixas etárias usadas). As idades de todos os isolados modernos foram fixadas em 0 anos antes do presente. Testamos uma série de árvores coalescentes prioritárias, como o tamanho constante coalescente, horizonte Bayesiano 88 e modelos de população exponencial, todos os quais foram usados ou testados em estudos genômicos de patógenos antigos anteriores 17 , 89 , 90 , 91. Também testamos a árvore do horizonte de nascimento-morte anterior, que ganhou força nos últimos anos 91 , 92 , 93 porque pode explicar variáveis epidemiológicas e modelar disparidades de amostragem ao longo do tempo 94 . Além disso, usamos jModelTest v.2.1.10 (ref. 95 ) para identificar o modelo de substituição de melhor ajuste para nosso conjunto de dados. O modelo de transversão indicado foi implementado em BEAUti usando um modelo GTR (4 categorias de taxa gama) e o parâmetro de taxa de substituição AG fixado em 1,0 (como indicado anteriormente 93 ). Todas as árvores prioritárias foram usadas em combinação com uma taxa de clock relaxada log-normal com uma distribuição a priori uniforme variando entre 1 × 10 −3 e 1 × 10−6 substituições por sítio por ano para o alinhamento SNP (1.405 sítios após uma deleção parcial de 95%), correspondendo a um intervalo de 3 × 10 −7 a 3 × 10 −10 em todo o genoma, que está dentro do intervalo de estimativas 17 . Como parte da configuração da topologia filogenética, todos os genomas dos ramos 1-4 (antigos e modernos), bem como a linhagem 0.ANT3, foram restringidos a serem clados monofiléticos independentes. Para o tamanho da população constante e os priors da árvore de população exponencial, todos os outros parâmetros foram definidos como padrão. Para a árvore do horizonte coalescente anterior, uma distribuição anterior de Jeffreys (1/ x) foi usado para os tamanhos de população e uma dimensão de 5 foi usada para permitir variações nos tamanhos de grupo e população ao longo do tempo, com um limite superior de 380.000 para o tamanho efetivo da população (padrão). Além disso, para a árvore do horizonte de nascimento-morte anterior, usamos uma prévia uniforme para a taxa de se tornar não infecciosa que variou entre 0,03 e 70, para considerar possíveis períodos infecciosos variando de 30 anos (infecções ao longo da vida em reservatórios de roedores 96 , 97 ) a 5 dias (período infeccioso médio para peste bubônica 98 ). Usamos uma distribuição beta anterior com média = 0,1 (alfa = 10,0, beta = 90,0) para a probabilidade de amostragem ρ no tempo 0 e uma distribuição uniforme variando entre 0 e 0,1 para a proporção de amostragem s. Para este último, foram permitidos dois turnos ao longo do tempo. Finalmente, permitiu-se que o número reprodutivo R variasse entre 0 e 4,0 usando uma longa distribuição normal anterior de mediana = 1,0 e dp = 0,7, que está dentro da faixa de estimativas anteriores para peste bubônica e pneumônica durante epidemias medievais 98 .
A adequação de todas as árvores prioritárias foi avaliada usando amostragem de caminho conforme implementado no pacote de seleção de modelo do BEAST2 v.6.6. A amostragem de caminho foi executada em 50 etapas, com 20 milhões de estados como o comprimento da cadeia para cada etapa. As probabilidades log-marginais resultantes favoreceram com 'forte suporte' 99 o modelo coalescente skyline para a presente análise (fator log Bayes = 8,35 quando comparado com o segundo melhor modelo) (Tabela Suplementar 20 ). Portanto, o modelo de skyline coalescente foi escolhido para análise posterior. Para avaliar o sinal temporal no presente conjunto de dados, usamos TempEst v.1.5.3 para estimar a distância raiz-ponta em relação às idades dos espécimes em uma análise de regressão linear 100 . Para TempEst, usamos uma árvore de máxima parcimônia calculada em MEGA7 (ref.86 ) no formato NEXUS. Além disso, usamos o ponto médio dos intervalos de datas arqueológicas ou de radiocarbono para todos os genomas antigos como datas de ponta. Todas as idades do genoma moderno foram definidas para 0 anos antes do presente. Os valores resultantes do coeficiente de correlação r (0,39) e R 2 (0,16) apoiaram a existência de um sinal temporal no presente conjunto de dados. Além disso, usamos a abordagem BETS 101para uma avaliação de sinal temporal que leva em consideração todos os parâmetros de análise. O BETS compara as estimativas de verossimilhança marginal (log) produzidas a partir de um modelo isócrono (todas as datas de amostragem definidas como 0 anos antes do presente) com um modelo heterócrono (incluindo tempos de amostragem reais). Como anteriormente, a amostragem de caminho foi executada em 50 etapas com 20 milhões de estados como o comprimento da cadeia para cada etapa. O fator (log)-Bayes estimado de 129,33 deu forte suporte ao modelo heterócrono; portanto, indicou a presença de um sinal temporal no presente conjunto de dados.
Para a análise de datação molecular usando uma configuração de modelo de horizonte coalescente, realizamos amostragem Monte Carlo de cadeia de Markov usando 2 cadeias independentes de 300 a 400 milhões de estados cada. Após a conclusão, as execuções foram combinadas usando LogCombiner v.2.6.7 e a convergência foi avaliada usando Tracer v.1.6 ( http://tree.bio.ed.ac.uk/software/tracer/ ) garantindo que os tamanhos de amostra efetivos fossem maiores de 200 para cada distribuição posterior estimada após um burn-in de 10%. Árvores de credibilidade máxima de clade foram construídas usando TreeAnnotator no pacote BEAST2 v.6.6 87com um burn-in de 10% e foram então visualizados no FigTree v.1.4.4. Em paralelo com a análise de datação molecular, realizamos uma amostragem da análise anterior para testar um possível overfitting da anterior aos dados. Realizamos amostragem Monte Carlo de cadeia de Markov para 2 cadeias independentes de 600 milhões de estados cada. Após a conclusão da corrida, as corridas foram combinadas e a convergência foi avaliada após um burn-in de 30%. Os resultados indicam que as distribuições posteriores do relógio relaxado log-normal não correlacionado e o tempo para as estimativas de ancestral comum mais recente não são concordantes com aqueles obtidos ao usar uma análise informada por dados (Fig. 13 suplementar ).
Como a maioria das estruturas filogenéticas Bayesianas (como BEAST2) são baseadas em árvores bifurcadas e, portanto, são pobres em resolver nós multifurcando, complementamos nossa abordagem usando TreeTime v.0.8.4 (ref. 35 ) para inferir uma filogenia calibrada no tempo usando um abordagem de máxima verossimilhança. O TreeTime demonstrou resolver politomias de maneira consistente com as datas de ponta da amostra. Geramos uma filogenia de máxima verossimilhança enraizada usando RAxML (Fig. 10 suplementar ) a partir do mesmo alinhamento de SNP usado para BEAST2 (95% de deleção parcial). A árvore de máxima verossimilhança foi então usada como entrada para TreeTime, que foi executado usando todas as datas de amostragem conhecidas para genomas modernos e o ponto médio da faixa etária para os genomas antigos (Tabela Suplementar 22). TreeTime foi executado usando a árvore coalescente Kingman antes com a configuração do horizonte. Um modelo de substituição apropriado foi escolhido para os dados usando a opção de inferência -gtr. A filogenia em escala de tempo foi inferida usando um relógio relaxado não correlacionado e com as opções de otimização de comprimento de ramificação, manter raiz e manter politomia. Além disso, os intervalos de tempo de divergência foram estimados a partir da árvore de maior verossimilhança usando a opção -confiança. As análises foram executadas usando um número máximo de 500 e 1.000 iterações (opção de número máximo de iterações) e produziram resultados consistentes. A árvore de tempo resultante pode ser encontrada na Fig. 11 Complementar .
Resumo do relatório
Mais informações sobre o projeto de pesquisa estão disponíveis no Nature Research Reporting Summary vinculado a este artigo.
