В настоящий момент вирус клещевого энцефалита (ВКЭ) принято разделять на 3 субтипа: дальневосточный (ВКЭ-ДВ), европейский (ВКЭ-Евр) и сибирский (ВКЭ-Сиб) [1]. Также недавно в литературе было опубликовано несколько работ с описанием двух новых субтипов: Гимайлайского и Байкальского [2—5]. В названии каждого варианта отражены основные места циркуляции первых его представителей. Однако впоследствии штаммы ВКЭ всех 3 субтипов были изолированы за пределами территорий, обозначенных в их названии. Так ВКЭ-Евр, помимо Центральной Европы, обнаружен в странах Балтии [6, 7], на территории Западной и Восточной Сибири [8—11], а также в Южной Корее [12—14]. ВКЭ-Сиб встречается наиболее широко и распространен по всему ареалу ВКЭ [15]. ВКЭ-ДВ в большей степени распространен в дальневосточном регионе Евразии, однако имеются находки на территории стран Балтии, Крыма, Урала, Сибири, а также в европейской части России [16, 17]. Филогеографические связи наиболее удаленных друг от друга очагов неясны.
Цель данной работы — филогенетический анализ штаммов ВКЭ-Евр, выделенных на территории Восточной и Западной Сибири, с установлением времени образования их общего предка при помощи Байесовского метода. Кроме того, основываясь на опубликованных ранее работах, посвященных географическому распространению ВКЭ, мы проанализировали филогеографические отношения между исследуемой группой штаммов ВКЭ-Евр с территории Сибири и ближайшими изолятами с территории Европы и европейской части России.
Материал и методы
В анализе использованы 53 нуклеотидные последовательности гена полипротеина (10 245 н.о.) ВКЭ-Евр из международной базы данных GenBank. Также в анализ включены нуклеотидные последовательности двух штаммов (№214, №172) ВКЭ-Евр из коллекции Иркутского научно-исследовательского противочумного института Сибири и Дальнего Востока, изолированные на территории Иркутской области (Эхирит-Булагатский и Усть-Удинский районы). Время изоляции штаммов №214 и №172 — 1967 и 1968 гг. соответственно. Оба штамма изолированы из ликвора больных клещевым энцефалитом. Реизоляцию штаммов проводили на 2—3-суточных беспородных белых мышах путем заражения в мозг в объеме 0,02 мл с последующим однократным закреплением в пассаже.
Выделение суммарной РНК осуществляли с помощью набора Рибо-преп, обратную транскрипцию — набором Реверта-L производства «АмплиСенс» (ФГУН ЦНИИ эпидемиологии Роспотребнадзора, Москва) в соответствии с инструкциями изготовителя. Амплификацию фрагментов вирусного генома проводили с помощью набора Синтол (Москва). Секвенирование нуклеотидных последовательностей из ПЦР-продукта осуществляли с использованием набора реактивов ABI PRISM BigDye Terminator v.1.1 Cycle Sequencing Kit на приборе Genetic Analyzer 3500 xL («Applied Biosystems»).
Расшифрованные нуклеотидные последовательности штаммов №214 (MK562430) и №172 (MK560446) депонированы в базу данных GenBank.
Выявление полиморфизма анализируемого набора нуклеотидных последовательностей проводили в программе IQTREE v.1.6.12 [18] и MEGA X [19]. Насыщение заменами в 1+2 и 3 позициях кодона тестировали в программе DAMBE v.6.4 методом X. Xia и соавт. [20], в основе которого лежит сравнение индекса насыщения заменами (an index of substitution saturation, Iss) и критического индекса насыщения (Iss.c). В случае Iss>Iss.c позиция кодона считается перенасыщенной нуклеотидными заменами и рассматривается как непригодная для проведения филогенетической реконструкции и наоборот. Статистическую значимость при расчете индексов оценивали при помощи «двухвостового» (two-tailed) теста в той же программе.
Для проведения филогенетического анализа и оценки времени дивергенции исследуемой группы применяли пакет программ BEAST v.1.8.4 (Байесовский филогенетический метод) [21]. Воспроизводимость результатов для каждого сочетания моделей оценивали в 4 независимых запусках алгоритма Монте-Карло по схеме Марковской цепи (MCMC), длина запуска при этом составляла от 50 до 100 млн итерации (запуск MCMC продолжался до тех пор, пока значение эффективного размера выборки (ESS) не достигало 200). Частоту сохранения данных выбирали таким образом, чтобы суммарное количество образцов (деревьев) в каждом запуске было равно 50 000 (величина отжига при этом составляла 10% от общей длины цепи).
Сравнение комбинаций моделей проводили на основании значений функции маргинального правдоподобия, вычисленных методами Path Sampling и Stepping-Stone Sampling (PS/SS) [22], количество шагов равнялось 100, длина MCMC для каждого шага — 1 млн итераций. Дополнительное сравнение параметров эволюционных моделей (учет параметра α для гамма-распределения неравномерности скорости накопления замен (+G4), а также доли инвариантных сайтов (+I) в выравнивании) проводили на основании значений Байесовского информационного критерия (BIC), рассчитанного с использованием пакета Phangorn, реализованного на языке программирования R [23]. В анализе рассмотрены 3 эволюционные модели (HKY, GTR, SRD06), модель строгих молекулярных часов, популяционные модели экспоненциального роста генетического разнообразия популяции (EG), а также непараметрические популяционные модели Bayesian Skyline (BSL) [24] и Bayesian Skygrid (BSG) [25], способные учитывать множественные изменения генетического разнообразия во времени. Расслабленные молекулярные часы не использовались в анализе по причине низкого (0,08) значения коэффициента вариации для скорости накопления нуклеотидных замен между ветвями филогенетического дерева. Статистическая неопределенность для времени дивергенции и эволюционной скорости оценивалась при помощи 95% интервала максимума плотности апостериорной вероятности (95% HPD).
Оценка филогеографических связей «сибирских» штаммов ВКЭ-Евр проводилась с использованием регрессионного анализа в соответствии с методическими указаниями из работы D. Heinze и соавт. [26]. Генетическую дистанцию вычисляли относительно штамма Sofjin-HO в программе MEGA X (модель замен — Kimura 2-parameters+G4+I). Географическую дистанцию рассчитывали от места изоляции штамма Sofjin при помощи сервиса Google Earth (https://www.google.com/earth/index.html). В анализе также использовались нуклеотидные последовательности штаммов ВКЭ и шотландского энцефаломиелита овец (ШЭО) из работы D. Heinze и соавт. [26] с учетом указанного в той же работе их географического расстояния от штамма Sofjin.
Результаты и обсуждение
Ранее в работах ряда исследователей была дана подробная характеристика «сибирских» штаммов ВКЭ-Евр, включающая описание вирулентных и инвазивных свойств выделенных штаммов, процент гомологии на нуклеотидном и аминокислотном уровнях, были указаны характерные аминокислотные замены в различных белках структуры полипротеина ВКЭ-Евр, а также показана относительно высокая генетическая стабильность ВКЭ-Евр в целом [8—11]. Исследователи рассмотрели филогенетические связи между «сибирскими» штаммами ВКЭ-Евр и их положение на общем филогенетическом дереве. Однако в упомянутых работах расчет времени образования кластера, сформированного штаммами, выделенными на территории Сибири, проведен не был. Для решения этой задачи мы применили Байесовский филогенетический подход с использованием в анализе двух новых нуклеотидных последовательностей ВКЭ-Евр, кодирующих ген полипротеина.
Филогенетический анализ. Анализ полиморфизма нуклеотидных последовательностей в исследуемом наборе показал следующие результаты: общее количество полиморфных сайтов составило 1309, количество информативных сайтов (parsimony informative sites) — 1109, среднее и максимальное количество несовпадающих нуклеотидов — 213 и 316 соответственно. Согласно полученным значениям, исследуемый набор геномных данных можно характеризовать как информативный. Анализ нуклеотидных последовательностей в программе DAMBE методом X. Xia и соавт. [20] показал отсутствие насыщения заменами в 1+2 и 3 позициях кодона — в обоих случаях значение Iss (0,016 и 0,095 для 1+2 и 3 позиций соответственно) оказалось значительно меньше Iss.c (0,57 и 0,55 для 1+2 и 3 позиций соответственно); p-value составило 0,00.
На основании предварительного теста эволюционных моделей, проведенного при помощи пакета Phangorn, в наш анализ добавлены расчеты неравномерности скорости накопления нуклеотидных замен (количество дискретных категорий — 4) и вычисление доли инвариантных сайтов (табл. 1).
Таблица 1. Сравнение эволюционных моделей GTR и HKY с различными сочетаниями параметров формы гамма-распределения (4 дискретные категории) и доли инвариантных сайтов
| Модель | Сравнение моделей | |||
| logLik | AIC | AICc | BIC | |
| GTR+G4+I | –32752,9 | 65731,8 | 65734,3 | 66549,3 |
| GTR+G4 | –32770,8 | 65765,7 | 65768,2 | 66576,0 |
| GTR+I | –32794,5 | 65813,1 | 65815,6 | 66623,4 |
| HKY+G4+I | –32959,5 | 66137,1 | 66139,5 | 66925,7 |
| HKY+G4 | –32985,0 | 66186,1 | 66188,4 | 66967,4 |
| HKY+I | –33002,4 | 66220,92 | 66223,24 | 67002,25 |
| GTR | –33453,6 | 67129,23 | 67131,69 | 67932,27 |
| HKY | –33764,6 | 67743,39 | 67745,67 | 68517,48 |
Примечание. logLik — значение логарифма правдоподобия; AIC — значение информационного критерия Акаике; AICc — значение скорректированного критерия Акаике; BIC — значение Байесовского информационного критерия (критерий Шварца). Полужирным начертанием выделены модели с наивысшим значением BIC, которые в дальнейшем были использованы для проведения филогенетического анализа в BEAST.
Четыре независимые запуска МСМС для каждого сочетания моделей показали схождение цепи в одном и том же оптимуме значений, показатели ESS в каждом запуске превышали необходимый порог (>200). Методы PS/SS показали наивысшее значение функции маргинального правдоподобия для сочетания строгих молекулярных часов, эволюционной модели SRD06 и популяционной модели BSG (MSC+SRD06+BSG).
Модель SRD06 подразумевает разделение массива данных на две партиции — 1+2 и 3 позиции кодонов. Для каждой партиции мы применили эволюционную модель HKY+G4 с учетом частот нуклеотидов (данное сочетание параметров показало себя наилучшим образом при анализе 177 локусов различных РНК-содержащих вирусов [27]).
Предварительный запуск с использованием расслабленных молекулярных часов показал значение коэффициента вариации <10% (μ=0,08; 95% HPD, 0—0,15), в связи с чем данную модель в работе не рассматривали [28]. Результаты сравнения моделей методами PS/SS представлены в табл. 2.
Таблица 2. Сравнение значений функций маргинального правдоподобия для всех рассматриваемых сочетаний моделей
| Строгие часы | |||||
| модель замен | популяционная модель | маргинальное правдоподобие | эволюционная скорость1 (95% HPD) | возраст корня2 (95% HPD) | |
| PS | SS | ||||
| SRD06 (G4) | BSG | –32306,8 | –32308,3 | 1,3E-5 (9,5E-6—1,8E-5) | 1173 (785—1557) |
| BSL | –32312,9 | –32313,7 | 1,2Е-5 (4,7Е-6—2,0E-5) | 1268 (706—2730) | |
| EG | –32339,7 | –32341,4 | 1,2Е-5 (2,2Е-6—2,2E-5) | 1319 (551—3958) | |
| GTR+G4+I | BSG | –33085,5 | –33082,6 | 1,3E-5 (9,2E-6—1,8E-5) | 1170 (796—1595) |
| BSL | –33082,8 | –33082,0 | 1,1Е-5 (4,4Е-6—1,9E-5) | 1315 (665—2755) | |
| EG | –33104,2 | –33035,2 | 1,1E-5 (2,5E-6—2,0E-5) | 1405 (598—4092) | |
| HKY+G4+I | BSG | –33266,8 | –33267,2 | 1,3E-5 (1,0E-5—1,8E-5) | 1154 (787—1530) |
| BSL | –33273,9 | –33274,9 | 1,22E-5 (4,7E-6—2,0E-5) | 1272 (673—2685) | |
| EG | –33296,7 | –33296,8 | 1,2E-5 (2,7E-6—2,1E-5) | 1332 (630—3621) | |
Примечание. Полужирным начертанием выделено наиболее предпочтительное сочетание моделей с наивысшим значением функции маргинального правдоподобия; значения эволюционной скорости и возраста корня даны по медиане апостериорного распределения, аппроксимированного в результате запуска МСМС; 1 — значение эволюционной скорости выражено в количестве нуклеотидных замен на сайт в год; 2 — значение возраста корня дерева выражено в годах. За ноль как точку отсчета принят 2017 г.
На рис. 1 изображено филогенетическое дерево, реконструированное в программе BEAST v.1.8.4 с использованием 53 нуклеотидных последовательностей (10 245 н.о.) ВКЭ-Евр, выделенных в период 1951—2017 гг. Скорость фиксации нуклеотидных замен была рассчитана с относительно высокой точностью и составила 1,3E-5 нуклеотидов на сайт (или на одну позицию в геноме) в год или 1 нуклеотидная замена в полипротеине за 7,5 года (95% HPD, 1,0E-5—1,8E-5 нуклеотидов на сайт в год или 1 нуклеотидная замена в промежутке 5,4—9,7 года). Применение модели SRD06 показало, что скорость фиксации нуклеотидных замен в 3 позиции кодона превышает скорости для 1+2 позиций в 7,9 раза. При пересчете в абсолютных значениях (с учетом средней эволюционной скорости, выраженной в годах, см. выше) скорость фиксации нуклеотидных замен для 1+2 позиций кодона составила 2,9E-6 нуклеотидов на сайт в год (95% HPD, 2,2E-6—4,0E-6), для 3 позиции кодона скорость составила 2,3E-5 нуклеотидных замен на сайт в год (1,7E-5—3.2E-5).
Рис. 1. Филогенетическое дерево, реконструированное на основании 53 белок-кодирующих нуклеотидных последовательностей (10245 н.о.) ВКЭ-Евр при помощи Байесовского метода, реализованного в программе BEAST.
Слева от основных узлов указаны значения апостериорной вероятности. Серым цветом выделен кластер S, сформированный штаммами ВКЭ-Евр с территории Сибири. Для основного узла кластера S серой горизонтальной полосой отображены значения 95% HPD. Ниже филогенетического дерева изображена временная шкала, выраженная в календарных годах. Геномы, расшифрованные в ходе данной работы, выделены полужирным начертанием.
Двенадцать штаммов с территории Западной и Восточной Сибири сформировали кластер S (апостериорная вероятность — 1). Время существования ближайшего общего предка для данного кластера оценивается приблизительно 1928 г.; доверительный интервал составляет 50 лет (95% HPD, 1900—1950 гг.). Вследствие небольшого количества информативных сайтов у 12 рассматриваемых штаммов (7 сайтов при длине последовательностей 10 245 н.о.) внутри кластера S наблюдаются политомия и низкие значения апостериорной вероятности. Однако при наличии в общем наборе нуклеотидных последовательностей достаточного количества информации для расчета эволюционной скорости время дивергенции для неразрешенных кластеров можно установить с высокой точностью.
Филогеография «сибирских» штаммов ВКЭ-Евр. Исследователями D. Heinze и соавт. [26] был проведен анализ пространственной и временной динамики флавивирусов, переносимых клещами (tick-borne flaviviruses (TBFV)). Авторами исследования был показан линейный характер зависимости между генетической дистанцией и географическим расстоянием в км (R2=0,91) для TBFV, в том числе для ВКЭ трех основных субтипов. На основании выявленной зависимости была создана предполагаемая модель глобального распространения TBFV, согласно которой отделение ВКЭ-Евр от основной генетической линии (иначе говоря, образование ВКЭ-Евр) произошло на территории Европы, откуда берут свое начало все его современные представители, в том числе штаммы ВКЭ-Евр, выделенные на территории Сибири.
Повторное проведение регрессионного анализа на основе работы D. Heinz и соавт. с использованием «сибирских» штаммов ВКЭ-Евр показало, что данные штаммы располагаются значительно выше выявленной линии регрессии (рис. 2), что явно указывает на их относительно недавнее появление на территории Сибири (в той же работе схожие выводы были сделаны относительно вируса Повассан, выделенного на территории Приморского края).
Рис. 2. Линейная регрессия, характеризующая взаимосвязь между генетической дистанцией и географическим расстоянием штаммов ВКЭ и ШЭО (R2=0,92 без учета «сибирских» штаммов ВКЭ-Евр).
Заключение
По результатам проведенного анализа эволюционная скорость для ВКЭ-Евр составила 1,3E-5 нуклеотидов на сайт в год или 1 нуклеотидная замена на 7,5 года (95% HPD, 1,0E-5—1,8E-5 нуклеотидов на сайт в год или 1 замена в промежутке 5,4—9,7 года). Рассчитанное время образования группы штаммов ВКЭ-Евр, выделенных на территории Сибири, датируется приблизительно 1928 г. (95% HPD, 1900—1950 гг.). Относительно молодой эволюционный возраст кластера, высокое генетическое сходство изолятов, а также результаты регрессионного анализа свидетельствуют о недавнем заносе и последующем увеличении генетического разнообразия ВКЭ-Евр на территории Сибири.
Права человека и животных. Экспериментальная работа на животных выполнена в соответствии с «Правилами лабораторной практики в Российской Федерации», утвержденными Приказом Министерства здравоохранения Российской Федерации от 01.04.16. №199н.
Финансирование. Дополнительного финансирования данной работы не было.
Авторы заявляют об отсутствии конфликта интересов.
The authors declare no conflicts of interest.