Sekwencjonowanie genomu i porównanie dwóch modeli zwierzęcych naczelnych, makaka kynomolgusa i chińskiego makaka rhesus

Makaki są najbardziej rozpowszechnionymi ssakami naczelnymi, obejmującymi ponad 20 gatunków, które różniły się od siebie do 5-6 milionów lat temu2. Rodzaj Macaca jest blisko spokrewniony z ludźmi, dzieląc ostatniego wspólnego przodka ∼ 25 milionów lat temu3. Bliskie pokrewieństwo między ludźmi a makakami sprawiło, że kilka gatunków stało się atrakcyjnymi modelami zwierzęcymi dla różnych analiz biomedycznych. Chociaż indyjski podgatunek makaka rhesus (Macaca mulatta mulatta) był początkowo z wyboru modelem badawczym, zakaz wywozu tego makaka znacznie ograniczył dostępność tych zwierząt, prowadząc do zwiększonego wykorzystania innych gatunków i/lub podgatunków makaków, w szczególności chińskiego makaka rhesus (CR) (Macaca mulatta lasiota) i makaka cynomolgus lub makaka krabożernego (CE) (Macaca fascicularis).

Sekwencjonowaliśmy genomy samicy makaka CR i samicy makaka CE przy użyciu strategii whole-genome shotgun na platformie sekwencjonowania następnej generacji. Krótko mówiąc, analiza sekwencji genomu mitochondrialnego zweryfikowała przewidywane pochodzenie obu osobników (Supplementary Section 1). Następnie skonstruowaliśmy 19 i 18 wielokrotnych sparowanych bibliotek genomowego DNA o stopniowo wzrastających rozmiarach insertów odpowiednio dla makaka CR i makaka CE. Całkowity rozmiar zmontowanych genomów makaka CR i makaka CE wynosił odpowiednio ∼ 2,84 Gb i 2,85 Gb, zapewniając średnio 47-krotne i 54-krotne pokrycie (Tabela 1 i Dodatkowa Sekcja 1). Rusztowania zostały przypisane do chromosomów zgodnie z wykazaną syntezą z sekwencjami genomów makaka indyjskiego (IR)4 i człowieka. Około 97% rusztowań makaka CR i 92% rusztowań makaka CE udało się umieścić na chromosomach. Zastosowaliśmy również RNA-seq do profilowania transkryptów w różnych tkankach od jednego makaka IR i dwóch makaków CE (Metody Online). Zintegrowana analiza łącząca dane genomowe i transkryptomowe została następnie wykorzystana do określenia struktury transkryptów i ustalenia profilu ekspresji każdego genu (Supplementary Section 2).

Różnorodność genetyczna makaków została oceniona przez porównanie całych genomów i wyrównanie krótkich odczytów przy użyciu genomu makaka IR jako odniesienia. W sumie wykryliśmy >20 milionów różnic pojedynczych nukleotydów i 740 827 zdarzeń indel w trzech gatunkach lub podgatunkach makaków (Supplementary Section 3), co zapewni obfitą heterogeniczność genetyczną do wykorzystania w przyszłych aplikacjach i analizach biomedycznych. Sklasyfikowaliśmy wszystkie zmienne miejsca pojedynczego nukleotydu do trzech klas (wspólne, stałe i unikalne warianty) w oparciu o ich obecność lub brak u trzech osobników (ryc. 1a). Warianty unikalne stanowiły >71,7% wszystkich wariantów, co nie jest zaskakujące, biorąc pod uwagę fakt, że nawet w populacji panmictic, 44% alleli powinno być singletonami w próbce trzech osobników. Warto zauważyć, że duża liczba różnic genetycznych była wspólna dla co najmniej dwóch makaków. Używając tylko stałych i unikalnych wariantów, oszacowaliśmy, że najwyższy współczynnik dywergencji, 0.40%, był pomiędzy makakiem CE i makakiem IR (Rys. 1b). Jednakże, dywergencja sekwencji między makakiem CE i makakiem CR (0.34%), chociaż nominalnie różne gatunki, była zbliżona do tej obserwowanej między podgatunkami, makakiem CR i makakiem IR (0.31%).

Poprzednie badania oparte na ograniczonych danych sekwencyjnych sugerowały, że starożytna introgresja mogła wystąpić z makaków CR do makaków CE żyjących w nakładającej się geograficznej strefie dystrybucji na półwyspie Indochińskim5,6,7. Dwa sekwencjonowane genomy pozwoliły nam na ilościowe określenie wpływu tej introgresji na poziomie całego genomu. W szczególności zbadaliśmy, czy sygnał DNA zgodny z hybrydyzacją międzygatunkową był widoczny w genomach makaka CR i makaka CE. Obliczyliśmy stosunek dywergencji między makakiem CE i makakiem CR i porównaliśmy go ze stosunkiem dywergencji między makakiem CR i makakiem IR dla 50-kb okien w obrębie wyrównanych genomów (Dodatkowa Sekcja 4). Do tych obliczeń zignorowaliśmy zmiany w miejscach CpG, ponieważ wiadomo, że ewoluują one szczególnie szybko. Ponad 27% okien wykazywało współczynnik dywergencji mniejszy od zera, co sugeruje, że makaki CE i CR są bliżej spokrewnione niż podgatunki CR i IR makaków w tych regionach (Dodatkowa sekcja 4). Ponadto, >93% 50-kb okien genomowych wykazywało niższy wskaźnik dywergencji między makakami CE i CR w porównaniu z makakami CE i IR. Dlatego też nieusortowane polimorfizmy przodków nie mogły całkowicie wyjaśnić wysokiego odsetka niespójnych regionów obserwowanych między makakiem CE i makakiem CR. Ponadto, łącząc wcześniejsze dane dotyczące polimorfizmu pojedynczego nukleotydu (SNP) z populacji makaków IR i CR z danymi z naszych własnych sekwencjonowanych osobników makaków CR i CE8, zauważyliśmy, że nasz osobnik z makaka CE grupował się w obrębie populacji makaka CR (Dodatkowa sekcja 4). Wspiera to występowanie silnego przepływu genów z genomu makaka CR do genomu makaka CE. Badając stopień asymetrii w dywergencji między makakiem CE i CR oraz między makakiem CE i IR, oszacowaliśmy, że ∼30% genomu makaka CE pochodzi od makaka CR (Dodatkowa sekcja 4).

Następnie staraliśmy się zidentyfikować domniemane regiony introgresji (PIRs) w genomie makaka CE, które mogły zostać wniesione przez przepływ genów. Użyliśmy symulowanych danych (pod neutralnym modelem braku migracji) jako kontroli (Metody Online i Uzupełniająca Sekcja 4) i zidentyfikowaliśmy 8,942 PIRs obejmujące 778 Mb ze znacznie niższym niż oczekiwano wskaźnikiem dywergencji między makakiem CE i CR (Fig. 2a,b). Po połączeniu nakładających się PIRs, stwierdziliśmy, że większość PIRs (>98%) była krótsza niż 500 kb. Ponieważ rozkład długości PIR jest funkcją czasu, w którym nastąpił przepływ genów9, przewaga krótkich PIR sugeruje, że przepływ genów nastąpił w dłuższym okresie czasu ewolucyjnego i jest mało prawdopodobne, aby był po prostu konsekwencją bardzo niedawnego przepływu genów spowodowanego przez człowieka. Zaobserwowaliśmy również wyraźną różnicę w zmienności między chromosomami X i autosomalnymi (Dodatkowa sekcja 4), która mogła wynikać z przepływu genów napędzanego przez mężczyzn. Jednym z prawdopodobnych czynników przyczyniających się do ograniczonego przepływu genów od samic makaka CR do samców makaka CE jest fakt, że samice makaka CR wykazują wyraźną sezonowość owulacyjną i kopulują tylko podczas owulacji, podczas gdy samice makaka CE nie wykazują wyraźnej sezonowości reprodukcyjnej i pozostają seksualnie otwarte przez cały rok10. Ponadto, biorąc pod uwagę, że rozproszenie jest u makaków napędzane głównie przez samce z powodu filopatii samic, mogłoby to również tłumaczyć przepływ genów z samców makaków CR do samic makaków CE i brak przepływu w odwrotnym kierunku. Populacje te mogą być zatem interesujące dla badania fizjologicznych i behawioralnych aspektów reprodukcji między różnymi gatunkami.

Silna selekcja na korzyść nowych korzystnych alleli powoduje „selektywne zamiatanie”, które zmniejsza różnorodność genetyczną w stosunku do nieselekcjonowanych regionów. Opracowaliśmy algorytm do identyfikacji przypuszczalnych regionów wymiatania zawierających zmniejszoną zmienność między trzema gatunkami/podgatunkami makaków i wygenerowaliśmy dane symulacyjne przy założeniu neutralnego modelu, aby ocenić ich znaczenie statystyczne (Metody Online i Uzupełniająca Sekcja 5). Zidentyfikowaliśmy 217 silnie selektywnych regionów, które wykazywały zmniejszony poziom zmienności między makakami i które znacznie odbiegały od neutralnych oczekiwań (Rys. 2c, Uzupełniająca sekcja 5 i Uzupełniająca tabela 17). Warto zauważyć, że jeden z dziesięciu największych regionów selektywnego wymiatania, zlokalizowany na chromosomie 14 makaka, zawiera tylko jeden gen, czynnik wiążący SET 2 (SBF2) (ryc. 2c). Jest więc prawdopodobne, że gen ten, kodujący obwodowe białko błonowe z rodziny fosfataz białkowo-tyrozynowych, był celem pozytywnej selekcji podczas wczesnej ewolucji makaków. Z potencjalnego zainteresowania biomedycznego, defekty w ludzkim ortologu SBF2 powodują autosomalną recesywną demielinizacyjną postać choroby Charcot-Marie-Tooth (CMT4B2).

Aby ujawnić potencjalne cele pozytywnej selekcji w każdej gałęzi makaków, przypisaliśmy 14 978 1:1 ortologów genów dla człowieka, szympansa i trzech gatunków/podgatunków makaków przez wyrównanie genomów (Dodatkowa sekcja 6). Porównanie tria ortologów makaków ujawniło, że mają one niezwykle wysoki poziom podobieństwa sekwencji nukleotydów w obrębie regionów genów. Warto zauważyć, że 20,7% ortologów wykazuje wyższy stopień podobieństwa między makakiem CR a makakiem CE niż między makakiem CR a makakiem IR, co może sugerować wpływ introgresji. Porównanie kategorii genów oparte na ontologii genu między liniami Macaca, Hominid i Murid wykazało, że procesy oparte na mikrotubulach i szlak sygnałowy receptora insuliny ewoluowały szczególnie szybko w linii Macaca (Supplementary Section 6). Testy współczynnika prawdopodobieństwa oparte na modelu branch site ujawniły 16 pozytywnie wyselekcjonowanych genów specyficznie w gałęzi IR makaka, 7 w gałęzi CR makaka i 13 w gałęzi CE makaka (Supplementary Section 6). Intrygujący jest fakt, że 31 z 36 pozytywnie wyselekcjonowanych genów u makaków koduje białka wiążące, które odgrywają główną rolę w regulacji ekspresji genów. Warto również zauważyć, że dwa geny kodujące białka dendrytu, CLCN2 w linii IR makaka i regulowane aktywnością białko związane z cytoszkieletem (ARC) w linii CE makaka, doświadczyły pozytywnej selekcji. Te dwa geny, wraz z innymi pięcioma pozytywnie wyselekcjonowanymi genami, są już znane jako istotne dla ludzkich chorób genetycznych, co wskazuje na prawdopodobne znaczenie ich funkcji biologicznych.

Dostępność sekwencji genomu makaka CR i makaka CE pozwoliła nam ocenić ich różnorodność genetyczną, jak również różnice genetyczne między makakami i ludźmi, co jest ważne, biorąc pod uwagę znaczące wykorzystanie makaków w badaniach biomedycznych. Porównanie genomów makaków ujawniło brak 25 ludzkich genów jednokopijnych (Online Methods i Supplementary Section 6), w tym genu receptora chemokiny, IL32, który może odgrywać rolę zarówno we wrodzonej, jak i adaptacyjnej odpowiedzi immunologicznej, a zatem jest ważny do rozważenia, gdy te makaki są wykorzystywane w badaniach nad chorobami zakaźnymi. Ponadto, w sumie 170 genów związanych z chorobami lub odpornością u jednego lub drugiego gatunku makaków zawiera mutacje typu frameshift lub przedwczesne kodony stop, które, jak można przypuszczać, spowodowały pseudogenizację tych genów (Dodatkowa sekcja 6). Autentyczność tych okrojonych mutacji jest poparta zarówno danymi transkryptomu, jak i odczytami sekwencjonowania o dużej głębokości, a także niezależną walidacją PCR. Trzydzieści dwa z tych genów funkcjonują w szlakach odpornościowych i wydaje się, że zostały utracone u makaków. Na przykład, ważny gen odporności wrodzonej, DEFA4, który koduje jeden z mikrobójczych i cytotoksycznych peptydów wytwarzanych przez neutrofile11, został pseudogenizowany we wszystkich trzech makakach z powodu utraty jego pierwszego eksonu. Ponadto, gen receptora Toll-podobnego 4 (TLR4) zawierał delecję 1-bp, która generuje przedwczesny kodon stop w jego trzecim eksonie we wszystkich trzech makakach (Dodatkowa sekcja 6). Wykazano, że TLR4 podlega pozytywnej selekcji u naczelnych Starego Świata12. Warto zauważyć, że niektóre ludzkie geny związane z chorobami również zawierają przesunięcia ramki w swoich homologach u makaków. Na przykład stwierdziliśmy, że wszystkie trzy makaki miały przedwczesny kodon stop w drugim eksonie genu receptora opioidowego mu1 (OPRM1), który koduje białko rozmieszczone w całej neuraksji i obwodowym układzie nerwowym, i który jest głównym celem opioidów13 (Dodatkowa sekcja 6).

Zbadaliśmy również różnice genetyczne w ortologach, które są szczególnie ważne w badaniach biomedycznych. Cytoplazmatyczne białko 5α tripartite-motif (kodowane przez gen TRIM5), które może ograniczać replikację szerokiej gamy retrowirusów, jest kluczowym biomarkerem stosowanym do wyboru zwierzęcych modeli infekcji HIV14. Aby zbadać różnorodność genetyczną TRIM5 w całej populacji, amplifikowaliśmy PCR i sekwencjonowaliśmy TRIM5 od 33 niespokrewnionych osobników makaków CE pochodzenia wietnamskiego i 28 osobników makaków CR (Metody Online). Nie wykryliśmy wcześniej zgłoszonej chimery15 Trim5-cyklofiliny A (TRIM-CypA2) u żadnego osobnika, co sugeruje, że ten genotyp jest rzadki w tych populacjach. Jednakże, 19 nonsynonimicznych polimorfizmów i jedna mikrodelecja zostały zidentyfikowane w genie TRIM5 w stosunku do makaka IR; prawie wszystkie te polimorfizmy wykazywały różne częstotliwości pomiędzy dwoma populacjami (Fig. 3 i Supplementary Section 7). Zidentyfikowaliśmy również 6-bp delecję w genie TRIM5 u makaka CE, która skutkuje utratą dwóch aminokwasów (Thr339 i Phe340). Ostatnie badania wskazują, że delecja tych reszt może prowadzić do zwiększonej patogenności HIV lub SIV16. Wysoka częstotliwość (97,5%) tej mutacji została wykryta w populacji makaków CE, co wskazuje, że delecja ta stała się praktycznie utrwalona u makaków CE. Natomiast w populacji makaków CR częstotliwość występowania tej mutacji wynosi około 50% i jest tylko nieznacznie wyższa niż w populacji makaków IR (36%)17. Różnice w częstości występowania tej delecji 6-bp i innych polimorfizmów między makakami o różnym pochodzeniu geograficznym mogą być odpowiedzialne za obserwowane różnice w oporności na HIV między tymi gatunkami/podgatunkami makaków16. Zbadaliśmy również zmienność genetyczną w innych genach związanych z chorobą w tej samej populacji makaków CE lub CR, obserwując, że mutacje często występują z różną częstotliwością u tych dwóch gatunków (Dodatkowa sekcja 7).

Aby zbadać ortologi ludzkich domen białkowych nadających się do leczenia w makakach i stworzyć zasób do terapeutycznego wykorzystania „genomu nadającego się do leczenia”, przesiewaliśmy ortologów makaków dla obecnie znanych domen leków. Prawie wszystkie ortologi lekowalne mogą być wykryte u trzech gatunków/podgatunków makaków, co wskazuje, że te modele zwierzęce są prawdopodobnie funkcjonalnie równoważne. Jednakże, w bardzo niewielu przypadkach, ortolog znaleziony u makaka różni się od swojego ludzkiego odpowiednika. Na przykład, mitochondrialna acylotransferaza (GLYATL2), która przenosi grupę acylową do glicyny, została całkowicie utracona we wszystkich trzech makakach. Ponadto, zidentyfikowaliśmy 19 ludzkich genów z domenami lekowalnymi, które stały się pseudogenami u makaków (Dodatkowa sekcja 7). Na przykład, gen receptora hormonu przytarczyc 1 (PTH1R), cel leku przeciwosteoporotycznego teriparatydu (Forteo)18,19, zawiera przedwczesny kodon stop u makaków. Jeden z celów rekombinowanego ludzkiego czynnika wzrostu keratynocytów (Palifermin20), receptor czynnika wzrostu fibroblastów 3, kodowany przez FGFR3, również uległ pseudogenizacji u makaków z powodu obecności przedwczesnego kodonu stop.

Dodatkowym zainteresowaniem biomedycznym cieszą się skompensowane odchylenia patogenetyczne. Reprezentują one ludzkie potencjalnie patologiczne allele missense, w których zastępujące je aminokwasy są identyczne z resztami aminokwasów typu dzikiego w ortologicznych pozycjach w innych organizmach. Zidentyfikowaliśmy 931 skompensowanych patogennych odchyleń w czterech blisko spokrewnionych gatunkach naczelnych (szympans i trzy makaki), z których 220 różniło się między naczelnymi, w tym 65, które różniły się między trzema gatunkami makaków (Uzupełniająca sekcja 8 i Uzupełniająca tabela 26). Na przykład, jedna mutacja (R40→H40) w genie transkarbamylazy ornityny (OTC) była widoczna u dwóch podgatunków makaka rhesus, ale nie u makaka CE. Na podstawie przykładów zidentyfikowanych różnic genetycznych przedstawionych powyżej, jasne jest, że potencjalne istnienie takich różnic międzygatunkowych należy rozważyć przy wyborze makaków do wykorzystania jako modele chorób.

Porównanie profili ekspresji genów (Supplementary Section 9) między makakiem CE i makakiem IR ujawniło, że ich ortologowie wyświetlali konserwowane profile ekspresji w tych samych tkankach. Jednakże zauważyliśmy, że jądro wykazywało poziomy ekspresji, które były bardziej rozbieżne pomiędzy tymi ortologami, które miały niższe współczynniki korelacji Pearsona (Supplementary Section 9). Obserwacja, że więcej genów wykazuje niespójne poziomy ekspresji w jądrze w porównaniu z innymi tkankami, może być związana z szybkim tempem ewolucji genów ekspresjonowanych przez plemniki naczelnych21. Dane transkryptomu posłużyły również do zidentyfikowania kilku nowych genów u makaków CE w odniesieniu do makaków rhesus.

Podsumowując, nasze sekwencjonowanie i analizy genomów dwóch makaków potwierdziły, że introgresywna hybrydyzacja prawdopodobnie odegrała ważną rolę w tworzeniu genomu istniejącego makaka CE pochodzącego z kontynentu. Tak więc, makak CE może być użytecznym modelem do badania wymiany genów między gatunkami naczelnych, a w konsekwencji roli tego procesu w ewolucji naczelnych i specjacji. Dwa nowe genomy makaków przedstawione tutaj podkreślają również stopień zróżnicowania istniejący pomiędzy tymi szeroko stosowanymi modelami zwierzęcymi naczelnych. Bogata różnorodność genetyczna widoczna u poszczególnych makaków z różnych populacji geograficznych jest przedmiotem bezpośredniego zainteresowania prymatologii, medycyny przedklinicznej, genetyki populacyjnej i badań filogeograficznych.