МикроРНК — семейство коротких некодирующих РНК, осуществляющих пост-транскрипционную регуляцию экспрессии генов [1]. Необходимым условием для связывания микроРНК с мРНК является комплементарность seed-региона микроРНК (2-7 нуклеотиды с 5′-конца) с последовательностью мРНК-мишени [2]. В результате такого связывания происходит остановка трансляции мРНК или же ее деградация, причем вероятность деградации напрямую связана с количеством комплементарных связей за пределами seed-региона [2]. Хорошо известно, что молекулы микроРНК могут играть роль как опухолевых супрессоров, так и онкогенов для множества видов рака [3–5].

В ходе созревания микроРНК ферменты Drosha и Dicer могут неточно осуществлять обрезку шпильки примикроРНК, в результате чего образуются изоформы микроРНК, отличающиеся от канонической микроРНК несколькими нуклеотидами на концах молекулы [6]. Вариация длины микроРНК с 5′-конца играет особую роль, так как при изменениях смещается seed-регион; таким образом, 5′-изоформы микроРНК обладают другим набором мишеней даже в случае изменения длины на один нуклеотид.

Колоректальный рак (КРР) занимает третье место по частоте заболеваемости и второе по частоте летальных исходов среди онкологических заболеваний в мире [7]. Широко известны примеры участия микроРНК в механизмах прогрессии и метастизрования КРР. Например, семейство микроРНК miR-200 подавляет экспрессию генов ZEB1, ZEB2, кодирующих ключевые транскрипционные факторы (ТФ) для эпителиально-мезенхимного перехода (ЭМП) [8]. Соответственно, транскрипционное и/или эпигенетическое подавление  экспрессии miR-200 способствует ЭМП и метастазированию рака [8]. Профили экспрессии молекул микроРНК активно используют также для поиска диагностических и прогностических маркеров КРР [9]. На сегодняшний день изучение роли изоформ микроРНК при КРР ограничено изучением их уровней экспрессии в опухолевых и здоровых тканях [10], при этом, насколько нам известно, функциональную активность 5′-изоформ микроРНК при КРР исследователи не анализировали.

Большое количество микроРНК и на порядки большее число порождаемых ими регуляторных взаимодействий (в среднем для одной микроРНК предсказано около 200 мишеней [1]) требуют применения биоинформатических подходов; один из наиболее распространенных методов — анализ регуляторных сетей. В рамках данного подхода молекулам микроРНК и генам ставят в соответствие вершины сети, а паре взаимодействующих молекул  — ребро, соединяющее микроРНК и мишень [11]. Для построения регуляторных сетей традиционно используют два подхода: литературные базы данных взаимодействий и корреляционный анализ по выборке образцов с известной экспрессией мРНК и микроРНК. Нами ранее был разработан и программно реализован алгоритм miRGTF-net, позволяющий объединить эти два подхода и добавить в сеть другой класс регуляторных молекул — транскрипционных факторов [12]. Использование данного алгоритма позволяет максимально полно и достоверно описать ландшафт внутриклеточных взаимодействий в интересующем типе клеток/тканей.

Целью работы было выяснить функциональную роль 5′-изоформ микроРНК в тканях КРР, проанализировав их профиль экспрессии, предсказав мишени на основе нуклеотидных последовательностей и интегрировав полученные данные с данными активности ТФ с помощью алгоритма miRGTF-net.

МАТЕРИАЛЫ И МЕТОДЫ

Предсказание мишеней изоформ микроРНК

Последовательности шпилек при-микроРНК и канонические позиции их разрезания ферментами Drosha и Dicer извлекали из базы данных miRBase версии 21 (https:// www.mirbase.org). Для обозначения 5′-изоформ микроРНК использовали стандартную номенклатуру: число после вертикальной черты обозначает сдвиг позиции разрезания относительно канонического в направлении от 5′- к 3′-концу. Например, hsa-miR-10a-5p|+1 соответствует последовательности микроРНК hsa-miR-10a-5p без первого нуклеотида на 5′-конце молекулы.

Нуклеотидные последовательности микроРНК и их изоформ вводили в программы miRDB версии 6.0 [13] и TargetScan версии 7.2 [2] для определения мРНК-мишеней. Согласно рекомендациям разработчиком miRDB, выбирали предсказания с качеством связывания не менее 80. Число предсказаний TargetScan уравнивали с числом предсказаний miRDB, выбирая соответствующее число наиболее сильных взаимодействий для каждой изоформ микроРНК.

Сбор и анализ данных секвенирования мРНК и изоформ микроРНК

Публично доступные исходные данные секвенирования мРНК и изоформ микроРНК проекта The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) получали с портала GDC (https://portal.gdc.cancer.gov/). Данные нормировали с помощью пакета edgeR версии 3.30.0 [14], использовали алгоритм нормализации Trimmed Mean of M-values (TMM), в результате получили таблицы Reads Per Kilobase of transcript per Million mapped reads (TMM-RPKM) для экспрессии мРНК и Reads Per Million mapped reads (TMM-RPM) — для экспрессии микроРНК.

Таблицу экспрессии 5′-изоформ микроРНК сортировали по суммарной экспрессии в рассматриваемых образцах, после чего считали кумулятивную функцию распределения. Наименьшее число 5′-изоформ микроРНК, покрывающих 95% всех прочтений секвенирования, обозначали высокоэкспрессированными 5′-изоформами микроРНК и использовали для дальнейшего анализа.

Построение регуляторной сети взаимодействий изоформ микроРНК, их мишеней и транскрипционных факторов

Алгоритм miRGTF-net [12] использовали для построения регуляторной сети взаимодействий изоформ микроРНК, их мишеней и транскрипционных факторов. Основным преимуществом алгоритма является возможность интеграции данных экспрессии мРНК и изоформ микроРНК (TCGA-COAD) с биологическими базами данных:

• TRRUST версии 2 (https://www.grnpedia.org/trrust/): взаимодействия ТФ и генов;
• TransmiR версии 2 (http://www.cuilab.cn/transmir): взаимодействия ТФ и микроРНК;
• miRDB, TargetScan: взаимодействия 5′-изоформ микроРНК и их мишеней (см. выше);
• miRIAD (https://www.miriad-database.org): коэкспрессия генов-хозяев и их интронных микроРНК.

Использовали каноническую последовательность шагов алгоритма miRGTF-net. Вкратце, сначала строили сеть на основе взаимодействий из баз данных. Затем для каждого ребра рассчитывали коэффициент корреляции Спирмена по соответствующим значениям экспрессии TCGA-COAD. Ребра, соответствующие молекулам со слабо коррелированными уровнями экспрессии, удаляли (отсечку на абсолютное значение корреляции Спирмена выбирали по 0,9-квантили распределения корреляций). Кроме того, удаляли ребра, соединяющие положительно коррелированные изоформы микроРНК и их мишени и соединяющие отрицательно коррелированных геновхозяев с их интронными микроРНК.

Далее оценивали силу линейной зависимости между экспрессией каждой вершины и ее прямыми регуляторами. Соответствующие линейные модели строили с помощью гребневой регрессии. Качество моделей оценивали с помощью коэффициента детерминации R², силу и направление регуляции оценивали с помощью стандартизованных β-коэффициентов регрессии. Для фильтрации вершин и ребер сети использовали пороговые значения, установленные в пакете miRGTF-net по умолчанию: модуль β-коэффициента не менее 0,3, 90% наибольших значений коэффициентов детерминации. Таким образом, полученная сеть содержала вершины, соответствующие как регуляторам экспрессии, так и регулируемым генам и изоформам микроРНК.

Поиск сильно связных компонент в сети проводили с помощью пакета NetworkX версии 2.8 (https://networkx.org). Регуляторные сети визуализировали с помощью программ Gephi (https://gephi.org) и yED Graph Editor (yWorks GmbH; Германия).

Анализ обогащения по функциональной принадлежности

Для функциональной аннотации списков генов (мишеней 5′-изоформ микроРНК) использовали веб-сервис DAVID версии декабря 2021 г. [15] и аннотацию биологических путей Gene Ontology (GO) [16].

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

Профиль экспрессии изоформ микроРНК в образцах колоректального рака

Анализируемая выборка TCGA-COAD состояла из профилей экспрессии мРНК и 5′-изоформ микроРНК 426 образцов первичных опухолей КРР. Нами было выделено 55 высокоэкспрессированных изоформ микроРНК, 10 из которых составляли неканонические изоформы микроРНК (рис. 1). Были выявлены две неканонические 5′-изоформы микроРНК, на каждую из которых приходилось более 1% от тотальной экспрессии микроРНК в рассматриваемых образцах: hsa-miR-1925p|+1 (2,4%) и hsa-miR-10a-5p|+1 (1,3%).

Вариация последовательности микроРНК на ее 5′-конце меняет seed-регион молекулы, вследствие чего может меняться множество потенциальных геновмишеней. Последовательности найденных канонических и неканонических изоформ микроРНК были использованы для биоинформатического предсказания их мишеней. Как и ожидалось, пары изоформ микроРНК, отличающихся на один нуклеотид с 5′-конца, имели слабо пересекающиеся множества мишеней (таблица). Например, каноническая форма микроРНК hsa-miR-10a-5p и ее 5′-изоформа без первого нуклеотида имели всего 11 общих мишеней из 267 мишеней в объединении (4,1%). Максимальная доля общих мишеней была достигнута для микроРНК hsa-miR-29a-3p и ее более длинной изоформы: 246 общих мишеней, 788 мишеней в объединении (31,2%).

Регуляторная сеть взаимодействий изоформ микроРНК, их мишеней и транскрипционных факторов

Следующим шагом биоинформатического анализа было построение регуляторной сети взаимодействий в клетках КРР. Алгоритм miRGTF-net позволяет строить такие сети, интегрируя два типа данных: биологически обоснованные взаимодействия из баз данных и профили экспрессии мРНК и изоформ микроРНК в выборке образцов. Регуляторая сеть содержала взаимодействия четырех типов:

• ТФ, регулирующие экспрессию генов;
• ТФ, регулярующие экспрессию микроРНК;
• 5′-изоформы микроРНК, регулирующие экспрессию генов;
• коэкспрессию генов-хозяев и их интронных микроРНК.

Данные экспрессии мРНК и 5′-изоформ микроРНК в выборке TCGA-COAD использовали для выбора взаимодействий, подкрепленных значимой корреляцией в рассматриваемых образцах.

Построенная сеть состояла из 333 молекул: 24 5′-изоформы микроРНК, 166 ТФ и 143 гена, не кодирующих ТФ (рис. 2). Из 456 взаимодействий 42 соответствовали подавлению экспрессии геновмишеней изоформами микроРНК, 413 регуляции экспрессии генов и микроРНК с помощью ТФ и лишь одна пара ген-микроРНК соответствовала коэкспрессии генахозяина и интронной микроРНК: HOXB3 и hsa-mir-10a.

Наибольшее число антикоррелированных мишеней (семь) было найдено для канонической микроРНК hsamiR-148a-3p. Данный список состоял из известных онкогенов, включая регуляторы и маркеры пролиферации (CSF1, ETS1, FLT1, MEIS1, MITF, RUNX2, категория GO:0008284 “positive regulation of cell proliferation”) и молекулу из семейства интегринов ITGA5, участвующую в регуляции клеточной пролиферации, инвазии и миграции путем передачи сигнала в клетки [17]. Данная микроРНК тоже присутствовала в наибольшой компоненте сильной связности регуляторной сети (т. е. подсети, в которой существует ориентированный путь между двумя любыми вершинами), напрямую связанной с ЭМП и эстрогеновым сигнальным путем (рис. 3). Таким образом, микроРНК hsa-miR-148a-3p потенциально подавляет экспрессию проопухолевых генов, играя роль возможного опухолевого супрессора при КРР.

Второй по количеству регулируемых генов оказалась пара, состоящая из канонической микроРНК hsa-miR-203a-3p|0 и ее 5′-изоформы hsa-miR-203a-3p|+1. В то время как у данных молекул не было общих антикоррелированных мишеней, обе изоформы микроРНК выполняли единую функцию, потенциально подавляя экспрессию онкогенов. Так, экспрессия канонической микроРНК hsa-miR203a-3p|0 отрицательно коррелировала с экспрессией ТФ SNAI2, являющегося одним из драйверов ЭМП [18], а неканоническая 5′-изоформа hsa-miR-203a-3p|+1 предположительно регулировала экспрессию белка внеклеточного матрикса, кодируемого геном TNC, который также играет ключевую роль в ЭМП при КРР [19].

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

В настоящем исследовании с помощью биоинформатического анализа была рассмотрена функциональная активность 5′-изоформ микроРНК в злокачественных колоректальных опухолях. Показано, что изоформы микроРНК, отличающиеся одним нуклеотидом на 5′-конце молекулы, обладают слабо пересекающимися множествами генов-мишеней (31,2% максимум). В число наиболее активных регуляторов вошли hsa-miR-148a-3p (каноническая микроРНК), hsa-miR-203a-3p|0 (каноническая микроРНК) и ее 5′-изоформа hsa-miR-203a-3p|+1. Интересно, что все три найденные микроРНК предположительно подавляли экспрессию проопухолевых генов, причем множества антикоррелированных мишеней канонической формы и 5′-изоформы miR-203a не пересекались.

Функциональную активность 5′-изоформ микроРНК ранее изучали в контексте рака молочной железы. Показано, что две 5′-изоформы микроРНК hsa-miR-183-5p оказывают различное влияние на транскриптом клеток MDA-MB-231, в частности, были найдены гены, опосредованно регулируемые изоформами в разных направлениях (EGFR, NRAS) [20]. По нашим данным, анализ мишени 5′-изоформ микроРНК при КРР проведен впервые.

Четыре из семи отобранных нами потенциальных мишеней hsa-miR-148a-3p были валидированы in vitro в ранее опубликованных исследованиях: CSF1, ITGA5 [21], MITF [22], RUNX2 [23]. Нами ранее было также обнаружено, что гипоксия клеточных линий КРР HT-29 и Caco-2 приводила к подавлению экспрессии микроРНК hsa-miR148a-3p, что влекло за собой повышение экспрессии гена-мишени ITGA5 [24]. В ряде других работ показано, что miR-148a оказывает проапоптотическое действие и ингибирует пролиферацию, миграцию и инвазию клеток КРР, подавляя экспрессию Bcl-2 [25], ErbB3 [26] и WNT10b [27]. Помимо КРР, роль hsa-miR-148a-3p в качестве супрессора пролиферации опухолевых клеток была показана в контексте раков молочной железы, простаты и уротелия [28]. Таким образом, полученные нами данные о противоопухолевой роли miR-148a хорошо согласуются с существующей литературой.

Аналогичную картину можно наблюдать для канонической формы микроРНК miR-203a: предсказанное нами взаимодействие miR-203a и SNAI2 было ранее валидировано in vitro [29], оверэкспрессия miR-203a в клеточных линиях КРР приводила к ингибированию инвазии и миграции клеток [30]. Таким образом, выявленные нами потенциальные мишени неканонической изоформы микроРНК hsa-miR-203a-3p|+1 согласуются с известными фактами о функциональной активности канонической микроРНК.

ВЫВОДЫ

Применение методов построения и анализа регуляторных сетей взаимодействий позволило нам определить роль функциональной активности некоторых 5′-изоформ микроРНК в клетках колоректального рака. Показано, что 5′-изоформы микроРНК hsa-miR-203a-3p могут регулировать различные гены-мишени, играя при этом схожую антиопухолевую роль. Дальнейшие экспериментальные исследования, например, оверэкспрессия 5′-изоформ hsa-miR-203a-3p и других микроРНК in vitro и in vivo, необходимы для понимания молекулярных механизмов развития и прогрессирования колоректальных опухолей.