El análisis del pangenoma revela tanto «especialistas» específicos del nicho como microbios que desempeñan funciones diversas en los microbiomas del cáncer colorrectal.

Cada vez hay más evidencia que relaciona el microbioma intestinal con un papel clave en la patogénesis del cáncer colorrectal (CCR),^[1-3], sin embargo, traducir estas asociaciones en herramientas de prevención del CCR basadas en el microbioma ha resultado ser un desafío. Aunque algunos taxones bacterianos específicos están relacionados de manera reproducible con el CCR, muchas especies y subespecies bacterianas intestinales a menudo difieren entre estudios, variando según la geografía, la demografía y el historial clínico.^[4-6] Esta variabilidad puede reflejar una propiedad fundamental de los microbiomas del CCR como ecosistemas definidos por un estado funcional convergente que puede ser alcanzado por taxones filogenéticamente diversos, es decir, una comunidad polifilética análoga a los microbiomas en las enfermedades inflamatorias intestinales.^[7]
,
^[8]

Históricamente, los esfuerzos para analizar el papel del microbioma en la enfermedad se han basado principalmente en la elaboración de perfiles taxonómicos.^[9] Si bien este enfoque ha identificado varios organismos asociados al CCR de manera reproducible, implícitamente trata la abundancia de especies como un indicador de la capacidad funcional relevante para la enfermedad. Esta suposición tiene limitaciones en dos aspectos importantes. Primero, las especies bacterianas exhiben una variación genómica sustancial a nivel de cepa, de modo que los genes que contribuyen a las interacciones con el huésped, la adquisición de nutrientes o la aptitud microbiana pueden estar presentes solo en un subconjunto de cepas.^[10] En segundo lugar, los organismos filogenéticamente diferentes pueden codificar funciones similares a través de la influencia de la evolución convergente facilitada por la transferencia horizontal de genes (THG). La elaboración de perfiles a nivel de especie puede, por lo tanto, pasar por alto tanto los rasgos asociados a la enfermedad restringidos a cepas particulares como las capacidades funcionales compartidas distribuidas en organismos taxonómicamente diversos.

Estas limitaciones motivan un marco evolutivo funcional que evalúa el contenido de genes asociados a la enfermedad de forma independiente de la abundancia de especies. Dicho marco puede determinar si un rasgo asociado se conserva dentro de un linaje específico o se dispersa entre cepas filogenéticamente distantes, distinguiendo así las asociaciones específicas del linaje de las diferencias a nivel de cepa y la convergencia funcional más amplia. Para aplicar este marco, analizamos los grupos de genes coasociados (CAG) identificados en nuestro trabajo anterior.^[3]
,
^[5] Los CAG son grupos de genes codificantes de proteínas que coexisten y se cuantifican juntos, análogos a una unidad taxonómica operativa a nivel de especie, pero resueltos al nivel del contenido de genes a nivel de subespecie. Este enfoque puede detectar asociaciones que son portadas solo por un pequeño subconjunto de cepas y puede identificar elementos funcionales relacionados que se repiten en diferentes contextos genómicos. Luego, medimos la coherencia filogenética de cada CAG en relación con la filogenia del genoma central correspondiente para evaluar si el contenido de genes asociado al CCR y a la salud se conservaba predominantemente dentro de un linaje o se dispersaba filogenéticamente.

Una hipótesis central de nuestro estudio es que el microbioma del CCR se define por rasgos funcionales convergentes. Al llevar a cabo esta investigación, reconocemos que la hipótesis alternativa, que el microbioma del CCR se define por CAG específicos, no es mutuamente excluyente y, de hecho, está respaldada por evidencia, incluida la nuestra.^[3]

Nuestra investigación se basa en dos objetivos con implicaciones translacionales complementarias. Primero, buscamos identificar la arquitectura funcional compartida del microbioma del CCR ("jugadores versátiles"), razonando que esto está integrado en el contenido de genes funcionales polifiléticos de los CAG asociados al CCR, lo que refleja una adaptación oportunista que confiere aptitud.^[11] En segundo lugar, al identificar las características monofiléticas de los CAG asociados al CCR específicos del nicho ("especialistas"), nuestro objetivo es identificar los rasgos específicos del linaje que pueden servir como objetivos para nuevas estrategias de detección o intervención terapéutica basadas en el microbioma del CCR. Prevemos una combinación de "especialistas" y "jugadores versátiles" en el CCR y en la salud, y proponemos que distinguir entre estos dos modos es importante desde una perspectiva translacional: los rasgos "especialistas" conservados en el linaje ofrecen principalmente objetivos reproducibles para el desarrollo de biomarcadores (y, si se pueden identificar rasgos metabólicos específicos del linaje asociados, también podrían motivar el diseño de inhibidores de moléculas pequeñas específicos del taxón), mientras que la arquitectura funcional compartida de los rasgos "jugadores versátiles" puede revelar objetivos de intervención farmacológica que son robustos a la variabilidad interespecífica entre estudios (y también podría ser útil para el diseño de biomarcadores que identifiquen un rasgo funcional). Al enmarcar la investigación del microbioma del CCR en términos de la distinción evolutiva y ecológica entre la asociación monofilética y polifilética, nuestro trabajo ofrece tanto una prueba de principio mecanicista para vincular la ecología microbiana intestinal con la carcinogénesis colorrectal como una estrategia sistemática para traducir la investigación del microbioma en herramientas clínicas.

Resultados

Un análisis del pan-genoma de las bacterias intestinales asociadas con el CCR, la salud o ambas en tres cohortes independientes

Con el fin de evaluar las contribuciones a nivel de gen al CCR, generamos bases de datos de referencia del pan-genoma para un conjunto de 12 taxones bacterianos que se habían identificado previamente como enriquecidos en el CCR (n = 8) o en la salud (n = 4) (Tabla S1).^[3]
,
^[12]
,
^[13] También estudiamos Faecalibacterium prausnitzii (Tabla S1), que se ha descrito en una relación compleja con el CCR según un metaanálisis de estudios publicados del microbioma.^[3] El pan-genoma de cada especie se construyó como un catálogo de genes codificantes de proteínas no redundantes agrupados en grupos según los patrones de coabundancia en los genomas de referencia (Figura 1A). Se utilizaron tres conjuntos de datos de estudios del microbioma publicados que comprenden un total de 379 individuos con CCR y 347 controles sanos, extraídos de tres contextos geográficos y dietéticos distintos (Japón, Estados Unidos y Europa), para cuantificar las abundancias de los grupos de genes en el CCR y en la salud (Yachida 2019;^[14] Hannigan 2018,^[15] Zeller 2014.^[16]) Las abundancias relativas de cada grupo de genes se calcularon de forma independiente para cada muestra de microbioma como el número de lecturas alineadas de forma única a los genes dentro del grupo, dividido por la longitud total del gen, normalizado al número de lecturas alineadas al pan-genoma total (lecturas por kilobase por millón de lecturas [RPKM]) (Figura 1A).

En el conjunto de datos de Yachida, 509 de un total de 3711 (13,7 %) grupos de genes diferían significativamente entre el CCR y los controles (300 enriquecidos en el CCR, 209 enriquecidos en la salud) (p < 0,05, prueba U de Mann-Whitney) en las 13 especies (Tabla S2A). En el conjunto de datos de Hannigan, 170 de un total de 3509 (4,8 %) grupos de genes diferían significativamente entre el CCR y los controles (85 enriquecidos en el CCR, 85 enriquecidos en la salud) (p < 0,05, prueba U de Mann-Whitney) en las 13 especies (Tabla S2A). En el conjunto de datos de Zeller, 885 de un total de 3683 (24,0 %) grupos de genes diferían significativamente entre el CCR y los controles (p < 0,05, prueba U de Mann-Whitney) en las 13 especies (Tabla S2A). En general, estas observaciones reflejaron nuestro informe anterior en el que encontramos que 2319 de 22295 CAG (10,4 %) identificados en 3 conjuntos de datos públicos estaban asociados con el CCR o la salud.^[3]

Después de aplicar el procedimiento de Benjamini-Hochberg al conjunto de datos de Yachida para corregir la prueba de hipótesis múltiple y minimizar la tasa de falsos descubrimientos (FDR-BH, con un umbral establecido en 0,20, prueba U de Mann-Whitney), 276 de un total de 3711 (7,4 %) grupos de genes del pan-genoma siguieron siendo significativamente diferentes entre el CCR y los controles (149 enriquecidos en el CCR, 127 enriquecidos en la salud) en las 13 especies (Figura 1B, Tabla S2B). Para las especies que contenían grupos de genes significativamente asociados, el número de grupos de genes significativos osciló entre 2 (B. animalis) y 83 (F. animalis), con 11 especies que tenían al menos dos grupos de genes significativos y 10 que tenían tanto grupos de genes enriquecidos en el CCR como en la salud, lo que indica una heterogeneidad importante a nivel de especie.

Entre los 3505 grupos de genes cuantificados en las tres cohortes, 92 grupos de genes estaban asociados positivamente con el CCR en exactamente dos cohortes, en comparación con una media derivada de una permutación de 77,5, lo que sitúa el recuento observado en el percentil 99 de la distribución nula (Tabla S2C). Dos estaban asociados positivamente con el CCR en las tres cohortes, en comparación con una media derivada de una permutación de 1,16 (percentil 88 de la distribución nula).

Interpretamos este patrón en el contexto de una heterogeneidad biológica y técnica sustancial entre las tres cohortes. Las tasas de detección por cohorte diferían en varios pliegues (el 4,8 % de los grupos de genes fueron significativos en Hannigan frente al 24,0 % en Zeller), y los estudios difieren en la plataforma de secuenciación, la profundidad de la secuenciación y la población de origen. La concordancia a nivel de grupos de genes individuales no es, por lo tanto, el resultado esperado de la adición de cohortes, y la afirmación central de este trabajo se refiere a la convergencia a nivel de la función codificada en lugar de la reproducibilidad de grupos de genes específicos en las poblaciones.

Para ilustrar los contextos genómicos de estas asociaciones, seleccionamos el grupo 81 de F. animalis, que fue uno de los grupos de genes más fuertemente enriquecidos en el CCR, teniendo en cuenta los tres estudios (Fig. S3, Tabla S2D). Los genomas de F. animalis que contenían este grupo se distribuyeron de forma parafilética en 24 de los 73 genomas de referencia utilizados para construir el pan-genoma (Figura 1C), lo que significa que no eran exclusivos de ningún clado a nivel de subespecie. Los genes contenidos en este grupo se encuentran adyacentes dentro de todos los genomas de referencia, lo que sugiere una transcripción coordinada como parte de un operón, con anotaciones tanto para el metabolismo como para el transporte de azúcares (Figura 1D). Estimamos qué proporción de copias del genoma de F. animalis contiene este operón normalizándolo con el grupo del genoma central, y observamos que el grupo 81 estaba presente en proporciones similares en cada uno de los grupos de CCR y de salud, respectivamente, entre los tres estudios (Figura 1E). La mayor aptitud relativa de las cepas de F. animalis que contienen este operón (como se manifiesta por su mayor proporción en la población) podría tener una conexión causal con el CCR, ya sea aumentando la probabilidad de su desarrollo o confiriendo una ventaja de crecimiento en el entorno de un microbioma del CCR.^[17]

Las distribuciones de las puntuaciones de monofilía son comparables en el CCR y la salud

A continuación, preguntamos si los grupos del pan-genoma asociados al CCR y a la salud eran, en general, representativos de clados monofiléticos (es decir, "especialistas") dentro de cada una de las especies incluidas en este análisis. Para evaluar el grado de monofilía observado para los grupos del pan-genoma asociados al CCR y a la salud, definimos una "puntuación de monofilía" que se podía aplicar a cualquier grupo del pan-genoma. Esta puntuación se calculó como la proporción de genomas que codifican el gen en relación con el número de genomas presentes en el clado monofilético más pequeño que abarca este gen (Figura 2A). Interpretamos esta puntuación como una medida de la coherencia filogenética. Una puntuación alta indica que un grupo de genes está en gran medida restringido a un solo linaje, lo que es consistente con la herencia vertical. Una puntuación más baja indica una amplia distribución que requiere ganancias o pérdidas repetidas de genes en la filogenia de la especie. Si bien este patrón es compatible con la transferencia horizontal de genes, también puede surgir de un polimorfismo antiguo o de la pérdida diferencial de genes y, por lo tanto, no establece la THG por sí solo. La distribución de los genes por puntuación de monofilía varió ampliamente entre las especies, con grados similares de variación observados en todas las especies analizadas (Figura 2B–D).

Los grupos asociados a la salud de F. prausnitzii presentaron puntuaciones de monofilia general más altas en comparación con los grupos no asociados (p corregida por FDR = 0,023, prueba U de Mann-Whitney), mientras que los grupos asociados al CRC de F. prausnitzii mostraron puntuaciones de monofilia más bajas en relación con los grupos no asociados (p corregida por FDR = 0,047; Figura 2E). Estos hallazgos son coherentes con la presencia de características asociadas a la salud dentro de clados "especialistas" conservados a nivel de linaje, mientras que la distribución de los grupos de genes asociados al CRC parece ser más consistente con elementos genéticos polifiléticos y potencialmente móviles que actúan como "actividades secundarias". Este fue un hallazgo único de F. prausnitzii, ya que las otras especies mostraron distribuciones de monofilia similares independientemente de la asociación con el CRC. Tras la corrección para múltiples hipótesis, las puntuaciones de monofilia de las especies distintas de F. prausnitzii no fueron significativamente diferentes entre los microbiomas de CRC y los microbiomas sanos, ni en comparación con los grupos no asociados (la distribución nula para esta prueba). El análisis de los valores de significación con signo en función de la puntuación de monofilia reveló una distribución bimodal, coherente con la coexistencia de grupos de genes "especialistas" y "actividades secundarias" dentro de cada especie, sin un enriquecimiento claro de ninguno de los patrones en el CRC en comparación con la salud (Figura 2E, F).

No obstante, el examen de grupos de genes individuales asociados al CRC reveló temas funcionales compartidos, coherentes con un comportamiento de "actividad secundaria", que incluyen la adquisición de nutrientes y la maquinaria de interacción con el huésped, distribuidos en especies filogenéticamente diversas. Basándonos en la significación y la relevancia de la anotación funcional, destacamos los grupos del pangenoma asociados al CRC de B. caccae (grupo 85, puntuación de monofilia = 0,27), B. fragilis (grupo 15, puntuación de monofilia = 0,032) y F. animalis (grupo 11, monofilia = 0,14), que estaban presentes en mayor abundancia relativa en las muestras de CRC en comparación con los controles sanos (Figura 3). Estos grupos de genes contienen genes anotados para RagB/SusD, receptor dependiente de TonB y secreción de tipo IV, respectivamente (Figura 3A-I). La distribución parafilética de estos grupos en cada especie es coherente con la adquisición o pérdida repetida, que podría incluir la adquisición horizontal, un tema que se repite en todo el pangenoma asociado al CRC y que se examina con más detalle a continuación.

Los grupos de genes del pangenoma de Fusobacterium que codifican la regulación de la energía y la secreción de tipo IV están enriquecidos en el CRC

A la luz de la literatura que relaciona múltiples especies de Fusobacterium con los cánceres del tracto gastrointestinal,[17, 18] analizamos los grupos del pangenoma de forma más amplia en todo el género Fusobacterium. Tras la corrección para múltiples pruebas de hipótesis (FDR-BH = 0,2), se encontró que 121 grupos estaban asociados al CRC y 3 estaban asociados a la salud (Figura 4A). La distribución de monofilia de estos grupos no mostró diferencias significativas con los grupos no asociados, lo que sugiere que el contenido de genes de Fusobacterium asociados al CRC no muestra un enriquecimiento a nivel de género para ninguno de los comportamientos de "especialista" o "actividad secundaria" en conjunto (Figura S4). Del grupo de grupos de genes enriquecidos en los microbiomas de CRC, seleccionamos para una evaluación adicional tres grupos de alta significación y relevancia de la anotación funcional (Figura 4B-D). Estos grupos tenían puntuaciones de monofilia que oscilaban entre 0,23 (grupo 200) y 0,44 (grupo 37) hasta 0,76 (grupo 58), lo que destaca la diversidad de las distribuciones filogenéticas de los genes asociados al CRC, incluso dentro de un género: algunos grupos fuertemente asociados al CRC están restringidos a regiones locales dentro del árbol filogenético de Fusobacterium, mientras que otros están dispersos de forma polifilética. Los genes funcionalmente relevantes codificados en estos grupos de genes incluían un regulador de respuesta de co-agregación (Figura 4E), un transductor de energía TonB (Figura 4F) y un sistema de secreción de tipo IV, respectivamente (revisado en Discusión). La secreción de tipo IV, un tema funcional recurrente en los grupos asociados al CRC que se examinan con más detalle a continuación, se codificó en un grupo de baja monofilia (grupo 200, monofilia = 0,23), lo que es coherente con un mecanismo de "actividad secundaria" asociado al cáncer que podría compartirse en los sublinajes de Fusobacterium.

Estabilidad temporal de los grupos de genes dentro de los microbiomas intestinales individuales

Una cuestión práctica para traducir los hallazgos de los estudios del microbioma del CRC al diseño de una prueba de detección basada en el microbioma es si las señales microbianas clave son estables en el tiempo. Por lo tanto, a continuación, buscamos caracterizar la estabilidad temporal de nuestro marco basado en grupos del pangenoma y definir las escalas de tiempo en las que nuestro marco podría aplicarse. Se evaluó la robustez de nuestros análisis utilizando muestras longitudinales del grupo de Fukuyama,^[19], realizando cálculos de autocorrelación de las densidades de los grupos dentro y entre los individuos. En todas las especies bacterianas, las correlaciones de abundancia de los grupos intraindividuales fueron más altas que las correlaciones entre los individuos. (Δ mediana intra- menos interindividual osciló entre +0,06 y +0,41, p < 10-10) (Figura 5A).

Para la mayoría de las especies que examinamos, la estructura a nivel de grupo también se mantuvo estable antes y después de la limpieza del colon, lo que sugiere que incluso las perturbaciones agudas importantes no inducen cambios sustanciales en las abundancias de los genes de los grupos microbianos (Δ mediana pre- menos post-preparación intestinal osciló entre -0,00 y +0,04, p < 0,04) (Figura 5B y S5). Además, en diferentes intervalos de tiempo, las correlaciones de densidad de los grupos del pangenoma fueron abrumadoramente invariantes, lo que apoya la suposición de que la estructura del grupo del pangenoma es temporalmente persistente (Figura 5B). Cabe destacar que dos especies bacterianas, Escherichia coli y Fusobacterium animalis, mostraron correlaciones intraindividuales dependientes del tiempo comparativamente más bajas, lo que sugiere un nicho más precario en el ecosistema microbiano intestinal y, por lo tanto, un mayor potencial de variabilidad temporal en sus perfiles del pangenoma en relación con otros taxones.

Asociaciones sinérgicas interespecies de grupos de genes en el CRC

Imaginando que las relaciones ecológicas y evolutivas entre las bacterias en los microbiomas del CRC están integradas en nuestro análisis de grupos de genes, a continuación, preguntamos si los grupos de genes de diferentes especies podrían actuar juntos, reforzando o oponiéndose a su asociación con el CRC. Para probar esto, evaluamos todas las combinaciones de pares de grupos utilizando modelos de regresión logística binaria que incorporan términos de interacción. Los modelos se clasificaron por la significación nominal del término de interacción, y las 200 asociaciones principales se conservaron para una evaluación exploratoria (Tabla S10). Dado el espacio de búsqueda de alta dimensión, estos resultados se utilizaron para priorizar las asociaciones candidatas no aditivas para la interpretación biológica en lugar de establecer interacciones interespecies directas. Entre las interacciones destacadas, seis de los ocho grupos del pangenoma tenían puntuaciones de monofilia bajas (