Авторитетный инструмент для диагностики подвидов европейских медоносных пчел на основе информативных SNP предков

Аннотация

Справочная информация

Благодаря многочисленным эндемичным подвидам, представляющим четыре из пяти эволюционных линий, в Европе сохраняется значительная часть генетического разнообразия Apis mellifera. Это разнообразие и естественный ареал распространения были изменены антропогенными факторами. Сохранение этого природного наследия зависит от наличия точных инструментов для диагностики подвидов. На основе данных о последовательностях 2145 рабочих пчел из 22 популяций, отобранных по всей Европе, мы использовали два высокодискриминантных подхода (PCA и FST) для отбора наиболее информативных SNP для определения происхождения.

Результаты

Используя метод контролируемого машинного обучения (ML) и набор из 3896 генотипированных особей, мы смогли показать, что 4094 отобранных однонуклеотидных полиморфизма (SNPs) обеспечивают точное предсказание родословной европейских медоносных пчел. Лучшей ML-моделью оказался линейный классификатор опорных векторов (Linear Support Vector Classifier), который правильно отнес большинство особей к одному из 14 подвидов или к различным генетическим происхождениям со средней точностью 96,2% ± 0,8 SD. Всего 3,8% особей были неправильно классифицированы, что, скорее всего, связано с ограниченной дифференциацией между подвидами, вызванной близким географическим соседством, или вмешательством человека в генетическую целостность референтных подвидов, или их комбинацией.

Выводы

Представленный здесь диагностический инструмент будет способствовать устойчивому сохранению и поддержке селекционной деятельности с целью сохранения генетического наследия европейских медоносных пчел.

Справочная информация

Медоносные пчелы (Apis mellifera L.) являются наиболее важными управляемыми опылителями и в настоящее время находятся под угрозой из-за множества факторов давления по всему миру [1, 2]. Этот вид демонстрирует значительную изменчивость в пределах своего естественного ареала и состоит, по меньшей мере, из 30 описанных подвидов, принадлежащих к различным эволюционным линиям [3,4,5,6]. В Европе сосредоточена большая часть этого разнообразия медоносной пчелы с многочисленными эндемичными подвидами, представляющими четыре эволюционные линии, а именно: африканскую линию (А), линию Центральной и Восточной Европы (С), линию Западной и Северной Европы (М) и линию Ближнего Востока и Центральной Азии (О) [7, 8]. Однако на это разнообразие и естественный ареал распространения европейских медоносных пчел повлияли антропогенные факторы в такой степени, что несколько локально адаптированных популяций оказались под угрозой из-за интрогрессии и скрещивания [9,10,11]. Крупномасштабное разведение маток, коммерческая торговля и миграция пчел на большие расстояния могут снизить генетическое разнообразие и привести к генетической гомогенизации смешанных популяций [9, 12] и потенциальной последующей утрате местных адаптаций. Фактически, было продемонстрировано, что местноадаптированные медоносные пчелы имеют более высокую выживаемость [13], из чего следует, что сохранение основной генотипической вариации должно быть приоритетом для долгосрочной устойчивости популяций [14]. Для сохранения природного наследия медоносных пчел и, следовательно, их адаптационного потенциала к будущим глобальным изменениям, необходимо содействовать устойчивому разведению сертифицированных местных подвидов.

Многочисленные усилия по сохранению местных медоносных пчел были инициированы в Европе [9, 10, 15, 16]. Успех таких усилий по сохранению, включая программы генетического улучшения [17, 18], зависит от спаривания в интересующей популяции, что осложняется системой спаривания медоносных пчел, в которой девственные королевы свободно спариваются с множеством трутней из окружающих колоний [19, 20]. Помимо использования изолированных пасек для спаривания или искусственного осеменения, успешные меры контроля спаривания могут включать различные методы управления матками и трутнями [21] и регулярный мониторинг генетического происхождения и родства. В некоторых странах и регионах Европы импорт маток ограничивается местными подвидами медоносных пчел [22, 23] или экотипами [24, 25]. В таких случаях при торговле матками или колониями через государственные границы необходимо проверять происхождение маток. Кроме того, подтверждение генетического происхождения продуктов пчеловодства в виде сертифицированной маркировки местных пчел может помочь пчеловодам лучше реализовывать свою продукцию [26]. Таким образом, для осуществления эффективного пограничного контроля, повышения экономической ценности продуктов пчеловодства и поддержки обоснованных решений по сохранению и селекции по всей Европе существует потребность в диагностических генетических тестах для надежного определения подвида происхождения.

С развитием технологии высокопроизводительного секвенирования и генотипирования в последнее десятилетие для многих видов теперь доступны эталонные геномы, данные о последовательности всего генома и тысячи индивидуальных генотипов. В этих зачастую огромных наборах данных можно найти высокоинформативные однонуклеотидные полиморфизмы (SNP), которые затем можно использовать для генотипирования большего числа особей [27, 28]. Такие панели генотипирования, основанные на выбранном наборе информативных SNP, были разработаны для многих видов, включая человека, и могут использоваться для интрогрессии, генетического происхождения, структуры популяции, идентификации генетического фонда и пищевой криминалистики [29,30,31].

Для отбора информативных SNP из больших панелей генотипирования или данных последовательности использовались различные подходы (обзор в [32, 33]). Наиболее распространенным и популярным методом отбора является популяционная дифференциация, оцениваемая по FST, которая основана на различиях в частоте аллелей между популяциями, выражающих вариации между популяциями по отношению к общей популяции [34, 35]. Анализ главных компонент (PCA) также использовался для выявления информативных SNP, поскольку он снижает размерность признаков при незначительной потере информации и особенно выгоден при сложной структуре популяций [28, 36]. Учитывая набор информативных SNP-маркеров, применяется контролируемая классификация и так называемые тесты присвоения, в результате которых особь относят к заранее определенным классам (т.е. подвидам или популяциям происхождения). Классическое применение тестов на отнесение в популяционной генетике впервые использовало контролируемые параметрические подходы на основе правдоподобия [37, 38]. Недавно в вычислительной популяционной геномике появились новые методы, которые в совокупности называются контролируемым машинным обучением (НК) [39]. Общий подход для любых контролируемых классификаторов ML заключается в разделении данных на эталонное (обучающее) множество для «обучения» функции, которая может различать классы данных [40]. Затем эта функция используется для прогнозирования вероятности принадлежности «неизвестного образца» (теста) к любому заданному классу (например, подвиду). Точность классификации, выраженная как доля тестовых особей, правильно отнесенных к популяции их происхождения, зависит от свойств обучающего набора данных (т.е. от количества образцов, генетического разнообразия, уровня дифференциации популяций, степени перекрытия в распределении данных и качества эталонных образцов) [41]. Классификаторы ML нацелены на оптимизацию предсказательной точности алгоритма, а не на выполнение оценки параметров вероятностной модели, и они обладают потенциалом быть агностичными к оценке данного набора данных, т.е. без предположений о процессах, ведущих к дифференциации, включая эволюционную историю [39].

Для медоносных пчел были разработаны различные панели SNP, например, для выявления и оценки интрогрессии С-линейного ряда в М-линейных подвидах A. m. iberiensis и A. m. mellifera [15, 42,43,44,45,46]. Последний подвид является родным для северной и западной Европы и когда-то занимал большую часть европейской территории, но сейчас находится под угрозой исчезновения и даже был полностью вытеснен на большей части своего ареала [10, 47, 48]. Кроме того, были разработаны панели SNP для определения уровня африканизации и происхождения медоносных пчел Нового Света и Австралии [46, 49, 50]. Однако для большинства подвидов A. mellifera, популяции которых были генетически исследованы в меньшей степени или вообще не исследованы, молекулярных знаний на таком уровне детализации все еще не хватает. Эти подвиды и локально адаптированные популяции или экотипы представляются более уязвимыми в связи с существующими многочисленными угрозами для медоносных пчел.

Проект SmartBees был инициирован с целью разработки новых инструментов для описания и сохранения разнообразия медоносных пчел в Европе. Мы разработали молекулярный инструмент, состоящий из высокоинформативных SNP-маркеров, пригодных для отнесения особей медоносных пчел к их родовому подвиду, на основе всесторонней выборки европейского разнообразия медоносных пчел. На основе данных последовательности 1995 рабочих пчел, представляющих 22 популяции, четыре эволюционные линии и 14 подвидов, мы отобрали 4400 информативных SNP, используя два мощных и широко используемых подхода (FST и PCA). Из них 4165 SNP, для которых можно было разработать зонды и которые прошли проверку качества декодирования BeadChip, были генотипированы в 3903 отдельных пчелах с использованием платформы Illumina Infinium. Окончательная фильтрация контроля качества оставила 4094 надежных SNP для построения статистической модели с использованием алгоритмов машинного обучения (ML) для отнесения европейских медоносных пчел к 14 различным генетическим происхождениям. Лучшей моделью оказался линейный классификатор опорных векторов (Linear Support Vector Classifier, Linear SVC), который смог правильно отнести 96,2% протестированных образцов к их генетическому происхождению. Таким образом, представленный здесь метод точно идентифицирует европейские подвиды, что имеет решающее значение для поддержки стратегий управления в программах устойчивого разведения и сохранения медоносных пчел.

Результаты

Образцы и секвенирование пула

В общей сложности 22 популяции, представляющие четыре европейские эволюционные линии и 14 подвидов, были отобраны из их родных ареалов по всей Европе и прилегающим регионам (Таблицы 1 и S1). Каждая отобранная популяция включала до 100 рабочих пчел из неродственных колоний, в общей сложности 2145 образцов, что представляет собой наиболее полную выборку для изучения европейских медоносных пчел на сегодняшний день. Образцы из каждой популяции были гомогенизированы, объединены в пул и из них была выделена ДНК. Секвенирование на Illumina HiSeq 2500 позволило получить 1,6 миллиарда парных фрагментов (3,2 миллиарда отдельных чтений) со средней длиной чтения 125 п.н. и общей глубиной покрытия генома 2800x. Статистику секвенирования и вариантов можно найти в таблице S2.

Отобранные SNPs

В то время как основные эволюционные линии легко дифференцировались с помощью всего нескольких SNP (рис. S1A), дифференцировать близкородственные подвиды с меньшим числом генетических маркеров оказалось сложнее. Учитывая сложную, иерархическую структуру популяций европейских медоносных пчел, мы использовали два мощных и широко используемых подхода, PCA (Рисунок S1) и FST, чтобы определить наиболее дискриминантные маркеры для дифференциации подвидов европейских медоносных пчел (см. подробности в Методах и дополнительных материалах и методах). На основе вариантов, полученных из данных пула последовательностей, мы отобрали 4400 информативных SNP, из которых в общей сложности 4165 SNP прошли метрику качества декодирования для генотипирования с использованием пользовательского BeadChip Illumina Infinium, что указывает на то, что 99% первоначально представленных зондов были пригодны для генотипирования. SNPs распределены по всем 16 хромосомам медоносной пчелы, а также в неразмещенных контигах (Таблица S3), со средним расстоянием между SNPs 64 кб. Информация о SNP и геномное положение 4165 SNP, отобранных для дифференциации европейских подвидов медоносной пчелы, представлены в Дополнительном файле 1.

Генотипирование и визуализация образцов

Из 4165 SNPs 4094 были успешно генотипированы в 3896 отдельных пчелах с помощью технологии Illumina Infinium BeadChip (Таблица 1). Поскольку только 71 SNP не дал никаких данных, коэффициент успешного генотипирования (валидации SNP) составил 98%. Средняя частота вызовов на особь составила 0,87, варьируя среди образцов каждого подвида от 0,84 у A. m. cypria до 0,89 у A. m. adami (Таблица S4). Более чем у трети образцов коэффициент срабатывания превысил 0,9.

Данные о генотипах особей из пула секвенирования визуализированы на графике t-SNE [51], который сводит высокоразмерные данные к двумерной карте, где каждая особь представлена точкой (рис. 1). Генотипированные образцы были сгруппированы в несколько отдельных кластеров в соответствии с их эволюционной линией или подвидом происхождения (рис. 1). В пределах каждой линии большинство особей одного географического происхождения были тесно сгруппированы вместе и в целом хорошо отделены от соседних групп. Единственный подвид A-линии в нашем исследовании, A. m. ruttneri, был помещен в центр, промежуточный по отношению к другим кластерам. В роде О пчелы A. m. cypria были хорошо отделены от A. m. anatoliaca, A. m. caucasia и A. m. remipes, которые выглядят менее хорошо дифференцированными. Два подвида M-линейки были хорошо дифференцированы, при этом популяции A. m. mellifera были сгруппированы в три субкластера, разделяющие отдаленные (Бурзянский регион, Россия, верхний кластер A. m. mellifera на рис. 1) или изолированные (остров Лезё, Дания, нижний кластер A. m. mellifera) регионы отбора образцов. Образцы С-линейки сгруппированы в три субкластера: (i) A. m. ligustica, (ii) пчелы A. m. carnica, включая часть образцов «A. m. carpatica», и (iii) гетерогенный субкластер из A. m. macedonica, A. m. cecropia, A. m. adami, «A. m. rodopica» и остальных пчел «A. m. carpatica». График t-SNE с маркировкой образцов в соответствии с их пулом происхождения представлен на рисунке S2.

Классификация образцов с помощью машинного обучения

Мы использовали методы машинного обучения (ML) для построения модели классификации и отнесения европейских медоносных пчел к подвидам их происхождения. Из всех протестированных алгоритмов ML наилучшие результаты показала линейная модель SVC (Таблица S5). Модель рассчитывает вероятность предсказания принадлежности образца к любой из 14 референтных популяций. Каждый образец был отнесен к тому подвиду, который показал наибольшую вероятность предсказания, варьирующуюся от 0,29 до 1,0 с медианой 0,98 (Рисунок S3).

Матрица смешения использовалась для обобщения, описания и визуализации производительности линейной модели классификации SVC на наборе тестовых данных (вневыборочные данные, N = 2505), для которых были известны истинные значения (подвиды). Для родословных модель способна предсказать все образцы со 100% точностью (Рисунок S4). Для подвидов матрица смешения показала, что для большинства из них модель точно предсказала родословную тестовых образцов (N = 2505), лишь с несколькими исключениями (рис. 2a). Точность колебалась от 65 до 100%, указывая на то, что некоторые подвиды легче различить, чем другие. В целом 96,2% тестовых образцов были предсказаны правильно, в то время как 95 особей (3,8%) были неправильно классифицированы, т.е, например: четыре пчелы A. m. ligustica были предсказаны как A. m. carnica, две пчелы «A. m. carpatica» — как A. m. carnica или A. m. macedonica, а 23 пчелы A. m. cecropia были предсказаны как A. m. macedonica.

Модель предсказывает вероятность того, что данный образец принадлежит к одному из 14 исследуемых подвидов. Исходя из этого, тестовые образцы были отнесены к определенному подвиду на основании наибольшей вероятности предсказания, даже если эта вероятность была низкой (см. выше). Поэтому для повышения достоверности классификации мы установили порог вероятности, чтобы гарантировать, что только образцы с высокой вероятностью принадлежат к любому из 14 подвидов, а тестовые образцы с низкой вероятностью предсказания считались неприсвоенными. На рис. 2б показан пример установки порога вероятности на уровне 90%. Установив этот порог, мы увеличили долю действительно отнесенных образцов с 96,1 до 99,6%, а уровень ошибочной классификации снизился с 3,9 до 0,4%. Однако 407 тестовых особей остались «неприсвоенными», например, 22 из 23 пчел A. m. cecropia, предсказанных как A. m. macedonica, больше не считались ошибочно классифицированными, а попали в категорию неприсвоенных.

Обсуждение

В этом исследовании мы провели крупномасштабную и всестороннюю выборку, следуя стандартизированной процедуре, и поставили цель охватить как можно больше генетического разнообразия медоносных пчел в Европе путем глубокого секвенирования объединенных популяций. Далее мы применили два мощных метода отбора SNP [32, 33] для изучения разнообразия на разных уровнях дифференциации (линии, подвиды, популяции). Впоследствии эти информативные маркеры предков были использованы для построения модели классификации образцов европейских медоносных пчел на подвиды.

Значительное разнообразие медоносных пчел представляет собой проблему, когда речь идет о создании дискриминационного инструмента, применимого во всей Европе. Четыре европейские линии были легко различимы генетически с помощью всего 200 SNPs из-за их древней дивергенции [52], но трудности возникли на более низком иерархическом уровне дифференциации. Подвиды из одной эволюционной линии разошлись совсем недавно [53] и, таким образом, генетически очень близки. Более того, в Европе есть районы, где подвидовая изменчивость A. mellifera еще не описана полностью, а в других районах интрогрессия, опосредованная человеком, способствует размыванию естественных границ между подвидами [42, 48, 54]. Национальные селекционные программы также могут нарушать естественный поток генов и способствовать изменению генетического фона исходного подвида [11, 12, 55, 56]. На самом деле, в нашем исследовании, применяя строгую фильтрацию, мы выявили лишь несколько уникальных SNP, которые относились исключительно к одной популяции. Аналогичным образом, другие исследования популяционной геномики обнаружили высокую степень обмена аллелями между эволюционными линиями и внутри них [7, 53]. Напротив, мы обнаружили различия в средней частоте вызовов на особь между подвидами, что частично может объясняться наличием нулевых аллелей (аллелей, не дающих сигнала), что свидетельствует об изменении последовательности или видоспецифических делециях в сайте зонда. Зонды, которые не сработали для определенных подвидов (т.е. отсутствующие данные), на самом деле содержат ценную информацию и даже обогащают нашу модель.

Мы использовали подход машинного обучения (ML) для построения модели классификации подвидов. ML использует преимущества высокой размерности исходных данных и обеспечивает повышение точности предсказания в подходе без модели [39, 40]. Таким образом, можно выявить тонкие различия, что было особенно актуально в нашем исследовании из-за большого количества близкородственных подвидов, которые мы хотели различать. Нашей моделью, показавшей наилучшие результаты, была линейная SVC, член семейства опорных векторных машин (SVM), которые, как известно, хорошо обобщают, поскольку они разработаны для максимизации разницы между любыми двумя классами (подвидами) [57]. Типичные биологические приложения SVM включают предсказание функций белков, предсказание сайтов инициации транскрипции и классификацию данных экспрессии генов (обзор приведен в 57). В области популяционной генетики тщательный ML-подход для выбора наилучшей модели, как правило, еще не получил широкого распространения, хотя для вывода происхождения были разработаны специальные модели [58, 59]. Здесь мы используем комплексный подход ML на основе данных о генотипах для диагностики подвидов медоносной пчелы.

Несмотря на всестороннюю выборку, тщательный отбор SNP и применение новейших методов классификации, в диагностической системе остаются некоторые ограничения. Например, в рамках C-линейки мы столкнулись с проблемами при дифференциации образцов в соответствии с предполагаемым подвидом. Такая неправильная классификация особей может быть объяснена сочетанием различных факторов: (i) эта линия имеет сравнительно недавнее происхождение [53] и (ii) состоит из множества сильно взаимосвязанных подвидов в непосредственной географической близости (см. Рисунок S1D); (iii) таксономический статус некоторых популяций еще не полностью определен [60,61,62]; и (iv) генетический фон некоторых популяций изменяется в результате интрогрессии из-за вмешательства человека [63]. Кроме того, нельзя исключать ошибки маркировки образцов, не содержащих данных, как дополнительный источник неправильной классификации, особенно если мы имеем в виду те образцы, для которых модель предсказала другой подвид с высокой вероятностью. Контролируемый ML полагается на качество исходных данных для классификации, поэтому в будущем мы стремимся улучшить обучающие данные, чтобы повысить точность предсказания модели и снизить процент ошибочной классификации.

Важно также отметить, что, установив порог вероятности для отнесения к тому или иному подвиду, уровень ошибочной классификации был значительно снижен для некоторых подвидов. Хотя такой порог повышает уверенность в предсказании подвидов, он также подразумевает, однако, что довольно много особей осталось «неприсвоенными». Какой порог используется в качестве отсечки для классификации подвидов, зависит от конкретных обстоятельств и области применения. Например, для сохранения небольшой исчезающей популяции порог может быть установлен ниже, чтобы сохранить генетическое разнообразие, чем, например, для чистой селекционной линии, отбираемой по определенным признакам.

В целом, более ранние методы, основанные на морфометрии, вариациях мтДНК, микросателлитных локусах или даже SNPs, были эффективны для дифференциации между эволюционными линиями и, в некоторой степени, между подвидами одной линии [22, 42, 45, 64,65,66,67]. Тем не менее, наш диагностический инструмент является наиболее полным на сегодняшний день инструментом, позволяющим надежно классифицировать европейских медоносных пчел на подвиды в рамках одного анализа. Более того, преимущество нашего подхода заключается в том, что это динамический инструмент, который может быть обновлен для включения большего количества подвидов путем генотипирования новых образцов и добавления их данных для перестройки классификационной модели с использованием ML с дополнительными подвидами. Текущие исследования показывают, что этот подход применим к A. m. siciliana из Сицилии. Более того, отдельные пчелы из Южной Африки, протестированные с помощью нашей системы, были отклонены как имеющие европейское происхождение (т.е. низкая вероятность предсказания принадлежности к любому из подвидов). Таким образом, этот динамический инструмент может легко включать новые популяции, подлежащие дискриминации, и даже может быть оптимизирован для дифференциации популяций/экотипов внутри подвидов или для оценки степени интрогрессии.

Выводы

Основной вывод исследования заключается в том, что наша модель может классифицировать пчел на каждый из европейских подвидов с высокой точностью. Следовательно, поскольку пчелы, включенные в этот проект, были собраны на обширной территории от России и Армении на востоке до Португалии на западе, и от Мальты на юге до Шотландии на севере, мы пришли к выводу, что большая часть природного разнообразия европейских медоносных пчел все еще может считаться сохранившейся, несмотря на вмешательство человека на протяжении более 150 лет. Сохранение in situ этого генетического наследия — наш долг [68], и мы считаем, что представленный инструмент диагностики подвидов медоносной пчелы внесет полезный вклад. Он будет полезен в самых разных областях: для пчеловодов, которые хотят знать подвид своих пчел; для менеджеров по охране природы в Европе, где подвидовая диагностика необходима для мониторинга уровня гибридизации колоний в консерваториях; для ветеринаров для контроля торговли матками; для пчеловодов для сертификации подвидового происхождения своих маток; и для пчеловодов для проверки подлинности их пчелопродукции.

Методы

Секвенирование образцов

Для данного исследования в общей сложности были отобраны образцы 22 популяций, все в пределах их родного ареала (Таблицы 1 и S1), которые относятся к различным подвидам и генетическим происхождениям в соответствии с классификацией Руттнера [8] и последним пересмотром рода Apis, проведенным Энгелем [62]: A. m. iberiensis Engel 1999, A. m. mellifera Linnaeus 1758, A. m. carnica Pollman 1879, A. m. caucasia Pollmann 1889, A. m. ligustica Spinola 1806, A. m. macedonica Ruttner 1988, A. m. cecropia Kiesenwetter 1860, A. m. cypria Pollman 1879, A. m. adami Ruttner 1975, A. m. anatoliaca Maa 1953, A. m. remipes Gerstaecker 1862; кроме того, мы включаем A. m. ruttneri Sheppard et al. 1997 [69], «A. m. carpatica» Foti 1965 [60] и «A. m. rodopica» Petrov 1991 [61]. Существует некоторая неопределенность и нерешенный таксономический статус некоторых популяций, а описания подвидов в литературе не всегда были выполнены в соответствии со стандартами, изложенными в Международном кодексе зоологической номенклатуры (ICZN) [62]. Таким образом, в литературе встречаются различные точки зрения на то, что следует считать подвидом или экотипом. В данной работе мы не ставим целью разрешить или обосновать какую-либо классификацию. Наконец, мы рассмотрели 14 подвидов/генетических истоков (перечисленных выше) для нашего диагностического инструмента, которые были использованы в качестве категорий в модели классификации машинного обучения.

Каждая отобранная популяция включала до 100 (от 86 до 100) рабочих пчел из неродственных колоний, которые использовались для последующего секвенирования пула. Были предприняты усилия для охвата всего ареала распространения любого подвида, при этом учитывалась внутривидовая изменчивость, когда это было необходимо. Мы сосредоточились на сборе репрезентативных образцов для каждого подвида, отбирая образцы в основном у пчеловодов, которые, как известно, не импортируют пчел, чтобы минимизировать риск включения гибридов. Кроме того, мы выбирали только одну рабочую пчелу на пасеку, чтобы избежать родственных особей и включить как можно больше разнообразия в популяцию. Также для того, чтобы убедиться в подвидовой принадлежности собранных образцов, в некоторых случаях (где это было возможно) проводился морфометрический анализ и/или мы использовали уже генотипированных пчел [55, 65, 66, 70, 71, 72]. Подробная информация о происхождении образцов и соответствующие ссылки представлены в таблице S1.

Выделение ДНК, подготовка библиотек и секвенирование пулов

Головы или тораксы до 100 пчел (Таблица S1) из каждого пула были гомогенизированы, ДНК была выделена из всех образцов с помощью метода очистки на основе магнитных бусинок (NucleoMag® Blood 100 мкл, Macherey-Nagel, Германия). Затем из каждого пула ДНК были созданы библиотеки секвенирования с помощью набора для подготовки библиотек TruSeq DNA PCR-Free и секвенированы на платформе Illumina HiSeq 2500. Биоинформационная обработка, включая обрезку, картирование и поиск вариантов сгенерированных данных о последовательности пула, проводилась с использованием передового опыта и стандартного программного обеспечения (подробности в дополнительных материалах и методах). Конвейер для анализа данных пула последовательностей доступен на сайте https://github.com/jlanga/smsk_popoolation.

Выбор информативных маркеров предков

В ряде исследований было отобрано ограниченное число SNP для дифференциации основных эволюционных линий [15, 45, 46], однако для близкородственных подвидов необходимо большее число маркеров и более тонкий подход к отбору. Таким образом, мы использовали два различных подхода (PCA и FST) [28, 34] для определения и отбора информативных SNP, чтобы выявить наиболее дискриминантные маркеры на разных уровнях: (i) SNP для дифференциации четырех основных эволюционных линий, (ii) SNP для дифференциации подвидов внутри эволюционных линий, и (iii) SNP для идентификации конкретных популяций внутри подвидов (например, экотипов).

Во-первых, мы создали матрицу с минорными частотами аллелей для каждого SNP и пула секвенирования, которая использовалась для проведения PCA для отбора SNP, различающих основные эволюционные линии (Рисунок S1A). Во-вторых, PCA был выполнен отдельно для подмножеств пулов из каждой линии, чтобы выбрать информативные SNP для дискриминации подвидов внутри каждой линии (Рисунок S1B-D). Мы использовали пакет FactoMineR R [73] и собственные скрипты R для отбора на каждом иерархическом уровне SNPs с наибольшим вкладом в значимые PCs. С помощью этой процедуры было отобрано 300 PCA-информативных SNP для различения четырех эволюционных линий, 200 SNP для M-линии, 600 SNP для O-линии и 1100 SNP для наиболее сложной C-линии (рис. S1D). Предварительное моделирование с использованием частот аллелей из пула секвенирования показало, что этот подход был высокоэффективным для различения линий и подвидов (Рисунок S1).

Для выбора дополнительных SNP, которые могут различать пулы, парные значения FST [74] между всеми популяциями были рассчитаны для каждого SNP с двумя настройками (свободный и строгий варианты) с использованием PoPoolation2 [75]. При свободном варианте настроек будет получено больше SNP с меньшей определенностью и более низким качеством, что, в свою очередь, потенциально снижает успешность генотипирования. Этот недостаток уравновешивается тем, что свободный вариант повышает вероятность выявления высокоинформативных популяционно-специфических (уникальных) SNP. При любом варианте настройки (свободном и строгом) парные значения FST каждого пула по отношению ко всем другим пулам суммировались для каждого SNP, а затем ранжировались по наибольшему суммарному значению FST. Ожидается, что фиксированный и уникальный SNP в одном пуле будет иметь максимальную сумму 21, что означает, что этот вариант присутствует только в данной конкретной популяции. Разумный компромисс между уникальными и надежными SNP был достигнут путем отбора 20 лучших SNP с наибольшей суммарной FST из свободного варианта и 80 лучших SNP из строгого варианта для каждого пула. В 22 пулах с помощью F_ST было отобрано 2200 информативных SNP, специфичных для популяции.

В целом, на основе PCA и FST было отобрано 4400 информативных SNP (Таблица S3). Эти высокоинформативные маркеры не только важны для отнесения особей к подвидам, как это представлено в данном исследовании, но, благодаря варьированию частот аллелей в различных популяциях, они могут быть использованы, например, для классификации новых подвидов и для дальнейших последующих исследований.

Дизайн зондов

Зонды для 4400 отобранных SNP были оценены для генотипирования на платформе Illumina Infinium с использованием программного обеспечения DesignStudio® компании Illumina, которое требует в качестве входных данных фланкирующую область в 50 п.н. вверх и вниз по течению от каждого SNP. SNP отбрасывались, если не удавалось сконструировать зонд во фланкирующей области или если зонды имели более одного совпадения при выравнивании с эталонным геномом медоносной пчелы. Окончательный список из 4197 SNP был передан в компанию Illumina для разработки и производства зондов. SNPs распределены по всем 16 хромосомам медоносной пчелы, а также в неразмещенных контигах (Таблица S3; Дополнительный файл 1), среднее расстояние между SNPs составляет 64 кб.

Валидационные образцы и генотипирование

В общей сложности 3958 особей пчел были генотипированы по выбранным SNP, включая 2050 тех же рабочих пчел, которые использовались для секвенирования пула, а также 1908 вновь собранных особей (Таблица 1). Эти новые дополнительные образцы были получены из нескольких различных источников и разного качества, включая целых пчел в этаноле, пчел, выжатых на FTA-карты, образцы тканей из летных мышц или очищенную ДНК. Они были получены с селекционных пасек проекта SmartBees [76] и из колоний, исследованных на чувствительное к варроа гигиеническое поведение в рамках проекта SmartBees [77]. Образцы были генотипированы с помощью изготовленного на заказ BeadChip массива Infinium iSelect XT 96. Результаты анализировались с помощью программного обеспечения GenomeStudio® компании Illumina, и генотипы каждого образца экспортировались для дальнейшего анализа. Для первоначальной визуализации результатов генотипирования мы построили графики многообразия t-distributed stochastic neighbor embedding (t-SNE). Эта техника визуализирует высокоразмерные данные, задавая каждой точке данных местоположение на двух- или трехмерной карте [78]. Выбросы и образцы, которые были помечены как один подвид, но были явно сгруппированы с другим кластером, были удалены, в общей сложности 62 образца, оставив N = 1988 эталонных образцов последовательности пула. Это было сделано с целью создания высококачественного и репрезентативного набора эталонных данных для выделения подвидов.

Классификация образцов с использованием алгоритмов машинного обучения (ML)

Чтобы построить модель для классификации и прогнозирования видовой принадлежности неизвестных образцов европейских медоносных пчел, мы использовали методы ML с помощью среды scikit-learn python [79]. Сначала 1988 генотипированных особей из пулов были перемешаны, затем 70% из них (N = 1391) были использованы в качестве обучающих данных. Оставшиеся 30% (N = 597), а также дополнительные вновь собранные особи (N = 1908) рассматривались как данные вне выборки (NTotal = 2505) для последующей валидации (Таблица 1) [40]. Были протестированы различные контролируемые алгоритмы ML, включая RandomForest, LogisticRegression, SupportVector Machine (SVM) и Linear SupportVectorClassifier (SVC) (Таблица S5; подробная информация о выборе модели в дополнительных материалах и методах). Вкратце, данные о генотипах были преобразованы в матрицу, совместимую с машинным обучением (кодирование one-hot) [80]. Информация о классе, такая как родословная и подвид каждого образца, была добавлена в матрицу, которая использовалась для обучения различных моделей машинного обучения для предсказания происхождения образца. Линейная SVC была одной из лучших моделей по средней точности, оцененной с помощью перекрестной валидации, и была выбрана (Таблица S5, Рисунок S5).

После обучения линейной модели SVC она была использована для классификации данных вне выборки (N = 2505). Образцы были классифицированы в соответствии с наибольшей вероятностью предсказания принадлежности к любому из подвидов. Для обобщения и визуализации результатов на вневыборочных данных, для которых известны истинные значения, была создана матрица смешения [81]. Каждая строка матрицы представляет истинный класс, а каждый столбец — предсказанный класс, основанный на наибольшей вероятности для каждого подвида. Полученные проценты сравнивают список ожидаемых значений со списком предсказаний модели.

Для того чтобы модель можно было применять в практической деятельности по сохранению и разведению, мы определили порог в 90% на основе наблюдаемого распределения вероятностей предсказания (Рисунок S3), которые соответствуют значениям, найденным в литературе по пчеловодству [43, 82]. Если вероятность предсказания для любого образца меньше порога в 90%, он считался «неприсвоенным», в то время как если она превышала порог, его относили к соответствующему подвиду.

Доступность данных и материалов

Все данные последовательности из пулов, проанализированных в ходе настоящего исследования, были переданы в архив коротких чтений NCBI (SRA) под номером доступа BioProject PRJNA666033: https://www.ncbi.nlm.nih.gov/sra/?term=PRJNA666033. Схема анализа данных последовательности пулов доступна на сайте https://github.com/jlanga/smsk_popoolation.

Сокращения

ДНК: дезоксирибонуклеиновая кислота

FTA: Flinders Technology Associates

МКЗН: Международный кодекс зоологической номенклатуры

ML: Машинное обучение

мтДНК: Митохондриальная дезоксирибонуклеиновая кислота

РСА: Анализ главных компонент

ПЦР: Полимеразная цепная реакция

SD: Стандартное отклонение

SNP: Однонуклеотидный полиморфизм

SVC: Классификатор опорных векторов

SVM: Машинка опорных векторов

t-SNE: t-распределенное стохастическое соседнее встраивание

Дополнительная информация

Дополнительный файл 1.
Информативные маркеры подвидов медоносной пчелы. Информация о SNP и геномное положение 4165 SNP, отобранных в данном исследовании для дифференциации европейских подвидов медоносной пчелы.

Дополнительный файл 2: Дополнительные материалы и методы.
Дополнительные материалы и методы, подробно описывающие используемые наборы данных, лабораторные методы, биоинформационный конвейер, подход к выбору SNP и классификацию образцов с помощью алгоритмов машинного обучения.

Дополнительный файл 3: Рисунок S1.
PCA-полосы с отобранными по PCA SNPs и 100 смоделированными особями на основе частот аллелей пулов. (A) Используя 300 SNPs, эволюционные линии M, C и O были хорошо разделены с помощью первых двух PC, в то время как линия A может быть дифференцирована с помощью третьего компонента (не показано). (B) В линии M, в которой представлены только два европейских подвида (A. m. mellifera и A. m. iberiensis), первый ПК с использованием 200 SNPs уже дает достаточно информации. (C) В родословной O представлены четыре подвида, которые разделяются с помощью 600 SNP, используя 4 ПК. (D) Для линии C, содержащей множество подвидов (комплекс подвидов), 1100 SNP были выбраны для получения лучшего разрешения.

Дополнительный файл 4: Рисунок S2.
Визуализация с помощью многоугольного графика t-SNE образцов 1988 медоносных пчел из пула секвенирования, индивидуально генотипированных по 4094 SNPs. Образцы выделены цветом в соответствии с названием пула, которое представляет подвид и страну происхождения, как указано в таблице 1.

Дополнительный файл 5: Рисунок S3.
Гистограмма вероятностей предсказания для данных вне выборки. Порог присвоения 90% включает 2098 из 2505 образцов (=84%).

Дополнительный файл 6: Рисунок S4.
Матрица путаницы для предсказания эволюционной линии по данным вне выборки: Африканская линия (A), Центрально- и Восточноевропейская линия (C), Западно- и Североевропейская линия и Ближневосточная и Центральноазиатская линия (O). Каждая строка матрицы представляет истинный класс (род), а каждый столбец — предсказанный класс, основанный на наибольшей вероятности предсказания. Полученные проценты сравнивают список ожидаемых значений со списком предсказаний модели.

Дополнительный файл 7: Рисунок S5.
Кривая обучения для наиболее эффективной модели, линейной SVC, со средним и стандартным отклонением оценок 10-кратной перекрестной валидации.

Дополнительный файл 8: Таблица S1.
Дополнительная информация о выборке, ссылки и происхождение образцов для секвенирования пула.

Дополнительный файл 9: Таблица S2.
Пул секвенирования и статистика вариантов.

Дополнительный файл 10: Таблица S3.
Распределение выбранных информативных SNP по геному медоносной пчелы.

Дополнительный файл 11: Таблица S4.
Средняя частота вызовов на образец для каждого подвида.

Дополнительный файл 12: Таблица S5.
Статистика точности для различных протестированных моделей, оцененная с помощью 10-кратной перекрестной валидации.

Ссылки

1.Potts SG, Biesmeijer JC, Kremen C, Neumann P, Schweiger O, Kunin WE. Global pollinator declines: trends, impacts and drivers. Trends Ecol Evol. 2010;25:345–53.
Article PubMed Google Scholar 

2.IPBES. Summary for policymakers of the assessment report of the intergovernmental science-policy platform on biodiversity and ecosystem services on pollinators, pollination and food production. Bonn: Secretariat of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services; 2016.
Google Scholar 

3.Dogantzis KA, Zayed A. Recent advances in population and quantitative genomics of honey bees. Curr Opin Insect Sci. 2019;31:93–8.
Article PubMed Google Scholar 

4.Chen C, Liu Z, Pan Q, Chen X, Wang H, Guo H, et al. Genomic analyses reveal demographic history and temperate adaptation of the newly discovered honey bee subspecies Apis mellifera sinisxinyuan n. ssp. Mol Biol Evol. 2016;33:1337–48.
Article CAS PubMed PubMed Central Google Scholar 

5.Sheppard WS, Meixner MD. Apis mellifera pomonella, a new honey bee subspecies from Central Asia. Apidologie. 2003;34:367–75.
Article Google Scholar 

6.Meixner MD, Leta MA, Koeniger N, Fuchs S. The honey bees of Ethiopia represent a new subspecies of Apis mellifera—Apis mellifera simensis n. ssp. Apidologie. 2011;42(3):425–37.
Article Google Scholar 

7.Cridland JM, Tsutsui ND, Ramírez SR. The complex demographic history and evolutionary origin of the western honey bee, Apis mellifera. Genome Biol Evol. 2017;9(2):457–72.
Article PubMed PubMed Central Google Scholar 

8.Ruttner F. Biogeography and taxonomy of honeybees. Berlin: Springer Verlag; 1988.
Book Google Scholar 

9.De La Rúa P, Jaffé R, Dall’olio R, Muñoz I, Serrano J. Biodiversity, conservation and current threats to European honeybees. Apidologie. 2009;40:263–84.
Article Google Scholar 

10.Pinto MA, Henriques D, Chávez-Galarza J, Kryger P, Garnery L, van der Zee R, et al. Genetic integrity of the dark European honey bee (Apis mellifera mellifera) from protected populations: a genome-wide assessment using SNPs and mtDNA sequence data. J Apic Res. 2014;53(2):269–78.
Article Google Scholar 

11.Bouga M, Harizanis PC, Kilias G, Alahiotis S. Genetic divergence and phylogenetic relationships of honey bee Apis mellifera (Hymenoptera: Apidae) populations from Greece and Cyprus using PCR–RFLP analysis of three mtDNA segments. Apidologie. 2005;36(3):335–44.
Article CAS Google Scholar 

12.Dall’Olio R, Marino A, Lodesani M, Moritz RF. Genetic characterization of Italian honeybees, Apis mellifera ligustica, based on microsatellite DNA polymorphisms. Apidologie. 2007;38(2):207–17.
Article CAS Google Scholar 

13.Büchler R, Costa C, Hatjina F, Andonov S, Meixner MD, Conte YL, et al. The influence of genetic origin and its interaction with environmental effects on the survival of Apis mellifera L. colonies in Europe. J Apic Res. 2014;53(2):205–14.
Article Google Scholar 

14.Frankham R, Ballou JD, Briscoe DA. Introduction to conservation genetics. Cambridge: Cambridge University Press; 2002.
Book Google Scholar 

15.Parejo M, Wragg D, Gauthier L, Vignal A, Neumann P, Neuditschko M. Using whole-genome sequence information to Foster conservation efforts for the European dark honey bee, Apis mellifera mellifera. Front Ecol Evol. 2016;4:140.
Article Google Scholar 

16.Bertrand B, Alburaki M, Legout H, Moulin S, Mougel F, Garnery L. MtDNA COI-COII marker and drone congregation area: an efficient method to establish and monitor honeybee (Apis mellifera L.) conservation centres. Mol Ecol Resour. 2015;15:673–83.
Article CAS PubMed Google Scholar 

17.Büchler R, Uzunov A. Selecting for Varroa resistance in German honey bees. Bee World. 2016;93(2):49–52.
Article Google Scholar 

18.Uzunov A, Brascamp EW, Büchler R. The basic concept of honey bee breeding programs. Bee World. 2017;94(3):84–7.
Article Google Scholar 

19.Baudry E, Solignac M, Garnery L, Gries M, Cornuet J, Koeniger N. Relatedness among honeybees (Apis mellifera) of a drone congregation. Proc R Soc Lond B Biol Sci. 1998;265(1409):2009–14.
Article Google Scholar 

20.Tarpy DR, Delaney DA, Seeley TD. Mating frequencies of honey bee queens (Apis mellifera L.) in a population of feral colonies in the northeastern United States. PLoS One. 2015;10(3):e0118734.
Article PubMed PubMed Central CAS Google Scholar 

21.Büchler R, Andonov S, Bienefeld K, Costa C, Hatjina F, Kezic N, et al. Standard methods for rearing and selection of Apis mellifera queens. J Apic Res. 2013;52(1):1–30.
Article Google Scholar 

22.Bouga M, Alaux C, Bienkowska M, Büchler R, Carreck NL, Cauia E, et al. A review of methods for discrimination of honey bee populations as applied to European beekeeping. J Apic Res. 2011;50(1):51–84.
Article Google Scholar 

23.Lodesani M, Costa C. Bee breeding and genetics in Europe. Bee World. 2003;84(2):69–85.
Article Google Scholar 

24.Muñoz I, Pinto MA, De la Rúa P. Temporal changes in mitochondrial diversity highlights contrasting population events in Macaronesian honey bees. Apidologie. 2013;44(3):295–305.
Article CAS Google Scholar 

25.Miguel I, Garnery L, Iriondo M, Baylac M, Manzano C, Steve Sheppard W, et al. Origin, evolution and conservation of the honey bees from La Palma Island (Canary Islands): molecular and morphological data. J Apic Res. 2015;54(5):427–40.
Article Google Scholar 

26.Soares S, Grazina L, Mafra I, Costa J, Pinto MA, Oliveira MBP, et al. Towards honey authentication: differentiation of Apis mellifera subspecies in European honeys based on mitochondrial DNA markers. Food Chem. 2019;283:294–301.
Article CAS PubMed Google Scholar 

27.Rosenberg NA, Li LM, Ward R, Pritchard JK. Informativeness of genetic markers for inference of ancestry. Am J Hum Genet. 2003;73:1402–22.
Article CAS PubMed PubMed Central Google Scholar 

28.Paschou P, Ziv E, Burchard EG, Choudhry S, Rodriguez-Cintron W, Mahoney MW, et al. PCA-correlated SNPs for structure identification in worldwide human populations. PLoS Genet. 2007;3:1672–86.
Article CAS PubMed Google Scholar 

29.Lewis J, Abas Z, Dadousis C, Lykidis D, Paschou P, Drineas P. Tracing cattle breeds with principal components analysis ancestry informative SNPs. PloS One. 2011;6(4):e18007.
Article CAS PubMed PubMed Central Google Scholar 

30.Heaton MP, Leymaster KA, Kalbfleisch TS, Kijas JW, Clarke SM, McEwan J, et al. SNPs for parentage testing and traceability in globally diverse breeds of sheep. PloS One. 2014;9(4):e94851.
Article PubMed PubMed Central CAS Google Scholar 

31.Montes I, Conklin D, Albaina A, Creer S, Carvalho GR, Santos M, et al. SNP discovery in European anchovy (Engraulis encrasicolus, L) by high-throughput transcriptome and genome sequencing. PLoS One. 2013;8(8):e70051.
Article CAS PubMed PubMed Central Google Scholar 

32.Wilkinson S, Wiener P, Archibald AL, Law A, Schnabel RD, McKay SD, et al. Evaluation of approaches for identifying population informative markers from high density SNP chips. BMC Genet. 2011;12(1):45.
Article PubMed PubMed Central Google Scholar 

33.Ding L, Wiener H, Abebe T, Altaye M, Go RC, Kercsmar C, et al. Comparison of measures of marker informativeness for ancestry and admixture mapping. BMC Genomics. 2011;12(1):622.
Article CAS PubMed PubMed Central Google Scholar 

34.Holsinger KE, Weir BS. Genetics in geographically structured populations: defining, estimating and interpreting F ST. Nat Rev Genet. 2009;10(9):639.
Article CAS PubMed PubMed Central Google Scholar 

35.Montes I, Laconcha U, Iriondo M, Manzano C, Arrizabalaga H, Estonba A. Reduced single nucleotide polymorphism panels for assigning Atlantic albacore and Bay of Biscay anchovy individuals to their geographic origin: toward sustainable fishery management. J Agric Food Chem. 2017;65(21):4351–8.
Article CAS PubMed Google Scholar 

36.Liu N, Zhao H. A non-parametric approach to population structure inference using multilocus genotypes. Hum Genomics. 2006;2:353–64.
Article CAS PubMed PubMed Central Google Scholar 

37.Paetkau D, Calvert W, Stirling I, Strobeck C. Microsatellite analysis of population structure in Canadian polar bears. Mol Ecol. 1995;4(3):347–54.
Article CAS PubMed Google Scholar 

38.Piry S, Alapetite A, Cornuet JM, Paetkau D, Baudouin L, Estoup A. GENECLASS2: a software for genetic assignment and first-generation migrant detection. J Hered. 2004;95(6):536–9.
Article CAS PubMed Google Scholar 

39.Schrider DR, Kern AD. Supervised machine learning for population genetics: a new paradigm. Trends Genet. 2018;34(4):301–12.
Article CAS PubMed PubMed Central Google Scholar 

40.Tarca AL, Carey VJ, Chen XW, Romero R, Drăghici S. Machine learning and its applications to biology. PLoS Comput Biol. 2007;3(6):116.
Article CAS Google Scholar 

41.Guinand B, Topchy A, Page KS, Burnham-Curtis MK, Punch WF, Scribner KT. Comparisons of likelihood and machine learning methods of individual classification. J Hered. 2002;93(4):260–9.
Article CAS PubMed Google Scholar 

42.Muñoz I, Henriques D, Johnston JS, Chávez-Galarza J, Kryger P, Pinto MA. Reduced SNP panels for genetic identification and introgression analysis in the dark honey bee (Apis mellifera mellifera). PLoS One. 2015;10:e0124365.
Article PubMed PubMed Central CAS Google Scholar 

43.Parejo M, Henriques D, Pinto MA, Soland-Reckeweg G, Neuditschko M. Empirical comparison of microsatellite and SNP markers to estimate introgression in Apis mellifera mellifera. J Apic Res. 2018;57(4):504–6.
Article Google Scholar 

44.Garnery L, Franck P, Baudry E, Vautrin D, Cornuet JM, Solignac M. Genetic diversity of the west European honey bee (Apis mellifera mellifera) and A. M Iber II Microsatellite. Loci Genet Sel Evol. 1998;30:49–74.
Article Google Scholar 

45.Henriques D, Parejo M, Vignal A, Wragg D, Wallberg A, Webster MT, et al. Developing reduced SNP assays from whole-genome sequence data to estimate introgression in an organism with complex genetic patterns, the Iberian honeybee (Apis mellifera iberiensis). Evol Appl. 2018;11(8):1270–82.
Article CAS PubMed PubMed Central Google Scholar 

46.Chapman NC, Harpur BA, Lim J, Rinderer TE, Allsopp MH, Zayed A, et al. A SNP test to identify Africanized honeybees via proportion of ‘African’ancestry. Mol Ecol Resour. 2015;15(6):1346–55.
Article CAS PubMed Google Scholar 

47.Rortais A, Arnold G, Alburaki M, Legout H, Garnery L. Review of the DraI COI-COII test for the conservation of the black honeybee (Apis mellifera mellifera). Conserv Genet Resour. 2011;3:383–91.
Article Google Scholar 

48.Jensen AB, Palmer KA, Boomsma JJ, Pedersen BV. Varying degrees of Apis mellifera ligustica introgression in protected populations of the black honeybee, Apis mellifera mellifera, in Northwest Europe. Mol Ecol. 2005;14:93–106.
Article PubMed Google Scholar 

49.Harpur BA, Chapman NC, Krimus L, Maciukiewicz P, Sandhu V, Sood K, et al. Assessing patterns of admixture and ancestry in Canadian honey bees. Insect Soc. 2015;62(4):479–89.
Article Google Scholar 

50.Chapman NC, Bourgeois AL, Beaman LD, Lim J, Harpur BA, Zayed A, et al. An abbreviated SNP panel for ancestry assignment of honeybees (Apis mellifera). Apidologie. 2017;48(6):776–83.
Article CAS Google Scholar 

51.Platzer A. Visualization of SNPs with t-SNE. PLoS One. 2013;8(2):56883.
Article CAS Google Scholar 

52.Whitfield CW, Behura SK, Berlocher SH, Clark AG, Johnston JS, Sheppard WS, et al. Thrice out of Africa: ancient and recent expansions of the honey bee, Apis mellifera. Science. 2006;314:642–5.
Article CAS PubMed Google Scholar 

53.Wallberg A, Han F, Wellhagen G, Dahle B, Kawata M, Haddad N, et al. A worldwide survey of genome sequence variation provides insight into the evolutionary history of the honeybee Apis mellifera. Nat Genet. 2014;46(10):1081–8.
Article CAS PubMed Google Scholar 

54.Franck P, Garnery L, Celebrano G, Solignac M, Cornuet JM. Hybrid origins of honeybees from Italy (Apis mellifera ligustica) and Sicily (A. m. sicula). Mol Ecol. 2000;9:907–21.
Article CAS PubMed Google Scholar 

55.Uzunov A, Meixner MD, Kiprijanovska H, Andonov S, Gregorc A, Ivanova E, et al. Genetic structure of Apis mellifera macedonica in the Balkan peninsula based on microsatellite DNA polymorphism. J Apic Res. 2014;53:288–95.
Article Google Scholar 

56.Nedić N, Francis RM, Stanisavljević L, Pihler I, Kezić N, Bendixen C, et al. Detecting population admixture in honey bees of Serbia. J Apic Res. 2014;53(2):303–13.
Article Google Scholar 

57.Yang ZR. Biological applications of support vector machines. Brief Bioinform. 2004;5(4):328–38.
Article CAS PubMed Google Scholar 

58.Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19:1655–64.
Article CAS PubMed PubMed Central Google Scholar 

59.Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Article CAS PubMed PubMed Central Google Scholar 

60.Foti N, Lungu M, Pelimon C, Barac I, Copaitici M, Marza E. Researches on morphological characteristics and biological features of the bee population in Romania. In: Proceedings of XXth Jubiliar International Congress of Beekeeping Apimondia; 1965. p. 171–176.
Google Scholar 

61.Petrov P. Systematics of Bulgarian bees. Pchelarstvo. 1991;9:15–7.
Google Scholar 

62.Engel MS. The taxonomy of recent and fossil honey bees (Hymenoptera: Apidae; Apis). J Hymenopt Res. 1999;8:165–96.
Google Scholar 

63.Ivanova E, Bouga M, Staykova T, Mladenovic M, Rasic S, Charistos L, et al. The genetic variability of honey bees from the southern Balkan Peninsula, based on alloenzymic data. J Apic Res. 2012;51(4):329–35.
Article CAS Google Scholar 

64.Garnery L, Cornuet JM, Solignac M. Evolutionary history of the honey bee Apis mellifera inferred from mitochondrial DNA analysis. Mol Ecol. 1992;1:145–54.
Article CAS PubMed Google Scholar 

65.Francis RM, Kryger P, Meixner M, Bouga M, Ivanova E, Andonov S, et al. The genetic origin of honey bee colonies used in the COLOSS genotype-environment interactions experiment: a comparison of methods. J Apic Res. 2014;53(2):188–204.
Article Google Scholar 

66.Meixner MD, Pinto MA, Bouga M, Kryger P, Ivanova E, Fuchs S. Standard methods for characterising subspecies and ecotypes of Apis mellifera. J Apic Res. 2013;52:1–28.
Article Google Scholar 

67.Muñoz I, Henriques D, Jara L, Johnston JS, Chávez-Galarza J, De La Rúa P, et al. SNP s selected by information content outperform randomly selected microsatellite loci for delineating genetic identification and introgression in the endangered dark European honeybee (Apis mellifera mellifera). Mol Ecol Resour. 2017;17(4):783–95.
Article PubMed CAS Google Scholar 

68.Jenks DT. The convention on biological diversity–an efficient framework for the preservation of life on earth. Nw J Intl Bus. 1994;15:636.
Google Scholar 

69.Sheppard WS, Arias MC, Grech A, Meixner MD. Apis mellifera ruttneri, a new honey bee subspecies from Malta. Apidologie. 1997;28(5):287–93.
Article Google Scholar 

70.Chávez-Galarza J, Henriques D, Johnston JS, Azevedo JC, Patton JC, Muñoz I, et al. Signatures of selection in the I berian honey bee (Apis mellifera iberiensis) revealed by a genome scan analysis of single nucleotide polymorphisms. Mol Ecol. 2013;22(23):5890–907.
Article PubMed CAS Google Scholar 

71.Miguel I, Iriondo M, Garnery L, Sheppard WS, Estonba A. Gene flow within the M evolutionary lineage of Apis mellifera: role of the Pyrenees, isolation by distance and post-glacial re-colonization routes in the western Europe. Apidologie. 2007;38(2):141–55.
Article CAS Google Scholar 

72.Ilyasov RA, Poskryakov AV, Petukhov AV, Nikolenko AG. Molecular genetic analysis of five extant reserves of black honeybee Apis mellifera mellifera in the Urals and the Volga region. Russ J Genet. 2016;52(8):828–39.
Article CAS Google Scholar 

73.Lê S, Josse J, Husson F. FactoMineR: an R package for multivariate analysis. J Stat Softw. 2008;25(1):1–18.
Article Google Scholar 

74.Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70.
CAS PubMed Google Scholar 

75.Kofler R, Pandey RV, Schlötterer C. PoPoolation2: identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq). Bioinformatics. 2011;27(24):3435–6.
Article CAS PubMed PubMed Central Google Scholar 

76.Büchler R, Uzunov A, Hoppe A, Bienefeld K. Field testing and selection on European honey bee populations (Smartbees project 2015-2018). In: Abstract book of the 8th EurBee congress of Apidology. Ghent: Ghent University; 2018.

77.Farajzadeh L, Wegener J, Momeni J, Nielsen O, Bienefeld K, Bendixen C. Whole-genome analysis of uncapping behaviour of individual honey bees towards Varroa destructor-parasitized brood. In: Proceedings of the 46th International Apicultural Congress. Montreal: Apimondia; 2019. p. 8–12.

78.Maaten LVD, Hinton G. Visualizing data using t-SNE. J Mach Learn Res. 2008;9:2579–605.
Google Scholar 

79.Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: machine learning in Python. J Mach Learn Res. 2011;12:2825–30.
Google Scholar 

80.Yang KK, Wu Z, Bedbrook CN, Arnold FH. Learned protein embeddings for machine learning. Bioinformatics. 2018;34(15):2642–8.
Article CAS PubMed PubMed Central Google Scholar 

81.Talbot J, Lee B, Kapoor A, Tan DS. EnsembleMatrix: interactive visualization to support machine learning with multiple classifiers. In: Proceedings of the SIGCHI Conference on Human Factors in Computing Systems. New York: ACM; 2009. p. 1283–1292.

82.Henriques D, Browne KA, Barnett MW, Parejo M, Kryger P, Freeman TC, et al. High sample throughput genotyping for estimating C-lineage introgression in the dark honeybee: an accurate and cost-effective SNP-based tool. Sci Rep. 2018;8:8552.
Article PubMed PubMed Central CAS Google Scholar