Śmiertelne epidemie dżumy wśród łowców-zbieraczy z jeziora Bajkał 5500 lat temu / Ruairidh Macleod, Frederik V. Seersholm, Bianca De Sanctis, Angela Lieverse, Adrian Timpson, Rick Schulting, Jesper T. Stenderup, Charleen Gaunitz ,Lasse Vinner, Olga Iwanowna Goriunova, Włodzimierz Iwanowicz Bazalijski, Siergiej W. Wasiljew, Erin Jessup, Yucheng Wang, Christopher Bronk Ramsey, Mark G. Thomas, Russell Corbett-Detig, Astrid KN Iversen, Andrzej W. Weber ,Martin Sikora i Eske Willerslev

Analiza genomów starożytnych patogenów znacząco poszerzyła naszą wiedzę na temat historii ewolucji ludzkich chorób zakaźnych (na przykład Salmonella enterica ⁴ i wirusowe zapalenie wątroby typu B ⁵ ), chociaż miało to miejsce głównie w kontekście społeczności rolniczych lub pasterskich. Y. pestis , czynnik etiologiczny dżumy, jest prawdopodobnie najlepiej zbadany pod tym względem i miał niszczycielskie konsekwencje dla populacji ludzkich przez tysiąclecia. Historyczne wybuchy dżumy odpowiadają za niektóre z najbardziej tragicznych wydarzeń w historii ludzkości ¹ . Odzyskanie starożytnego DNA od ofiar dżumy dostarczyło niezwykłych spostrzeżeń na temat pochodzenia i ewolucji dżumy w czasie tych wydarzeń ^6 , 7 , i, co niezwykłe, ujawniło zakażenia u prehistorycznych osobników w całej Europie ⁸ . Historycznie i obecnie dżuma jest związana z transmisją za pośrednictwem pcheł od gryzoni, które z powodzeniem przystosowały się do ludzkiej niszy komensalnej w neolicie ⁹ . Analiza genomiczna prehistorycznej dżumy wskazuje, że w przypadku wczesnych, różnicujących się szczepów, kluczowe adaptacje genetyczne niezbędne do transmisji choroby za pośrednictwem pcheł oraz zakażenia dymieniczego nie występowały ^2 , 3 , co prowadzi do niepewności co do drogi transmisji i ciężkości tych szczepów.

Wykrycie wczesnych przypadków dżumy u wielu pokoleń rolników z późnego neolitu zostało wykorzystane do powiązania wybuchów choroby z przedłużającym się spadkiem demograficznym między około 5300–4900 kalibrowanymi latami przed teraźniejszością (cal  bp ) ^10 , 11, chociaż alternatywne wyjaśnienie przypisuje spadek kryzysowi rolniczemu12 ^, 13. Pierwsza interpretacja była kontrowersyjna, a inni sugerowali, że infekcje bardziej przypominają łagodne zapalenie jelit przenoszone przez żywność14. Podobieństwo lub brak podobieństwa tych wczesnych szczepów do Y. pseudotuberculosis — najbliższego krewnego Y. pestis — było ważnym punktem zainteresowania podczas takich dyskusji, a na podstawie istniejących starożytnych genomów szacuje się , że Y. pestis oddzielił się od Y. pseudotuberculosis w pewnym momencie w ciągu ostatnich 50 000 lat (odnośniki ^{8 , 11 , 15} ).

Badania prehistorycznych genomów dżumy ze szczepów z późnego neolitu i epoki brązu (LNBA) datowane są przeważnie na 4700–2400  pz (odnośniki ^{3 , 8 , 16} ) i są zazwyczaj definiowane jako jedna z dwóch linii, w zależności od obecności (LNBA+) lub braku (LNBA−) genu ymt3 . ymt koduje toksynę myszy Yersinia , która zwiększa przeżywalność bakterii w przewodzie pokarmowym pcheł w okresie przejściowym między gryzoniami a żywicielami ludzkimi, a tym samym przenoszoną przez ukąszenia pcheł dymieniczą postać dżumy u ludzi¹⁷ . Linie Y. pestis , które rozeszły się przed tymi kladami LNBA, zidentyfikowano również u garstki neolitycznych osobników szwedzkich (5200–4850  pb kal .) ^10 , 11 i łotewskiego osobnika o zachodnich przodkach łowców-zbieraczy (5300–5050  pb kal. ) ¹⁵ . Tym genomomom brakuje klasycznych genów wirulencji ( profaga YpfΦ i ymt ), chociaż analiza pangenomiczna wykazała obecność locus kodującego mitogen pochodzący z Y. pseudotuberculosis (YPM), superantygenową toksynę związaną z Y. pseudotuberculosis (ale nie późniejszymi szczepami Y. pestis ). Rodzi to intrygujące pytania o możliwą ciężkość wczesnych szczepów dżumy; późniejsze szczepy LNBA− wykazują znaczną utratę genów, chociaż ich potencjał wirulencji jest nieznany ³ . W tych badaniach jak dotąd brakowało dowodów dotyczących wpływu demograficznego zakażeń dżumą na populacje prehistoryczne.

Środkowoholoceńscy łowcy-zbieracze z okolic jeziora Bajkal, w południowo-wschodniej Syberii, stali się przedmiotem intensywnych badań archeologicznych prowadzonych przez Baikal Archaeology Project, dostarczając ważnych zbiorów danych do określania prehistorycznych sposobów życia łowców-zbieraczy18 ^, 19. Grupy te wykazują niezwykłą ciągłość sposobu życia i utrzymania łowców-zbieraczy, o czym świadczy obszerny zapis archeologiczny miejsc pochówku z okresu 8500–3500 lat  p.n.e. (ref. ²⁰ ). Genomy badanych łowców-zbieraczy wskazują na długoterminową ciągłość starożytnego północnego Eurazji i północno-wschodniej Azji aż do ok. 4500–4000 lat p.n.e. ( ref. ^21 , 22 ) (rozszerzone dane, rys. 1 i 2 ). W tym okresie przypadki dżumy ze szczątków ludzkich odpowiadających szczepowi LNBA− są dokumentowane sporadycznie wśród pochówków z wczesnej epoki brązu22 ^{, 23.} Zdarzenia przenoszenia się chorób odzwierzęcych powodujące zakażenia dżumą w tym regionie do dziś stanowią poważny problem zdrowotny24 ^. Są one głównie związane ze świstakami, głównym odzwierzęcym rezerwuarem dżumy w regionie25 ^, 26.  Aby zbadać zdrowie i strukturę społeczności w prehistorycznych grupach łowców-zbieraczy, przeanalizowaliśmy starożytne DNA ludzi i patogenów z czterech miejsc cmentarzysk w Cis-Bajkale (zachodni i północny region jeziora) w dwóch oddzielnych epidemiach datowanych na 5520–5265 kJ  i 5315–4235 kJ  ( 95,4% przedziały ufności dla modelowanych zakresów dat na podstawie osób z wykrytymi przypadkami dżumy, skorygowane o efekty zbiorników słodkowodnych; Uwaga dodatkowa  4 ). Długi ogon drugiego zakresu dat wybuchu epidemii (rys. 1b ) wynika z faktu, że obejmuje on tylko dwie bezpośrednie daty, chociaż najbardziej prawdopodobny zakres dat wynosi około  5050–4850 p.n.e.

Obraz w pełnym rozmiarze

a , Lokalizacje dotkniętych cmentarzy nad rzeką Angara na północny zachód od jeziora Bajkał oraz współdzielenie IBD między objętymi badaniem mieszkańcami cmentarzy (parzyste linie podziału między stanowiskami; szara skala wskazuje całkowite współdzielenie IBD w segmentach większych niż 3 centymorgany (cM) łącznie ponad 10 cM w relacjach parami między stanowiskami; Notatka uzupełniająca  2 ). Wstawka przedstawia objęte badaniem osoby w każdym stanowisku (31 Ust’-Ida I, 8 Bratskij Kamień, 2 Serowo i 5 Szumilcha), ze wskazanymi wykryciami dżumy i wspólnymi grobami wskazanymi zacienionymi obszarami wokół osób. Mapy utworzone przy użyciu Natural Earth Data. b , Oszacowania gęstości jądra naniesione w modelach bayesowskich zakresów dat dla wczesnej (czerwonej) i późnej (ciemnożółtej) fazy wybuchów dżumy nad Bajkałem. Jaśniejsze zacieniowane obszary odpowiadają zsumowanym rozkładom prawdopodobieństwa przed modelowaniem; wykorzystane daty pochodzą od osób zidentyfikowanych wyłącznie z dżumą. Wstawka przedstawia 95,4% przedziały ufności dla modelowanych zakresów dat z Bajkału w porównaniu z tymi z innych przypadków dżumy sprzed epidemii LNBA: RV2039 z Łotwy ^[15] , Warberg 1 i Warberg 2 z Niemiec ^[61] oraz FRA102 ze Szwecji ^[11] . c . Oszacowanie gęstości jądra na podstawie modelowanych rozkładów datowania radiowęglowego dla czterech cmentarzysk, dla wszystkich ludzi w wieku po odstawieniu od piersi (lub związanych z nimi zawieszek z zębami jelenia) datowanych radiowęglowo, niezależnie od pobranego DNA. Modelowanie datowania radiowęglowego przeprowadzono za pomocą OxCal w wersji 4.4.4 [ ^65 , 66] .

Wygenerowaliśmy sekwencjonowane metodą shotgun starożytne DNA od 46 osób z późnego neolitu i przeanalizowaliśmy te dane pod kątem obecności patogenów ( Metody ). Wykazało to wyjątkowo wysokie występowanie Y. pestis wśród tych osób, bardziej niż jakiegokolwiek innego patogenu. Y. pestis wykryto u 18 osób, co wskazuje na 2 odrębne fazy wybuchów zakażeń dżumą — oddzielone od 4 do 6 wieków — na 4 cmentarzyskach (Ryc. 1a,b ). Występują one w dwóch fazach w Szumilichie, Ust’-Idzie I, Brackim Kamieniu i Serowie (patrz Ryc. 1c ), przy czym przypadki z Brackiego Kamienia występują zarówno w pierwszej, jak i drugiej fazie. Wszystkie te miejsca znajdują się na brzegach rzeki Angara, głównego cieku wodnego wypływającego z jeziora Bajkał, bogatego w ^{rybołówstwo27} . Dane dotyczące stabilnych izotopów węgla i azotu u osobników z Ust’-Idy I dowodzą spożycia zarówno lokalnych ryb, jak i zwierzyny lądowej ²⁷ . Pochówki w Ust’-Idzie I i Szumilichej odpowiadają tradycji pogrzebowej z Isakowa (charakteryzującej się ciałami zazwyczaj zorientowanymi równolegle do rzeki oraz obecnością darów grobowych, takich jak gliniane naczynia w kształcie mitry, kamienne groty strzał oraz kościane lub porożne groty), podczas gdy te z Brackiego Kamienia i Serowa odpowiadają tradycji pogrzebowej z Serowa (z ciałami często zorientowanymi prostopadle do rzeki; głównymi cechami darów grobowych są bifasady i jajowate garnki; patrz również Notatka uzupełniająca  1 ). W Ust’-Idzie I wykryliśmy również odczyty powiązane z patogenem odzwierzęcym Brucella , przyczyną brucelozy, u jednego osobnika (#26.04; Notatka uzupełniająca  3 ). Dwie epidemie dżumy pogrupowano według dominujących praktyk pogrzebowych na każdym cmentarzysku: groby w stylu Isakovo w pierwszej epidemii i groby w stylu Serovo w drugiej fazie (Uwaga uzupełniająca  1 ), które są współczesne nad jeziorem Bajkał między około 6000–5000 kal bp (ref. ²⁷ ). Okres ten jest definiowany lokalnie jako późny neolit, zgodnie z syberyjską terminologią archeologiczną, gdzie neolit ​​jest definiowany na podstawie kryteriów technologicznych, takich jak wprowadzenie łuku i strzał, naczyń glinianych i technik mielenia kamienia (nie ma roślin i zwierząt domowych innych niż psy), chociaż społeczności te pozostają łowcami-zbieraczami aż do wkroczenia pasterstwa w późnej epoce brązu. Stanowiska grobowe stanowią zdecydowaną większość zapisu archeologicznego Cis-Bajkału, a określenia takie jak późny neolit ​​są późniejszymi kategoriami, stosowanymi do odrębnych zestawów cech pochówku i darów grobowych, które ogólnie odpowiadają różnym okresom. Wszystkie cztery cmentarzyska były również wykorzystywane we wczesnym neolicie (7650–6660 r.  p.n.e. ) i wczesnej epoce brązu (4970–3470 r. p.n.e.) bp ) ²⁷ , chociaż tutaj brane są pod uwagę tylko ich późnoneolityczne komponenty.

Dzielenie się segmentami tożsamości przez pochodzenie (IBD) między osobnikami na tych cmentarzyskach wskazuje na niedawne wspólne pochodzenie. Chociaż odległość między nimi sięgała 340 km, rzeka Angara z pewnością ułatwiała podróżowanie. Stwierdzono bardzo niskie wskaźniki chowu wsobnego, a wysoką efektywną wielkość populacji, opartą na seriach homozygotyczności, wywnioskowano za pomocą hapROH (estymacja metodą maksymalnego prawdopodobieństwa: 18 219 osobników, 95% przedział ufności 9445–42 062). Jest to zgodne ze scenariuszem wysoce mobilnych, egzogamicznych grup łowców-zbieraczy.

Wśród analizowanych tu osobników łowców-zbieraczy, najwięcej wykrytych zakażeń dżumą odnotowano w Ust’-Idzie I, będącym jednocześnie największym cmentarzyskiem w Isakowie na obszarze Cis-Bajkału. W tym przypadku wskaźnik detekcji wyniósł 35% (11 z 31 zsekwencjonowanych osobników), w tym pochówki nr 14 i 56.01, dla których dane dotyczące genomu ludzkiego były wcześniej raportowane ^[21] . W innych miejscach zidentyfikowaliśmy jeden genom dżumy o wysokim pokryciu w Szumilichie, cztery genomy o niższym pokryciu z Brackiego Kamienia oraz jeden genom o średnim pokryciu z Serowa. Ogólnie rzecz biorąc, obserwujemy 39% wskaźnik detekcji wśród osobników z późnego neolitu na tych cmentarzyskach (z cementu zębowego). Dla porównania, ilościowe badanie PCR znanych ofiar zarazy w średniowieczu w Smithfield w Londynie w Wielkiej Brytanii28 wykazało wskaźnik detekcji na poziomie 5,7% w przypadku kości i 37% w przypadku miazgi zębowej (łącznie 20%), co wskazuje na wysoki odsetek fałszywie ujemnych wyników wykrywania zarazy przy użyciu starożytnego DNA. Aby zapobiec błędnej interpretacji danych, wszystkie starożytne osoby, których dane przesiewowe pochodzą z dotkniętych miejsc, zostały tutaj uwzględnione (pokrycie ludzkiego genomu autosomalnego waha się od 0,001× do 1,9×, średnio 0,65×). Bezpośrednie datowanie radiowęglowe uzyskano dla prawie wszystkich osób w późnoneolitycznych komponentach tych cmentarzysk (łącznie 58, wliczając te wcześniej zgłoszone z Ust’-Idy I29 ^; Dane uzupełniające  7 ).

Genomy Y. pestis zidentyfikowane pomiędzy dwiema fazami epidemii nad Bajkałem okazały się rozbieżne przodkowo do obecnie znanego kladu starożytnych i współczesnych szczepów dżumy (Rys. 2 ). Z pewnością przypisujemy je do Y. pestis na podstawie ich pozycji filogenetycznej, a także obecności genów wirulencji i plazmidów charakterystycznych dla Y. pestis (Rozszerzone dane Rys. 3 i Notatka uzupełniająca  3 ). Ta filogeneza została zbudowana przy użyciu genomów uzyskanych z pochówku Shumilikha nr 34 (pokrycie 6,4×) z pierwszej fazy oraz z pochówku Bratskii Kamen nr 22 (1,6×) i pochówku Serovo nr 10 (1,0×) z drugiej fazy. Osiem genomów o niższym pokryciu zostało umieszczonych filogenetycznie przy użyciu UShER ³⁰. Algorytm UShER znajduje najbardziej oszczędne umiejscowienie na drzewie, wybierając węzeł z największą liczbą potomków, jeśli wiele z nich jest równie oszczędnych, i ignoruje brakujące genotypy. Umieszczenie wszystkich genomów o niskim pokryciu w tym samym węźle bazalnym wynika częściowo z braków danych, chociaż jest zgodne z pozycją trzech genomów Bajkału o wyższym pokryciu. Wnioskowanie bayesowskie dotyczące dat węzłów przeprowadzono zgodnie z podejściem uwzględniającym efekty rekombinacji w filogenezach bakterii ^31 , 32 ( Metody i Notatka uzupełniająca  3 ). Pojawienie się Y. pestis jako gatunku klonalnego Y. pseudotuberculosis następuje w pewnym czasie między rozejściem się linii, która dała początek Y. pestis (oznaczonej jako węzeł A na ryc. 2 ) a najnowszym wspólnym przodkiem dostępnych genomów Y. pestis (węzeł B na ryc. 2 ). Górna granica podana przez pierwszą prawdopodobnie będzie znacząco zmieniona przez niedobór sekwencji genomów Y. pseudotuberculosis i mogłaby być nowsza, gdyby zidentyfikowano filogenetycznie bliższe serowary. Niemniej jednak ta dolna granica (ze średnią datą 5709 lat temu) rewiduje poprzednie szacunki rozbieżności na 4810–5122 lat temu33 ^, czego można by się spodziewać po uwzględnieniu genomów Y. pestis starszych niż ten zakres (inne szacunki wahały się od 6000 do 50 000 lat temu⁸ i 7400 lat temu15 ⁾ . Filogeneza wspiera wniosek, że Y. pestis pierwotnie wyewoluował z wariantu szczepu O:1 Y. pseudotuberculosis (reprezentowanego przez genom z serotypu O:1c, akcesja Europejskiego Archiwum Nukleotydów (ENA): SAMEA7160327), co jest zgodne z wcześniejszymi ustaleniami³⁴donosząc o inaktywacji klastra genów antygenu O jako kroku w kierunku ewolucji Y. pestis . Pomiędzy tymi dwiema fazami obserwujemy niewielkie różnice genetyczne między szczepami w odrębnych mutacjach prywatnych w szczepach pierwszej i drugiej fazy (z rygorystycznymi filtrami do określania genotypów; Metody i Notatka uzupełniająca  3 ); jest to również wyraźnie widoczne w położeniu węzłów na ryc. 2. Chociaż wiadomo, że wskaźniki mutacji w Y. pestis są bardzo zmienne w różnych liniach ³³ , wynik ten jest zgodny ze scenariuszem pokrewnych szczepów powstałych w wyniku oddzielnych zdarzeń odzwierzęcych z lokalnego rezerwuaru zwierzęcego.

Obraz w pełnym rozmiarze

Po prawej, ogólna topologia kompleksu gatunków Y. pseudotuberculosis jest pokazana na drzewie z adnotacjami o mutacjach, opartym na 448 genomach (długości gałęzi wskazują odległość mutacji). Wstawka, uproszczona wersja tej filogenezy, ze szczególnym uwzględnieniem prehistorycznych szczepów dżumy (niektóre długości gałęzi zostały skrócone). Trzy próbki z jeziora Bajkalskiego o większym pokryciu zostały bezpośrednio włączone do konstrukcji filogenezy RAxML, podczas gdy osiem próbek o mniejszym pokryciu zostało filogenetycznie umieszczonych później, a wszystkie dzieliły najbardziej oszczędne położenie w najbardziej podstawowym węźle Y. pestis . Ich położenie w tym węźle nie stanowi pozycji rozgałęzienia, dlatego ich włączenie jest w sąsiedztwie tego węzła. Lewy górny róg, daty węzłów wewnętrznych oszacowane za pomocą BactDating. BRK, Bratskii Kamen; SER, Serovo; SHU, Shumilikha; UID, Ust’-Ida I.

Aby umieścić te wybuchy dżumy w odpowiednim kontekście, wzięliśmy pod uwagę wzorce pokrewieństwa biologicznego, obróbkę pochówku i wiek w chwili zgonu na dotkniętych cmentarzyskach łowców-zbieraczy. W miejscu z najwyższym pozytywnym wykryciem dżumy (i największą próbką), Ust’-Ida I, datowania radiowęglowe dla późnoneolitycznego komponentu Isakovo są wyjątkowo ściśle skupione jak na stosunkowo duży cmentarz ²⁹ (Rozszerzone dane Rys. 4 ). Modelowane zakresy dat dla wszystkich ofiar dżumy we wczesnej fazie wskazują na bardzo wąski przedział czasowy, rzędu kilku dekad (Uwaga uzupełniająca  4 ), co wspiera scenariusz, że pochówki te były równoczesne. Potwierdza to dodatkowo duże podobieństwo genomów dżumy, co jest zgodne z zakażeniami dżumą występującymi podczas pojedynczego wybuchu lub w bardzo krótkim okresie czasu. Rekonstruując najbardziej prawdopodobne rodzinne rodowody, stwierdzamy, że relacje i wiek członków rodziny są zgodne ze zdarzeniem śmiertelnym w okresie krótszym niż jedno pokolenie (Rys. 3 ). Żadne z par wieku w chwili śmierci nie wskazuje na przykład na dzieci, które osiągnęły podobny wiek co rodzice, ani na rodzeństwo przyrodnie i przyrodnie w bardzo różnym wieku (największa różnica wieku między rodzeństwem wynosi dziewięć lat, a dzieli je rodzeństwo w środku). W przypadku obecności kilku pokoleń, ich wnioskowane zakresy wieku w chwili śmierci są generalnie zgodne z przewidywanymi, gdyby wszyscy krewni zmarli w tym samym czasie (na przykład 12–15-latek ma ojca w wieku 35–50 lat).

Obraz w pełnym rozmiarze

Osoby, u których wykryto dżumę, oznaczono sylwetkami pałeczek. Rodowody sporządzono na podstawie 30 osób pobranych w próbkach z Ust’-Idy I i 8 z Brackiego Kamienia; przedstawiono tylko bliskie pokrewieństwa rodzinne; chociaż w Ust’-Idzie I wykryto wiele innych pokrewieństwa trzeciego lub czwartego stopnia (patrz Notatka uzupełniająca  2 ). Rodowody rekonstruowano na podstawie zgonów, które nastąpiły w tym samym czasie. Dwie próbki z Serowa (niepokazane) okazały się krewnymi czwartego stopnia. yo, lat.

Grupa grobowa Isakowo na cmentarzysku Ust’-Ida I jest niezwykła z kilku innych powodów wśród cmentarzysk łowców-zbieraczy w Cis-Bajkalsku. Oprócz ściśle skupionych dat radiowęglowych, śmiertelność dzieci jest nieproporcjonalnie wysoka (obserwowana również na Brackim Kamieniu, patrz ryc. 4 ), a także występuje wysoka częstość grobów z wieloma pochówkami (ponad połowa na tym stanowisku) bez śladów późniejszego otwierania grobów i dodawania nowych pochówków. Sugeruje to współwystępowanie zgonów we wspólnych grobach, co jest zgodne z katastrofalnym zdarzeniem śmiertelnym. Za pomocą KIN ³⁵ wykryto również pokrewieństwo wujkowskie między cmentarzyskami Ust’-Ida I i Szumilcha, ale nie zostało to potwierdzone oczekiwanym wzorcem wspólnego IBD (Uwaga uzupełniająca  2 ). Mimo to wysoki stopień dzielenia się IBD między poszczególnymi osobami (ryc. 1 ) na odcinku zaledwie 37 km wzdłuż rzeki Angara sugeruje, że jednoczesne wybuchy dżumy mogą być powiązane z bliskim kontaktem grup w tym czasie.

Obraz w pełnym rozmiarze

Wykres gęstości jądra modelowanych prawdopodobieństw wieku w chwili śmierci, oparty na modelu zerowym ciągłego prawdopodobieństwa zgonu w dowolnym wieku (patrz Notatka uzupełniająca  5 ). Pokazano wszystkie istotne zespoły szczątków kostnych człowieka z regionu Cis-Bajkał badanego przez AL, obejmujące ponad 20 osób. Rozmiary próby dla przedstawionych stanowisk: Bracki Kamień (późny neolit ​​(LN)), n  = 20; Chużyr-Nuge XIV, n  = 81; Lokomotiw, n  = 101; Szamanka II, n  = 156; Szumilicha, n  = 36 (wczesna epoka brązu); Ust’-Ida I, n  = 48; Wiercholensk, n  = 27. Osoby uwzględniono tylko z dominującego okresu użytkowania grobów w każdym stanowisku. Wykluczono osoby odstające od normy z późniejszych lub wcześniejszych tradycji grobowych.

Jeśli chodzi o wykrywanie dżumy w grupach grobowych, nie znajdujemy statystycznie istotnego wzorca współwystępowania dżumy wśród krewnych (Notatka uzupełniająca  3 ), chociaż osoby dotknięte chorobą wydają się być powiązane w wielu przypadkach. Pochówek w Brackim Kamieniu obejmuje wspólny grób 3 młodych dziewcząt w wieku od 4 do 9 lat (Ryc. 3 , po lewej), z podobnymi datami radiowęglowymi (Notatka uzupełniająca  4 ). Dwie z nich (#19.01 i #19.03) zostały uznane za spokrewnione trzeciego stopnia (najpewniej kuzynki); trzecia miała niewystarczające zachowanie DNA, aby z całą pewnością wnioskować o pokrewieństwie, ale wszystkie trzy miały wspólny haplotyp mitochondrialny z trzema rzadkimi mutacjami prywatnymi i dlatego prawdopodobnie były bliskimi krewnymi ze strony matki. Dane genomu dla Y. pestis zidentyfikowano u wszystkich trzech, co sugeruje wybuch zakażenia dżumą w rodzinie, z synchronicznymi zgonami trojga dzieci. Podobnie w Ust’-Idzie I siostrzeniec i ciotka (nr 20.01 i nr 20.02) zostali pochowani we wspólnym grobie, przy czym Y. pestis zidentyfikowano u obu (ryc. 3 , pomarańczowy rodowód). Nastoletnia siostrzenica ciotki została jednak pochowana w innym, wspólnym grobie z niespokrewnionym nastoletnim mężczyzną (co może sugerować pokrewieństwo niebiologiczne); jego ojciec z kolei (zielony rodowód) został pochowany w zupełnie oddzielnym grobie.

Ponadto niektóre pary rodzeństwa pochowane razem we wspólnych grobach wykazują tylko jedną osobę wykrytą jako dodatnią pod kątem dżumy, jak ma to miejsce w przypadku rodzeństwa w grobie nr 25 (ryc. 3 , czerwony rodowód). W innym przykładzie, w przypadku siostry (#26.01) i brata (#26.04), siostra jest uznawana za dodatnią, podczas gdy brat nie (chociaż odczyty Y. pestis są wykrywane tuż poniżej progu pewnej identyfikacji; Dane uzupełniające ). Obserwacje te są zgodne z wysokim wskaźnikiem fałszywie ujemnych wyników w analizie paleogenomicznej dżumy²⁸. Brat był również zakażony prawdopodobnie nieśmiertelną brucelozą (Uwaga uzupełniająca  3 ). W kilku przypadkach bliscy członkowie rodziny znajdują się w różnych grobach na cmentarzu, na przykład w grobie nr 8 trzecie rodzeństwo pary w grobie 26. Widoczny jest wzór, w którym dwóch blisko spokrewnionych członków rodziny jest pochowanych razem, a jeden lub więcej innych jest pochowanych dalej. Może to być zgodne z bardziej rozciągniętą sekwencją zgonów zamiast pojedynczego zdarzenia śmiertelnego, jeśli wspólne groby wskazują na współistniejące zgony, odzwierciedlając scenariusz opóźnionego przeniesienia choroby z człowieka na człowieka. Nie stwierdzono żadnych innych przyczyn zgonu poza genetycznie wykrytą infekcją dżumą (chociaż inne wykryte drobnoustroje mogą wskazywać na koinfekcje bakteryjne w momencie zgonu; uwaga uzupełniająca  3 ). Co istotne, musieli istnieć ocaleni, aby pochować zmarłych, a typowe dla Isakowa traktowanie i wyposażenie grobów, a także uznanie pokrewieństwa biologicznego sugerują bardziej rozciągniętą sekwencję zdarzeń śmiertelności.

Na Bajkale głównym współczesnym odzwierzęcym rezerwuarem dżumy jest świstak syberyjski ( Marmota sibirica ), a polowanie na świstaki dla mięsa i futra historycznie skutkowało przewlekłymi zakażeniami dżumą, zwłaszcza u młodych mężczyzn, którzy są narażeni podczas skórowania i rzeźnictwa ³⁶ . Od XIX wieku świstaki były najczęściej wybieranym gatunkiem łownym przez rdzennych myśliwych w tym regionie, początkowo poprzez łapanie w pułapki ³⁷ , i istnieją obszerne historyczne opisy „zarazy tarbagan” w wyniku spożycia zakażonych świstaków w okolicach jeziora Bajkał ³⁸ . Zaopatrzenie się w świstaki przez prehistorycznych łowców-zbieraczy jest wyraźnie udokumentowane obecnością licznych zębów świstaków jako darów grobowych w grobach Kitoi z wczesnego neolitu ^19 , 39 , chociaż nie zostały one znalezione w grobach z późnego neolitu. Spożycie surowych lub niedogotowanych organów świstaka skutkuje zakażeniem posocznicowym drogą fekalno-oralną, podczas gdy bliski kontakt ze świstakami zakażonymi współczesnymi szczepami Y. pestis wywołuje dżumę dymieniczą lub płucną (lub często obie), przy czym ta druga często występuje wtórnie do zakażenia posocznicowego40 lub wdychania zakaźnych kropelek krwi podczas, na przykład, obdzierania ze skory⁴¹ . Opisana powyżej częstość wykrytych zakażeń wśród wspólnie pochowanych krewnych byłaby zgodna z przenoszeniem dżumy między ludźmi, szczególnie poprzez transmisję płucną w przypadku równoczesnych zgonów.

Uderzającym aspektem danych osteologicznych dotyczących wieku w chwili śmierci na cmentarzach Ust’-Ida I i Bratskij Kamień — dwóch cmentarzach, na których wykryto wiele przypadków dżumy — jest to, że ich profile demograficzne są silnie przechylone w kierunku śmiertelności dzieci. Oba wykazują szczyt śmiertelności w przedziale wiekowym 7,5–11 lat — to znaczy u dzieci przed okresem dojrzewania (Uwaga uzupełniająca  1 ). W analizie profili śmiertelności na cmentarzach łowców-zbieraczy w cis-bajkalskim regionie z połowy holocenu, te dwa cmentarze wyraźnie odstają od normy pod względem proporcji zgonów dzieci (Rys. 4 ). Wynik ten okazał się wysoce istotny statystycznie, biorąc pod uwagę zerowy model profili śmiertelności (Uwaga uzupełniająca  5 ). Z kolei przedział wiekowy 20–25 lat wykazuje najniższą śmiertelność na Ust’-Idzie I, a zgony w wieku 20–35 lat są całkowicie nieobecne na Bratskim Kamieniu (Uwaga uzupełniająca 4 ). Rodzice są również wyraźnie nieobecni w grupach rodowodowych; chociaż istnieje wiele relacji między rodzeństwem i kuzynostwem, istnieje tylko jeden przypadek relacji rodzic–potomstwo. Jednak stosunek płci u tych osobników wydaje się niezmieniony (22 XY i 24 XX).

W kontekście powszechnego zakażenia szczepem dżumy o nieznanej zjadliwości, tę różnicę w śmiertelności między dziećmi a dorosłymi można interpretować na wiele sposobów, biorąc pod uwagę dostępne dane bioarcheologiczne i obecną wiedzę na temat odporności człowieka. Po pierwsze, dorośli mogliby w dużej mierze składać się z osób, które były już narażone na dżumę i wyzdrowiały z niej w dzieciństwie, a tym samym nabyły odporność ochronną, zapobiegającą ponownemu zakażeniu lub śmierci. Oznaczałoby to regularne nawracanie epidemii, czego nie jesteśmy w stanie potwierdzić w naszych badaniach, a także sugerowałoby, że osoby starsze z większym prawdopodobieństwem nabyły odporność, jednak śmiertelność w rzeczywistości nieznacznie wzrasta po 20–35 roku życia (po pierwotnym szczycie zachorowań około 10. roku życia). Alternatywnie, nie można wykluczyć zmienności śmiertelności spowodowanej różnicami w zachowaniu między grupami wiekowymi (na przykład podział zadań grupowych lub ról według wieku, skutkujący większym narażeniem dzieciństwa na świstaki), chociaż istnieje niewiele analogicznych precedensów w odniesieniu do świstaków konkretnie, a nie potwierdza tego brak podwyższonej śmiertelności dzieci na żadnych innych cmentarzyskach łowców-zbieraczy na Bajkalsku (ryc. 4 ). Wreszcie możliwe jest, że dzieci mogą być bardziej narażone na śmierć ze względu na wrodzone różnice w odpowiedziach immunologicznych między dorosłymi i dziećmi przed okresem dojrzewania. Wiadomo, że dzieci są bardziej podatne na zakażenia bakteriami Gram-ujemnymi42 ^, o czym świadczy dzisiejszy profil epidemiologiczny zakażeń Yersinia enterocolitica i Y. pseudotuberculosis43.

Ewolucja linii Y. pestis jest w znacznym stopniu kształtowana przez procesy utraty genów44 ^, wzorzec typowy dla bakterii patogennych w procesie przejściowym do obligatoryjnego pasożytnictwa, który został również zidentyfikowany w szczepach dżumy LNBA− ³ . Z analizy zasięgu klasycznych genów wirulencji dżumy odkryliśmy, że geny wirulencji nieobecne w opublikowanych szczepach LNBA− i pre-LNBA z Riņņukalns (RV 2039) i Falbygden ¹¹ są również nieobecne w Bajkale ( ymt i profag YpfΦ ; ryc. 5 ), uniemożliwiając wystąpienie dżumy dymieniczej (rozszerzone dane, ryc. 3 ). Jednakże, ponieważ ta analiza genów wirulencji opiera się na tradycyjnym mapowaniu pojedynczego odniesienia, jest ograniczona do zawartości genetycznej obecnej w nowoczesnym odniesieniu. Aby scharakteryzować możliwą zmienność przodków Y. pseudotuberculosis w szczepach Cis-Baikal, która może przyczynić się do naszej interpretacji ich patogenności, zmapowaliśmy odczyty sekwencjonowania na wykres zmienności pan-genomowej reprezentujący różnorodność genetyczną w 82 kompletnych zestawach gatunków kompleksu Y. pseudotuberculosis (56 Y. pestis, 24 Y. pseudotuberculosis i 1 Yersinia similis , w oparciu o ref. ¹¹ ). Odkryliśmy, że dwa szczepy dżumy z jeziora Bajkał miały podobne poziomy różnorodności przodków Yersinia , które stwierdzono tylko w Y. pseudotuberculosis i Y. similis, jak w innych szczepach pre-LNBA (rys. 5e ). Na przykład wykryliśmy obecność ypm, genu kodującego superantygen YPM znanego ze współczesnych szczepów Y. pseudotuberculosis45 ⁱ ostatnio zaobserwowanego w szczepach pre-LNBA i LNBA−dżumy¹¹. W nowoczesnym Y. pseudotuberculosis występują trzy allele tego genu : ypmA , ypmB i ypmC , przy czym ypmA jest uważany za najbardziej zjadliwą formę genu46 ^.

Obraz w pełnym rozmiarze

a , Polimorfizmy pojedynczego nukleotydu (SNP) w genie ypm (brakujące genotypy pokazane jako puste). Ref, reference. b , Pozycje SNP (czerwone strzałki) w strukturze białkowej YPMa (Protein Data Bank: 1PM4 ). Rysunek wykonano przy użyciu danych z bazy danych AlphaFold, numer akcesyjny A0A0U1QV71 . c , Niezakorzenione drzewo sąsiadujących genów z sekwencji wariantów ypm , w tym opisanej tu formy dżumy sprzed LNBA. d , Obecność otwartych ramek odczytu wokół locus ypm . e , Porównanie zawartości genów przodków między prehistorycznymi szczepami dżumy na podstawie wyrównań grafów referencyjnych, z wykorzystaniem znormalizowanej szerokości pokrycia.

YPM wiąże się z regionem niezmiennym cząsteczek antygenu leukocytarnego człowieka (HLA) klasy II i oddziałuje z domeną zmienną łańcucha β receptora limfocytów T. Poprzez połączenie HLA klasy II z receptorem limfocytów T, YPM promuje aktywację limfocytów T i uwalnianie szeregu cytokin prozapalnych, dodatkowo wzmacniając odpowiedź immunologiczną47 ^, 48 .

Sugerowano, że te odpowiedzi immunologiczne związane z YPM są przyczyną różnych powikłań zapalnych, w tym encefalopatii, szkarlatyny dalekowschodniej (FESLF; znanej również jako gorączka Izumi; szczególnie związanej z ypmA ) i zespołu podobnego do Kawasaki ^{47 , 49 , 50 , 51} . Obecnie FESLF występuje głównie u dzieci w wieku poniżej 14 lat, a choroba Kawasaki występuje głównie u dzieci w wieku 5 lat lub młodszych. Jednak choroba Kawasaki po zakażeniu Y. pseudotuberculosis może również dotyczyć starszych dzieci ⁵² . Te powikłania zapalne związane z YPM prawdopodobnie również dotyczyły głównie dzieci przed okresem dojrzewania, dodatkowo zaostrzając wczesną zachorowalność i śmiertelność wywołaną przez Y. pestis u młodych.

Co godne uwagi, odkryliśmy, że gen ypm z dwóch szczepów dżumy z Cis-Bajkal jest najbardziej podobny pod względem sekwencji do ypmA , różniąc się jedynie w 3 pozycjach zasad: 4653 (T>G (izoleucyna>arginina)), 4711 (C>T (synonim treonina)) i 4770 (G>A (glicyna>glutaminian); rys. 5a ). Te trzy pojedyncze polimorfizmy nukleotydu (SNP) wydają się być ustalone we wszystkich szczepach dżumy, w których obecny jest ten gen (szczepy pre-LNBA i LNBA−; rys. 5a ). Ponieważ dwa z tych trzech wariantów są mutacjami niesynonimicznymi (I54R i G93E znajdują się w odrębnych beta-harmonijkach struktury YPM), mogą potencjalnie wpływać na drugorzędową, trzeciorzędową lub czwartorzędową strukturę białka, interakcje białko-białko i rozpoznawanie przez układ odpornościowy (rys. 5b ). Ponadto, rekonstruując najbardziej prawdopodobną filogenezę trzech znanych wariantów ypm wraz z naszymi danymi, odkryliśmy, że wersja genu ypmB jest wysoce rozbieżna z ypmA , ypmC  i starożytnej dżumy ypm . Zakładając, że korzeń drzewa znajduje się między ypmB  a pozostałą różnorodnością, ypm ze starożytnej dżumy wydaje się rozbieżny przodkowo zarówno do ypmA, jak i ypmC (rys. 5c ).

Ponadto zidentyfikowaliśmy dziesięć otwartych ramek odczytu (ORF) wokół locus ypm , które są obecne w formie przodków dżumy, ale są nieobecne w formach późniejszych. Region ORF jest podobny do niestabilnego regionu genomu Y. pseudotuberculosis o wyraźnie niskiej zawartości GC. Odkryliśmy, że te ORF w genomach Baikal Y. pestis były podobne do tych otaczających wariant ypmB , jak zgłoszono w innych szczepach przed LNBA ¹¹ (rys. 5d ). Ten wzór z genem ypm podobnym do ypmA w połączeniu z locus ypm podobnym do ypmB — według naszej wiedzy — nie był wcześniej obserwowany. Możliwym wyjaśnieniem może być to, że ta różnorodność w szczepach jeziora Bajkał przed LNBA, najbardziej bazowych ze zbadanych szczepów dżumy, odzwierciedla trwającą lokalną adaptację do świstaków i innych żywicieli gryzoni w większym stopniu niż do ludzi, ponieważ regionalny rezerwuar żywicieli zwierzęcych prawdopodobnie w tym czasie znacznie przewyższał rezerwuar ludzi. Ta jak dotąd unikalna kombinacja może mieć wpływ na przykład na metylację genów i poziom transkrypcji ypm .

Te cechy genetyczne szczepów Baikal pre-LNBA, wraz z zależnymi od wieku różnicami w układzie odpornościowym, mogą częściowo wyjaśnić, dlaczego dzieci przed okresem dojrzewania stanowią większość wśród ofiar dżumy, chociaż ocena ich rzeczywistego wpływu wymaga badań funkcjonalnych.

Nasze odkrycia pokazują, że najwcześniejsze znane wybuchy dżumy miały miejsce u prehistorycznych łowców-zbieraczy na wieki przed zaobserwowaniem zakażeń u neolitycznych rolników. Wybuchy te były prawdopodobnie wynikiem odzwierzęcego przeniesienia zarazy z populacji dzikich świstaków nad jeziorem Bajkał. Wyniki te potwierdzają centralno- lub północno-wschodnioazjatyckie pochodzenie dżumy, podczas gdy wcześniej najwcześniejsze próbki odnotowano tylko w północnej Europie11 ^, 15. Jest to zgodne z szacunkami opartymi na analizie różnorodności współczesnego Y. pestis53. Nasza analiza filogenetyczna ujawnia, że ​​te zjadliwe szczepy dżumy są czasowo stosunkowo bliskie ostatniemu wspólnemu przodkowi Y. pestis i Y. pseudotuberculosis , co może wskazywać na szybką dywersyfikację wraz z przeniesieniem na żywicieli gryzoni z jednego lub kilku innych żywicieli zwierzęcych Y. pseudotuberculosis . Dodatkowo, podnosi to pytania dotyczące różnicowania taksonów w obrębie kompleksu gatunków Y. pseudotuberculosis (który obejmuje Y. pestis i Y. similis ), na które same dane dotyczące starożytnego genomu mogą nie być wystarczające, aby odpowiedzieć (biorąc pod uwagę, że konwencjonalne rozróżnienia mogą być również oparte na potencjale patogenicznym i zakresie żywicieli). Co więcej, wnioskowana wysoka śmiertelność tych epidemii jest bezpośrednio potwierdzona profilami śmiertelności i zbieżnymi datami radiowęglowymi w dotkniętych miejscach pochówku, co wskazuje, że dzieci i młodzież były szczególnie narażone; są to spostrzeżenia, których wcześniej brakowało w przypadku prehistorycznych zakażeń dżumą. Do tej pory najwcześniej wykryte szczepy dżumy charakteryzowały się niepewną patogennością; ich wirulencja była przedmiotem znacznej debaty, opartej wyłącznie na danych genetycznych ^{3 , 10 , 15 , 16 , 54} . W niniejszym artykule zintegrowaliśmy liczne dowody z miejsc pochówku dotkniętych dżumą (w tym genomy dżumy, wzorce pokrewieństwa biologicznego, profile śmiertelności i modelowane zakresy datowania radiowęglowego), aby scharakteryzować to, co naszym zdaniem jest śmiertelnymi konsekwencjami infekcji w tym okresie.

Kontekst tych ognisk jest ważny dla interpretacji zdrowia i epidemiologii w przeszłości. Fakt, że ogniska te występują w stosunkowo małych, mobilnych prehistorycznych grupach łowców-zbieraczy, podkreśla, że ​​zwiększona gęstość zaludnienia, udomowienie zwierząt i zmiany stylu życia wynikające z przejścia neolitycznego nie są koniecznymi warunkami dla znaczących ognisk odzwierzęcych. To dodatkowo rewiduje interpretacje dżumy jako unikalnego czynnika przyczyniającego się do spadku demograficznego w późnym neolicie w Europie, jak wcześniej sugerowano10 ^, 11, zwłaszcza biorąc pod uwagę widoczną dotkliwość ognisk zidentyfikowanych tutaj. Profil śmiertelności w Bajkale również ostro kontrastuje z oczekiwaniami zwolenników teorii neolitycznej transformacji epidemiologicznej, że największy ciężar chorób odzwierzęcych wśród prehistorycznych łowców-zbieraczy spadnie na producentów (20-40-latków) ⁵⁵ . Nasze odkrycia dodatkowo ujawniają wgląd w społeczny wymiar tych społeczności podczas epidemii, świadcząc o trosce o zmarłych (poprzez współpochówek bliskich krewnych i najwyraźniej współczesnych ofiar) oraz o równoczesnym wykrywaniu zakażeń wśród krewnych, którzy mieli kontakt interpersonalny za życia. Te dowody transmisji międzyludzkiej kontrastują z wcześniejszymi oczekiwaniami dotyczącymi podstawowych szczepów dżumy¹⁵.

Warto zauważyć, że dzieci wydają się ponosić główny ciężar śmiertelności z powodu zakażeń dżumą w Cis-Bajkale. W historycznych zapisach epidemii dżumy zaobserwowano różne wskaźniki śmiertelności w różnych grupach wiekowych. Podobnie do naszych ustaleń, zapisy parafialne z epidemii dżumy dymieniczej w Londynie (Wielka Brytania) w 1603 r. wykazały znacznie wyższy wskaźnik śmiertelności dzieci (około 5 razy wyższy) ^56 , 57. Jedną wyraźną różnicą między epidemiami w Cis-Bajkale a średniowieczną epidemią dżumy dymieniczej jest prawdopodobna droga transmisji. Chociaż w obu przypadkach mogła wystąpić transmisja drogą powietrzną i fekalno-oralną, transmisja przez ukąszenia pcheł związana z dżumą dymieniczą jest mało prawdopodobna w przypadku epidemii w Cis-Bajkale (biorąc pod uwagę brak ymt ). Rozprzestrzenianie się zakaźnych kropelek lub aerozoli poprzez kaszel jest udokumentowane jako główny sposób transmisji dżumy płucnej58 ^, co jest zgodne z naszymi ustaleniami dotyczącymi rozprzestrzeniania się z człowieka na człowieka, wywnioskowanymi na podstawie pokrewieństwa biologicznego i danych archeologicznych. Co ciekawe, nasze wyniki są zgodne z wcześniejszymi interpretacjami, zgodnie z którymi wczesne szczepy Y. pestis mogły stać się śmiertelnymi patogenami układu oddechowego ⁵⁴ .

Nasze wyniki wskazują, że najwcześniejszy zaobserwowany rozlew odzwierzęcy nie był zdarzeniem jednorazowym, lecz powtórzył się kilka wieków później, podkreślając znaczenie zakażeń odzwierzęcych w prehistorycznych społecznościach w wielu różnych środowiskach kulturowych i środowiskowych. Ponadto, niski zasięg identyfikacji brucelozy sugeruje dowody na transmisję odzwierzęcą ze zwierząt na ludzi w tych grupach (zakażenie następuje poprzez bezpośredni kontakt z zakażonymi zwierzętami59 , Notatka uzupełniająca  3 ). Powtarzające się wybuchy epidemii przodków rozbieżnych szczepów dżumy w grupach z regionu Cis-Bajkal w wieku 5500–5000 lat  ^p.n.e. dodatkowo sugerują długą historię dzikich gryzoni jako wieloletniego rezerwuaru rozlewu dżumy. Kolejne szczepy dżumy genetycznie najbliższe tym z regionu Cis-Bajkal pochodzą z około 5000 km na zachód, z północnej Europy. Biorąc pod uwagę tę odległość i fakt, że istnieje niewiele dowodów na zewnętrzny kontakt z grupami niebędącymi myśliwymi-zbieraczami w tym czasie, wspiera to hipotezę, że znaczny, obejmujący cały kontynent rezerwuar gryzoni Y. pestis może odpowiadać za częste, izolowane zdarzenia rozprzestrzeniania się w kolejnych tysiącleciach, zamiast ciągłej transmisji z człowieka na człowieka. Co więcej, potencjalny związek prehistorycznego rozprzestrzeniania się z zaopatrzeniem przez ludzi świstaków w Cis-Bajkale podkreśla prawdopodobną kluczową rolę innych gatunków gryzoni w składzie rezerwuarów Y. pestis . Około 352 gatunki rezerwuarowe zostały zidentyfikowane na podstawie współczesnych obserwacji, z których wiele jest ekologicznie od dawna zadomowionych (na przykład wiewiórki naziemne i myszoskoczki) ⁶⁰ . Scenariusz trwałego prehistorycznego rezerwuaru jest zgodny zarówno z wcześniejszymi ustaleniami szybkich, powtarzających się zakażeń rozbieżnymi szczepami dżumy w obrębie tej samej linii rodzinnej z neolitycznej Szwecji ¹¹ , jak i w Cis-Bajkale. Zauważamy, że ostatnie odkrycia wskazują, iż wczesne przypadki dżumy u zwierząt udomowionych ^61 , 62 mogą być zgodne z odwróconą zoonozą, chociaż biorąc pod uwagę długą historię gatunków gryzoni uznawanych za rezerwuary dżumy ⁶⁰, uważamy, że obecnie jest to najbardziej oszczędne źródło tych epidemii.

Łącznie nasze odkrycia podkreślają powszechność zakażeń odzwierzęcych, biorąc pod uwagę znacząco odmienny sposób życia prehistorycznych łowców-zbieraczy i europejskich rolników z epoki neolitu. Te spostrzeżenia są równie istotne dla wyzwań stojących przed dzisiejszym światem, jak 5500 lat temu, kiedy 75% nowych ludzkich patogenów pojawia się w wyniku transmisji od zwierząt ^[63]. Wgląd w historię ewolucji patogenów w okresach znaczących zmian demograficznych i technologicznych (w tym przypadku przed wpływem neolitu) może dostarczyć danych pozwalających na kontekstualizację głównych wyzwań, przed którymi stoi obecnie ludzkość, takich jak spowodowane zmianami klimatu zaburzenie nisz ekologicznych na całym świecie ^[64] .

Praca laboratoryjna

Starożytne DNA wyekstrahowano z cementu zębowego zębów trzonowych lub przedtrzonowych z archeologicznych pozostałości szkieletowych badanych w ramach Baikal Archaeology Project. Pobieranie próbek starożytnego DNA (aDNA) przeprowadzono w dedykowanych, czystych laboratoriach w Centrum GeoGenetyki Fundacji Lundbecka (Kopenhaga) oraz w Instytucie Archeologii University College London (Londyn). Cement wyizolowano specyficznie z korzeni zębów ⁶⁷ za pomocą wysterylizowanej ręcznej piły obrotowej, a następnie sproszkowano przed demineralizacją i trawieniem enzymatycznym. Pobrane alikwoty stanowiły około 50–100 mg materiału. Ekstrakcję, oczyszczanie i przygotowanie biblioteki aDNA do sekwencjonowania metodą shotgun przeprowadzono zgodnie z podejściem opisanym w pracy Allentoft i in. ⁶⁸ , stosując protokół biblioteki dwuniciowej według Margaryan i in. ⁶⁹ w pierwszym przypadku oraz protokół biblioteki jednoniciowej „Santa Cruz Reaction” ⁷⁰ dla próbek o niskiej zawartości matrycowego DNA. Stężenia uzyskanych bibliotek uzyskano za pomocą analizatora FragmentAnalyzer firmy Agilent i połączono w stężeniach ekwimolarnych do sekwencjonowania na przepływowych komórkach Illumina NovaSeq 6000 S4 (odczyty parzyste 100 pb) w GeoGenetics Sequencing Core (Kopenhaga). Wszystkie próbki przesiewano bez częściowego traktowania glikozylazą uracyl-DNA (UDG), a w niektórych przypadkach kolejne biblioteki zbudowano z traktowaniem UDG (zgodnie z ⁷¹ ).

W przypadku bibliotek, w których sekwencjonowanie przesiewowe wykazało obecność DNA Y. pestis, przeprowadzono wzbogacanie w roztworze. Wychwyt hybrydyzacyjny przeprowadzono za pomocą zestawu Arbor Sciences myBaits zgodnie z protokołem producenta o wysokiej czułości (High Sensitivity), ale tylko z jedną rundą wzbogacania. Połączone biblioteki z reakcji wychwytu poddano następnie ponownej amplifikacji przez 16 cykli i sekwencjonowaniu na tej samej platformie, co powyżej ^.

Wstępna bioinformatyka

Po przeprowadzeniu analizy baz danych Illumina za pomocą programu CASAVA (wersja 1.8.2) ⁷³ , sekwencje adapterowe i ogony poliN zostały przycięte z demultipleksowanych plików fastq za pomocą programu AdapterRemoval (wersja 2.0). Odczyty wyrównano z ludzkim genomem referencyjnym GRCh38 za pomocą programu bwa aln (wersja 0.7.18) ⁷⁴ (genom referencyjny hg19 został również użyty do analizy hapRoH, patrz poniżej). Wyrównane odczyty zostały przekonwertowane do plików BAM, scalone między bibliotekami na poziomie próbki, posortowane, przefiltrowane i zindeksowane za pomocą narzędzia Samtools (wersja 1.21) ⁷⁵ , a następnie duplikaty zidentyfikowano za pomocą narzędzia MarkDuplicates z programu Picard (wersja 2.18.7) z następującymi opcjami: „OPTICAL_DUPLICATE_PIXEL_DISTANCE = 12000 REMOVE_DUPLICATES = false TAGGING_POLICY = All VALIDATION_STRINGENCY = LENIENT”. Duplikaty odczytów zostały następnie odfiltrowane za pomocą narzędzia Samtools wraz z odczytami o jakości mapowania <30. Statystyki podsumowujące głębokość i zasięg sekwencjonowania wygenerowano za pomocą narzędzi BEDtools (wersja 2.23.0) ⁷⁶ i pysam ( https://github.com/pysam-developers/pysam ). Oszacowanie wzorców skażenia i uszkodzeń DNA człowieka przeprowadzono na poziomie biblioteki, przy użyciu programów contamMix ⁷⁷ , ANGSD (v.0.940) ⁷⁸ i mapDamage2.0 ⁷⁹ .

Analiza ludzkiego DNA

Płeć chromosomową wnioskowano na podstawie stosunku odczytów odpowiadających chromosomom Y i X, zgodnie z istniejącymi przedziałami ufności ^[80] . Nie wykryto aneuploidii chromosomowych. Haplogrupy mitochondrialne przypisano za pomocą haplogrepu (wersja 2.4.0) ^[81] po adnotacji wariantów za pomocą mutserve (wersja 1.3.0) ^[82 ]. Haplogrupy chromosomów Y przypisano zgodnie z podejściem opisanym w pozycji ¹¹ .

Do analizy eksploracyjnej pochodzenia poprzez analizę głównych składowych (PCA), pseudohaploidalne genotypy zostały wywołane poprzez losowy wybór wariantu z pileupu wygenerowanego przez Samtools. Próbki zostały następnie rzutowane na przestrzeń zmienności uzyskaną za pomocą smartpca ⁸³ do przeprowadzenia PCA na 2 086 279 SNP (filtrowanych tylko pod kątem transwersji i z niewielką częstością alleli >0,1%) z panelu referencyjnego starożytnych populacji euroazjatyckich ⁶⁸ . Ten ostatni został przeniesiony z hg19 za pomocą hgLiftOver ( https://genome.ucsc.edu/cgi-bin/hgLiftOver ), a efekty przeniesienia oceniono (Extended Data Fig. 2 ). Wszystkie zgłoszone próbki zostały włączone do PCA poprzez projekcję.

Genotypy diploidalne wywołano przy użyciu narzędzia bcftools (v.1.21), a w celu analizy współdzielenia segmentów IBD brakujące genotypy diploidalne wprowadzono przy użyciu narzędzia GLIMPSE ⁸⁴ (dla próbek z minimalnym pokryciem genomu autosomalnego 0,1×) zgodnie z podejściem opisanym w pozycji ⁶⁸, a segmenty IBD wywołano przy użyciu narzędzia IBDseq (v.r1206) ⁸⁵ , a następnie przeprowadzono klasterowanie genetyczne przy użyciu narzędzia IBD ⁶⁸ .

Przebiegi homozygotyczności wykryto z segmentów homozygotycznych ze względu na pochodzenie uzyskanych z IBDseq oraz z podzbioru danych pseudohaploidalnych do 1240 tys. pozycji SNP przy użyciu hapRoH (v.1) ⁸⁶. Pokrewieństwo biologiczne wnioskowano przy użyciu KIN ³⁵, początkowo uruchamiając KINgaroo na filtrowanych plikach BAM ukierunkowanych na 2 086 279 SNP opisanych powyżej, a następnie walidując na podstawie współdzielenia IBD wywnioskowanego z IBDseq. Następnie zrekonstruowano rodowody, biorąc pod uwagę uzyskane szacunki logarytmu prawdopodobieństwa dla scenariuszy pokrewieństwa, haplotypów jednorodzicielskich, płci i wieku w chwili śmierci (Uwaga uzupełniająca  2 ).

Badanie przesiewowe w kierunku taksonów patogenów

Dane sekwencjonowania shotgun wygenerowane z cementu zębowego zostały przebadane pod kątem obecności znanych ludzkich patogenów za pomocą procesu pracy Pathopipe ( https://github.com/martinsikora/pathopipe/ ), opisanego szczegółowo w pracy Sikora i in. ⁸⁷ . Odczyty klasyfikowano za pomocą szybkiego podejścia k-merowego, KrakenUniq ⁸⁸ (wersja 0.5.8), w oparciu o niestandardową bazę danych ludzkich patogenów i mikroorganizmów środowiskowych. Dla każdego rodzaju zidentyfikowanego w każdej próbce wykonano porównania parami za pomocą Bowtie2 ⁸⁹ (wersja 2.5.4) dla wszystkich odczytów sklasyfikowanych do tego rodzaju, porównując je ze wszystkimi dostępnymi genomami referencyjnymi gatunków dla tego samego rodzaju. Następnie przyporządkowanie obecności taksonów patogenów odbywa się na podstawie następujących progów wykrywania: liczba unikalnych odczytów >30, ranga k-merów = 1, skorygowany współczynnik pokrycia >0,5, średnia identyczność nukleotydów >0,97, średnia liczba miękkich przyciętych zasad <8, w oparciu o progi zastosowane przez Seersholm et al. ¹¹ .

Analiza DNA Y. pestis

W przypadku próbek, w których zakażenie dżumą zidentyfikowano za pomocą procesu przesiewowego patogenów, przeprowadziliśmy tradycyjne mapowanie pojedynczego odniesienia za pomocą Bowtie2⁸⁹ na referencyjnym genomie dżumowym (CO92; GCA_000009065.1 ), z parametrami ’-D 20 -R 3 -N 1 -L 20 -i S,1,0.50–end-to-end–no-unal’. Następnie duplikaty odczytów i odczyty o niskiej jakości mapowania (MQ < 30) zostały usunięte za pomocą samtools ⁷⁵, po czym obliczono średnią głębokość pokrycia w każdej próbce za pomocą BEDtools genomecov⁷⁶. Scharakteryzowaliśmy przypadki dżumowe na podstawie ich pokrycia jako: wstępne wykrycia (<0,01×), przypadki dżumowe o niższym pokryciu (0,01-1×) lub genomy dżumowe o wyższym pokryciu (>1×).

Dla trzech genomów o wyższym pokryciu wywołaliśmy genotypy w sposób próbkowy, używając HaplotypeCaller z GATK⁹⁰, a następnie przeprowadziliśmy kolejny krok wspólnego wywołania haplotypów przy użyciu GenotypeGVCF w scalonym zbiorze danych. Używając VariantFiltration (GATK) usunęliśmy wywołania o niskim poziomie ufności: niskiej jakości genotypu (<50), równowagi alleli mniejszej niż 0,9, głębokości odczytu mniejszej niż 3 lub głębokości odczytu większej niż 1000. Następnie przekonwertowaliśmy zbiór danych do formatu multifasta, używając konsensusu bcftools. Robiąc to, zastosowaliśmy maskę w regionach zawierających najwyższy odsetek odczytów o zerowej jakości mapowania, zazwyczaj zlokalizowanych w regionach powtarzalnych, zgodnie z podejściem opisanym w Seersholm i in. ¹¹ . Aby zachować współrzędne genomu referencyjnego i zapewnić spójność z innymi dopasowanymi próbkami, dane sekwencjonowania dla opublikowanych genomów referencyjnych dla Y. pseudotuberculosis i Y. similis pobrano z ENA, a odczyty dopasowano do genomu referencyjnego Y. pestis GCA_000009065. Następnie genotypy wywołano i przefiltrowano zgodnie z powyższym opisem, uzyskując wielokrotne dopasowanie sekwencji dla całego chromosomu Yersinia z 448 sekwencjami na współrzędnych genomu referencyjnego. Drzewo filogenetyczne wywnioskowano z pełnego pliku dopasowania obejmującego wszystkie sekwencje referencyjne i trzy próbki o wysokim pokryciu z tego badania, wykorzystując RAxML-NG ⁹¹ z modelem substytucji GTR + G i używając genomu referencyjnego Y. similis (SAMEA5779183) jako grupy zewnętrznej (patrz rozszerzone dane, rys. 3 ). Przekonwertowaliśmy wielosekwencyjne wyrównanie do haploidalnego modelu VCF za pomocą programu faToVcf³⁰, używając genomu referencyjnego Y. pestis NC_003143.1.fa jako odniesienia, a następnie zbudowaliśmy obiekt drzewa z adnotacjami mutacji na podstawie tego modelu VCF i filogenezy RAxML za pomocą programu UShER ^30. Ta nowa filogeneza UShER zachowuje oryginalną topologię, ale bezpośrednio przypisuje podstawienia w modelu VCF do gałęzi na drzewie za pomocą algorytmu Fitcha-Sankoffa^92 , 93, dzięki czemu długości krawędzi są wyrażone w jednostkach rzeczywistych podstawień. Użyliśmy programu matUtils ⁹⁴ do wyodrębnienia pliku .json z protobuf drzewa z adnotacjami mutacji, dostępnego interaktywnie online pod adresem https://bit.ly/Ypestis_MAT .

Dla 8 próbek o niższym pokryciu wywołaliśmy vcf tylko dla SNP za pomocą narzędzia bcftools, filtrując dla minimalnej jakości odwzorowania 30 i baz, minimalnej jakości odwzorowania 30 baz i maksymalnej głębokości 1000. Zachowaliśmy tylko miejsca, które były zmienne w panelu referencyjnym lub w których więcej niż jedna próbka o niższym pokryciu zawierała wywołany wariant, a następnie usunęliśmy wszystkie miejsca wariantów w masce niskiej jakości odwzorowania opisanej powyżej, uzyskując przefiltrowany vcf dla próbek o niższym pokryciu. Ten vcf o niskim pokryciu został użyty do filogenetycznego umieszczenia próbek o niskim pokryciu w drzewie z adnotacjami mutacji za pomocą UShER. Wszystkie 8 próbek o niższym pokryciu miało jedno, maksymalnie oszczędne położenie w węźle korzeniowym kladu Y. pestis .

Na koniec uruchomiliśmy Gubbins (v.3.4.3) ³¹ na multifasta wszystkich genomów LNBA− i pre-LNBA, a także genomu serowaru O:1c (SAMEA7160327) i genomu SAMN03121000 jako grupy zewnętrznej (patrz Supplemental Fig. 14 ). Aby określić datę filogenezy, biorąc pod uwagę rekombinację, uruchomiliśmy BactDating ³² na drzewie wyjściowym z Gubbins w oparciu o zmienność nierekombinacyjną. Użyliśmy 100 000 iteracji i zrelaksowanego modelu gamma, jak zasugerowano w ref. ^33. Konwergencja została potwierdzona za pomocą pliku śledzenia, jak pokazano w Supplement. Przedstawiamy szacunki wieku mediany i 95% przedziały ufności dla węzłów będących przedmiotem zainteresowania na rys. 2 .

Analiza wykresu wariacji

Aby scharakteryzować pełną różnorodność starożytnej dżumy, zbudowaliśmy pangenomowy graf zmienności całej znanej różnorodności w obrębie kompleksu gatunków Y. pseudotuberculosis ( Y. pestis , Y. pseudotuberculosis i Y. similis ). Użyliśmy Pangenome Graph Builder (pggb) ⁹⁵ na wszystkich dostępnych zespołach kompleksu Y. pseudotuberculosis z NCBI z poziomem zespołu scharakteryzowanym jako „chromosomowy” lub „kompletny”. Aby zapewnić poprawną konstrukcję grafu wokół plazmidów, zbudowaliśmy oddzielne grafy dla chromosomu i plazmidów, a następnie połączyliśmy je za pomocą narzędzi vg ⁹⁶ . Następnie zindeksowaliśmy graf i przeprowadziliśmy mapowanie krótkich odczytów Giraffe ⁹⁷ na grafie zmienności danych z tego badania i wszystkich publicznie dostępnych starożytnych danych shotgun. Na koniec zidentyfikowaliśmy węzły grafowe obecne w szczepach dżumy z jeziora Bajkalskiego, ale nieobecne we wszystkich współczesnych zespołach dżumy, i sklasyfikowaliśmy je na podstawie wzorca obecności/nieobecności w Y. pseudotuberculosis i Y. similis . Każdy węzeł został sklasyfikowany jako przodek (obecny zarówno w Y. pseudotuberculosis , jak i Y. similis ) lub pochodzący z Y. pseudotuberculosis lub Y. similis .

Oszacowanie wieku w chwili śmierci

Oszacowanie wieku w chwili zgonu opierało się na różnych ustalonych metodach antropologicznych. W przypadku osób niebędących dorosłymi (na ogół <20 lat) oceniano go poprzez formowanie się i wyrzynanie zębów, pomiary nasad i trzonów kości długich oraz zrost nasadowy, jak podsumowano w pozycjach piśmiennictwa ^{49 , 50 , 51 , 98 , 99.} Oszacowanie wieku osoby dorosłej koncentrowało się na zmianach morfologicznych szkieletu, mianowicie spojenia łonowego ^100 , 101 i powierzchni małżowiny usznej biodrowej ^{102 , 103 , 104} , ale także zamknięcia szwów podniebiennych i ektokranialnych ^{105 , 106 , 107.} W przypadku wszystkich osób rozważono jak najwięcej metod w oparciu o stan zachowania szkieletu i/lub uzębienia.

Datowanie radiowęglowe

Datowanie radiowęglowe przeprowadzono w Oxford Radiocarbon Accelerator Unit zgodnie z ustalonym protokołem w tej placówce108 ^. Nowe ustalenia z Bratskii Kamen, Serovo i Shumilkha przedstawiono tutaj po raz pierwszy, obok wcześniej opublikowanych dat radiowęglowych27 ^, 109 (Extended Data Fig. ⁴ iSupplementary Table 1 ) . Wszystkie daty ludzkie skorygowano o efekt zbiornika słodkowodnego (FRE), stosując równanie regresji dla południowo-zachodniego Bajkału/Angary110 ^, 111 ( Supplementary Note 4 ). Niewielka liczba dat na wisiorkach z zębów jelenia szlachetnego jest preferowana nad datami ludzkimi, tam gdzie jest dostępna, ponieważ unikają FRE. Bayesowskie modelowanie dat radiowęglowych przeprowadzono w OxCal 4.465 ^, wykorzystując nieinformacyjne modele jednofazowe z jednolitymi granicami. Aby zwizualizować zsumowane wiele dat, zastosowano modele estymacji gęstości jądra (KDE) i wykresy w modelach bayesowskich¹¹² .

Podsumowanie raportu

Więcej informacji na temat projektu badawczego można znaleźć w  podsumowaniu sprawozdawczości Nature Portfolio Reporting Summary, do którego link znajduje się w tym artykule.